This is the mail archive of the
newlib@sources.redhat.com
mailing list for the newlib project.
Re: bug in fdlibm e_pow.c
- From: "J. Johnston" <jjohnstn at redhat dot com>
- To: Stephen L Moshier <steve at moshier dot net>
- Cc: newlib at sources dot redhat dot com
- Date: Thu, 13 Jun 2002 19:14:51 -0400
- Subject: Re: bug in fdlibm e_pow.c
- Organization: Red Hat Inc.
- References: <Pine.LNX.4.33.0206131532130.25815-100000@moshier.net>
Stephen L Moshier wrote:
>
> 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:
Patch checked in and duplicated for ef_pow.c. Thanks.
-- Jeff J.
>
> *** 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];