http://sources.redhat.com/ml/newlib/2002/msg00230.html 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.

--- a/src/video/e_pow.h	Mon Sep 08 07:38:41 2008 +0000
+++ b/src/video/e_pow.h	Mon Sep 15 04:31:30 2008 +0000
@@ -208,7 +208,7 @@
             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 */
+        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;