Date: Thu, 13 Jun 2002 12:30:49 -0700 (PDT) From: Stephen L Moshier <steve@moshier.net> To: freebsd-gnats-submit@FreeBSD.org Subject: misc/39236: error in e_pow.c Message-ID: <200206131930.g5DJUn0x084577@www.freebsd.org>
next in thread | raw e-mail | index | archive | help
>Number: 39236
>Category: misc
>Synopsis: error in e_pow.c
>Confidential: no
>Severity: non-critical
>Priority: medium
>Responsible: freebsd-bugs
>State: open
>Quarter:
>Keywords:
>Date-Required:
>Class: sw-bug
>Submitter-Id: current-users
>Arrival-Date: Thu Jun 13 12:40:01 PDT 2002
>Closed-Date:
>Last-Modified:
>Originator: Stephen L Moshier
>Release: any
>Organization:
>Environment:
>Description:
pow(x,y) returns 0 when x is very close to -1.0 and y is very large.
>How-To-Repeat:
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);
}
>Fix:
freebsd: src/lib/msun/src/e_pow.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];
>Release-Note:
>Audit-Trail:
>Unformatted:
To Unsubscribe: send mail to majordomo@FreeBSD.org
with "unsubscribe freebsd-bugs" in the body of the message
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?200206131930.g5DJUn0x084577>