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, 28 insertions, 4 deletions

diff --git a/sysdeps/x86/fpu/powl_helper.c b/sysdeps/x86/fpu/powl_helper.c
index c9c92e1..b4b5bd6 100644
--- a/sysdeps/x86/fpu/powl_helper.c
+++ b/sysdeps/x86/fpu/powl_helper.c
@@ -18,6 +18,7 @@
 
 #include <math.h>
 #include <math_private.h>
+#include <stdbool.h>
 
 /* High parts and low parts of -log (k/16), for integer k from 12 to
    24.  */
@@ -63,15 +64,32 @@ acc_split (long double *rhi, long double *rlo, long double hi, long double lo,
 extern long double __powl_helper (long double x, long double y);
 libm_hidden_proto (__powl_helper)
 
-/* Given X a value that is finite and nonzero, or a NaN, and only
-   negative if Y is not an integer, and Y a finite nonzero value with
-   0x1p-79 <= |Y| <= 0x1p78, compute X to the power Y.  */
+/* Given X a value that is finite and nonzero, or a NaN, and Y a
+   finite nonzero value with 0x1p-79 <= |Y| <= 0x1p78, compute X to
+   the power Y.  */
 
 long double
 __powl_helper (long double x, long double y)
 {
-  if (isnan (x) || x < 0)
+  if (isnan (x))
     return __ieee754_expl (y * __ieee754_logl (x));
+  bool negate;
+  if (x < 0)
+    {
+      long double absy = fabsl (y);
+      if (absy >= 0x1p64L)
+	negate = false;
+      else
+	{
+	  unsigned long long yll = absy;
+	  if (yll != absy)
+	    return __ieee754_expl (y * __ieee754_logl (x));
+	  negate = (yll & 1) != 0;
+	}
+      x = fabsl (x);
+    }
+  else
+    negate = false;
 
   /* We need to compute Y * log2 (X) to at least 64 bits after the
      point for normal results (that is, to at least 78 bits
@@ -199,11 +217,17 @@ __powl_helper (long double x, long double y)
      fractional parts.  */
   long double log2_res_int = __roundl (log2_res_hi);
   long double log2_res_frac = log2_res_hi - log2_res_int + log2_res_lo;
+  /* If the integer part is very large, the computed fractional part
+     may be outside the valid range for f2xm1.  */
+  if (fabsl (log2_res_int) > 16500)
+    log2_res_frac = 0;
 
   /* Compute the final result.  */
   long double res;
   asm ("f2xm1" : "=t" (res) : "0" (log2_res_frac));
   res += 1.0L;
+  if (negate)
+    res = -res;
   asm ("fscale" : "=t" (res) : "0" (res), "u" (log2_res_int));
   return res;
 }