Fix exp10 spurious underflows (bug 16560)
Commit Message
This patch fixes spurious underflows from exp10 for arguments near 0
(part of bug 16560; that bug also includes spurious underflows from
exp2, which are not fixed by this patch). The problem is underflows
in the internal computation converting the exp10 argument to arguments
for exp (with extra precision), and the fix is simply to return 1
early for arguments near enough to 0 (just as arguments with large
enough magnitude have their own overflow / underflow logic at the
start of the function).
Tested x86_64 and x86 and ulps updated accordingly; also tested for
powerpc32 and mips64 to validate the ldbl-128ibm and ldbl-128 changes.
(auto-libm-test-out diffs omitted below.)
2014-06-20 Joseph Myers <joseph@codesourcery.com>
[BZ #16560]
* sysdeps/ieee754/dbl-64/e_exp10.c (__ieee754_exp10): Return 1 for
arguments close to 0.
* sysdeps/ieee754/ldbl-128/e_exp10l.c (__ieee754_exp10l):
Likewise.
* sysdeps/ieee754/ldbl-128ibm/e_exp10l.c (__ieee754_exp10l):
Likewise.
* math/auto-libm-test-in: Add more tests of exp10.
* math/auto-libm-test-out: Regenerated.
* sysdeps/x86_64/fpu/libm-test-ulps: Update.
Comments
Ping. This patch
<https://sourceware.org/ml/libc-alpha/2014-06/msg00549.html> is pending
review.
@@ -889,6 +889,24 @@ exp10 -max
exp10 0.75
# GCC bug 59666: results on directed rounding may be incorrect.
exp10 0x1.348e45573a1dd72cp+8 xfail-rounding:ldbl-128ibm
+exp10 0x1p-10
+exp10 -0x1p-10
+exp10 0x1p-20
+exp10 -0x1p-20
+exp10 0x1p-30
+exp10 -0x1p-30
+exp10 0x1p-40
+exp10 -0x1p-40
+exp10 0x1p-50
+exp10 -0x1p-50
+exp10 0x1p-60
+exp10 -0x1p-60
+exp10 0x1p-100
+exp10 -0x1p-100
+exp10 min
+exp10 -min
+exp10 min_subnorm
+exp10 -min_subnorm
exp2 0
exp2 -0
@@ -35,6 +35,8 @@ __ieee754_exp10 (double arg)
return DBL_MIN * DBL_MIN;
else if (arg > DBL_MAX_10_EXP + 1)
return DBL_MAX * DBL_MAX;
+ else if (fabs (arg) < 0x1p-56)
+ return 1.0;
GET_LOW_WORD (lx, arg);
lx &= 0xf8000000;
@@ -35,6 +35,8 @@ __ieee754_exp10l (long double arg)
return LDBL_MIN * LDBL_MIN;
else if (arg > LDBL_MAX_10_EXP + 1)
return LDBL_MAX * LDBL_MAX;
+ else if (fabsl (arg) < 0x1p-116L)
+ return 1.0L;
u.value = arg;
u.parts64.lsw &= 0xfe00000000000000LL;
@@ -35,6 +35,8 @@ __ieee754_exp10l (long double arg)
return LDBL_MIN * LDBL_MIN;
else if (arg > LDBL_MAX_10_EXP + 1)
return LDBL_MAX * LDBL_MAX;
+ else if (fabsl (arg) < 0x1p-109L)
+ return 1.0L;
u.ld = arg;
arg_high = u.d[0].d;
@@ -1359,7 +1359,9 @@ ldouble: 1
Function: "exp10_downward":
double: 1
+float: 1
idouble: 1
+ifloat: 1
ildouble: 2
ldouble: 2
@@ -1371,7 +1373,9 @@ ldouble: 1
Function: "exp10_towardzero":
double: 1
+float: 1
idouble: 1
+ifloat: 1
ildouble: 1
ldouble: 1
@@ -1717,13 +1721,17 @@ ldouble: 1
Function: "pow10_downward":
double: 1
+float: 1
idouble: 1
+ifloat: 1
ildouble: 2
ldouble: 2
Function: "pow10_towardzero":
double: 1
+float: 1
idouble: 1
+ifloat: 1
ildouble: 1
ldouble: 1