Fix expm1 spurious underflow exceptions (bug 6778).

9 files changed, 140 insertions, 45 deletions

diff --git a/ChangeLog b/ChangeLog
index b006d3f48a..2ba450800c 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,24 @@
+2012-07-06  Joseph Myers  <joseph@codesourcery.com>
+
+	[BZ #6778]
+	* sysdeps/i386/fpu/s_expm1.S (__expm1): Check for large negative
+	inputs and return -1 for them.  Do not check for +Inf in case not
+	reachable for +Inf.
+	* sysdeps/i386/fpu/s_expm1f.S (__expm1f): Likewise.
+	* sysdeps/i386/fpu/e_expl.S [USE_AS_EXPM1L] (csat): Do not define.
+	(IEEE754_EXPL) [USE_AS_EXPM1L]: Check for large negative inputs
+	and return -1 for them.  Do not check for +Inf in case not
+	reachable for +Inf.
+	* sysdeps/x86_64/fpu/e_expl.S [USE_AS_EXPM1L] (csat): Do not
+	define.
+	(IEEE754_EXPL) [USE_AS_EXPM1L]: Check for large negative inputs
+	and return -1 for them.  Do not check for +Inf in case not
+	reachable for +Inf.
+	* math/libm-test.inc (expm1_test): Add more tests.  Do not allow
+	spurious underflow.
+	* sysdeps/i386/fpu/libm-test-ulps: Update.
+	* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
+
 2012-07-06  Mike Frysinger  <vapier@gentoo.org>
 
 	* sunrpc/rpc_clntout.c: Change <rpc/types.h> to "rpc/types.h".
diff --git a/NEWS b/NEWS
index 06df26fef6..e07105748b 100644
--- a/NEWS
+++ b/NEWS
@@ -9,7 +9,7 @@ Version 2.17
 
 * The following bugs are resolved with this release:
 
-  14157, 14283, 14328, 14331
+  6778, 14157, 14283, 14328, 14331
 
 
 Version 2.16
diff --git a/math/libm-test.inc b/math/libm-test.inc
index 6adbb61dbc..b87a40d8d9 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -4039,12 +4039,32 @@ expm1_test (void)
   TEST_f_f (expm1, 11356.25L, 9.05128237311923300051376115753226014206e+4931L);
 #endif
 
+  TEST_f_f (expm1, -10.0, -0.9999546000702375151484644084844394493898L);
+  TEST_f_f (expm1, -16.0, -0.9999998874648252807408854862248209398728L);
+  TEST_f_f (expm1, -17.0, -0.9999999586006228121483334034897228104472L);
+  TEST_f_f (expm1, -18.0, -0.9999999847700202552873715638633707664826L);
+  TEST_f_f (expm1, -36.0, -0.9999999999999997680477169756430611687736L);
+  TEST_f_f (expm1, -37.0, -0.9999999999999999146695237425593420572195L);
+  TEST_f_f (expm1, -38.0, -0.9999999999999999686086720795197037129104L);
+  TEST_f_f (expm1, -44.0, -0.9999999999999999999221886775886620348429L);
+  TEST_f_f (expm1, -45.0, -0.9999999999999999999713748141945060635553L);
+  TEST_f_f (expm1, -46.0, -0.9999999999999999999894693826424461876212L);
+  TEST_f_f (expm1, -73.0, -0.9999999999999999999999999999999802074012L);
+  TEST_f_f (expm1, -74.0, -0.9999999999999999999999999999999927187098L);
+  TEST_f_f (expm1, -75.0, -0.9999999999999999999999999999999973213630L);
+  TEST_f_f (expm1, -78.0, -0.9999999999999999999999999999999998666385L);
+  TEST_f_f (expm1, -79.0, -0.9999999999999999999999999999999999509391L);
+  TEST_f_f (expm1, -80.0, -0.9999999999999999999999999999999999819515L);
+  TEST_f_f (expm1, -100.0, -1.0);
+  TEST_f_f (expm1, -1000.0, -1.0);
+  TEST_f_f (expm1, -10000.0, -1.0);
+  TEST_f_f (expm1, -100000.0, -1.0);
+
   errno = 0;
   TEST_f_f (expm1, 100000.0, plus_infty, OVERFLOW_EXCEPTION);
   check_int ("errno for expm1(large) == ERANGE", errno, ERANGE, 0, 0, 0);
   TEST_f_f (expm1, max_value, plus_infty, OVERFLOW_EXCEPTION);
-  /* Bug 6778: spurious underflow exception.  */
-  TEST_f_f (expm1, -max_value, -1, UNDERFLOW_EXCEPTION_OK);
+  TEST_f_f (expm1, -max_value, -1);
 
   END (expm1);
 }
