Fix pow of zero and infinity to large powers.

5 files changed, 145 insertions, 33 deletions

diff --git a/sysdeps/i386/fpu/e_pow.S b/sysdeps/i386/fpu/e_pow.S
index 63c44f1357..1abedf6284 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 cc8456d280..aa58ed2b67 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 5d850897c5..c0aa194c62 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 f936a72de7..26ffaaddc6 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 bd6d828027..562791d302 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