*[musl] catan(z)@ 2024-08-11 14:01 Damian McGuckin2024-08-11 20:08 ` Szabolcs Nagy 0 siblings, 1 reply; 13+ messages in thread 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 ^ permalink raw reply [flat|nested] 13+ messages in thread

*Re: [musl] catan(z)2024-08-11 14:01 [musl] catan(z) Damian McGuckin@ 2024-08-11 20:08 ` Szabolcs Nagy2024-08-12 3:39 ` Damian McGuckin 0 siblings, 1 reply; 13+ messages in thread 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 ^ permalink raw reply [flat|nested] 13+ messages in thread

*Re: [musl] catan(z)2024-08-11 20:08 ` Szabolcs Nagy@ 2024-08-12 3:39 ` Damian McGuckin2024-08-12 3:56 ` Damian McGuckin 0 siblings, 1 reply; 13+ messages in thread 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 ^ permalink raw reply [flat|nested] 13+ messages in thread

*Re: [musl] catan(z)2024-08-12 3:39 ` Damian McGuckin@ 2024-08-12 3:56 ` Damian McGuckin2024-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 ^ permalink raw reply [flat|nested] 13+ messages in thread

*Re: [musl] catan(z)2024-08-12 3:56 ` Damian McGuckin@ 2024-08-12 9:00 ` Szabolcs Nagy2024-08-15 1:16 ` Rich Felker 2024-08-15 13:18 ` Morten Welinder 1 sibling, 1 reply; 13+ messages in thread 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. ^ permalink raw reply [flat|nested] 13+ messages in thread

*Re: [musl] catan(z)2024-08-12 9:00 ` Szabolcs Nagy@ 2024-08-15 1:16 ` Rich Felker2024-08-15 2:12 ` Damian McGuckin 0 siblings, 1 reply; 13+ messages in thread 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 ^ permalink raw reply [flat|nested] 13+ messages in thread

*Re: [musl] catan(z)2024-08-15 1:16 ` Rich Felker@ 2024-08-15 2:12 ` Damian McGuckin2024-08-15 14:28 ` Rich Felker 0 siblings, 1 reply; 13+ messages in thread 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 ^ permalink raw reply [flat|nested] 13+ messages in thread

*Re: [musl] catan(z)2024-08-12 3:56 ` Damian McGuckin 2024-08-12 9:00 ` Szabolcs Nagy@ 2024-08-15 13:18 ` Morten Welinder2024-08-15 13:44 ` Rich Felker 1 sibling, 1 reply; 13+ messages in thread 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 ^ permalink raw reply [flat|nested] 13+ messages in thread

*Re: [musl] catan(z)2024-08-15 13:18 ` Morten Welinder@ 2024-08-15 13:44 ` Rich Felker2024-08-16 1:47 ` Damian McGuckin 0 siblings, 1 reply; 13+ messages in thread 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 ^ permalink raw reply [flat|nested] 13+ messages in thread

*Re: [musl] catan(z)2024-08-15 2:12 ` Damian McGuckin@ 2024-08-15 14:28 ` Rich Felker2024-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 ^ permalink raw reply [flat|nested] 13+ messages in thread

*2024-08-15 13:44 ` Rich FelkerRe: [musl] catan(z)@ 2024-08-16 1:47 ` Damian McGuckin0 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 ^ permalink raw reply [flat|nested] 13+ messages in thread

*Re: [musl] catan(z)2024-08-15 14:28 ` Rich Felker@ 2024-08-16 7:26 ` Szabolcs Nagy2024-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 ^ permalink raw reply [flat|nested] 13+ messages in thread

*Re: [musl] catan(z)2024-08-15 14:28 ` Rich Felker 2024-08-16 7:26 ` Szabolcs Nagy@ 2024-08-16 8:18 ` Eric Pruitt1 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 ^ permalink raw reply [flat|nested] 13+ messages in thread

end of thread, other threads:[~2024-08-16 8:18 UTC | newest]Thread overview:13+ messages (download: mbox.gz / follow: Atom feed) -- links below jump to the message on this page -- 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).