Fix and improve some magic numbers for the "medium size" case.

e_rem_pio2.c:
This case goes up to about 2**20pi/2, but the comment about it said that
it goes up to about 2**19pi/2.

It went too far above 2**pi/2, giving a multiplier fn with 21 significant
bits in some cases.  This would be harmful except for a numerical
accident.  It happens that the terms of the approximation to pi/2,
when rounded to 33 bits so that multiplications by 20-bit fn's are
exact, happen to be rounded to 32 bits so multiplications by 21-bit
fn's are exact too, so the bug only complicates the error analysis (we
might lose a bit of accuracy but have bits to spare).

e_rem_pio2f.c:
The bogus comment in e_rem_pio2.c was copied and the code was changed
to be bug-for-bug compatible with it, except the limit was made 90
ulps smaller than necessary.  The approximation to pi/2 was not
modified except for discarding some of it.

The same rough error analysis that justifies the limit of 2**20pi/2
for double precision only justifies a limit of 2**18pi/2 for float
precision.  We depended on exhaustive testing to check the magic numbers
for float precision.  More exaustive testing shows that we can go up
to 2**28pi/2 using a 53+25 bit approximation to pi/2 for float precision,
with a the maximum error for cosf() and sinf() unchanged at 0.5009
ulps despite the maximum error in rem_pio2f being ~0.25 ulps.  Implement
this.

2 files changed, 4 insertions, 4 deletions

diff --git a/lib/msun/src/e_rem_pio2.c b/lib/msun/src/e_rem_pio2.c
index 87e4c8c..c0e90b0 100644
--- a/lib/msun/src/e_rem_pio2.c
+++ b/lib/msun/src/e_rem_pio2.c
@@ -126,7 +126,7 @@ __ieee754_rem_pio2(double x, double *y)
 		}
 	    }
 	}
-	if(ix<=0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */
+	if(ix<0x413921fb) {	/* |x| ~< 2^20*(pi/2), medium size */
 medium:
 	    /* Use a specialized rint() to get fn.  Assume round-to-nearest. */
 	    STRICT_ASSIGN(double,fn,x*invpio2+0x1.8p52);
diff --git a/lib/msun/src/e_rem_pio2f.c b/lib/msun/src/e_rem_pio2f.c
index 2cdf78a..4439dbc 100644
--- a/lib/msun/src/e_rem_pio2f.c
+++ b/lib/msun/src/e_rem_pio2f.c
@@ -38,8 +38,8 @@ __FBSDID("$FreeBSD$");
 static const double
 half =  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
 invpio2 =  6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
-pio2_1  =  1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
-pio2_1t =  6.07710050650619224932e-11; /* 0x3DD0B461, 0x1A626331 */
+pio2_1  =  1.57079631090164184570e+00, /* 0x3FF921FB, 0x50000000 */
+pio2_1t =  1.58932547735281966916e-08; /* 0x3E5110b4, 0x611A6263 */
 
 #ifdef INLINE_REM_PIO2F
 extern inline
@@ -55,7 +55,7 @@ __ieee754_rem_pio2f(float x, double *y)
 	GET_FLOAT_WORD(hx,x);
 	ix = hx&0x7fffffff;
     /* 33+53 bit pi is good enough for medium size */
-	if(ix<=0x49490f80) {		/* |x| ~<= 2^19*(pi/2), medium size */
+	if(ix<0x4dc90fdb) {		/* |x| ~< 2^28*(pi/2), medium size */
 	    /* Use a specialized rint() to get fn.  Assume round-to-nearest. */
 	    STRICT_ASSIGN(double,fn,x*invpio2+0x1.8p52);
 	    fn = fn-0x1.8p52;