From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: X-Spam-Checker-Version: SpamAssassin 3.4.4 (2020-01-24) on inbox.vuxu.org X-Spam-Level: X-Spam-Status: No, score=-3.3 required=5.0 tests=HEADER_FROM_DIFFERENT_DOMAINS, MAILING_LIST_MULTI,RCVD_IN_DNSWL_MED,RCVD_IN_MSPIKE_H4, RCVD_IN_MSPIKE_WL autolearn=ham autolearn_force=no version=3.4.4 Received: from second.openwall.net (second.openwall.net [193.110.157.125]) by inbox.vuxu.org (Postfix) with SMTP id E359722F9C for ; Thu, 15 Aug 2024 16:28:39 +0200 (CEST) Received: (qmail 3657 invoked by uid 550); 15 Aug 2024 14:28:35 -0000 Mailing-List: contact musl-help@lists.openwall.com; run by ezmlm Precedence: bulk List-Post: List-Help: List-Unsubscribe: List-Subscribe: List-ID: Reply-To: musl@lists.openwall.com Received: (qmail 3616 invoked from network); 15 Aug 2024 14:28:34 -0000 Date: Thu, 15 Aug 2024 10:28:25 -0400 From: Rich Felker To: Damian McGuckin Cc: MUSL Message-ID: <20240815142825.GP10433@brightrain.aerifal.cx> References: <21ca5c9-b1e-71b5-87b-a37f81f691ab@esi.com.au> <20240811200812.GZ3766212@port70.net> <7e9463-98a3-4b2-c10-e3fbf79a6b8@esi.com.au> <5918d2a7-7b3-932b-2b4-b24390832244@esi.com.au> <20240812090022.GA3766212@port70.net> <20240815011655.GN10433@brightrain.aerifal.cx> MIME-Version: 1.0 Content-Type: multipart/mixed; boundary="XStn23h1fwudRqtG" Content-Disposition: inline In-Reply-To: User-Agent: Mutt/1.5.21 (2010-09-15) Subject: Re: [musl] catan(z) --XStn23h1fwudRqtG Content-Type: text/plain; charset=us-ascii Content-Disposition: inline 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 --XStn23h1fwudRqtG Content-Type: text/plain; charset=us-ascii Content-Disposition: attachment; filename="0001-catan-remove-no-up-reduction-mod-pi-and-unused-code.patch" >From 392529fb945841d96981235d16f8c1b926f134e4 Mon Sep 17 00:00:00 2001 From: Rich Felker 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 --XStn23h1fwudRqtG--