Fix exp2 spurious underflows (bug 16560) [committed]
Commit Message
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.
@@ -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
@@ -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
@@ -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
@@ -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);
@@ -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);
@@ -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