This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Fix pow of zero and infinity to large powers
- From: "Joseph S. Myers" <joseph at codesourcery dot com>
- To: libc-alpha at sourceware dot org
- Date: Tue, 20 Mar 2012 15:05:17 +0000 (UTC)
- Subject: Fix pow of zero and infinity to large powers
This patch fixes various problems with zero and infinity to large
powers (roughly, cases (4b), (5), (6) and (7) of bug 3866). (I was
originally intending to look at the patch in bug 369, but before
adding comprehensive tests for powers of negative finite values it
seemed appropriate to cover powers of negative zero and infinity -
resulting in the testcases in this patch, and the code changes to make
them pass.)
The general issue is x86 floating-to-integer conversion instructions
requiring inputs in a certain range, so the exponents need checking
first to make sure they are in that range (and in the powl case,
there's a further issue that values in the range (-2^64, 2^64) could
be odd integers but are too large to convert directly to integer, so
fprem is used to reduce to 63 bits in that case). pow and powf had
bogus logic "OK, the value is an integer, but is the number of bits
small enough so that all are coming from the mantissa?", which might
be relevant if you were checking bits of the floating-point
representation to see if it's odd - but is not relevant to this code
which is checking the low bit of the result of converting to integer,
rather than looking at the representation. This bogus logic caused
one test failure (pow (-0, -0x1.fffffffffffffp+52) returned +infinity
instead of -infinity); I removed it in all cases. One bug in the
generic C pow, where it failed to raise divide-by-zero for zero to
large negative powers, is also fixed.
One bug found in the course of testing, bug 13879 (a spurious overflow
exception), was filed in Bugzilla and the exception marked as allowed
with a comment pointing to the bug instead of fixing it (this looks
like a bug in the wrapper calling __kernel_standard in certain cases,
where it passes original long double arguments to __kernel_standard
which takes double, so causing overflow exceptions from the
conversion; probably all __kernel_standard calls need reviewing to
work out how to address this issue and how many places it affects).
Tested x86 and x86_64.
2012-03-20 Joseph Myers <joseph@codesourcery.com>
[BZ #3866]
* sysdeps/i386/fpu/e_pow.S (__ieee754_pow): Test for y outside the
range of signed 64-bit integers before using fistpll. Remove
checks for whether integers fit in mantissa bits.
* sysdeps/i386/fpu/e_powf.S (__ieee754_powf): Test for y outside
the range of signed 32-bit integers before using fistpl. Remove
checks for whether integers fit in mantissa bits.
* sysdeps/i386/fpu/e_powl.S (p64): New object.
(__ieee754_powl): Test for y outside the range of signed 64-bit
integers before using fistpll. Reduce 64-bit values to 63-bit
ones as needed.
* sysdeps/ieee754/dbl-64/e_pow.c (__ieee754_pow): Ensure
divide-by-zero is raised for zero to large negative powers.
* sysdeps/x86_64/fpu/e_powl.S (p64): New object.
(__ieee754_powl): Test for y outside the range of signed 64-bit
integers before using fistpll. Reduce 64-bit values to 63-bit
ones as needed.
* math/libm-test.inc (pow_test): Add more tests.
diff --git a/math/libm-test.inc b/math/libm-test.inc
index af3d645..64d12a5 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -5431,11 +5431,75 @@ pow_test (void)
TEST_ff_f (pow, 0, -11, plus_infty, DIVIDE_BY_ZERO_EXCEPTION);
check_int ("errno for pow(0,-odd) == ERANGE", errno, ERANGE, 0, 0, 0);
errno = 0;
+ TEST_ff_f (pow, 0, -0xffffff, plus_infty, DIVIDE_BY_ZERO_EXCEPTION);
+ check_int ("errno for pow(0,-odd) == ERANGE", errno, ERANGE, 0, 0, 0);
+#ifndef TEST_FLOAT
+ errno = 0;
+ TEST_ff_f (pow, 0, -0x1.fffffffffffffp+52L, plus_infty, DIVIDE_BY_ZERO_EXCEPTION);
+ check_int ("errno for pow(0,-odd) == ERANGE", errno, ERANGE, 0, 0, 0);
+#endif
+#ifdef TEST_LDOUBLE
+# if LDBL_MANT_DIG >= 64
+ errno = 0;
+ TEST_ff_f (pow, 0, -0x1.fffffffffffffffep+63L, plus_infty, DIVIDE_BY_ZERO_EXCEPTION);
+ check_int ("errno for pow(0,-odd) == ERANGE", errno, ERANGE, 0, 0, 0);
+# endif
+# if LDBL_MANT_DIG >= 106
+ errno = 0;
+ TEST_ff_f (pow, 0, -0x1.ffffffffffffffffffffffffff8p+105L, plus_infty, DIVIDE_BY_ZERO_EXCEPTION);
+ check_int ("errno for pow(0,-odd) == ERANGE", errno, ERANGE, 0, 0, 0);
+# endif
+# if LDBL_MANT_DIG >= 113
+ errno = 0;
+ TEST_ff_f (pow, 0, -0x1.ffffffffffffffffffffffffffffp+112L, plus_infty, DIVIDE_BY_ZERO_EXCEPTION);
+ check_int ("errno for pow(0,-odd) == ERANGE", errno, ERANGE, 0, 0, 0);
+# endif
+#endif
TEST_ff_f (pow, minus_zero, -1, minus_infty, DIVIDE_BY_ZERO_EXCEPTION);
check_int ("errno for pow(-0,-odd) == ERANGE", errno, ERANGE, 0, 0, 0);
errno = 0;
TEST_ff_f (pow, minus_zero, -11L, minus_infty, DIVIDE_BY_ZERO_EXCEPTION);
check_int ("errno for pow(-0,-odd) == ERANGE", errno, ERANGE, 0, 0, 0);
+ errno = 0;
+ TEST_ff_f (pow, minus_zero, -0xffffff, minus_infty, DIVIDE_BY_ZERO_EXCEPTION);
+ check_int ("errno for pow(-0,-odd) == ERANGE", errno, ERANGE, 0, 0, 0);
+ errno = 0;
+ TEST_ff_f (pow, minus_zero, -0x1fffffe, plus_infty, DIVIDE_BY_ZERO_EXCEPTION);
+ check_int ("errno for pow(-0,-even) == ERANGE", errno, ERANGE, 0, 0, 0);
+#ifndef TEST_FLOAT
+ errno = 0;
+ TEST_ff_f (pow, minus_zero, -0x1.fffffffffffffp+52L, minus_infty, DIVIDE_BY_ZERO_EXCEPTION);
+ check_int ("errno for pow(-0,-odd) == ERANGE", errno, ERANGE, 0, 0, 0);
+ errno = 0;
+ TEST_ff_f (pow, minus_zero, -0x1.fffffffffffffp+53L, plus_infty, DIVIDE_BY_ZERO_EXCEPTION);
+ check_int ("errno for pow(-0,-even) == ERANGE", errno, ERANGE, 0, 0, 0);
+#endif
+#ifdef TEST_LDOUBLE
+# if LDBL_MANT_DIG >= 64
+ errno = 0;
+ TEST_ff_f (pow, minus_zero, -0x1.fffffffffffffffep+63L, minus_infty, DIVIDE_BY_ZERO_EXCEPTION);
+ check_int ("errno for pow(-0,-odd) == ERANGE", errno, ERANGE, 0, 0, 0);
+ errno = 0;
+ TEST_ff_f (pow, minus_zero, -0x1.fffffffffffffffep+64L, plus_infty, DIVIDE_BY_ZERO_EXCEPTION);
+ check_int ("errno for pow(-0,-even) == ERANGE", errno, ERANGE, 0, 0, 0);
+# endif
+# if LDBL_MANT_DIG >= 106
+ errno = 0;
+ TEST_ff_f (pow, minus_zero, -0x1.ffffffffffffffffffffffffff8p+105L, minus_infty, DIVIDE_BY_ZERO_EXCEPTION);
+ check_int ("errno for pow(-0,-odd) == ERANGE", errno, ERANGE, 0, 0, 0);
+ errno = 0;
+ TEST_ff_f (pow, minus_zero, -0x1.ffffffffffffffffffffffffff8p+106L, plus_infty, DIVIDE_BY_ZERO_EXCEPTION);
+ check_int ("errno for pow(-0,-even) == ERANGE", errno, ERANGE, 0, 0, 0);
+# endif
+# if LDBL_MANT_DIG >= 113
+ errno = 0;
+ TEST_ff_f (pow, minus_zero, -0x1.ffffffffffffffffffffffffffffp+112L, minus_infty, DIVIDE_BY_ZERO_EXCEPTION);
+ check_int ("errno for pow(-0,-odd) == ERANGE", errno, ERANGE, 0, 0, 0);
+ errno = 0;
+ TEST_ff_f (pow, minus_zero, -0x1.ffffffffffffffffffffffffffffp+113L, plus_infty, DIVIDE_BY_ZERO_EXCEPTION);
+ check_int ("errno for pow(-0,-even) == ERANGE", errno, ERANGE, 0, 0, 0);
+# endif
+#endif
errno = 0;
TEST_ff_f (pow, 0, -2, plus_infty, DIVIDE_BY_ZERO_EXCEPTION);
@@ -5444,11 +5508,31 @@ pow_test (void)
TEST_ff_f (pow, 0, -11.1L, plus_infty, DIVIDE_BY_ZERO_EXCEPTION);
check_int ("errno for pow(0,-num) == ERANGE", errno, ERANGE, 0, 0, 0);
errno = 0;
+ TEST_ff_f (pow, 0, -0x1p24, plus_infty, DIVIDE_BY_ZERO_EXCEPTION);
+ check_int ("errno for pow(0,-num) == ERANGE", errno, ERANGE, 0, 0, 0);
+ errno = 0;
+ TEST_ff_f (pow, 0, -0x1p127, plus_infty, DIVIDE_BY_ZERO_EXCEPTION);
+ check_int ("errno for pow(0,-num) == ERANGE", errno, ERANGE, 0, 0, 0);
+ errno = 0;
+ /* Bug 13879: spurious OVERFLOW exception may be present. */
+ TEST_ff_f (pow, 0, -max_value, plus_infty, DIVIDE_BY_ZERO_EXCEPTION|OVERFLOW_EXCEPTION_OK);
+ check_int ("errno for pow(0,-num) == ERANGE", errno, ERANGE, 0, 0, 0);
+ errno = 0;
TEST_ff_f (pow, minus_zero, -2, plus_infty, DIVIDE_BY_ZERO_EXCEPTION);
check_int ("errno for pow(-0,-even) == ERANGE", errno, ERANGE, 0, 0, 0);
errno = 0;
TEST_ff_f (pow, minus_zero, -11.1L, plus_infty, DIVIDE_BY_ZERO_EXCEPTION);
check_int ("errno for pow(-0,-num) == ERANGE", errno, ERANGE, 0, 0, 0);
+ errno = 0;
+ TEST_ff_f (pow, minus_zero, -0x1p24, plus_infty, DIVIDE_BY_ZERO_EXCEPTION);
+ check_int ("errno for pow(-0,-num) == ERANGE", errno, ERANGE, 0, 0, 0);
+ errno = 0;
+ TEST_ff_f (pow, minus_zero, -0x1p127, plus_infty, DIVIDE_BY_ZERO_EXCEPTION);
+ check_int ("errno for pow(-0,-num) == ERANGE", errno, ERANGE, 0, 0, 0);
+ errno = 0;
+ /* Bug 13879: spurious OVERFLOW exception may be present. */
+ TEST_ff_f (pow, minus_zero, -max_value, plus_infty, DIVIDE_BY_ZERO_EXCEPTION|OVERFLOW_EXCEPTION_OK);
+ check_int ("errno for pow(-0,-num) == ERANGE", errno, ERANGE, 0, 0, 0);
TEST_ff_f (pow, 0x1p72L, 0x1p72L, plus_infty, OVERFLOW_EXCEPTION);
TEST_ff_f (pow, 10, -0x1p72L, 0);
@@ -5487,32 +5571,155 @@ pow_test (void)
/* pow (+inf, y) == +inf for y > 0. */
TEST_ff_f (pow, plus_infty, 2, plus_infty);
+ TEST_ff_f (pow, plus_infty, 0xffffff, plus_infty);
+#ifndef TEST_FLOAT
+ TEST_ff_f (pow, plus_infty, 0x1.fffffffffffffp+52L, plus_infty);
+#endif
+#ifdef TEST_LDOUBLE
+# if LDBL_MANT_DIG >= 64
+ TEST_ff_f (pow, plus_infty, 0x1.fffffffffffffffep+63L, plus_infty);
+# endif
+# if LDBL_MANT_DIG >= 106
+ TEST_ff_f (pow, plus_infty, 0x1.ffffffffffffffffffffffffff8p+105L, plus_infty);
+# endif
+# if LDBL_MANT_DIG >= 113
+ TEST_ff_f (pow, plus_infty, 0x1.ffffffffffffffffffffffffffffp+112L, plus_infty);
+# endif
+#endif
+ TEST_ff_f (pow, plus_infty, 0x1p24, plus_infty);
+ TEST_ff_f (pow, plus_infty, 0x1p127, plus_infty);
+ TEST_ff_f (pow, plus_infty, max_value, plus_infty);
/* pow (+inf, y) == +0 for y < 0. */
TEST_ff_f (pow, plus_infty, -1, 0.0);
+ TEST_ff_f (pow, plus_infty, -0xffffff, 0.0);
+#ifndef TEST_FLOAT
+ TEST_ff_f (pow, plus_infty, -0x1.fffffffffffffp+52L, 0.0);
+#endif
+#ifdef TEST_LDOUBLE
+# if LDBL_MANT_DIG >= 64
+ TEST_ff_f (pow, plus_infty, -0x1.fffffffffffffffep+63L, 0.0);
+# endif
+# if LDBL_MANT_DIG >= 106
+ TEST_ff_f (pow, plus_infty, -0x1.ffffffffffffffffffffffffff8p+105L, 0.0);
+# endif
+# if LDBL_MANT_DIG >= 113
+ TEST_ff_f (pow, plus_infty, -0x1.ffffffffffffffffffffffffffffp+112L, 0.0);
+# endif
+#endif
+ TEST_ff_f (pow, plus_infty, -0x1p24, 0.0);
+ TEST_ff_f (pow, plus_infty, -0x1p127, 0.0);
+ TEST_ff_f (pow, plus_infty, -max_value, 0.0);
/* pow (-inf, y) == -inf for y an odd integer > 0. */
TEST_ff_f (pow, minus_infty, 27, minus_infty);
+ TEST_ff_f (pow, minus_infty, 0xffffff, minus_infty);
+ TEST_ff_f (pow, minus_infty, 0x1fffffe, plus_infty);
+#ifndef TEST_FLOAT
+ TEST_ff_f (pow, minus_infty, 0x1.fffffffffffffp+52L, minus_infty);
+ TEST_ff_f (pow, minus_infty, 0x1.fffffffffffffp+53L, plus_infty);
+#endif
+#ifdef TEST_LDOUBLE
+# if LDBL_MANT_DIG >= 64
+ TEST_ff_f (pow, minus_infty, 0x1.fffffffffffffffep+63L, minus_infty);
+ TEST_ff_f (pow, minus_infty, 0x1.fffffffffffffffep+64L, plus_infty);
+# endif
+# if LDBL_MANT_DIG >= 106
+ TEST_ff_f (pow, minus_infty, 0x1.ffffffffffffffffffffffffff8p+105L, minus_infty);
+ TEST_ff_f (pow, minus_infty, 0x1.ffffffffffffffffffffffffff8p+106L, plus_infty);
+# endif
+# if LDBL_MANT_DIG >= 113
+ TEST_ff_f (pow, minus_infty, 0x1.ffffffffffffffffffffffffffffp+112L, minus_infty);
+ TEST_ff_f (pow, minus_infty, 0x1.ffffffffffffffffffffffffffffp+113L, plus_infty);
+# endif
+#endif
/* pow (-inf, y) == +inf for y > 0 and not an odd integer. */
TEST_ff_f (pow, minus_infty, 28, plus_infty);
+ TEST_ff_f (pow, minus_infty, 0x1p24, plus_infty);
+ TEST_ff_f (pow, minus_infty, 0x1p127, plus_infty);
+ TEST_ff_f (pow, minus_infty, max_value, plus_infty);
/* pow (-inf, y) == -0 for y an odd integer < 0. */
TEST_ff_f (pow, minus_infty, -3, minus_zero);
+ TEST_ff_f (pow, minus_infty, -0xffffff, minus_zero);
+ TEST_ff_f (pow, minus_infty, -0x1fffffe, plus_zero);
+#ifndef TEST_FLOAT
+ TEST_ff_f (pow, minus_infty, -0x1.fffffffffffffp+52L, minus_zero);
+ TEST_ff_f (pow, minus_infty, -0x1.fffffffffffffp+53L, plus_zero);
+#endif
+#ifdef TEST_LDOUBLE
+# if LDBL_MANT_DIG >= 64
+ TEST_ff_f (pow, minus_infty, -0x1.fffffffffffffffep+63L, minus_zero);
+ TEST_ff_f (pow, minus_infty, -0x1.fffffffffffffffep+64L, plus_zero);
+# endif
+# if LDBL_MANT_DIG >= 106
+ TEST_ff_f (pow, minus_infty, -0x1.ffffffffffffffffffffffffff8p+105L, minus_zero);
+ TEST_ff_f (pow, minus_infty, -0x1.ffffffffffffffffffffffffff8p+106L, plus_zero);
+# endif
+# if LDBL_MANT_DIG >= 113
+ TEST_ff_f (pow, minus_infty, -0x1.ffffffffffffffffffffffffffffp+112L, minus_zero);
+ TEST_ff_f (pow, minus_infty, -0x1.ffffffffffffffffffffffffffffp+113L, plus_zero);
+# endif
+#endif
/* pow (-inf, y) == +0 for y < 0 and not an odd integer. */
TEST_ff_f (pow, minus_infty, -2.0, 0.0);
+ TEST_ff_f (pow, minus_infty, -0x1p24, 0.0);
+ TEST_ff_f (pow, minus_infty, -0x1p127, 0.0);
+ TEST_ff_f (pow, minus_infty, -max_value, 0.0);
/* pow (+0, y) == +0 for y an odd integer > 0. */
TEST_ff_f (pow, 0.0, 27, 0.0);
+ TEST_ff_f (pow, 0.0, 0xffffff, 0.0);
+#ifndef TEST_FLOAT
+ TEST_ff_f (pow, 0.0, 0x1.fffffffffffffp+52L, 0.0);
+#endif
+#ifdef TEST_LDOUBLE
+# if LDBL_MANT_DIG >= 64
+ TEST_ff_f (pow, 0.0, 0x1.fffffffffffffffep+63L, 0.0);
+# endif
+# if LDBL_MANT_DIG >= 106
+ TEST_ff_f (pow, 0.0, 0x1.ffffffffffffffffffffffffff8p+105L, 0.0);
+# endif
+# if LDBL_MANT_DIG >= 113
+ TEST_ff_f (pow, 0.0, 0x1.ffffffffffffffffffffffffffffp+112L, 0.0);
+# endif
+#endif
/* pow (-0, y) == -0 for y an odd integer > 0. */
TEST_ff_f (pow, minus_zero, 27, minus_zero);
+ TEST_ff_f (pow, minus_zero, 0xffffff, minus_zero);
+ TEST_ff_f (pow, minus_zero, 0x1fffffe, plus_zero);
+#ifndef TEST_FLOAT
+ TEST_ff_f (pow, minus_zero, 0x1.fffffffffffffp+52L, minus_zero);
+ TEST_ff_f (pow, minus_zero, 0x1.fffffffffffffp+53L, plus_zero);
+#endif
+#ifdef TEST_LDOUBLE
+# if LDBL_MANT_DIG >= 64
+ TEST_ff_f (pow, minus_zero, 0x1.fffffffffffffffep+63L, minus_zero);
+ TEST_ff_f (pow, minus_zero, 0x1.fffffffffffffffep+64L, plus_zero);
+# endif
+# if LDBL_MANT_DIG >= 106
+ TEST_ff_f (pow, minus_zero, 0x1.ffffffffffffffffffffffffff8p+105L, minus_zero);
+ TEST_ff_f (pow, minus_zero, 0x1.ffffffffffffffffffffffffff8p+106L, plus_zero);
+# endif
+# if LDBL_MANT_DIG >= 113
+ TEST_ff_f (pow, minus_zero, 0x1.ffffffffffffffffffffffffffffp+112L, minus_zero);
+ TEST_ff_f (pow, minus_zero, 0x1.ffffffffffffffffffffffffffffp+113L, plus_zero);
+# endif
+#endif
/* pow (+0, y) == +0 for y > 0 and not an odd integer. */
TEST_ff_f (pow, 0.0, 4, 0.0);
+ TEST_ff_f (pow, 0.0, 0x1p24, 0.0);
+ TEST_ff_f (pow, 0.0, 0x1p127, 0.0);
+ TEST_ff_f (pow, 0.0, max_value, 0.0);
/* pow (-0, y) == +0 for y > 0 and not an odd integer. */
TEST_ff_f (pow, minus_zero, 4, 0.0);
+ TEST_ff_f (pow, minus_zero, 0x1p24, 0.0);
+ TEST_ff_f (pow, minus_zero, 0x1p127, 0.0);
+ TEST_ff_f (pow, minus_zero, max_value, 0.0);
TEST_ff_f (pow, 16, 0.25L, 2);
TEST_ff_f (pow, 0x1p64L, 0.125L, 256);
diff --git a/sysdeps/i386/fpu/e_pow.S b/sysdeps/i386/fpu/e_pow.S
index 63c44f1..1abedf6 100644
--- a/sysdeps/i386/fpu/e_pow.S
+++ b/sysdeps/i386/fpu/e_pow.S
@@ -230,6 +230,16 @@ ENTRY(__ieee754_pow)
testb $2, %dh
jz 16f // jump if x == +inf
+ // fistpll raises invalid exception for |y| >= 1L<<63, so test
+ // that (in which case y is certainly even) before testing
+ // whether y is odd.
+ fld %st // y : y
+ fabs // |y| : y
+ fcompl MO(p63) // y
+ fnstsw
+ sahf
+ jnc 16f
+
// We must find out whether y is an odd integer.
fld %st // y : y
fistpll (%esp) // y
@@ -239,20 +249,13 @@ ENTRY(__ieee754_pow)
sahf
jne 17f
- // OK, the value is an integer, but is the number of bits small
- // enough so that all are coming from the mantissa?
+ // OK, the value is an integer.
popl %eax
cfi_adjust_cfa_offset (-4)
popl %edx
cfi_adjust_cfa_offset (-4)
andb $1, %al
jz 18f // jump if not odd
- movl %edx, %eax
- orl %edx, %edx
- jns 155f
- negl %eax
-155: cmpl $0x00200000, %eax
- ja 18f // does not fit in mantissa bits
// It's an odd integer.
shrl $31, %edx
fldl MOX(minf_mzero, %edx, 8)
@@ -289,6 +292,16 @@ ENTRY(__ieee754_pow)
testb $2, %dh
jz 25f
+ // fistpll raises invalid exception for |y| >= 1L<<63, so test
+ // that (in which case y is certainly even) before testing
+ // whether y is odd.
+ fld %st // y : y
+ fabs // |y| : y
+ fcompl MO(p63) // y
+ fnstsw
+ sahf
+ jnc 25f
+
fld %st // y : y
fistpll (%esp) // y
fildll (%esp) // int(y) : y
@@ -297,16 +310,13 @@ ENTRY(__ieee754_pow)
sahf
jne 26f
- // OK, the value is an integer, but is the number of bits small
- // enough so that all are coming from the mantissa?
+ // OK, the value is an integer.
popl %eax
cfi_adjust_cfa_offset (-4)
popl %edx
cfi_adjust_cfa_offset (-4)
andb $1, %al
jz 27f // jump if not odd
- cmpl $0xffe00000, %edx
- jbe 27f // does not fit in mantissa bits
// It's an odd integer.
// Raise divide-by-zero exception and get minus infinity value.
fldl MO(one)
@@ -329,6 +339,14 @@ ENTRY(__ieee754_pow)
21: testb $2, %dh
jz 22f
+ // fistpll raises invalid exception for |y| >= 1L<<63, so test
+ // that (in which case y is certainly even) before testing
+ // whether y is odd.
+ fcoml MO(p63) // y
+ fnstsw
+ sahf
+ jnc 22f
+
fld %st // y : y
fistpll (%esp) // y
fildll (%esp) // int(y) : y
@@ -337,16 +355,13 @@ ENTRY(__ieee754_pow)
sahf
jne 23f
- // OK, the value is an integer, but is the number of bits small
- // enough so that all are coming from the mantissa?
+ // OK, the value is an integer.
popl %eax
cfi_adjust_cfa_offset (-4)
popl %edx
cfi_adjust_cfa_offset (-4)
andb $1, %al
jz 24f // jump if not odd
- cmpl $0xffe00000, %edx
- jae 24f // does not fit in mantissa bits
// It's an odd integer.
fldl MO(mzero)
ret
diff --git a/sysdeps/i386/fpu/e_powf.S b/sysdeps/i386/fpu/e_powf.S
index cc8456d..aa58ed2 100644
--- a/sysdeps/i386/fpu/e_powf.S
+++ b/sysdeps/i386/fpu/e_powf.S
@@ -224,6 +224,16 @@ ENTRY(__ieee754_powf)
testb $2, %dh
jz 16f // jump if x == +inf
+ // fistpl raises invalid exception for |y| >= 1L<<31, so test
+ // that (in which case y is certainly even) before testing
+ // whether y is odd.
+ fld %st // y : y
+ fabs // |y| : y
+ fcompl MO(p31) // y
+ fnstsw
+ sahf
+ jnc 16f
+
// We must find out whether y is an odd integer.
fld %st // y : y
fistpl (%esp) // y
@@ -233,18 +243,11 @@ ENTRY(__ieee754_powf)
sahf
jne 17f
- // OK, the value is an integer, but is the number of bits small
- // enough so that all are coming from the mantissa?
+ // OK, the value is an integer.
popl %edx
cfi_adjust_cfa_offset (-4)
testb $1, %dl
jz 18f // jump if not odd
- movl %edx, %eax
- orl %edx, %edx
- jns 155f
- negl %eax
-155: cmpl $0x01000000, %eax
- ja 18f // does not fit in mantissa bits
// It's an odd integer.
shrl $31, %edx
fldl MOX(minf_mzero, %edx, 8)
@@ -281,6 +284,16 @@ ENTRY(__ieee754_powf)
testb $2, %dh
jz 25f
+ // fistpl raises invalid exception for |y| >= 1L<<31, so test
+ // that (in which case y is certainly even) before testing
+ // whether y is odd.
+ fld %st // y : y
+ fabs // |y| : y
+ fcompl MO(p31) // y
+ fnstsw
+ sahf
+ jnc 25f
+
fld %st // y : y
fistpl (%esp) // y
fildl (%esp) // int(y) : y
@@ -289,14 +302,11 @@ ENTRY(__ieee754_powf)
sahf
jne 26f
- // OK, the value is an integer, but is the number of bits small
- // enough so that all are coming from the mantissa?
+ // OK, the value is an integer.
popl %edx
cfi_adjust_cfa_offset (-4)
testb $1, %dl
jz 27f // jump if not odd
- cmpl $0xff000000, %edx
- jbe 27f // does not fit in mantissa bits
// It's an odd integer.
// Raise divide-by-zero exception and get minus infinity value.
fldl MO(one)
@@ -319,6 +329,14 @@ ENTRY(__ieee754_powf)
21: testb $2, %dh
jz 22f
+ // fistpl raises invalid exception for |y| >= 1L<<31, so test
+ // that (in which case y is certainly even) before testing
+ // whether y is odd.
+ fcoml MO(p31) // y
+ fnstsw
+ sahf
+ jnc 22f
+
fld %st // y : y
fistpl (%esp) // y
fildl (%esp) // int(y) : y
@@ -327,14 +345,11 @@ ENTRY(__ieee754_powf)
sahf
jne 23f
- // OK, the value is an integer, but is the number of bits small
- // enough so that all are coming from the mantissa?
+ // OK, the value is an integer.
popl %edx
cfi_adjust_cfa_offset (-4)
testb $1, %dl
jz 24f // jump if not odd
- cmpl $0xff000000, %edx
- jae 24f // does not fit in mantissa bits
// It's an odd integer.
fldl MO(mzero)
ret
diff --git a/sysdeps/i386/fpu/e_powl.S b/sysdeps/i386/fpu/e_powl.S
index 5d85089..c0aa194c 100644
--- a/sysdeps/i386/fpu/e_powl.S
+++ b/sysdeps/i386/fpu/e_powl.S
@@ -32,6 +32,9 @@ limit: .double 0.29
ASM_TYPE_DIRECTIVE(p63,@object)
p63: .byte 0, 0, 0, 0, 0, 0, 0xe0, 0x43
ASM_SIZE_DIRECTIVE(p63)
+ ASM_TYPE_DIRECTIVE(p64,@object)
+p64: .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x43
+ ASM_SIZE_DIRECTIVE(p64)
.section .rodata.cst16,"aM",@progbits,16
@@ -243,6 +246,19 @@ ENTRY(__ieee754_powl)
testb $2, %dh
jz 16f // jump if x == +inf
+ // fistpll raises invalid exception for |y| >= 1L<<63, but y
+ // may be odd unless we know |y| >= 1L<<64.
+ fld %st // y : y
+ fabs // |y| : y
+ fcompl MO(p64) // y
+ fnstsw
+ sahf
+ jnc 16f
+ fldl MO(p63) // p63 : y
+ fxch // y : p63
+ fprem // y%p63 : p63
+ fstp %st(1) // y%p63
+
// We must find out whether y is an odd integer.
fld %st // y : y
fistpll (%esp) // y
@@ -295,6 +311,19 @@ ENTRY(__ieee754_powl)
testb $2, %dh
jz 25f
+ // fistpll raises invalid exception for |y| >= 1L<<63, but y
+ // may be odd unless we know |y| >= 1L<<64.
+ fld %st // y : y
+ fabs // |y| : y
+ fcompl MO(p64) // y
+ fnstsw
+ sahf
+ jnc 25f
+ fldl MO(p63) // p63 : y
+ fxch // y : p63
+ fprem // y%p63 : p63
+ fstp %st(1) // y%p63
+
fld %st // y : y
fistpll (%esp) // y
fildll (%esp) // int(y) : y
@@ -332,6 +361,18 @@ ENTRY(__ieee754_powl)
21: testb $2, %dh
jz 22f
+ // fistpll raises invalid exception for |y| >= 1L<<63, but y
+ // may be odd unless we know |y| >= 1L<<64.
+ fld %st // y : y
+ fcompl MO(p64) // y
+ fnstsw
+ sahf
+ jnc 22f
+ fldl MO(p63) // p63 : y
+ fxch // y : p63
+ fprem // y%p63 : p63
+ fstp %st(1) // y%p63
+
fld %st // y : y
fistpll (%esp) // y
fildll (%esp) // int(y) : y
diff --git a/sysdeps/ieee754/dbl-64/e_pow.c b/sysdeps/ieee754/dbl-64/e_pow.c
index f936a72..26ffaad 100644
--- a/sysdeps/ieee754/dbl-64/e_pow.c
+++ b/sysdeps/ieee754/dbl-64/e_pow.c
@@ -111,7 +111,7 @@ __ieee754_pow(double x, double y) {
if (((v.i[HIGH_HALF] & 0x7fffffff) == 0x7ff00000 && v.i[LOW_HALF] != 0)
|| (v.i[HIGH_HALF] & 0x7fffffff) > 0x7ff00000)
return y;
- if (ABS(y) > 1.0e20) return (y>0)?0:INF.x;
+ if (ABS(y) > 1.0e20) return (y>0)?0:1.0/ABS(x);
k = checkint(y);
if (k == -1)
return y < 0 ? 1.0/x : x;
diff --git a/sysdeps/x86_64/fpu/e_powl.S b/sysdeps/x86_64/fpu/e_powl.S
index bd6d828..7e6a7c9 100644
--- a/sysdeps/x86_64/fpu/e_powl.S
+++ b/sysdeps/x86_64/fpu/e_powl.S
@@ -32,6 +32,9 @@ limit: .double 0.29
ASM_TYPE_DIRECTIVE(p63,@object)
p63: .byte 0, 0, 0, 0, 0, 0, 0xe0, 0x43
ASM_SIZE_DIRECTIVE(p63)
+ ASM_TYPE_DIRECTIVE(p64,@object)
+p64: .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x43
+ ASM_SIZE_DIRECTIVE(p64)
.section .rodata.cst16,"aM",@progbits,16
@@ -227,6 +230,19 @@ ENTRY(__ieee754_powl)
testb $2, %dh
jz 16f // jump if x == +inf
+ // fistpll raises invalid exception for |y| >= 1L<<63, but y
+ // may be odd unless we know |y| >= 1L<<64.
+ fldl MO(p64) // 1L<<64 : y
+ fld %st(1) // y : 1L<<64 : y
+ fabs // |y| : 1L<<64 : y
+ fcomip %st(1), %st // 1L<<64 : y
+ fstp %st(0) // y
+ jnc 16f
+ fldl MO(p63) // p63 : y
+ fxch // y : p63
+ fprem // y%p63 : p63
+ fstp %st(1) // y%p63
+
// We must find out whether y is an odd integer.
fld %st // y : y
fistpll -8(%rsp) // y
@@ -284,6 +300,19 @@ ENTRY(__ieee754_powl)
testb $2, %dh
jz 25f
+ // fistpll raises invalid exception for |y| >= 1L<<63, but y
+ // may be odd unless we know |y| >= 1L<<64.
+ fldl MO(p64) // 1L<<64 : y
+ fld %st(1) // y : 1L<<64 : y
+ fabs // |y| : 1L<<64 : y
+ fcomip %st(1), %st // 1L<<64 : y
+ fstp %st(0) // y
+ jnc 25f
+ fldl MO(p63) // p63 : y
+ fxch // y : p63
+ fprem // y%p63 : p63
+ fstp %st(1) // y%p63
+
fld %st // y : y
fistpll -8(%rsp) // y
fildll -8(%rsp) // int(y) : y
@@ -315,6 +344,18 @@ ENTRY(__ieee754_powl)
21: testb $2, %dh
jz 22f
+ // fistpll raises invalid exception for |y| >= 1L<<63, but y
+ // may be odd unless we know |y| >= 1L<<64.
+ fldl MO(p64) // 1L<<64 : y
+ fxch // y : 1L<<64
+ fcomi %st(1), %st // y : 1L<<64
+ fstp %st(1) // y
+ jnc 22f
+ fldl MO(p63) // p63 : y
+ fxch // y : p63
+ fprem // y%p63 : p63
+ fstp %st(1) // y%p63
+
fld %st // y : y
fistpll -8(%rsp) // y
fildll -8(%rsp) // int(y) : y
--
Joseph S. Myers
joseph@codesourcery.com