This is the mail archive of the
java-patches@gcc.gnu.org
mailing list for the Java project.
e_pow.c bug
- From: Steve Moshier <steve at moshier dot net>
- To: java-patches at gcc dot gnu dot org
- Date: Sun, 11 Apr 2004 09:44:21 -0400 (EDT)
- Subject: e_pow.c bug
Math libraries copied from Sun Microsystem's fdlibm (www.netlib.org/fdlibm)
have a bug in the real power function, e_pow.c. The bug is that
pow(x,y) returns 0 when x is very close to -1.0 and y is very large.
The following test program (in C) 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 gcc/libjava/java/lang/e_pow.c:
*** e_pow.c 2003/07/08 21:27:37 1.1
--- e_pow.c 2004/04/11 13:30:15
*************** ivln2_l = 1.92596299112661746887e-08;
*** 196,202 ****
if(ix>0x3ff00000) return (hy>0)? huge*huge:tiny*tiny;
/* 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 = ivln2_h*t; /* ivln2_h has 21 sig. bits */
v = t*ivln2_l-w*ivln2;
--- 196,202 ----
if(ix>0x3ff00000) return (hy>0)? huge*huge:tiny*tiny;
/* 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 = ivln2_h*t; /* ivln2_h has 21 sig. bits */
v = t*ivln2_l-w*ivln2;