Skip to content
Permalink
Browse files

Fix x86_64 / x86 powl inaccuracy for integer exponents (bug 19848).

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.

	[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.
  • Loading branch information...
Joseph Myers
Joseph Myers committed Mar 24, 2016
1 parent 7e1ff08 commit c898991d8bcfacc825097ba389ffccc5367c2b2d
@@ -1,3 +1,17 @@
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.

2016-03-23 Aurelien Jarno <aurelien@aurel32.net>

* sysdeps/unix/sysv/linux/futimens.c (futimens) [__NR_utimensat]:
@@ -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

Large diffs are not rendered by default.

Oops, something went wrong.
@@ -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
@@ -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
@@ -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
@@ -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

0 comments on commit c898991

Please sign in to comment.
You can’t perform that action at this time.