Fix expm1 spurious underflow exceptions (bug 6778).

6 files changed, 96 insertions, 42 deletions

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