From 4da6db51880289f0bf41b39e05cf9bb1c4769c47 Mon Sep 17 00:00:00 2001 From: Joseph Myers Date: Mon, 23 Jun 2014 20:12:33 +0000 Subject: Fix pow overflow in non-default rounding modes (bug 16315). This patch fixes bug 16315, bad pow handling of overflow/underflow in non-default rounding modes. Tests of pow are duly converted to ALL_RM_TEST to run all tests in all rounding modes. There are two main issues here. First, various implementations compute a negative result by negating a positive result, but this yields inappropriate overflow / underflow values for directed rounding, so either overflow / underflow results need recomputing in the correct sign, or the relevant overflowing / underflowing operation needs to be made to have a result of the correct sign. Second, the dbl-64 implementation sets FE_TONEAREST internally; in the overflow / underflow case, the result needs recomputing in the original rounding mode. Tested x86_64 and x86 and ulps updated accordingly. [BZ #16315] * sysdeps/i386/fpu/e_pow.S (__ieee754_pow): Ensure possibly overflowing or underflowing operations take place with sign of result. * sysdeps/i386/fpu/e_powf.S (__ieee754_powf): Likewise. * sysdeps/i386/fpu/e_powl.S (__ieee754_powl): Likewise. * sysdeps/ieee754/dbl-64/e_pow.c: Include . (__ieee754_pow): Recompute overflowing and underflowing results in original rounding mode. * sysdeps/x86/fpu/powl_helper.c: Include . (__powl_helper): Allow negative argument X and scale negated value as needed. Avoid passing value outside [-1, 1] to f2xm1. * sysdeps/x86_64/fpu/e_powl.S (__ieee754_powl): Ensure possibly overflowing or underflowing operations take place with sign of result. * sysdeps/x86_64/fpu/multiarch/e_pow.c [HAVE_FMA4_SUPPORT]: Include . * math/auto-libm-test-in: Add more tests of pow. * math/auto-libm-test-out: Regenerated. * math/libm-test.inc (pow_test): Use ALL_RM_TEST. (pow_tonearest_test_data): Remove. (pow_test_tonearest): Likewise. (pow_towardzero_test_data): Likewise. (pow_test_towardzero): Likewise. (pow_downward_test_data): Likewise. (pow_test_downward): Likewise. (pow_upward_test_data): Likewise. (pow_test_upward): Likewise. (main): Don't call removed functions. * sysdeps/i386/fpu/libm-test-ulps: Update. * sysdeps/x86_64/fpu/libm-test-ulps: Likewise. --- sysdeps/ieee754/dbl-64/e_pow.c | 61 ++++++++++++++++++++++++++++-------------- 1 file changed, 41 insertions(+), 20 deletions(-) (limited to 'sysdeps/ieee754/dbl-64/e_pow.c') 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 #include "endian.h" #include "upow.h" #include @@ -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 */ -- cgit 1.4.1