diff --git a/sysdeps/i386/fpu/e_expl.S b/sysdeps/i386/fpu/e_expl.S
index bab0a081b8..e42c9a1213 100644
--- a/sysdeps/i386/fpu/e_expl.S
+++ b/sysdeps/i386/fpu/e_expl.S
@@ -60,10 +60,12 @@ c1:	.byte 0x20, 0xfa, 0xee, 0xc2, 0x5f, 0x70, 0xa5, 0xec, 0xed, 0x3f
 	.byte 0, 0, 0, 0, 0, 0
 	ASM_SIZE_DIRECTIVE(c1)
 #endif
+#ifndef USE_AS_EXPM1L
 	ASM_TYPE_DIRECTIVE(csat,@object)
 csat:	.byte 0, 0, 0, 0, 0, 0, 0, 0x80, 0x0e, 0x40
 	.byte 0, 0, 0, 0, 0, 0
 	ASM_SIZE_DIRECTIVE(csat)
+#endif
 
 #ifdef PIC
 # define MO(op) op##@GOTOFF(%ecx)
@@ -88,9 +90,26 @@ ENTRY(IEEE754_EXPL)
 #ifdef PIC
 	LOAD_PIC_REG (cx)
 #endif
-#ifndef USE_AS_EXPM1L
+#ifdef USE_AS_EXPM1L
+	xorb	$0x80, %ah
+	cmpl	$0xc006, %eax
+	fstsw	%ax
+	movb	$0x45, %dh
+	jb	4f
+
+	/* Below -64.0 (may be -NaN or -Inf). */
+	andb	%ah, %dh
+	cmpb	$0x01, %dh
+	je	2f		/* Is +-NaN, jump.  */
+	jmp	1f		/* -large, possibly -Inf.  */
+
+4:	/* In range -64.0 to 64.0 (may be +-0 but not NaN or +-Inf).  */
+	/* Test for +-0 as argument.  */
+	andb	%ah, %dh
+	cmpb	$0x40, %dh
+	je	2f
+#else
 	movzwl	4+8(%esp), %eax
-#endif
 	andl	$0x7fff, %eax
 	cmpl	$0x400d, %eax
 	jle	3f
@@ -108,16 +127,8 @@ ENTRY(IEEE754_EXPL)
 	andb	$2, %ah
 	jz	3f
 	fchs
-3:
-#ifdef USE_AS_EXPM1L
-	/* Test for +-0 as argument.  */
-	fstsw	%ax
-	movb	$0x45, %dh
-	andb	%ah, %dh
-	cmpb	$0x40, %dh
-	je	2f
 #endif
-	FLDLOG			/* 1  log2(base)      */
+3:	FLDLOG			/* 1  log2(base)      */
 	fmul	%st(1), %st	/* 1  x log2(base)    */
 	frndint			/* 1  i               */
 	fld	%st(1)		/* 2  x               */
@@ -154,13 +165,16 @@ ENTRY(IEEE754_EXPL)
 #endif
 	fstp	%st(1)		/* 0  */
 	jmp	2f
-1:	testl	$0x200, %eax	/* Test sign.  */
-	jz	2f		/* If positive, jump.  */
-	fstp	%st
+1:
 #ifdef USE_AS_EXPM1L
+	/* For expm1l, only negative sign gets here.  */
+	fstp	%st
 	fld1
 	fchs
 #else
+	testl	$0x200, %eax	/* Test sign.  */
+	jz	2f		/* If positive, jump.  */
+	fstp	%st
 	fldz			/* Set result to 0.  */
 #endif
 2:	ret
diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
index 97249194c7..ab02bde5ff 100644
--- a/sysdeps/i386/fpu/libm-test-ulps
+++ b/sysdeps/i386/fpu/libm-test-ulps
@@ -1779,6 +1779,9 @@ idouble: 1
 ifloat: 1
 
 # expm1
