[musl] catan(z)

[musl] catan(z)

[Damian McGuckin]

In this routine, there are 2 lines of code

     t = 0.5 * atan2(2.0 * x, a);
     w = _redupi(t);

The first computes atan2() which returns a number in the range [-pi,+pi]
which means that t is a number in the range [-pi/2,+pi/2].

As far as I understand, the routine _redupi(t) accepts a argument and
reduces it into the range [-pi, +pi]. Am I mistaken?

If I am not wrong, and 't' is already in the range that _redupi(t) will 
return, then _redupi() just returns its argument (within a rounding 
error). So it seems like a no-op to me.

Can anybody comment?

The routine is originally from Moshier's Cephes library. I have Moshier's 
book but I cannot find any hints there.

Thanks - Damian

Re: [musl] catan(z)

[Szabolcs Nagy]

* Damian McGuckin <damianm@esi.com.au> [2024-08-12 00:01:01 +1000]:
> 
> In this routine, there are 2 lines of code
> 
>     t = 0.5 * atan2(2.0 * x, a);
>     w = _redupi(t);
> 
> The first computes atan2() which returns a number in the range [-pi,+pi]
> which means that t is a number in the range [-pi/2,+pi/2].
> 
> As far as I understand, the routine _redupi(t) accepts a argument and
> reduces it into the range [-pi, +pi]. Am I mistaken?

*reduces into [-pi/2, pi/2]

> 
> If I am not wrong, and 't' is already in the range that _redupi(t) will
> return, then _redupi() just returns its argument (within a rounding error).
> So it seems like a no-op to me.
> 
> Can anybody comment?

it only changes the endpoints (t=+-pi/2 case) otherwise no-op

i dont know if we care about the endpoints but even if we do
there are better ways to fix them up

i think _redupi can be removed, nice catch.

> 
> The routine is originally from Moshier's Cephes library. I have Moshier's
> book but I cannot find any hints there.
> 
> Thanks - Damian

Re: [musl] catan(z)

[Damian McGuckin]

On Sun, 11 Aug 2024, Szabolcs Nagy wrote:

> * Damian McGuckin <damianm@esi.com.au> [2024-08-12 00:01:01 +1000]:
>>
>> In this routine, there are 2 lines of code
>>
>>     t = 0.5 * atan2(2.0 * x, a);
>>     w = _redupi(t);
>>
>> The first computes atan2() which returns a number in the range [-pi,+pi]
>> which means that t is a number in the range [-pi/2,+pi/2].
>>
>> As far as I understand, the routine _redupi(t) accepts a argument and
>> reduces it into the range [-pi, +pi]. Am I mistaken?
>
> *reduces into [-pi/2, pi/2]

Yes. Silly me. Why? Because ...

 	atan2() returns a number in [-pi, +pi], 't' is in [-pi/2,+pi/2],

 	hence, at least in this case, _redupi(t) just maps that 't' into
 	that same range.

There is some argument that if you handle the special cases at infinity 
separately (which I think MUSL should do but I do not have time at the 
moment), then one can assume that because pi/2 is irrational, then one 
should never have to deal with the end points in the chunk of code where 
those two lines of code seen above should appear. I will have a chat 
sometime with the guy who wrote that logic in a WG14 paper when I get a 
really clear head and can line him up.

Thanks - Damian

Re: [musl] catan(z)

[Damian McGuckin]

On Mon, 12 Aug 2024, Damian McGuckin wrote:

> There is some argument that if you handle the special cases at infinity 
> separately (which I think MUSL should do but I do not have time at the 
> moment), then one can assume that because pi/2 is irrational, then one 
> should never have to deal with the end points in the chunk of code where 
> those two lines of code seen above should appear. I will have a chat 
> sometime with the guy who wrote that logic in a WG14 paper when I get a 
> really clear head and can line him up.

Consider

 	atan2(y, x)

For any finite y and finite non-zero x floating point number arguments, 
i.e. rational numbers, the result of atan2(y, x) must be rational and so 
is never +/- pi (which is irrational and only occurs when the ration y/x 
is a mathematical infinity, not an overflowing infinity). So, we can 
ignore the endpoints as long as our special case handling takes care of 
the case of zero x.

