This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Fix exp2, exp2f spurious underflows (bug 18219) [committed]
- From: Joseph Myers <joseph at codesourcery dot com>
- To: <libc-alpha at sourceware dot org>
- Date: Tue, 23 Jun 2015 14:36:03 +0000
- Subject: Fix exp2, exp2f spurious underflows (bug 18219) [committed]
- Authentication-results: sourceware.org; auth=none
The dbl-64 and flt-32 implementations of exp2 functions produce
spurious underflow exceptions. The underlying reason is the same in
both cases: the computation works as (2^a - 1)*2^b + 2^b for suitably
chosen a and b, where a has small magnitude so 2^a - 1 can be computed
with a low-degree polynomial approximation, and (2^a - 1)*2^b can
underflow even when the final result does not. This patch fixes this
by adjusting the threshold for when scaling is used to avoid
intermediate underflow so it works for any possible value of a where
the final result would not underflow.
Tested for x86_64 and x86. Committed.
(auto-libm-test-out diffs omitted below.)
2015-06-23 Joseph Myers <joseph@codesourcery.com>
[BZ #18219]
* sysdeps/ieee754/dbl-64/e_exp2.c (__ieee754_exp2): Reduce
threshold on absolute value of exponent for which scaling is used.
* sysdeps/ieee754/flt-32/e_exp2f.c (__ieee754_exp2f): Likewise.
* math/auto-libm-test-in: Add more tests of exp2.
* math/auto-libm-test-out: Regenerated.
diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index 9bae9d5..cf7fc35 100644
--- a/math/auto-libm-test-in
+++ b/math/auto-libm-test-in
@@ -1271,9 +1271,16 @@ exp2 max
exp2 -max
exp2 0.75
exp2 100.5
+exp2 -116.5
+exp2 -123.5
+exp2 -124.5
+exp2 -125.5
exp2 127
exp2 -149
exp2 1000.25
+exp2 -1019.5
+exp2 -1020.5
+exp2 -1021.5
exp2 1023
exp2 -1074
exp2 16383
diff --git a/sysdeps/ieee754/dbl-64/e_exp2.c b/sysdeps/ieee754/dbl-64/e_exp2.c
index 948a756..30f0a8f 100644
--- a/sysdeps/ieee754/dbl-64/e_exp2.c
+++ b/sysdeps/ieee754/dbl-64/e_exp2.c
@@ -96,7 +96,9 @@ __ieee754_exp2 (double x)
/* 3. Compute ex2 = 2^(t/512+e+ex). */
ex2_u.d = exp2_accuratetable[tval & 511];
tval >>= 9;
- unsafe = abs (tval) >= -DBL_MIN_EXP - 1;
+ /* x2 is an integer multiple of 2^-54; avoid intermediate
+ underflow from the calculation of x22 * x. */
+ unsafe = abs (tval) >= -DBL_MIN_EXP - 56;
ex2_u.ieee.exponent += tval >> unsafe;
scale_u.d = 1.0;
scale_u.ieee.exponent += tval - (tval >> unsafe);
diff --git a/sysdeps/ieee754/flt-32/e_exp2f.c b/sysdeps/ieee754/flt-32/e_exp2f.c
index f3e3a8e..0b75a7e 100644
--- a/sysdeps/ieee754/flt-32/e_exp2f.c
+++ b/sysdeps/ieee754/flt-32/e_exp2f.c
@@ -89,7 +89,9 @@ __ieee754_exp2f (float x)
/* 3. Compute ex2 = 2^(t/255+e+ex). */
ex2_u.f = __exp2f_atable[tval & 255];
tval >>= 8;
- unsafe = abs(tval) >= -FLT_MIN_EXP - 1;
+ /* x2 is an integer multiple of 2^-30; avoid intermediate
+ underflow from the calculation of x22 * x. */
+ unsafe = abs(tval) >= -FLT_MIN_EXP - 32;
ex2_u.ieee.exponent += tval >> unsafe;
scale_u.f = 1.0;
scale_u.ieee.exponent += tval - (tval >> unsafe);
--
Joseph S. Myers
joseph@codesourcery.com