bug in fdlibm e_pow.c

[Stephen L Moshier <steve at moshier dot net>]

pow(x,y) returns 0 when x is very close to -1.0 and y is very large.
The following test program prints

pow(1.0000000000000002e+00 4.5035996273704970e+15) = 2.7182818284590455e+00
pow(-1.0000000000000002e+00 4.5035996273704970e+15) =0.0000000000000000e+00
pow(9.9999999999999978e-01 4.5035996273704970e+15) = 3.6787944117144222e-01
pow(-9.9999999999999978e-01 4.5035996273704970e+15) = 0.0000000000000000e+00

which is incorrect for the negative arguments raised to an odd integer
power.

-----
double pow (double, double);

int
main ()
{
  double x, y, z;

  x = 1.0 + pow (2.0, -52.0);
  y = 1.0 + pow (2.0, 52.0);
  z = pow (x, y);
  printf ("pow(%.16e %.16e) = %.16e\n", x, y, z);
  x = -x;
  z = pow (x, y);
  printf ("pow(%.16e %.16e) = %.16e\n", x, y, z);
  x = 1.0 - pow (2.0, -52.0);
  z = pow (x, y);
  printf ("pow(%.16e %.16e) = %.16e\n", x, y, z);
  x = -x;
  z = pow (x, y);
  printf ("pow(%.16e %.16e) = %.16e\n", x, y, z);
}
-----

Here is a patch for newlib/libm/math/epow.c:

*** e_pow.c     1999/07/13 23:49:20     1.1
--- e_pow.c     2002/06/13 17:06:57
***************
*** 224,230 ****
            if(ix>0x3ff00000) return (hy>0)? C[4]*C[4]:C[5]*C[5];
        /* now |1-x| is tiny <= 2**-20, suffice to compute
           log(x) by x-x^2/2+x^3/3-x^4/4 */
!           t = x-1;            /* t has 20 trailing zeros */
            w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25));
            u = C[25]*t;        /* ivln2_h has 21 sig. bits */
            v = t*C[26]-w*C[24];
--- 224,230 ----
            if(ix>0x3ff00000) return (hy>0)? C[4]*C[4]:C[5]*C[5];
        /* now |1-x| is tiny <= 2**-20, suffice to compute
           log(x) by x-x^2/2+x^3/3-x^4/4 */
!           t = ax-1;           /* t has 20 trailing zeros */
            w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25));
            u = C[25]*t;        /* ivln2_h has 21 sig. bits */
            v = t*C[26]-w*C[24];