Fix infinity loop with certain atan2 calls on 64-bit platforms

[Andreas Jaeger <aj at suse dot de>]

We got a report by Willus (Thanks!) that
 atan2(-.00756827042671106339, -.01792735857538728036); 
causes an infinity loop on amd64.

Analysis by Michael showed that this is a 64-bit portability problem
in mpsqrt.c.

I'm appending a patch together with a testcase.  glibc does not pass
the testsuite due to rounding errors on amd64 for me (with float
tests):

Failure: Test: atan2 (-0.00756827042671106339, -0.001792735857538728036) == -1.80338464113663849327153994380
Result:
 is:         -1.80338394641876220703e+00  -0x1.cdaa9200000000000000p+0
 should be:  -1.80338466167449951172e+00  -0x1.cdaa9e00000000000000p+0
 difference:  7.15255737304687500000e-07   0x1.80000000000000000000p-21
 ulp       :  6.0000
 max.ulp   :  0.0000
Maximal error of `atan2'
 is      : 6 ulp
 accepted: 3 ulp

So, I've changed the ulp file also (no changed needed on i686).

Ok to commit?

Andreas


2003-11-27  Andreas Jaeger  <aj@suse.de>

	* sysdeps/x86_64/fpu/libm-test-ulps: Add ulps for new atan2 test.

	* math/libm-test.inc (atan2_test): Add test that run infinitly.
	Reported by "Willus" <etc231etc231@willus.com>.

2003-11-27  Michael Matz  <matz@suse.de>

	* sysdeps/ieee754/dbl-64/mpsqrt.c (fastiroot): Fix 64-bit problem
	with wrong types.

============================================================
Index: math/libm-test.inc
--- math/libm-test.inc	24 Nov 2003 22:58:44 -0000	1.53
+++ math/libm-test.inc	27 Nov 2003 18:37:18 -0000
@@ -943,6 +943,8 @@ atan2_test (void)
   TEST_ff_f (atan2, 0.390625L, .00029L, 1.57005392693128974780151246612928941L);
   TEST_ff_f (atan2, 1.390625L, 0.9296875L, 0.981498387184244311516296577615519772L);
 
+  TEST_ff_f (atan2, -0.00756827042671106339L, -.001792735857538728036L, -1.80338464113663849327153994380L);
+  
   END (atan2);
 }
 
============================================================
Index: sysdeps/ieee754/dbl-64/mpsqrt.c
--- sysdeps/ieee754/dbl-64/mpsqrt.c	26 Aug 2002 22:40:37 -0000	1.8
+++ sysdeps/ieee754/dbl-64/mpsqrt.c	27 Nov 2003 18:37:18 -0000
@@ -83,9 +83,9 @@ void __mpsqrt(mp_no *x, mp_no *y, int p)
 /* with the relative error bounded by 2**-51.              */
 /***********************************************************/
 double fastiroot(double x) {
-  union {long i[2]; double d;} p,q;
+  union {int i[2]; double d;} p,q;
   double y,z, t;
-  long n;
+  int n;
   static const double c0 = 0.99674, c1 = -0.53380, c2 = 0.45472, c3 = -0.21553;
 
   p.d = x;
============================================================
Index: sysdeps/x86_64/fpu/libm-test-ulps
--- sysdeps/x86_64/fpu/libm-test-ulps	23 Mar 2003 00:52:10 -0000	1.7
+++ sysdeps/x86_64/fpu/libm-test-ulps	27 Nov 2003 18:37:18 -0000
@@ -27,6 +27,9 @@ ifloat: 3
 Test "atan2 (1.390625, 0.9296875) == 0.981498387184244311516296577615519772":
 float: 1
 ifloat: 1
+Test "atan2 (-0.00756827042671106339, -0.001792735857538728036) == -1.80338464113663849327153994380":
+float: 6
+ifloat: 6
 
 # atanh
 Test "atanh (0.75) == 0.972955074527656652552676371721589865":
@@ -921,8 +924,8 @@ ildouble: 1
 ldouble: 1
 
 Function: "atan2":
-float: 3
-ifloat: 3
+float: 6
+ifloat: 6
 
 Function: "atanh":
 float: 1

-- 
 Andreas Jaeger, aj@suse.de, http://www.suse.de/~aj
  SuSE Linux AG, Deutschherrnstr. 15-19, 90429 Nürnberg, Germany
   GPG fingerprint = 93A3 365E CE47 B889 DF7F  FED1 389A 563C C272 A126

Attachment: pgp00000.pgp
Description: PGP signature