From 2460d3aa21f04cdf28497683bd3e29183189f779 Mon Sep 17 00:00:00 2001
From: Joseph Myers <joseph@codesourcery.com>
Date: Wed, 21 Mar 2012 12:14:57 +0000
Subject: Fix pow of zero and infinity to large powers.

---
 sysdeps/x86_64/fpu/e_powl.S | 41 +++++++++++++++++++++++++++++++++++++++++
 1 file changed, 41 insertions(+)

(limited to 'sysdeps/x86_64/fpu/e_powl.S')

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
-- 
cgit v1.2.3