This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Fix spurious overflow exceptions from x86/x86_64 powl (bug 13872)
- From: "Joseph S. Myers" <joseph at codesourcery dot com>
- To: libc-alpha at sourceware dot org
- Date: Mon, 9 Apr 2012 17:32:00 +0000 (UTC)
- Subject: Fix spurious overflow exceptions from x86/x86_64 powl (bug 13872)
Bug 13872 is a problem with overflow exceptions for certain underflow
cases of x86 and x86_64 powl with very large exponents. The problem
is that the y*log2(x) computation overflows (to minus infinity), and
while there is code to check for such overflow (to avoid NaNs
appearing in the computation of the fractional part), this is too late
to avoid overflow exceptions from underflow cases. I propose this
patch to fix this problem by saturating the exponents to +/- 2**78
(powl(1-2**-64, 2**64) is about 1/e, powl(1+2**-63, 2**64) is about
e^2, and in either case raising to the power 2**14 results in
underflow to 0 / overflow to +infinity). This issue does not apply to
pow or powf because DBL_MAX * log2(DBL_TRUE_MIN) does not overflow the
range of long double and all these intermediate computations are in
long double.
Tested x86 and x86_64.
2012-04-09 Joseph Myers <joseph@codesourcery.com>
[BZ #13872]
* sysdeps/i386/fpu/e_powl.S (p78): New object.
(__ieee754_powl): Saturate large exponents rather than testing for
overflow of y*log2(x).
* sysdeps/x86_64/fpu/e_powl.S: Likewise.
* math/libm-test.inc (pow_test): Do not permit spurious overflow
exceptions.
diff --git a/math/libm-test.inc b/math/libm-test.inc
index 0533483..16efcb1 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -5622,8 +5622,7 @@ pow_test (void)
TEST_ff_f (pow, 0x1p72L, 0x1p72L, plus_infty, OVERFLOW_EXCEPTION);
TEST_ff_f (pow, 10, -0x1p72L, 0);
TEST_ff_f (pow, max_value, max_value, plus_infty, OVERFLOW_EXCEPTION);
- /* Bug 13872: spurious OVERFLOW exception may be present. */
- TEST_ff_f (pow, 10, -max_value, 0, OVERFLOW_EXCEPTION_OK);
+ TEST_ff_f (pow, 10, -max_value, 0);
TEST_ff_f (pow, 0, 1, 0);
TEST_ff_f (pow, 0, 11, 0);
@@ -5937,8 +5936,7 @@ pow_test (void)
TEST_ff_f (pow, -max_value, -0x1.ffffffffffffffffffffffffffffp+113L, plus_zero);
# endif
#endif
- /* Bug 13872: spurious OVERFLOW exception may be present. */
- TEST_ff_f (pow, -max_value, -max_value, plus_zero, OVERFLOW_EXCEPTION_OK);
+ TEST_ff_f (pow, -max_value, -max_value, plus_zero);
TEST_ff_f (pow, -max_value, 0xffffff, minus_infty, OVERFLOW_EXCEPTION);
TEST_ff_f (pow, -max_value, 0x1fffffe, plus_infty, OVERFLOW_EXCEPTION);
@@ -6062,8 +6060,7 @@ pow_test (void)
TEST_ff_f (pow, -min_value, 0x1.ffffffffffffffffffffffffffffp+113L, plus_zero);
# endif
#endif
- /* Bug 13872: spurious OVERFLOW exception may be present. */
- TEST_ff_f (pow, -min_value, max_value, plus_zero, OVERFLOW_EXCEPTION_OK);
+ TEST_ff_f (pow, -min_value, max_value, plus_zero);
#ifndef TEST_LDOUBLE /* Bug 13881. */
TEST_ff_f (pow, 0x0.ffffffp0, 10, 0.999999403953712118183885036774764444747L);
diff --git a/sysdeps/i386/fpu/e_powl.S b/sysdeps/i386/fpu/e_powl.S
index 0e7c05b..5b166ea 100644
--- a/sysdeps/i386/fpu/e_powl.S
+++ b/sysdeps/i386/fpu/e_powl.S
@@ -35,6 +35,9 @@ p63: .byte 0, 0, 0, 0, 0, 0, 0xe0, 0x43
ASM_TYPE_DIRECTIVE(p64,@object)
p64: .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x43
ASM_SIZE_DIRECTIVE(p64)
+ ASM_TYPE_DIRECTIVE(p78,@object)
+p78: .byte 0, 0, 0, 0, 0, 0, 0xd0, 0x44
+ ASM_SIZE_DIRECTIVE(p78)
.section .rodata.cst16,"aM",@progbits,16
@@ -166,6 +169,21 @@ ENTRY(__ieee754_powl)
fxch // x : y
fabs // |x| : y
fxch // y : |x|
+ // If y has absolute value at least 1L<<78, then any finite
+ // nonzero x will result in 0 (underflow), 1 or infinity (overflow).
+ // Saturate y to those bounds to avoid overflow in the calculation
+ // of y*log2(x).
+ fld %st // y : y : |x|
+ fabs // |y| : y : |x|
+ fcompl MO(p78) // y : |x|
+ fnstsw
+ sahf
+ jc 3f
+ fstp %st(0) // pop y
+ fldl MO(p78) // 1L<<78 : |x|
+ testb $2, %dl
+ jz 3f // y > 0
+ fchs // -(1L<<78) : |x|
.align ALIGNARG(4)
3: /* y is a real number. */
fxch // x : y
@@ -185,11 +203,6 @@ ENTRY(__ieee754_powl)
7: fyl2x // log2(x) : y
8: fmul %st(1) // y*log2(x) : y
- fxam
- fnstsw
- andb $0x45, %ah
- cmpb $0x05, %ah // is y*log2(x) == ħinf ?
- je 28f
fst %st(1) // y*log2(x) : y*log2(x)
frndint // int(y*log2(x)) : y*log2(x)
fsubr %st, %st(1) // int(y*log2(x)) : fract(y*log2(x))
@@ -198,13 +211,7 @@ ENTRY(__ieee754_powl)
faddl MO(one) // 2^fract(y*log2(x)) : int(y*log2(x))
fscale // 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x))
fstp %st(1) // 2^fract(y*log2(x))*2^int(y*log2(x))
- jmp 29f
-
-28: fstp %st(1) // y*log2(x)
- fldl MO(one) // 1 : y*log2(x)
- fscale // 2^(y*log2(x)) : y*log2(x)
- fstp %st(1) // 2^(y*log2(x))
-29: testb $2, %dh
+ testb $2, %dh
jz 292f
// x is negative. If y is an odd integer, negate the result.
fldt 24(%esp) // y : abs(result)
diff --git a/sysdeps/x86_64/fpu/e_powl.S b/sysdeps/x86_64/fpu/e_powl.S
index 0626ce4..10ede22 100644
--- a/sysdeps/x86_64/fpu/e_powl.S
+++ b/sysdeps/x86_64/fpu/e_powl.S
@@ -35,6 +35,9 @@ p63: .byte 0, 0, 0, 0, 0, 0, 0xe0, 0x43
ASM_TYPE_DIRECTIVE(p64,@object)
p64: .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x43
ASM_SIZE_DIRECTIVE(p64)
+ ASM_TYPE_DIRECTIVE(p78,@object)
+p78: .byte 0, 0, 0, 0, 0, 0, 0xd0, 0x44
+ ASM_SIZE_DIRECTIVE(p78)
.section .rodata.cst16,"aM",@progbits,16
@@ -151,6 +154,21 @@ ENTRY(__ieee754_powl)
fxch // x : y
fabs // |x| : y
fxch // y : |x|
+ // If y has absolute value at least 1L<<78, then any finite
+ // nonzero x will result in 0 (underflow), 1 or infinity (overflow).
+ // Saturate y to those bounds to avoid overflow in the calculation
+ // of y*log2(x).
+ fldl MO(p78) // 1L<<78 : y : |x|
+ fld %st(1) // y : 1L<<78 : y : |x|
+ fabs // |y| : 1L<<78 : y : |x|
+ fcomip %st(1), %st // 1L<<78 : y : |x|
+ fstp %st(0) // y : |x|
+ jc 3f
+ fstp %st(0) // pop y
+ fldl MO(p78) // 1L<<78 : |x|
+ testb $2, %dl
+ jz 3f // y > 0
+ fchs // -(1L<<78) : |x|
.align ALIGNARG(4)
3: /* y is a real number. */
fxch // x : y
@@ -170,11 +188,6 @@ ENTRY(__ieee754_powl)
7: fyl2x // log2(x) : y
8: fmul %st(1) // y*log2(x) : y
- fxam
- fnstsw
- andb $0x45, %ah
- cmpb $0x05, %ah // is y*log2(x) == ħinf ?
- je 28f
fst %st(1) // y*log2(x) : y*log2(x)
frndint // int(y*log2(x)) : y*log2(x)
fsubr %st, %st(1) // int(y*log2(x)) : fract(y*log2(x))
@@ -183,13 +196,7 @@ ENTRY(__ieee754_powl)
faddl MO(one) // 2^fract(y*log2(x)) : int(y*log2(x))
fscale // 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x))
fstp %st(1) // 2^fract(y*log2(x))*2^int(y*log2(x))
- jmp 29f
-
-28: fstp %st(1) // y*log2(x)
- fldl MO(one) // 1 : y*log2(x)
- fscale // 2^(y*log2(x)) : y*log2(x)
- fstp %st(1) // 2^(y*log2(x))
-29: testb $2, %dh
+ testb $2, %dh
jz 292f
// x is negative. If y is an odd integer, negate the result.
fldt 24(%rsp) // y : abs(result)
--
Joseph S. Myers
joseph@codesourcery.com