diff options
Diffstat (limited to 'sysdeps/ieee754')
| -rw-r--r-- | sysdeps/ieee754/dbl-64/e_pow.c | 61 |
1 files changed, 41 insertions, 20 deletions
diff --git a/sysdeps/ieee754/dbl-64/e_pow.c b/sysdeps/ieee754/dbl-64/e_pow.c index 1c5f4b311e..2dd66d87b7 100644 --- a/sysdeps/ieee754/dbl-64/e_pow.c +++ b/sysdeps/ieee754/dbl-64/e_pow.c @@ -34,6 +34,7 @@ /* round to nearest mode of IEEE 754 standard. */ /* */ /***************************************************************************/ +#include <math.h> #include "endian.h" #include "upow.h" #include <dla.h> @@ -91,27 +92,33 @@ __ieee754_pow (double x, double y) { /* if y<-1 or y>1 */ double retval; - SET_RESTORE_ROUND (FE_TONEAREST); + { + SET_RESTORE_ROUND (FE_TONEAREST); - /* Avoid internal underflow for tiny y. The exact value of y does - not matter if |y| <= 2**-64. */ - if (ABS (y) < 0x1p-64) - y = y < 0 ? -0x1p-64 : 0x1p-64; - z = log1 (x, &aa, &error); /* x^y =e^(y log (X)) */ - t = y * CN; - y1 = t - (t - y); - y2 = y - y1; - t = z * CN; - a1 = t - (t - z); - a2 = (z - a1) + aa; - a = y1 * a1; - aa = y2 * a1 + y * a2; - a1 = a + aa; - a2 = (a - a1) + aa; - error = error * ABS (y); - t = __exp1 (a1, a2, 1.9e16 * error); /* return -10 or 0 if wasn't computed exactly */ - retval = (t > 0) ? t : power1 (x, y); + /* Avoid internal underflow for tiny y. The exact value of y does + not matter if |y| <= 2**-64. */ + if (ABS (y) < 0x1p-64) + y = y < 0 ? -0x1p-64 : 0x1p-64; + z = log1 (x, &aa, &error); /* x^y =e^(y log (X)) */ + t = y * CN; + y1 = t - (t - y); + y2 = y - y1; + t = z * CN; + a1 = t - (t - z); + a2 = (z - a1) + aa; + a = y1 * a1; + aa = y2 * a1 + y * a2; + a1 = a + aa; + a2 = (a - a1) + aa; + error = error * ABS (y); + t = __exp1 (a1, a2, 1.9e16 * error); /* return -10 or 0 if wasn't computed exactly */ + retval = (t > 0) ? t : power1 (x, y); + } + if (__isinf (retval)) + retval = huge * huge; + else if (retval == 0) + retval = tiny * tiny; return retval; } @@ -164,7 +171,21 @@ __ieee754_pow (double x, double y) return y < 0 ? 0.0 : INF.x; } /* if y even or odd */ - return (k == 1) ? __ieee754_pow (-x, y) : -__ieee754_pow (-x, y); + if (k == 1) + return __ieee754_pow (-x, y); + else + { + double retval; + { + SET_RESTORE_ROUND (FE_TONEAREST); + retval = -__ieee754_pow (-x, y); + } + if (__isinf (retval)) + retval = -huge * huge; + else if (retval == 0) + retval = -tiny * tiny; + return retval; + } } /* x>0 */ |