This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Fix powf inaccuracy (bug 18956) [committed]
- From: Joseph Myers <joseph at codesourcery dot com>
- To: <libc-alpha at sourceware dot org>
- Date: Sat, 26 Sep 2015 00:27:47 +0000
- Subject: Fix powf inaccuracy (bug 18956) [committed]
- Authentication-results: sourceware.org; auth=none
The flt-32 version of powf can be inaccurate because of bugs in the
extra-precision calculation of (x-1)/(x+1) or (x-1.5)/(x+1.5) as part
of calculating log(x) with extra precision: a constant used (as part
of adding 1 or 1.5 through integer arithmetic) is incorrect, and then
the code fails to mask a computed high part before using it in
arithmetic that relies on s_h*t_h being exactly representable. This
patch fixes these bugs.
Tested for x86_64 and x86. x86_64 ulps for powf removed and
regenerated to reflect reduced ulps from the increased accuracy for
existing tests. Committed.
(auto-libm-test-out diffs omitted below.)
2015-09-26 Joseph Myers <joseph@codesourcery.com>
[BZ #18956]
* sysdeps/ieee754/flt-32/e_powf.c (__ieee754_powf): Add 0x00400000
not 0x0040000 for high bit of mantissa. Mask with 0xfffff000 when
extracting high part.
* math/auto-libm-test-in: Add another test of pow.
* math/auto-libm-test-out: Regenerated.
* sysdeps/x86_64/fpu/libm-test-ulps: Update.
diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index 4b93857..e86be23 100644
--- a/math/auto-libm-test-in
+++ b/math/auto-libm-test-in
@@ -3219,6 +3219,7 @@ pow 0x1.ce78f2p+0 -0x2.7f1f78p+4
pow 0xf.fffffp+124 -0x5.b5b648p+0
pow 0x1.430d4cp+0 0x5.0e462p+4
pow 0x9.8b82ap-4 -0x1.99907ap+12
+pow 0xd.73035p-4 -0x1.47bb8p+8
sin 0
sin -0
diff --git a/sysdeps/ieee754/flt-32/e_powf.c b/sysdeps/ieee754/flt-32/e_powf.c
index 1804296..c72fe37 100644
--- a/sysdeps/ieee754/flt-32/e_powf.c
+++ b/sysdeps/ieee754/flt-32/e_powf.c
@@ -166,7 +166,9 @@ __ieee754_powf(float x, float y)
GET_FLOAT_WORD(is,s_h);
SET_FLOAT_WORD(s_h,is&0xfffff000);
/* t_h=ax+bp[k] High */
- SET_FLOAT_WORD(t_h,((ix>>1)|0x20000000)+0x0040000+(k<<21));
+ SET_FLOAT_WORD (t_h,
+ ((((ix>>1)|0x20000000)+0x00400000+(k<<21))
+ & 0xfffff000));
t_l = ax - (t_h-bp[k]);
s_l = v*((u-s_h*t_h)-s_h*t_l);
/* compute log(ax) */
diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
index 127a8e1..12c3dd1 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -1966,8 +1966,8 @@ Function: "log_vlen8_avx2":
float: 2
Function: "pow":
-float: 3
-ifloat: 3
+float: 1
+ifloat: 1
ildouble: 1
ldouble: 1
@@ -2003,25 +2003,25 @@ ldouble: 2
Function: "pow_downward":
double: 1
-float: 4
+float: 1
idouble: 1
-ifloat: 4
+ifloat: 1
ildouble: 4
ldouble: 4
Function: "pow_towardzero":
double: 1
-float: 8
+float: 1
idouble: 1
-ifloat: 8
+ifloat: 1
ildouble: 1
ldouble: 1
Function: "pow_upward":
double: 1
-float: 8
+float: 1
idouble: 1
-ifloat: 8
+ifloat: 1
ildouble: 2
ldouble: 2
--
Joseph S. Myers
joseph@codesourcery.com