This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Fix exp2 spurious underflows (bug 16560) [committed]
- From: Joseph Myers <joseph at codesourcery dot com>
- To: <libc-alpha at sourceware dot org>
- Date: Thu, 12 Feb 2015 19:03:38 +0000
- Subject: Fix exp2 spurious underflows (bug 16560) [committed]
- Authentication-results: sourceware.org; auth=none
This patch fixes the remaining part of bug 16560, spurious underflows
from exp2 of arguments close to 0 (when the result is close to 1, so
should not underflow), by just using 1+x instead of a more complicated
calculation when the argument is sufficiently small.
Tested for x86_64, x86 and mips64. Committed.
(auto-libm-test-out diffs omitted below.)
2015-02-12 Joseph Myers <joseph@codesourcery.com>
[BZ #16560]
* math/e_exp2l.c [LDBL_MANT_DIG == 106] (LDBL_EPSILON): Undefine
and redefine.
(__ieee754_exp2l): Do not multiply small fractional parts by
M_LN2l.
* sysdeps/i386/fpu/e_exp2l.S (__ieee754_exp2l): Just add 1 to
small argument.
* sysdeps/ieee754/dbl-64/e_exp2.c (__ieee754_exp2): Likewise.
* sysdeps/ieee754/flt-32/e_exp2f.c (__ieee754_exp2f): Likewise.
* sysdeps/x86_64/fpu/e_exp2l.S (__ieee754_exp2l): 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 ad43113..6bcfd54 100644
--- a/math/auto-libm-test-in
+++ b/math/auto-libm-test-in
@@ -938,6 +938,24 @@ exp2 1023
exp2 -1074
exp2 16383
exp2 -16400
+exp2 0x1p-10
+exp2 -0x1p-10
+exp2 0x1p-20
+exp2 -0x1p-20
+exp2 0x1p-30
+exp2 -0x1p-30
+exp2 0x1p-40
+exp2 -0x1p-40
+exp2 0x1p-50
+exp2 -0x1p-50
+exp2 0x1p-60
+exp2 -0x1p-60
+exp2 0x1p-100
+exp2 -0x1p-100
+exp2 min
+exp2 -min
+exp2 min_subnorm
+exp2 -min_subnorm
expm1 0
expm1 -0
diff --git a/math/e_exp2l.c b/math/e_exp2l.c
index bb7feef..9eb7baf 100644
--- a/math/e_exp2l.c
+++ b/math/e_exp2l.c
@@ -20,6 +20,13 @@
#include <math_private.h>
#include <float.h>
+/* To avoid spurious underflows, use this definition to treat IBM long
+ double as approximating an IEEE-style format. */
+#if LDBL_MANT_DIG == 106
+# undef LDBL_EPSILON
+# define LDBL_EPSILON 0x1p-106L
+#endif
+
long double
__ieee754_exp2l (long double x)
{
@@ -31,6 +38,8 @@ __ieee754_exp2l (long double x)
{
int intx = (int) x;
long double fractx = x - intx;
+ if (fabsl (fractx) < LDBL_EPSILON / 4.0L)
+ return __scalbnl (1.0L + fractx, intx);
return __scalbnl (__ieee754_expl (M_LN2l * fractx), intx);
}
else
diff --git a/sysdeps/i386/fpu/e_exp2l.S b/sysdeps/i386/fpu/e_exp2l.S
index 203dd00..2bf9a25 100644
--- a/sysdeps/i386/fpu/e_exp2l.S
+++ b/sysdeps/i386/fpu/e_exp2l.S
@@ -18,7 +18,15 @@ ENTRY(__ieee754_exp2l)
andb %ah, %dh
cmpb $0x05, %dh
je 1f /* Is +-Inf, jump. */
- fld %st
+ movzwl 4+8(%esp), %eax
+ andl $0x7fff, %eax
+ cmpl $0x3fbe, %eax
+ jge 3f
+ /* Argument's exponent below -65, result rounds to 1. */
+ fld1
+ faddp
+ ret
+3: fld %st
frndint /* int(x) */
fsubr %st,%st(1) /* fract(x) */
fxch
diff --git a/sysdeps/ieee754/dbl-64/e_exp2.c b/sysdeps/ieee754/dbl-64/e_exp2.c
index 3666c6a..f964a5a 100644
--- a/sysdeps/ieee754/dbl-64/e_exp2.c
+++ b/sysdeps/ieee754/dbl-64/e_exp2.c
@@ -61,6 +61,9 @@ __ieee754_exp2 (double x)
double rx, x22, result;
union ieee754_double ex2_u, scale_u;
+ if (fabs (x) < DBL_EPSILON / 4.0)
+ return 1.0 + x;
+
{
SET_RESTORE_ROUND_NOEX (FE_TONEAREST);
diff --git a/sysdeps/ieee754/flt-32/e_exp2f.c b/sysdeps/ieee754/flt-32/e_exp2f.c
index 01cd444..f3e3a8e 100644
--- a/sysdeps/ieee754/flt-32/e_exp2f.c
+++ b/sysdeps/ieee754/flt-32/e_exp2f.c
@@ -54,6 +54,9 @@ __ieee754_exp2f (float x)
float rx, x22, result;
union ieee754_float ex2_u, scale_u;
+ if (fabsf (x) < FLT_EPSILON / 4.0f)
+ return 1.0f + x;
+
{
SET_RESTORE_ROUND_NOEXF (FE_TONEAREST);
diff --git a/sysdeps/x86_64/fpu/e_exp2l.S b/sysdeps/x86_64/fpu/e_exp2l.S
index 7abf425..7d42a93 100644
--- a/sysdeps/x86_64/fpu/e_exp2l.S
+++ b/sysdeps/x86_64/fpu/e_exp2l.S
@@ -19,7 +19,15 @@ ENTRY(__ieee754_exp2l)
andb %ah, %dh
cmpb $0x05, %dh
je 1f /* Is +-Inf, jump. */
- fld %st
+ movzwl 8+8(%rsp), %eax
+ andl $0x7fff, %eax
+ cmpl $0x3fbe, %eax
+ jge 3f
+ /* Argument's exponent below -65, result rounds to 1. */
+ fld1
+ faddp
+ ret
+3: fld %st
frndint /* int(x) */
fsubr %st,%st(1) /* fract(x) */
fxch
--
Joseph S. Myers
joseph@codesourcery.com