aboutsummaryrefslogtreecommitdiff
path: root/sysdeps/i386
diff options
context:
space:
mode:
authorJoseph Myers <joseph@codesourcery.com>2012-11-28 13:40:54 +0000
committerJoseph Myers <joseph@codesourcery.com>2012-11-28 13:40:54 +0000
commit1bead169c32a3a688de863709b863207b7aafddd (patch)
tree7c3dcf66e7b4d92a9bcc5e3bb67f2cbbb19172f3 /sysdeps/i386
parent0817d63dd1f8e165f8ef6590bf4feddf06705381 (diff)
downloadglibc-1bead169c32a3a688de863709b863207b7aafddd.zip
glibc-1bead169c32a3a688de863709b863207b7aafddd.tar.gz
glibc-1bead169c32a3a688de863709b863207b7aafddd.tar.bz2
Fix powl inaccuracy for x86_64 and x86 (bug 13881).
Diffstat (limited to 'sysdeps/i386')
-rw-r--r--sysdeps/i386/fpu/e_powl.S55
-rw-r--r--sysdeps/i386/fpu/libm-test-ulps9
2 files changed, 34 insertions, 30 deletions
diff --git a/sysdeps/i386/fpu/e_powl.S b/sysdeps/i386/fpu/e_powl.S
index ac4842c..7e29775 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 limit,@object
-limit: .double 0.29
- ASM_SIZE_DIRECTIVE(limit)
+ .type p3,@object
+p3: .byte 0, 0, 0, 0, 0, 0, 0x20, 0x40
+ ASM_SIZE_DIRECTIVE(p3)
.type p63,@object
p63: .byte 0, 0, 0, 0, 0, 0, 0xe0, 0x43
ASM_SIZE_DIRECTIVE(p63)
@@ -141,7 +141,15 @@ ENTRY(__ieee754_powl)
fchs // -0x1p-79 : x
jmp 3f
-9: /* OK, we have an integer value for y. */
+9: /* OK, we have an integer value for y. Unless very small
+ (we use < 8), 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
+ fnstsw
+ sahf
+ jnc 2f
popl %eax
cfi_adjust_cfa_offset (-4)
popl %edx
@@ -182,7 +190,7 @@ ENTRY(__ieee754_powl)
cfi_adjust_cfa_offset (8)
.align ALIGNARG(4)
-2: // y is a large integer (absolute value at least 1L<<63), but
+2: // y is a large integer (absolute value at least 8), but
// may be odd unless at least 1L<<64. So it may be necessary
// to adjust the sign of a negative result afterwards.
fxch // x : y
@@ -205,34 +213,21 @@ ENTRY(__ieee754_powl)
fchs // -(1L<<78) : |x|
.align ALIGNARG(4)
3: /* y is a real number. */
- fxch // x : y
- fldl MO(one) // 1.0 : x : y
- fldl MO(limit) // 0.29 : 1.0 : x : y
- fld %st(2) // x : 0.29 : 1.0 : x : y
- fsub %st(2) // x-1 : 0.29 : 1.0 : x : y
- fabs // |x-1| : 0.29 : 1.0 : x : y
- fucompp // 1.0 : x : y
- fnstsw
- fxch // x : 1.0 : y
- sahf
- ja 7f
- fsub %st(1) // x-1 : 1.0 : y
- fyl2xp1 // log2(x) : y
- jmp 8f
-
-7: fyl2x // log2(x) : y
-8: fmul %st(1) // y*log2(x) : y
- 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))
- fxch // fract(y*log2(x)) : int(y*log2(x))
- f2xm1 // 2^fract(y*log2(x))-1 : int(y*log2(x))
- 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))
+ subl $28, %esp
+ cfi_adjust_cfa_offset (28)
+ fstpt 12(%esp) // x
+ fstpt (%esp) // <empty>
+ mov %edx, 24(%esp)
+ call HIDDEN_JUMPTARGET (__powl_helper) // <result>
+ mov 24(%esp), %edx
+ addl $28, %esp
+ cfi_adjust_cfa_offset (-28)
testb $2, %dh
jz 292f
// x is negative. If y is an odd integer, negate the result.
+#ifdef PIC
+ LOAD_PIC_REG (cx)
+#endif
fldt 24(%esp) // y : abs(result)
fld %st // y : y : abs(result)
fabs // |y| : y : abs(result)
diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
index 8be0086..5b595bc 100644
--- a/sysdeps/i386/fpu/libm-test-ulps
+++ b/sysdeps/i386/fpu/libm-test-ulps
@@ -2480,6 +2480,11 @@ ifloat: 1
ildouble: 1
ldouble: 1
+# pow
+Test "pow (0x0.ffffffp0, -0x1p24) == 2.7182819094701610539628664526874952929416":
+ildouble: 1
+ldouble: 1
+
# pow_downward
Test "pow_downward (1.0625, 1.125) == 1.070582293028761362162622578677070098674":
double: 1
@@ -3782,6 +3787,10 @@ ifloat: 1
ildouble: 1
ldouble: 1
+Function: "pow":
+ildouble: 1
+ldouble: 1
+
Function: "pow_downward":
double: 1
float: 1