```mailing list of musl libc
help / color / mirror / code / Atom feed```
```* [musl] catan(z)
@ 2024-08-11 14:01 Damian McGuckin
2024-08-11 20:08 ` Szabolcs Nagy
From: Damian McGuckin @ 2024-08-11 14:01 UTC (permalink / raw)
To: MUSL

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)
2024-08-11 14:01 [musl] catan(z) Damian McGuckin
@ 2024-08-11 20:08 ` Szabolcs Nagy
2024-08-12  3:39   ` Damian McGuckin
From: Szabolcs Nagy @ 2024-08-11 20:08 UTC (permalink / raw)
To: Damian McGuckin; +Cc: MUSL

* 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)
2024-08-11 20:08 ` Szabolcs Nagy
@ 2024-08-12  3:39   ` Damian McGuckin
2024-08-12  3:56     ` Damian McGuckin
From: Damian McGuckin @ 2024-08-12  3:39 UTC (permalink / raw)
To: MUSL

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)
2024-08-12  3:39   ` Damian McGuckin
@ 2024-08-12  3:56     ` Damian McGuckin
2024-08-12  9:00       ` Szabolcs Nagy
2024-08-15 13:18       ` Morten Welinder
0 siblings, 2 replies; 13+ messages in thread
From: Damian McGuckin @ 2024-08-12  3:56 UTC (permalink / raw)
To: MUSL

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)
2024-08-12  3:56     ` Damian McGuckin
@ 2024-08-12  9:00       ` Szabolcs Nagy
2024-08-15  1:16         ` Rich Felker
2024-08-15 13:18       ` Morten Welinder
From: Szabolcs Nagy @ 2024-08-12  9:00 UTC (permalink / raw)
To: Damian McGuckin; +Cc: MUSL

* 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)
2024-08-12  9:00       ` Szabolcs Nagy
@ 2024-08-15  1:16         ` Rich Felker
2024-08-15  2:12           ` Damian McGuckin
From: Rich Felker @ 2024-08-15  1:16 UTC (permalink / raw)
To: Damian McGuckin, MUSL

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)
2024-08-15  1:16         ` Rich Felker
@ 2024-08-15  2:12           ` Damian McGuckin
2024-08-15 14:28             ` Rich Felker
From: Damian McGuckin @ 2024-08-15  2:12 UTC (permalink / raw)
To: MUSL

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)
2024-08-12  3:56     ` Damian McGuckin
2024-08-12  9:00       ` Szabolcs Nagy
@ 2024-08-15 13:18       ` Morten Welinder
2024-08-15 13:44         ` Rich Felker
From: Morten Welinder @ 2024-08-15 13:18 UTC (permalink / raw)
To: musl

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)
2024-08-15 13:18       ` Morten Welinder
@ 2024-08-15 13:44         ` Rich Felker
2024-08-16  1:47           ` Damian McGuckin
From: Rich Felker @ 2024-08-15 13:44 UTC (permalink / raw)
To: Morten Welinder; +Cc: musl

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)
2024-08-15  2:12           ` Damian McGuckin
@ 2024-08-15 14:28             ` Rich Felker
2024-08-16  7:26               ` Szabolcs Nagy
2024-08-16  8:18               ` Eric Pruitt
0 siblings, 2 replies; 13+ messages in thread
From: Rich Felker @ 2024-08-15 14:28 UTC (permalink / raw)
To: Damian McGuckin; +Cc: MUSL

[-- 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)
2024-08-15 13:44         ` Rich Felker
@ 2024-08-16  1:47           ` Damian McGuckin
0 siblings, 0 replies; 13+ messages in thread
From: Damian McGuckin @ 2024-08-16  1:47 UTC (permalink / raw)
To: musl; +Cc: Morten Welinder

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

```* Re: [musl] catan(z)
2024-08-15 14:28             ` Rich Felker
@ 2024-08-16  7:26               ` Szabolcs Nagy
2024-08-16  8:18               ` Eric Pruitt
1 sibling, 0 replies; 13+ messages in thread
From: Szabolcs Nagy @ 2024-08-16  7:26 UTC (permalink / raw)
To: Rich Felker; +Cc: Damian McGuckin, MUSL

* 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)
2024-08-15 14:28             ` Rich Felker
2024-08-16  7:26               ` Szabolcs Nagy
@ 2024-08-16  8:18               ` Eric Pruitt
1 sibling, 0 replies; 13+ messages in thread
From: Eric Pruitt @ 2024-08-16  8:18 UTC (permalink / raw)
To: musl

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

```end of thread, other threads:[~2024-08-16  8:18 UTC | newest]

2024-08-11 14:01 [musl] catan(z) Damian McGuckin
2024-08-11 20:08 ` Szabolcs Nagy
2024-08-12  3:39   ` Damian McGuckin
2024-08-12  3:56     ` Damian McGuckin
2024-08-12  9:00       ` Szabolcs Nagy
2024-08-15  1:16         ` Rich Felker
2024-08-15  2:12           ` Damian McGuckin
2024-08-15 14:28             ` Rich Felker
2024-08-16  7:26               ` Szabolcs Nagy
2024-08-16  8:18               ` Eric Pruitt
2024-08-15 13:18       ` Morten Welinder
2024-08-15 13:44         ` Rich Felker
2024-08-16  1:47           ` Damian McGuckin
```

```Code repositories for project(s) associated with this public inbox

https://git.vuxu.org/mirror/musl/

This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for NNTP newsgroup(s).```