Fix ldbl-128ibm powl overflow handling (bug 19674) [committed]

[Joseph Myers <joseph at codesourcery dot com>]

The ldbl-128ibm implementation of powl has some problems in the case
of overflow or underflow, which are mainly visible in non-default
rounding modes.

* When overflow or underflow is detected early, the correct sign of an
  overflowing or underflowing result is not allowed for.  This is
  mostly hidden in the default rounding mode by the errno-setting
  wrappers recomputing the result (except in non-default
  error-handling modes such as -lieee), but visible in other rounding
  modes where a result that is not zero or infinity causes the
  wrappers not to do the recomputation.

* The final scaling is done before the sign is incorporated in the
  result, but should be done afterwards for correct overflowing and
  underflowing results in directed rounding modes.

This patch fixes those problems.  Tested for powerpc.  Committed.

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

	[BZ #19674]
	* sysdeps/ieee754/ldbl-128ibm/e_powl.c (__ieee754_powl): Include
	sign in overflowing and underflowing results when overflow or
	underflow is detected early.  Include sign in result before rather
	than after scaling.

diff --git a/sysdeps/ieee754/ldbl-128ibm/e_powl.c b/sysdeps/ieee754/ldbl-128ibm/e_powl.c
index 90340e8..861b44a 100644
--- a/sysdeps/ieee754/ldbl-128ibm/e_powl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/e_powl.c
@@ -148,7 +148,7 @@ long double
 __ieee754_powl (long double x, long double y)
 {
   long double z, ax, z_h, z_l, p_h, p_l;
-  long double y1, t1, t2, r, s, t, u, v, w;
+  long double y1, t1, t2, r, s, sgn, t, u, v, w;
   long double s2, s_h, s_l, t_h, t_l, ay;
   int32_t i, j, k, yisint, n;
   uint32_t ix, iy;
@@ -263,6 +263,11 @@ __ieee754_powl (long double x, long double y)
   if (((((u_int32_t) hx >> 31) - 1) | yisint) == 0)
     return (x - x) / (x - x);
 
+  /* sgn (sign of result -ve**odd) = -1 else = 1 */
+  sgn = one;
+  if (((((u_int32_t) hx >> 31) - 1) | (yisint - 1)) == 0)
+    sgn = -one;			/* (-ve)**(odd int) */
+
   /* |y| is huge.
      2^-16495 = 1/2 of smallest representable value.
      If (1 - 1/131072)^y underflows, y > 1.4986e9 */
@@ -272,15 +277,15 @@ __ieee754_powl (long double x, long double y)
       if (iy > 0x47d654b0)
 	{
 	  if (ix <= 0x3fefffff)
-	    return (hy < 0) ? huge * huge : tiny * tiny;
+	    return (hy < 0) ? sgn * huge * huge : sgn * tiny * tiny;
 	  if (ix >= 0x3ff00000)
-	    return (hy > 0) ? huge * huge : tiny * tiny;
+	    return (hy > 0) ? sgn * huge * huge : sgn * tiny * tiny;
 	}
       /* over/underflow if x is not close to one */
       if (ix < 0x3fefffff)
-	return (hy < 0) ? huge * huge : tiny * tiny;
+	return (hy < 0) ? sgn * huge * huge : sgn * tiny * tiny;
       if (ix > 0x3ff00000)
-	return (hy > 0) ? huge * huge : tiny * tiny;
+	return (hy > 0) ? sgn * huge * huge : sgn * tiny * tiny;
     }
 
   ay = y > 0 ? y : -y;
@@ -351,11 +356,6 @@ __ieee754_powl (long double x, long double y)
   t1 = ldbl_high (t1);
   t2 = z_l - (((t1 - t) - dp_h[k]) - z_h);
 
-  /* s (sign of result -ve**odd) = -1 else = 1 */
-  s = one;
-  if (((((u_int32_t) hx >> 31) - 1) | (yisint - 1)) == 0)
-    s = -one;			/* (-ve)**(odd int) */
-
   /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
   y1 = ldbl_high (y);
   p_l = (y - y1) * t1 + y * t2;
@@ -367,22 +367,22 @@ __ieee754_powl (long double x, long double y)
     {
       /* if z > 16384 */
       if (((j - 0x40d00000) | lj) != 0)
-	return s * huge * huge;	/* overflow */
+	return sgn * huge * huge;	/* overflow */
       else
 	{
 	  if (p_l + ovt > z - p_h)
-	    return s * huge * huge;	/* overflow */
+	    return sgn * huge * huge;	/* overflow */
 	}
     }
   else if ((j & 0x7fffffff) >= 0x40d01b90)	/* z <= -16495 */
     {
       /* z < -16495 */
       if (((j - 0xc0d01bc0) | lj) != 0)
-	return s * tiny * tiny;	/* underflow */
+	return sgn * tiny * tiny;	/* underflow */
       else
 	{
 	  if (p_l <= z - p_h)
-	    return s * tiny * tiny;	/* underflow */
+	    return sgn * tiny * tiny;	/* underflow */
 	}
     }
   /* compute 2**(p_h+p_l) */
@@ -408,8 +408,8 @@ __ieee754_powl (long double x, long double y)
   t1 = z - t * u / v;
   r = (z * t1) / (t1 - two) - (w + z * w);
   z = one - (r - z);
-  z = __scalbnl (z, n);
-  math_check_force_underflow_nonneg (z);
-  return s * z;
+  z = __scalbnl (sgn * z, n);
+  math_check_force_underflow (z);
+  return z;
 }
 strong_alias (__ieee754_powl, __powl_finite)

-- 
Joseph S. Myers
joseph@codesourcery.com