+Test "expm1 (-45.0) == -0.9999999999999999999713748141945060635553":
+ildouble: 1
+ldouble: 1
 Test "expm1 (1) == M_El - 1.0":
 ildouble: 1
 Test "expm1 (11356.25) == 9.05128237311923300051376115753226014206e+4931":
diff --git a/sysdeps/i386/fpu/s_expm1.S b/sysdeps/i386/fpu/s_expm1.S
index 9883f9b353..d2754de911 100644
--- a/sysdeps/i386/fpu/s_expm1.S
+++ b/sysdeps/i386/fpu/s_expm1.S
@@ -51,19 +51,31 @@ ENTRY(__expm1)
 	jae	HIDDEN_JUMPTARGET (__exp)
 
 	fldl	4(%esp)		// x
-	fxam			// Is NaN or +-Inf?
+	fxam			// Is NaN, +-Inf or +-0?
+	xorb	$0x80, %ah
+	cmpl	$0xc043, %eax	// is num <= -38.0?
 	fstsw	%ax
 	movb	$0x45, %ch
+	jb	4f
+
+	// Below -38.0 (may be -NaN or -Inf).
+	andb	%ah, %ch
+#ifdef	PIC
+	LOAD_PIC_REG (dx)
+#endif
+	cmpb	$0x01, %ch
+	je	5f		// If -NaN, jump.
+	jmp	2f		// -large, possibly -Inf.
+
+4:	// In range -38.0 to 704.0 (may be +-0 but not NaN or +-Inf).
 	andb	%ah, %ch
 	cmpb	$0x40, %ch
 	je	3f		// If +-0, jump.
 #ifdef	PIC
 	LOAD_PIC_REG (dx)
 #endif
-	cmpb	$0x05, %ch
-	je	2f		// If +-Inf, jump.
 
-	fldt	MO(l2e)		// log2(e) : x
+5:	fldt	MO(l2e)		// log2(e) : x
 	fmulp			// log2(e)*x
 	fld	%st		// log2(e)*x : log2(e)*x
 	frndint			// int(log2(e)*x) : log2(e)*x
@@ -79,9 +91,7 @@ ENTRY(__expm1)
 	fsubrp	%st, %st(1)	// 2^(log2(e)*x)
 	ret
 
-2:	testl	$0x200, %eax	// Test sign.
-	jz	3f		// If positive, jump.
-	fstp	%st
+2:	fstp	%st
 	fldl	MO(minus1)	// Set result to -1.0.
 3:	ret
 END(__expm1)
diff --git a/sysdeps/i386/fpu/s_expm1f.S b/sysdeps/i386/fpu/s_expm1f.S
index 45257d7524..fc82b92341 100644
--- a/sysdeps/i386/fpu/s_expm1f.S
+++ b/sysdeps/i386/fpu/s_expm1f.S
@@ -51,19 +51,31 @@ ENTRY(__expm1f)
 	jae	HIDDEN_JUMPTARGET (__expf)
 
 	flds	4(%esp)		// x
-	fxam			// Is NaN or +-Inf?
+	fxam			// Is NaN, +-Inf or +-0?
+	xorb	$0x80, %ah
+	cmpl	$0xc190, %eax	// is num <= -18.0?
 	fstsw	%ax
 	movb	$0x45, %ch
+	jb	4f
+
+	// Below -18.0 (may be -NaN or -Inf).
+	andb	%ah, %ch
+#ifdef	PIC
+	LOAD_PIC_REG (dx)
+#endif
+	cmpb	$0x01, %ch
+	je	5f		// If -NaN, jump.
+	jmp	2f		// -large, possibly -Inf.
+
+4:	// In range -18.0 to 88.5 (may be +-0 but not NaN or +-Inf).
 	andb	%ah, %ch
 	cmpb	$0x40, %ch
 	je	3f		// If +-0, jump.
 #ifdef	PIC
 	LOAD_PIC_REG (dx)
 #endif
-	cmpb	$0x05, %ch
-	je	2f		// If +-Inf, jump.
 
-	fldt	MO(l2e)		// log2(e) : x
+5:	fldt	MO(l2e)		// log2(e) : x
 	fmulp			// log2(e)*x
 	fld	%st		// log2(e)*x : log2(e)*x
 	frndint			// int(log2(e)*x) : log2(e)*x
@@ -79,9 +91,7 @@ ENTRY(__expm1f)
 	fsubrp	%st, %st(1)	// 2^(log2(e)*x)
 	ret
 
