Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
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.
  • Loading branch information
Joseph Myers committed Sep 26, 2015
1 parent 6ace393 commit fa752c6
Show file tree
Hide file tree
Showing 6 changed files with 49 additions and 11 deletions.
10 changes: 10 additions & 0 deletions ChangeLog
View file
@@ -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]
Expand Down
4 changes: 2 additions & 2 deletions NEWS
View file
Expand Up @@ -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.
Expand Down
1 change: 1 addition & 0 deletions math/auto-libm-test-in
View file
Expand Up @@ -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
Expand Down
25 changes: 25 additions & 0 deletions math/auto-libm-test-out
View file
Expand Up @@ -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
Expand Down
4 changes: 3 additions & 1 deletion sysdeps/ieee754/flt-32/e_powf.c
View file
Expand Up @@ -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) */
Expand Down
16 changes: 8 additions & 8 deletions sysdeps/x86_64/fpu/libm-test-ulps
View file
Expand Up @@ -1966,8 +1966,8 @@ Function: "log_vlen8_avx2":
float: 2

Function: "pow":
float: 3
ifloat: 3
float: 1
ifloat: 1
ildouble: 1
ldouble: 1

Expand Down Expand Up @@ -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

Expand Down

0 comments on commit fa752c6

Please sign in to comment.