Fixed lots of 1 ULP errors caused by a broken approximation for pi/2.

We approximate pi with more than float precision using pi_hi+pi_lo in
the usual way (pi_hi is actually spelled pi in the source code), and
expect (float)0.5*pi_lo to give the low part of the corresponding
approximation for pi/2.  However, the high part for pi/2 (pi_o_2) is
rounded to nearest, which happens to round up, while the high part for
pi was rounded down.  Thus pi_o_2+(float)0.5*pi (in infinite precision)
was a very bad approximation for pi/2 -- the low term has the wrong
sign and increases the error drom less than half an ULP to a full ULP.

This fix rounds up instead of down for pi_hi.  Consistently rounding
down instead of up should work, and is the method used in e_acosf.c
and e_asinf.c.  The reason for the difference is that we sometimes
want to return precisely pi/2 in e_atan2f.c, so it is convenient to
have a correctly rounded (to nearest) value for pi/2 in a variable.
a_acosf.c and e_asinf.c also differ in directly approximating pi/2
instead pi; they multiply by 2.0 instead of dividing by 0.5 to convert
the approximation.

These complications are not directly visible in the double precision
versions because rounding to nearest happens to round down.
This commit is contained in:
Bruce Evans 2004-06-02 17:09:05 +00:00
parent 3c751c1b6c
commit 4f8f819975

View File

@ -25,8 +25,8 @@ tiny = 1.0e-30,
zero = 0.0,
pi_o_4 = 7.8539818525e-01, /* 0x3f490fdb */
pi_o_2 = 1.5707963705e+00, /* 0x3fc90fdb */
pi = 3.1415925026e+00, /* 0x40490fda */
pi_lo = 1.5099578832e-07; /* 0x34222168 */
pi = 3.1415927410e+00, /* 0x40490fdb */
pi_lo = -8.7422776573e-08; /* 0xb3bbbd2e */
float
__ieee754_atan2f(float y, float x)