d/freebsd-dev
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

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

Browse Source 
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.
This commit is contained in:
Bruce Evans 2008-02-28 16:22:36 +00:00
parent dbf772aca8
commit a278d99026
Notes: svn2git 2020-12-20 02:59:44 +00:00 
svn path=/head/; revision=176640
2 changed files with 4 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
lib/msun/src/e_rem_pio2.c
View File

@ -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);

6
lib/msun/src/e_rem_pio2f.c
View File

@ -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;