diff mbox

Fix exp2 spurious underflows (bug 16560) [committed]

Message ID alpine.DEB.2.10.1502121903170.5331@digraph.polyomino.org.uk
State New
Headers show

Commit Message

Joseph Myers Feb. 12, 2015, 7:03 p.m. UTC 
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 mbox

Patch

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