about summary refs log tree commit diff
path: root/sysdeps
diff options
context:
space:
mode:
authorJoseph Myers <joseph@codesourcery.com>2015-09-26 00:27:06 +0000
committerJoseph Myers <joseph@codesourcery.com>2015-09-26 00:27:06 +0000
commitfa752c698146ca3e9f7747d33059fbef9bb02b0e (patch)
treea41295f8a99770d272a53c714effd0f5749b2735 /sysdeps
parent6ace39382120dd0668b2168591611f71b04ac8f3 (diff)
downloadglibc-fa752c698146ca3e9f7747d33059fbef9bb02b0e.tar.gz
glibc-fa752c698146ca3e9f7747d33059fbef9bb02b0e.tar.xz
glibc-fa752c698146ca3e9f7747d33059fbef9bb02b0e.zip
Fix powf inaccuracy (bug 18956).
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.

	[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.
Diffstat (limited to 'sysdeps')
-rw-r--r--sysdeps/ieee754/flt-32/e_powf.c4
-rw-r--r--sysdeps/x86_64/fpu/libm-test-ulps16
2 files changed, 11 insertions, 9 deletions
diff --git a/sysdeps/ieee754/flt-32/e_powf.c b/sysdeps/ieee754/flt-32/e_powf.c
index 18042960c2..c72fe37d3b 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 127a8e1b4d..12c3dd1e01 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