This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Fix expm1 spurious underflow exceptions (bug 6778)
- From: "Joseph S. Myers" <joseph at codesourcery dot com>
- To: <libc-alpha at sourceware dot org>
- Date: Thu, 5 Jul 2012 22:08:12 +0000
- Subject: Fix expm1 spurious underflow exceptions (bug 6778)
Bug 6778 is spurious underflow from expm1 (x86 versions, including
x86_64 for expm1l only) for large negative arguments (where exp would
underflow, but expm1 should return a value very close to -1).
This patch fixes this by adding appropriate checks for large negative
arguments to each function. (The new tests of expm1 for negative
arguments are intended to cover ranges around where the result starts
to round to -1 exactly for float, double and the supported long double
formats.)
Tested x86 and x86_64 and ulps updated accordingly.
2012-07-05 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.
diff --git a/math/libm-test.inc b/math/libm-test.inc
index 6adbb61..b87a40d 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 bab0a08..e42c9a1 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 9724919..ab02bde 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 9883f9b..d2754de 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 45257d7..fc82b92 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 e6b842b..1c37c86 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 b64e52d..be28024 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
--
Joseph S. Myers
joseph@codesourcery.com