Re: [musl] catan(z)

[Rich Felker <dalias@libc.org>]

[-- 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