Fix expm1 spurious underflow exceptions (bug 6778).

This commit is contained in:
Joseph Myers 2012-07-06 11:17:41 +00:00
parent fb21f89b75
commit f17ac40d7c
9 changed files with 140 additions and 45 deletions

View file

@ -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".

2
NEWS
View file

@ -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

View file

@ -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);
}

View file

@ -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

View file

@ -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":

View file

@ -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)

View file

@ -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)

View file

@ -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

View file

@ -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