Fix spurious overflow exceptions from x86/x86_64 powl (bug 13872).

5 files changed, 53 insertions, 33 deletions

diff --git a/ChangeLog b/ChangeLog
index 3686491a2d..a85a1c06ef 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,13 @@
 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.
+
 	[BZ #11521]
 	* math/s_ctan.c: Include <float.h>.
 	(__ctan): Avoid internal overflow or cancellation in calculating
diff --git a/NEWS b/NEWS
index bb303f8b93..357f3b4e69 100644
--- a/NEWS
+++ b/NEWS
@@ -18,9 +18,10 @@ Version 2.16
   13530, 13531, 13532, 13533, 13547, 13551, 13552, 13553, 13555, 13559,
   13566, 13583, 13592, 13618, 13637, 13656, 13658, 13673, 13691, 13695,
   13704, 13705, 13706, 13726, 13738, 13760, 13761, 13786, 13792, 13806,
-  13824, 13840, 13841, 13844, 13846, 13851, 13852, 13854, 13871, 13873,
-  13879, 13883, 13892, 13895, 13908, 13910, 13911, 13912, 13913, 13915,
-  13916, 13917, 13918, 13919, 13920, 13921, 13926, 13928, 13938, 13963
+  13824, 13840, 13841, 13844, 13846, 13851, 13852, 13854, 13871, 13872,
+  13873, 13879, 13883, 13892, 13895, 13908, 13910, 13911, 13912, 13913,
+  13915, 13916, 13917, 13918, 13919, 13920, 13921, 13926, 13928, 13938,
+  13963
 
 * ISO C11 support:
 
diff --git a/math/libm-test.inc b/math/libm-test.inc
index a551b9f3f4..2809d57d09 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -5682,8 +5682,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);
@@ -5997,8 +5996,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);
@@ -6122,8 +6120,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 0e7c05bb82..5b166eab4b 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 0626ce4172..10ede22648 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)