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 <math.h>.
	(__ieee754_pow): Recompute overflowing and underflowing results in
	original rounding mode.
	* sysdeps/x86/fpu/powl_helper.c: Include <stdbool.h>.
	(__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.h>.
	* 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.

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 */