From a278d990261cf677d2a942d4f16fe6be2806668e Mon Sep 17 00:00:00 2001 From: Bruce Evans Date: Thu, 28 Feb 2008 16:22:36 +0000 Subject: [PATCH] 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. --- lib/msun/src/e_rem_pio2.c | 2 +- lib/msun/src/e_rem_pio2f.c | 6 +++--- 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 87e4c8c17062..c0e90b0e0cdb 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 2cdf78a05ab0..4439dbcf9833 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;