diff mbox

Fix x86_64 / x86 powl inaccuracy for integer exponents (bug 19848) [committed]

Message ID alpine.DEB.2.10.1603240133160.9621@digraph.polyomino.org.uk
State New
Headers show

Commit Message

Joseph Myers March 24, 2016, 1:33 a.m. UTC 
Bug 19848 reports cases where powl on x86 / x86_64 has error
accumulation, for small integer exponents, larger than permitted by
glibc's accuracy goals, at least in some rounding modes.  This patch
further restricts the exponent range for which the
small-integer-exponent logic is used to limit the possible error
accumulation.

Tested for x86_64 and x86 and ulps updated accordingly.  Committed.

(auto-libm-test-out diffs omitted below.)

2016-03-24  Joseph Myers  <joseph@codesourcery.com>

	[BZ #19848]
	* sysdeps/i386/fpu/e_powl.S (p3): Rename to p2 and change value
	from 8 to 4.
	(__ieee754_powl): Compare integer exponent against 4 not 8.
	* sysdeps/x86_64/fpu/e_powl.S (p3): Rename to p2 and change value
	from 8 to 4.
	(__ieee754_powl): Compare integer exponent against 4 not 8.
	* math/auto-libm-test-in: Add more tests of pow.
	* math/auto-libm-test-out: Regenerated.
	* sysdeps/i386/i686/fpu/multiarch/libm-test-ulps: Update.
	* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
diff mbox

Patch

diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index 4b753de..34a6323 100644
--- a/math/auto-libm-test-in
+++ b/math/auto-libm-test-in
@@ -3717,6 +3717,44 @@  pow 0x2.000582p0 -1022
 pow 2 -0x3.fe513p+8
 pow 2 -0x3.fe4e8p+8
 
+pow 10 -1
+pow 10 -2
+pow 10 -3
+pow 10 -4
+pow 10 -5
+pow 10 -6
+pow 10 -7
+
+pow 0x0.ffffffffffffffffp0 1
+pow 0x0.ffffffffffffffffp0 2
+pow 0x0.ffffffffffffffffp0 3
+pow 0x0.ffffffffffffffffp0 4
+pow 0x0.ffffffffffffffffp0 5
+pow 0x0.ffffffffffffffffp0 6
+pow 0x0.ffffffffffffffffp0 7
+pow 0x0.ffffffffffffffffp0 -1
+pow 0x0.ffffffffffffffffp0 -2
+pow 0x0.ffffffffffffffffp0 -3
+pow 0x0.ffffffffffffffffp0 -4
+pow 0x0.ffffffffffffffffp0 -5
+pow 0x0.ffffffffffffffffp0 -6
+pow 0x0.ffffffffffffffffp0 -7
+
+pow 0x1.0000000000000002p0 1
+pow 0x1.0000000000000002p0 2
+pow 0x1.0000000000000002p0 3
+pow 0x1.0000000000000002p0 4
+pow 0x1.0000000000000002p0 5
+pow 0x1.0000000000000002p0 6
+pow 0x1.0000000000000002p0 7
+pow 0x1.0000000000000002p0 -1
+pow 0x1.0000000000000002p0 -2
+pow 0x1.0000000000000002p0 -3
+pow 0x1.0000000000000002p0 -4
+pow 0x1.0000000000000002p0 -5
+pow 0x1.0000000000000002p0 -6
+pow 0x1.0000000000000002p0 -7
+
 pow 1.0625 1.125
 pow 1.5 1.03125
 pow 0x1.7d1a0a6f2p+681 1.5
diff --git a/sysdeps/i386/fpu/e_powl.S b/sysdeps/i386/fpu/e_powl.S
index 77d2abf..923ee37 100644
--- a/sysdeps/i386/fpu/e_powl.S
+++ b/sysdeps/i386/fpu/e_powl.S
@@ -26,9 +26,9 @@ 
 	.type one,@object
 one:	.double 1.0
 	ASM_SIZE_DIRECTIVE(one)
-	.type p3,@object
-p3:	.byte 0, 0, 0, 0, 0, 0, 0x20, 0x40
-	ASM_SIZE_DIRECTIVE(p3)
+	.type p2,@object
+p2:	.byte 0, 0, 0, 0, 0, 0, 0x10, 0x40
+	ASM_SIZE_DIRECTIVE(p2)
 	.type p63,@object
 p63:	.byte 0, 0, 0, 0, 0, 0, 0xe0, 0x43
 	ASM_SIZE_DIRECTIVE(p63)
@@ -146,11 +146,11 @@  ENTRY(__ieee754_powl)
 	jmp	3f
 
 9:	/* OK, we have an integer value for y.  Unless very small
-	   (we use < 8), use the algorithm for real exponent to avoid
+	   (we use < 4), use the algorithm for real exponent to avoid
 	   accumulation of errors.  */
 	fld	%st		// y : y : x
 	fabs			// |y| : y : x
-	fcompl	MO(p3)		// y : x
+	fcompl	MO(p2)		// y : x
 	fnstsw
 	sahf
 	jnc	3f
diff --git a/sysdeps/i386/i686/fpu/multiarch/libm-test-ulps b/sysdeps/i386/i686/fpu/multiarch/libm-test-ulps
index bbb644a..585be78 100644
--- a/sysdeps/i386/i686/fpu/multiarch/libm-test-ulps
+++ b/sysdeps/i386/i686/fpu/multiarch/libm-test-ulps
@@ -1903,14 +1903,14 @@  ldouble: 4
 Function: "pow_towardzero":
 double: 1
 idouble: 1
-ildouble: 1
-ldouble: 1
+ildouble: 4
+ldouble: 4
 
 Function: "pow_upward":
 double: 1
 idouble: 1
-ildouble: 2
-ldouble: 2
+ildouble: 4
+ldouble: 4
 
 Function: "sin":
 ildouble: 1
diff --git a/sysdeps/x86_64/fpu/e_powl.S b/sysdeps/x86_64/fpu/e_powl.S
index 1f68cf0..4a7f3a1 100644
--- a/sysdeps/x86_64/fpu/e_powl.S
+++ b/sysdeps/x86_64/fpu/e_powl.S
@@ -26,9 +26,9 @@ 
 	.type one,@object
 one:	.double 1.0
 	ASM_SIZE_DIRECTIVE(one)
-	.type p3,@object
-p3:	.byte 0, 0, 0, 0, 0, 0, 0x20, 0x40
-	ASM_SIZE_DIRECTIVE(p3)
+	.type p2,@object
+p2:	.byte 0, 0, 0, 0, 0, 0, 0x10, 0x40
+	ASM_SIZE_DIRECTIVE(p2)
 	.type p63,@object
 p63:	.byte 0, 0, 0, 0, 0, 0, 0xe0, 0x43
 	ASM_SIZE_DIRECTIVE(p63)
@@ -136,12 +136,12 @@  ENTRY(__ieee754_powl)
 	jmp	3f
 
 9:	/* OK, we have an integer value for y.  Unless very small
-	   (we use < 8), use the algorithm for real exponent to avoid
+	   (we use < 4), use the algorithm for real exponent to avoid
 	   accumulation of errors.  */
-	fldl	MO(p3)		// 8 : y : x
-	fld	%st(1)		// y : 8 : y : x
-	fabs			// |y| : 8 : y : x
-	fcomip	%st(1), %st	// 8 : y : x
+	fldl	MO(p2)		// 4 : y : x
+	fld	%st(1)		// y : 4 : y : x
+	fabs			// |y| : 4 : y : x
+	fcomip	%st(1), %st	// 4 : y : x
 	fstp	%st(0)		// y : x
 	jnc	3f
 	mov	-8(%rsp),%eax
diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
index 445b475..7e7707b 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -2016,16 +2016,16 @@  double: 1
 float: 1
 idouble: 1
 ifloat: 1
-ildouble: 1
-ldouble: 1
+ildouble: 4
+ldouble: 4
 
 Function: "pow_upward":
 double: 1
 float: 1
 idouble: 1
 ifloat: 1
-ildouble: 2
-ldouble: 2
+ildouble: 4
+ldouble: 4
 
 Function: "pow_vlen16":
 float: 3