-2:	testl	$0x200, %eax	// Test sign.
-	jz	3f		// If positive, jump.
-	fstp	%st
+2:	fstp	%st
 	fldl	MO(minus1)	// Set result to -1.0.
 3:	ret
 END(__expm1f)
diff --git a/sysdeps/x86_64/fpu/e_expl.S b/sysdeps/x86_64/fpu/e_expl.S
index e6b842bf26..1c37c86f7f 100644
--- a/sysdeps/x86_64/fpu/e_expl.S
+++ b/sysdeps/x86_64/fpu/e_expl.S
@@ -60,10 +60,12 @@ c1:	.byte 0x20, 0xfa, 0xee, 0xc2, 0x5f, 0x70, 0xa5, 0xec, 0xed, 0x3f
 	.byte 0, 0, 0, 0, 0, 0
 	ASM_SIZE_DIRECTIVE(c1)
 #endif
+#ifndef USE_AS_EXPM1L
 	ASM_TYPE_DIRECTIVE(csat,@object)
 csat:	.byte 0, 0, 0, 0, 0, 0, 0, 0x80, 0x0e, 0x40
 	.byte 0, 0, 0, 0, 0, 0
 	ASM_SIZE_DIRECTIVE(csat)
+#endif
 
 #ifdef PIC
 # define MO(op) op##(%rip)
@@ -85,9 +87,26 @@ ENTRY(IEEE754_EXPL)
    For the i686 the code can be written better.
    -- drepper@cygnus.com.  */
 	fxam			/* Is NaN or +-Inf?  */
-#ifndef USE_AS_EXPM1L
+#ifdef USE_AS_EXPM1L
+	xorb	$0x80, %ah
+	cmpl	$0xc006, %eax
+	fstsw	%ax
+	movb	$0x45, %dh
+	jb	4f
+
+	/* Below -64.0 (may be -NaN or -Inf). */
+	andb	%ah, %dh
+	cmpb	$0x01, %dh
+	je	2f		/* Is +-NaN, jump.  */
+	jmp	1f		/* -large, possibly -Inf.  */
+
+4:	/* In range -64.0 to 64.0 (may be +-0 but not NaN or +-Inf).  */
+	/* Test for +-0 as argument.  */
+	andb	%ah, %dh
+	cmpb	$0x40, %dh
+	je	2f
+#else
 	movzwl	8+8(%rsp), %eax
-#endif
 	andl	$0x7fff, %eax
 	cmpl	$0x400d, %eax
 	jle	3f
@@ -105,16 +124,8 @@ ENTRY(IEEE754_EXPL)
 	andb	$2, %ah
 	jz	3f
 	fchs
-3:
-#ifdef USE_AS_EXPM1L
-	/* Test for +-0 as argument.  */
-	fstsw	%ax
-	movb	$0x45, %dh
-	andb	%ah, %dh
-	cmpb	$0x40, %dh
-	je	2f
 #endif
-	FLDLOG			/* 1  log2(base)      */
+3:	FLDLOG			/* 1  log2(base)      */
 	fmul	%st(1), %st	/* 1  x log2(base)    */
 	frndint			/* 1  i               */
 	fld	%st(1)		/* 2  x               */
@@ -151,13 +162,16 @@ ENTRY(IEEE754_EXPL)
 #endif
 	fstp	%st(1)		/* 0  */
 	jmp	2f
-1:	testl	$0x200, %eax	/* Test sign.  */
-	jz	2f		/* If positive, jump.  */
-	fstp	%st
+1:
 #ifdef USE_AS_EXPM1L
+	/* For expm1l, only negative sign gets here.  */
+	fstp	%st
 	fld1
 	fchs
 #else
+	testl	$0x200, %eax	/* Test sign.  */
+	jz	2f		/* If positive, jump.  */
+	fstp	%st
 	fldz			/* Set result to 0.  */
 #endif
 2:	ret
diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
index b64e52d7d8..be28024b51 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -1657,6 +1657,9 @@ float: 1
 ifloat: 1
 
 # expm1
+Test "expm1 (-45.0) == -0.9999999999999999999713748141945060635553":
+ildouble: 1
+ldouble: 1
 Test "expm1 (0.75) == 1.11700001661267466854536981983709561":
 double: 1
 idouble: 1