I think that is correct .... or is my brain still not working properly 
after too many late nights watching the Olympics.

Thanks - Damian

Re: [musl] catan(z)

[Szabolcs Nagy]

* Damian McGuckin <damianm@esi.com.au> [2024-08-12 13:56:18 +1000]:
> On Mon, 12 Aug 2024, Damian McGuckin wrote:
> 
> > There is some argument that if you handle the special cases at infinity
> > separately (which I think MUSL should do but I do not have time at the
> > moment), then one can assume that because pi/2 is irrational, then one
> > should never have to deal with the end points in the chunk of code where
> > those two lines of code seen above should appear. I will have a chat
> > sometime with the guy who wrote that logic in a WG14 paper when I get a
> > really clear head and can line him up.
> 
> Consider
> 
> 	atan2(y, x)
> 
> For any finite y and finite non-zero x floating point number arguments, i.e.
> rational numbers, the result of atan2(y, x) must be rational and so is never
> +/- pi (which is irrational and only occurs when the ration y/x is a
> mathematical infinity, not an overflowing infinity). So, we can ignore the
> endpoints as long as our special case handling takes care of the case of
> zero x.
> 
> I think that is correct .... or is my brain still not working properly after
> too many late nights watching the Olympics.

well atan2 cannot return exact +-pi, but one would still
expect the right sign for the quadrant (y=+-0, x<0 case
or finite y, x=-inf case)

and then _redupi uses (int)(t/pi+-0.5) where that /pi is
not exact the same way the +-pi returned from atan2 was
not exact.

i.e currently in the atan2==+-pi case the sign gets flipped
by _redupi (even though double pi < exact pi so no need to
reduce)

i doubt the sign flip is needed, but even if that's the
right result there are better ways to handle that case.

Re: [musl] catan(z)

[Rich Felker]

On Mon, Aug 12, 2024 at 11:00:22AM +0200, Szabolcs Nagy wrote:
> * Damian McGuckin <damianm@esi.com.au> [2024-08-12 13:56:18 +1000]:
> > On Mon, 12 Aug 2024, Damian McGuckin wrote:
> > 
> > > There is some argument that if you handle the special cases at infinity
> > > separately (which I think MUSL should do but I do not have time at the
> > > moment), then one can assume that because pi/2 is irrational, then one
> > > should never have to deal with the end points in the chunk of code where
> > > those two lines of code seen above should appear. I will have a chat
> > > sometime with the guy who wrote that logic in a WG14 paper when I get a
> > > really clear head and can line him up.
> > 
> > Consider
> > 
> > 	atan2(y, x)
> > 
> > For any finite y and finite non-zero x floating point number arguments, i.e.
> > rational numbers, the result of atan2(y, x) must be rational and so is never
> > +/- pi (which is irrational and only occurs when the ration y/x is a
> > mathematical infinity, not an overflowing infinity). So, we can ignore the
> > endpoints as long as our special case handling takes care of the case of
> > zero x.
> > 
> > I think that is correct .... or is my brain still not working properly after
> > too many late nights watching the Olympics.
> 
> well atan2 cannot return exact +-pi, but one would still
> expect the right sign for the quadrant (y=+-0, x<0 case
> or finite y, x=-inf case)
> 
> and then _redupi uses (int)(t/pi+-0.5) where that /pi is
> not exact the same way the +-pi returned from atan2 was
> not exact.
> 
> i.e currently in the atan2==+-pi case the sign gets flipped
> by _redupi (even though double pi < exact pi so no need to
> reduce)
> 
> i doubt the sign flip is needed, but even if that's the
> right result there are better ways to handle that case.

If I'm following, it sounds like the current behavior is either a
no-op or wrong and removing the _redupi step entirely would make it
correct or at least closer to correct. Do you agree?

Any suggestions for a test to run? I did some basic smoke testing for
large positive real inputs and got expected results.

Rich

Re: [musl] catan(z)

[Damian McGuckin]

On Wed, 14 Aug 2024, Rich Felker wrote:

