Fix sysdeps/ieee754 pow handling of sNaN arguments (bug 20916) [committed]

[Joseph Myers <joseph at codesourcery dot com>]

Various pow function implementations mishandle sNaN arguments in
various ways.  This includes returning sNaN instead of qNaN for sNaN
arguments.  For arguments (1, sNaN) and (sNaN, 0), TS 18661-1
semantics are also that the result should be qNaN, whereas with a qNaN
argument there the result should be 1, but for the dbl-64
implementation of pow there are issues with sNaN arguments beyond not
implementing the TS 18661-1 semantics in those special cases.

This patch makes the implementations in sysdeps/ieee754 follow the TS
18661-1 semantics consistently.  Because x86 / x86_64 implementations
still need fixing, testcases are not included with this patch; they
will be included with the fix for the x86 / x86_64 versions.

Tested for x86_64, x86, mips64 and powerpc (with such testcases, which
pass in the mips64 and powerpc cases).  Committed.

2016-12-02  Joseph Myers  <joseph@codesourcery.com>

	[BZ #20916]
	* sysdeps/ieee754/dbl-64/e_pow.c (__ieee754_pow): Do not return 1
	for arguments (sNaN, 0) or (1, sNaN).  Do arithmetic on NaN
	arguments to compute result.
	* sysdeps/ieee754/flt-32/e_powf.c (__ieee754_powf): Do not return
	1 for arguments (sNaN, 0) or (1, sNaN).
	* sysdeps/ieee754/ldbl-128/e_powl.c (__ieee754_powl): Likewise.
	* sysdeps/ieee754/ldbl-128ibm/e_powl.c (__ieee754_powl): Likewise.

diff --git a/sysdeps/ieee754/dbl-64/e_pow.c b/sysdeps/ieee754/dbl-64/e_pow.c
index db6ecf7..8f9b1c0 100644
--- a/sysdeps/ieee754/dbl-64/e_pow.c
+++ b/sysdeps/ieee754/dbl-64/e_pow.c
@@ -74,8 +74,8 @@ __ieee754_pow (double x, double y)
       qx = u.i[HIGH_HALF] & 0x7fffffff;
       /* Is x a NaN?  */
       if ((((qx == 0x7ff00000) && (u.i[LOW_HALF] != 0)) || (qx > 0x7ff00000))
-	  && y != 0)
-	return x;
+	  && (y != 0 || issignaling (x)))
+	return x + x;
       if (y == 1.0)
 	return x;
       if (y == 2.0)
@@ -129,7 +129,7 @@ __ieee754_pow (double x, double y)
     {
       if (((v.i[HIGH_HALF] & 0x7fffffff) == 0x7ff00000 && v.i[LOW_HALF] != 0)
 	  || (v.i[HIGH_HALF] & 0x7fffffff) > 0x7ff00000)	/* NaN */
-	return y;
+	return y + y;
       if (fabs (y) > 1.0e20)
 	return (y > 0) ? 0 : 1.0 / 0.0;
       k = checkint (y);
@@ -143,9 +143,9 @@ __ieee754_pow (double x, double y)
   qy = v.i[HIGH_HALF] & 0x7fffffff;	/*   no sign   */
 
   if (qx >= 0x7ff00000 && (qx > 0x7ff00000 || u.i[LOW_HALF] != 0))	/* NaN */
-    return x;
+    return x + y;
   if (qy >= 0x7ff00000 && (qy > 0x7ff00000 || v.i[LOW_HALF] != 0))	/* NaN */
-    return x == 1.0 ? 1.0 : y;
+    return x == 1.0 && !issignaling (y) ? 1.0 : y + y;
 
   /* if x<0 */
   if (u.i[HIGH_HALF] < 0)
diff --git a/sysdeps/ieee754/flt-32/e_powf.c b/sysdeps/ieee754/flt-32/e_powf.c
index c72fe37..d9470f1 100644
--- a/sysdeps/ieee754/flt-32/e_powf.c
+++ b/sysdeps/ieee754/flt-32/e_powf.c
@@ -62,10 +62,10 @@ __ieee754_powf(float x, float y)
 	ix = hx&0x7fffffff;  iy = hy&0x7fffffff;
 
     /* y==zero: x**0 = 1 */
-	if(iy==0) return one;
+	if(iy==0 && !issignaling (x)) return one;
 
     /* x==+-1 */
-	if(x == 1.0) return one;
+	if(x == 1.0 && !issignaling (y)) return one;
 	if(x == -1.0 && isinf(y)) return one;
 
     /* +-NaN return x+y */
diff --git a/sysdeps/ieee754/ldbl-128/e_powl.c b/sysdeps/ieee754/ldbl-128/e_powl.c
index e6cd975..a344840 100644
--- a/sysdeps/ieee754/ldbl-128/e_powl.c
+++ b/sysdeps/ieee754/ldbl-128/e_powl.c
@@ -165,11 +165,12 @@ __ieee754_powl (_Float128 x, _Float128 y)
 
 
   /* y==zero: x**0 = 1 */
-  if ((iy | q.parts32.w1 | q.parts32.w2 | q.parts32.w3) == 0)
+  if ((iy | q.parts32.w1 | q.parts32.w2 | q.parts32.w3) == 0
+      && !issignaling (x))
     return one;
 
   /* 1.0**y = 1; -1.0**+-Inf = 1 */
-  if (x == one)
+  if (x == one && !issignaling (y))
     return one;
   if (x == -1 && iy == 0x7fff0000
       && (q.parts32.w1 | q.parts32.w2 | q.parts32.w3) == 0)
diff --git a/sysdeps/ieee754/ldbl-128ibm/e_powl.c b/sysdeps/ieee754/ldbl-128ibm/e_powl.c
index 861b44a..d6fbef6 100644
--- a/sysdeps/ieee754/ldbl-128ibm/e_powl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/e_powl.c
@@ -165,11 +165,11 @@ __ieee754_powl (long double x, long double y)
   iy = hy & 0x7fffffff;
 
   /* y==zero: x**0 = 1 */
-  if ((iy | ly) == 0)
+  if ((iy | ly) == 0 && !issignaling (x))
     return one;
 
   /* 1.0**y = 1; -1.0**+-Inf = 1 */
-  if (x == one)
+  if (x == one && !issignaling (y))
     return one;
   if (x == -1.0L && ((iy - 0x7ff00000) | ly) == 0)
     return one;

-- 
Joseph S. Myers
joseph@codesourcery.com