This is the mail archive of the libc-alpha@sourceware.org mailing list for the glibc project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

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


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]