> If I'm following, it sounds like the current behavior is either a
> no-op or wrong and removing the _redupi step entirely would make it
> correct or at least closer to correct. Do you agree?

Just a no-op.  Unless there is some subtle reason why Moshier did this why 
I am too dumb to understand.  I think the issue was that back in the 80s 
when he wrote it, there were some unstandardised atan2() routines which 
(Moshier alluded to in his book and which) returned results in [0,2*pi] 
rather than the now standardized routines with return in [-pi,+pi] and he 
was handling those issues.

> Any suggestions for a test to run? I did some basic smoke testing for
> large positive real inputs and got expected results.

Not off hand. Other issues taking my brain cycles right nnow. Young Mr 
Nagy might have a better idea.

Thanks - Damian

Re: [musl] catan(z)

[Rich Felker]

[-- Attachment #1: Type: text/plain, Size: 1298 bytes --]

On Thu, Aug 15, 2024 at 12:12:35PM +1000, Damian McGuckin wrote:
> On Wed, 14 Aug 2024, Rich Felker wrote:
> 
> >If I'm following, it sounds like the current behavior is either a
> >no-op or wrong and removing the _redupi step entirely would make it
> >correct or at least closer to correct. Do you agree?
> 
> Just a no-op.  Unless there is some subtle reason why Moshier did
> this why I am too dumb to understand.  I think the issue was that
> back in the 80s when he wrote it, there were some unstandardised
> atan2() routines which (Moshier alluded to in his book and which)
> returned results in [0,2*pi] rather than the now standardized
> routines with return in [-pi,+pi] and he was handling those issues.
> 
> >Any suggestions for a test to run? I did some basic smoke testing for
> >large positive real inputs and got expected results.
> 
> Not off hand. Other issues taking my brain cycles right nnow. Young
> Mr Nagy might have a better idea.

OK, that makes sense. Proposed patch attached. I think the code could
(maybe should) actually be reordered and simplified at some point
without confusing reuse of temporaries and without temp storage of the
real result in a complex var, but I made the patch a minimal
bad-code-removal for now. Further cleanup could be done separately.

Rich

[-- Attachment #2: 0001-catan-remove-no-up-reduction-mod-pi-and-unused-code.patch --]
[-- Type: text/plain, Size: 3400 bytes --]

From 392529fb945841d96981235d16f8c1b926f134e4 Mon Sep 17 00:00:00 2001
From: Rich Felker <dalias@aerifal.cx>
Date: Thu, 15 Aug 2024 10:07:36 -0400
Subject: [PATCH] catan: remove no-up reduction mod pi and unused code

the output of atan2 is already in the correct range and does not need
further reduction. the MAXNUM macros were both unused and incorrect.
---
 src/complex/catan.c  | 25 +------------------------
 src/complex/catanf.c | 28 +---------------------------
 src/complex/catanl.c | 24 +-----------------------
 3 files changed, 3 insertions(+), 74 deletions(-)

diff --git a/src/complex/catan.c b/src/complex/catan.c
index ccc2fb53..b4fe552a 100644
--- a/src/complex/catan.c
+++ b/src/complex/catan.c
@@ -60,29 +60,6 @@
 
 #include "complex_impl.h"
 
-#define MAXNUM 1.0e308
-
-static const double DP1 = 3.14159265160560607910E0;
-static const double DP2 = 1.98418714791870343106E-9;
-static const double DP3 = 1.14423774522196636802E-17;
-
-static double _redupi(double x)
-{
-	double t;
-	long i;
-
-	t = x/M_PI;
-	if (t >= 0.0)
-		t += 0.5;
-	else
-		t -= 0.5;
-
-	i = t;  /* the multiple */
-	t = i;
-	t = ((x - t * DP1) - t * DP2) - t * DP3;
-	return t;
-}
-
 double complex catan(double complex z)
 {
 	double complex w;
@@ -95,7 +72,7 @@ double complex catan(double complex z)
 	a = 1.0 - x2 - (y * y);
 
 	t = 0.5 * atan2(2.0 * x, a);
-	w = _redupi(t);
+	w = t;
 
 	t = y - 1.0;
 	a = x2 + (t * t);
diff --git a/src/complex/catanf.c b/src/complex/catanf.c
index 1d569f2d..faaa907a 100644
--- a/src/complex/catanf.c
+++ b/src/complex/catanf.c
@@ -55,32 +55,6 @@
 
 #include "complex_impl.h"
 
-#define MAXNUMF 1.0e38F
-
-static const double DP1 = 3.140625;
-static const double DP2 = 9.67502593994140625E-4;
-static const double DP3 = 1.509957990978376432E-7;
-
-static const float float_pi = M_PI;
-
-static float _redupif(float xx)
-{
-	float x, t;
-	long i;
-
-	x = xx;
-	t = x/float_pi;
-	if (t >= 0.0f)
-		t += 0.5f;
-	else
-		t -= 0.5f;
-
-	i = t;  /* the multiple */
-	t = i;
-	t = ((x - t * DP1) - t * DP2) - t * DP3;
-	return t;
-}
-
 float complex catanf(float complex z)
 {
 	float complex w;
@@ -93,7 +67,7 @@ float complex catanf(float complex z)
 	a = 1.0f - x2 - (y * y);
 
 	t = 0.5f * atan2f(2.0f * x, a);
-	w = _redupif(t);
+	w = t;
 
 	t = y - 1.0f;
 	a = x2 + (t * t);
diff --git a/src/complex/catanl.c b/src/complex/catanl.c
index e62526c0..cd2d2b00 100644
--- a/src/complex/catanl.c
+++ b/src/complex/catanl.c
@@ -67,28 +67,6 @@ long double complex catanl(long double complex z)
 	return catan(z);
 }
 #else
-static const long double PIL = 3.141592653589793238462643383279502884197169L;
-static const long double DP1 = 3.14159265358979323829596852490908531763125L;
-static const long double DP2 = 1.6667485837041756656403424829301998703007e-19L;
-static const long double DP3 = 1.8830410776607851167459095484560349402753e-39L;
-
-static long double redupil(long double x)
-{
-	long double t;
-	long i;
-
-	t = x / PIL;
-	if (t >= 0.0L)
-		t += 0.5L;
-	else
-		t -= 0.5L;
-
-	i = t;  /* the multiple */
-	t = i;
-	t = ((x - t * DP1) - t * DP2) - t * DP3;
-	return t;
-}
-
 long double complex catanl(long double complex z)
 {
 	long double complex w;
@@ -101,7 +79,7 @@ long double complex catanl(long double complex z)
 	a = 1.0L - x2 - (y * y);
 
 	t = atan2l(2.0L * x, a) * 0.5L;
-	w = redupil(t);
+	w = t;
 
 	t = y - 1.0L;
 	a = x2 + (t * t);
-- 
2.21.0

Re: [musl] catan(z)

[Szabolcs Nagy]

* Rich Felker <dalias@libc.org> [2024-08-15 10:28:25 -0400]:
> OK, that makes sense. Proposed patch attached. I think the code could
> (maybe should) actually be reordered and simplified at some point
> without confusing reuse of temporaries and without temp storage of the
> real result in a complex var, but I made the patch a minimal
> bad-code-removal for now. Further cleanup could be done separately.

patch looks good

Re: [musl] catan(z)

[Eric Pruitt]

On Thu, Aug 15, 2024 at 10:28:25AM -0400, Rich Felker wrote:
> OK, that makes sense. Proposed patch attached.

You misspelled "no-op" as "no-up" in the patch summary.

Eric

Re: [musl] catan(z)

[Morten Welinder]

atan2 definitely isn't supposed to always output rational numbers on
rational input.

atan2(+0,-1) is Pi.
atan2(-0,-1) is -Pi
atan2(1,1) is Pi/4 -- clearly not a rational number.

These aren't exactly rational results (unless you mean their
floating-point approximations).

M.

On Sun, Aug 11, 2024 at 11:56 PM Damian McGuckin <damianm@esi.com.au> wrote:
>
> On Mon, 12 Aug 2024, Damian McGuckin wrote:
>
> > There is some argument that if you handle the special cases at infinity
> > separately (which I think MUSL should do but I do not have time at the
> > moment), then one can assume that because pi/2 is irrational, then one
> > should never have to deal with the end points in the chunk of code where
> > those two lines of code seen above should appear. I will have a chat
> > sometime with the guy who wrote that logic in a WG14 paper when I get a
> > really clear head and can line him up.
>
> Consider
>
>         atan2(y, x)
>
> For any finite y and finite non-zero x floating point number arguments,
> i.e. rational numbers, the result of atan2(y, x) must be rational and so
> is never +/- pi (which is irrational and only occurs when the ration y/x
> is a mathematical infinity, not an overflowing infinity). So, we can
> ignore the endpoints as long as our special case handling takes care of
> the case of zero x.
>
> I think that is correct .... or is my brain still not working properly
> after too many late nights watching the Olympics.
>
> Thanks - Damian

Re: [musl] catan(z)

[Rich Felker]

On Thu, Aug 15, 2024 at 09:18:19AM -0400, Morten Welinder wrote:
> atan2 definitely isn't supposed to always output rational numbers on
> rational input.
> 
> atan2(+0,-1) is Pi.
> atan2(-0,-1) is -Pi
> atan2(1,1) is Pi/4 -- clearly not a rational number.
> 
> These aren't exactly rational results (unless you mean their
> floating-point approximations).

The way I read it, the claim was not that the exact mathematical value
of atan2 for some argument is rational, but that the floating point
number returned by the C function is necessarily rational (because all
floating point numbers are diadic rationals) and thus never actually
equal to ±pi. This means you don't have any issue with whatever
happens exactly at the endpoints.

> On Sun, Aug 11, 2024 at 11:56 PM Damian McGuckin <damianm@esi.com.au> wrote:
> >
> > On Mon, 12 Aug 2024, Damian McGuckin wrote:
> >
> > > There is some argument that if you handle the special cases at infinity
> > > separately (which I think MUSL should do but I do not have time at the
> > > moment), then one can assume that because pi/2 is irrational, then one
> > > should never have to deal with the end points in the chunk of code where
> > > those two lines of code seen above should appear. I will have a chat
> > > sometime with the guy who wrote that logic in a WG14 paper when I get a
> > > really clear head and can line him up.
> >
> > Consider
> >
> >         atan2(y, x)
> >
> > For any finite y and finite non-zero x floating point number arguments,
> > i.e. rational numbers, the result of atan2(y, x) must be rational and so
> > is never +/- pi (which is irrational and only occurs when the ration y/x
> > is a mathematical infinity, not an overflowing infinity). So, we can
> > ignore the endpoints as long as our special case handling takes care of
> > the case of zero x.
> >
> > I think that is correct .... or is my brain still not working properly
> > after too many late nights watching the Olympics.
> >
> > Thanks - Damian

Re: [musl] catan(z)

[Damian McGuckin]

On Thu, 15 Aug 2024, Rich Felker wrote:

> On Thu, Aug 15, 2024 at 09:18:19AM -0400, Morten Welinder wrote:
>> atan2 definitely isn't supposed to always output rational numbers on
>> rational input.
>>
>> atan2(+0,-1) is Pi.
>> atan2(-0,-1) is -Pi
>> atan2(1,1) is Pi/4 -- clearly not a rational number.
>>
>> These aren't exactly rational results (unless you mean their
>> floating-point approximations).
>
> The way I read it, the claim was not that the exact mathematical value
> of atan2 for some argument is rational, but that the floating point
> number returned by the C function is necessarily rational (because all
> floating point numbers are diadic rationals) and thus never actually
> equal to ?pi. This means you don't have any issue with whatever
> happens exactly at the endpoints.

Precisely. Sorry if I was unclear Morten.

Actually, I think we need to ensure we use a value of Pi above which is 
rounded towards zero (rather than one rounded to nearest tied to even) so 
that it never lies outside the exact range of [-PI,+PI] where PI in this 
case is an irrational number that we cannot represent with an IEEE 754 
floating point number.

If we can enhance catan(z) to explicitly handle the special cases of Annex 
G of the C standard, some of this discussion becomes a lot simpler. But 
at the time he wrote this routine, Moshier did not have the luxury of a
standard to consult.

- Damian