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.

1 files changed, 2 insertions, 2 deletions

diff --git a/lib/msun/src/e_atan2f.c b/lib/msun/src/e_atan2f.c
index d521942..af011318 100644
--- a/lib/msun/src/e_atan2f.c
+++ b/lib/msun/src/e_atan2f.c
@@ -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)