aboutsummaryrefslogtreecommitdiff
path: root/sysdeps/ieee754/ldbl-128
diff options
context:
space:
mode:
authorJoseph Myers <joseph@codesourcery.com>2014-06-29 11:49:08 +0000
committerJoseph Myers <joseph@codesourcery.com>2014-06-29 11:49:08 +0000
commitedea402804bce917cfd7cd1af76212e6364c23db (patch)
treed7688292b5b41de2f5707affccea4b992db30df2 /sysdeps/ieee754/ldbl-128
parentdd0ba018122e88937a5f14b6594b9a40693b2e58 (diff)
downloadglibc-edea402804bce917cfd7cd1af76212e6364c23db.zip
glibc-edea402804bce917cfd7cd1af76212e6364c23db.tar.gz
glibc-edea402804bce917cfd7cd1af76212e6364c23db.tar.bz2
Fix ldbl-128 powl sign of result in overflow / underflow cases (bug 17097).
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. [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.
Diffstat (limited to 'sysdeps/ieee754/ldbl-128')
-rw-r--r--sysdeps/ieee754/ldbl-128/e_powl.c26
1 files changed, 13 insertions, 13 deletions
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)