diff mbox

Fix ldbl-128 powl sign of result in overflow / underflow cases (bug 17097)

Message ID Pine.LNX.4.64.1406272150530.8357@digraph.polyomino.org.uk
State Committed
Headers show

Commit Message

Joseph Myers June 27, 2014, 9:52 p.m. UTC
This patch fixes bug 17097, ldbl-128 powl producing overflowing /
underflowing results with positive sign when the result should have
been negative.  This was shown up by the tests in non-default rounding
modes added by my patch for bug 16315, but isn't actually limited to
non-default rounding modes: rather, when rounding to nearest the
wrappers produced a result with the correct sign and so always hid the
bug unless -lieee was used to disable the wrappers.  The problem is
that in the cases where Y is large enough that the result overflows or
underflows for X not very close to 1, but not large enough to overflow
or underflow for all X != +/- 1 (in the latter case Y is always an
even integer), a positive overflowing / underflowing result is always
returned, rather than one with the correct sign.  This patch moves the
relevant part of computation of the sign earlier and returns a result
of the correct sign.

Tested for mips64.

2014-06-27  Joseph Myers  <joseph@codesourcery.com>

	[BZ #17097]
	* sysdeps/ieee754/ldbl-128/e_powl.c (__ieee754_powl): Return
	result with correct sign in case of exponents that produce
	overflow except for X very close to 1.

Comments

Andreas Jaeger June 29, 2014, 7:59 a.m. UTC | #1
On 06/27/2014 11:52 PM, Joseph S. Myers wrote:
> This patch fixes bug 17097, ldbl-128 powl producing overflowing /
> underflowing results with positive sign when the result should have
> been negative.  This was shown up by the tests in non-default rounding
> modes added by my patch for bug 16315, but isn't actually limited to
> non-default rounding modes: rather, when rounding to nearest the
> wrappers produced a result with the correct sign and so always hid the
> bug unless -lieee was used to disable the wrappers.  The problem is
> that in the cases where Y is large enough that the result overflows or
> underflows for X not very close to 1, but not large enough to overflow
> or underflow for all X != ± 1 (in the latter case Y is always an
> even integer), a positive overflowing / underflowing result is always
> returned, rather than one with the correct sign.  This patch moves the
> relevant part of computation of the sign earlier and returns a result
> of the correct sign.

Thanks,
Andreas
diff mbox

Patch

diff --git a/sysdeps/ieee754/ldbl-128/e_powl.c b/sysdeps/ieee754/ldbl-128/e_powl.c
index d131750..f531385 100644
--- a/sysdeps/ieee754/ldbl-128/e_powl.c
+++ b/sysdeps/ieee754/ldbl-128/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;
   u_int32_t ix, iy;
@@ -262,6 +262,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 */
@@ -277,9 +282,9 @@  __ieee754_powl (long double x, long double y)
 	}
       /* over/underflow if x is not close to one */
       if (ix < 0x3ffeffff)
-	return (hy < 0) ? huge * huge : tiny * tiny;
+	return (hy < 0) ? sgn * huge * huge : sgn * tiny * tiny;
       if (ix > 0x3fff0000)
-	return (hy > 0) ? huge * huge : tiny * tiny;
+	return (hy > 0) ? sgn * huge * huge : sgn * tiny * tiny;
     }
 
   ay = y > 0 ? y : -y;
@@ -366,11 +371,6 @@  __ieee754_powl (long double x, long double y)
   t1 = o.value;
   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 = y;
   o.value = y1;
@@ -386,11 +386,11 @@  __ieee754_powl (long double x, long double y)
     {
       /* if z > 16384 */
       if (((j - 0x400d0000) | o.parts32.w1 | o.parts32.w2 | o.parts32.w3) != 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) >= 0x400d01b9)	/* z <= -16495 */
@@ -398,11 +398,11 @@  __ieee754_powl (long double x, long double y)
       /* z < -16495 */
       if (((j - 0xc00d01bc) | o.parts32.w1 | o.parts32.w2 | o.parts32.w3)
 	  != 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) */
@@ -441,6 +441,6 @@  __ieee754_powl (long double x, long double y)
       o.parts32.w0 = j;
       z = o.value;
     }
-  return s * z;
+  return sgn * z;
 }
 strong_alias (__ieee754_powl, __powl_finite)