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.

diff --git a/ChangeLog b/ChangeLog
index 951b3f5a4f28b5ee161a0428f7f42d578ef1ecc7..e70cd2fbdc99caa310d3df4fefa37008924fb3d1 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,13 @@
+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.
+
 2015-09-25  Joseph Myers  <joseph@codesourcery.com>
 
        [BZ #18825]

diff --git a/NEWS b/NEWS
index 168e6411201dc243b23769175397ee23bf026f3a..aa9ca4bc9dcd93e6188dd00dd23255bf2d6c2af3 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -15,8 +15,8 @@ Version 2.23
   18084, 18086, 18240, 18265, 18370, 18421, 18480, 18525, 18595, 18610,
   18618, 18647, 18661, 18674, 18675, 18681, 18757, 18778, 18781, 18787,
   18789, 18790, 18795, 18796, 18803, 18820, 18823, 18824, 18825, 18857,
-  18863, 18870, 18872, 18873, 18875, 18887, 18921, 18951, 18952, 18961,
-  18966, 18967, 18970, 18977, 18980, 18981, 19003.
+  18863, 18870, 18872, 18873, 18875, 18887, 18921, 18951, 18952, 18956,
+  18961, 18966, 18967, 18970, 18977, 18980, 18981, 19003.
 
 * The obsolete header <regexp.h> has been removed.  Programs that require
   this header must be updated to use <regex.h> instead.

diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index 4b93857efb36dba4f0c4e700c882708a0cd4ddb2..e86be23fc856686495653c3848d59390ae749eee 100644 (file)
--- 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/math/auto-libm-test-out b/math/auto-libm-test-out
index 4126a7cd07cf72f618118135828d12e78fdd0949..bcfa832b7d69f2f1465860cf9a6e51ee1e4f0db5 100644 (file)
--- a/math/auto-libm-test-out
+++ b/math/auto-libm-test-out
@@ -226489,6 +226489,31 @@ pow 0x9.8b82ap-4 -0x1.99907ap+12
 = pow tonearest ldbl-128ibm 0x9.8b82ap-4L -0x1.99907ap+12L : plus_infty : inexact-ok overflow errno-erange
 = pow towardzero ldbl-128ibm 0x9.8b82ap-4L -0x1.99907ap+12L : 0xf.ffffffffffffbffffffffffffcp+1020L : inexact-ok overflow errno-erange-ok
 = pow upward ldbl-128ibm 0x9.8b82ap-4L -0x1.99907ap+12L : plus_infty : inexact-ok overflow errno-erange
+pow 0xd.73035p-4 -0x1.47bb8p+8
+= pow downward flt-32 0xd.73035p-4f -0x1.47bb8p+8f : 0x4.52398p+80f : inexact-ok
+= pow tonearest flt-32 0xd.73035p-4f -0x1.47bb8p+8f : 0x4.523988p+80f : inexact-ok
+= pow towardzero flt-32 0xd.73035p-4f -0x1.47bb8p+8f : 0x4.52398p+80f : inexact-ok
+= pow upward flt-32 0xd.73035p-4f -0x1.47bb8p+8f : 0x4.523988p+80f : inexact-ok
+= pow downward dbl-64 0xd.73035p-4 -0x1.47bb8p+8 : 0x4.523987c590d3p+80 : inexact-ok
+= pow tonearest dbl-64 0xd.73035p-4 -0x1.47bb8p+8 : 0x4.523987c590d3p+80 : inexact-ok
+= pow towardzero dbl-64 0xd.73035p-4 -0x1.47bb8p+8 : 0x4.523987c590d3p+80 : inexact-ok
+= pow upward dbl-64 0xd.73035p-4 -0x1.47bb8p+8 : 0x4.523987c590d34p+80 : inexact-ok
+= pow downward ldbl-96-intel 0xd.73035p-4L -0x1.47bb8p+8L : 0x4.523987c590d3192p+80L : inexact-ok
+= pow tonearest ldbl-96-intel 0xd.73035p-4L -0x1.47bb8p+8L : 0x4.523987c590d31928p+80L : inexact-ok
+= pow towardzero ldbl-96-intel 0xd.73035p-4L -0x1.47bb8p+8L : 0x4.523987c590d3192p+80L : inexact-ok
+= pow upward ldbl-96-intel 0xd.73035p-4L -0x1.47bb8p+8L : 0x4.523987c590d31928p+80L : inexact-ok
+= pow downward ldbl-96-m68k 0xd.73035p-4L -0x1.47bb8p+8L : 0x4.523987c590d3192p+80L : inexact-ok
+= pow tonearest ldbl-96-m68k 0xd.73035p-4L -0x1.47bb8p+8L : 0x4.523987c590d31928p+80L : inexact-ok
+= pow towardzero ldbl-96-m68k 0xd.73035p-4L -0x1.47bb8p+8L : 0x4.523987c590d3192p+80L : inexact-ok
+= pow upward ldbl-96-m68k 0xd.73035p-4L -0x1.47bb8p+8L : 0x4.523987c590d31928p+80L : inexact-ok
+= pow downward ldbl-128 0xd.73035p-4L -0x1.47bb8p+8L : 0x4.523987c590d3192757b32fb92c7p+80L : inexact-ok
+= pow tonearest ldbl-128 0xd.73035p-4L -0x1.47bb8p+8L : 0x4.523987c590d3192757b32fb92c74p+80L : inexact-ok
+= pow towardzero ldbl-128 0xd.73035p-4L -0x1.47bb8p+8L : 0x4.523987c590d3192757b32fb92c7p+80L : inexact-ok
+= pow upward ldbl-128 0xd.73035p-4L -0x1.47bb8p+8L : 0x4.523987c590d3192757b32fb92c74p+80L : inexact-ok
+= pow downward ldbl-128ibm 0xd.73035p-4L -0x1.47bb8p+8L : 0x4.523987c590d3192757b32fb92cp+80L : inexact-ok
+= pow tonearest ldbl-128ibm 0xd.73035p-4L -0x1.47bb8p+8L : 0x4.523987c590d3192757b32fb92cp+80L : inexact-ok
+= pow towardzero ldbl-128ibm 0xd.73035p-4L -0x1.47bb8p+8L : 0x4.523987c590d3192757b32fb92cp+80L : inexact-ok
+= pow upward ldbl-128ibm 0xd.73035p-4L -0x1.47bb8p+8L : 0x4.523987c590d3192757b32fb92ep+80L : inexact-ok
 sin 0
 = sin downward flt-32 0x0p+0f : 0x0p+0f : inexact-ok
 = sin tonearest flt-32 0x0p+0f : 0x0p+0f : inexact-ok

diff --git a/sysdeps/ieee754/flt-32/e_powf.c b/sysdeps/ieee754/flt-32/e_powf.c
index 18042960c22cd847fbb072662d0fbb4e53362302..c72fe37d3be72953dfa0df88cb1fd08de322199d 100644 (file)
--- 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 127a8e1b4dcb76c325ed73764d89d9426498d7fe..12c3dd1e01c1cfe844a14f0d9768a6814064ba5a 100644 (file)
--- 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