diff mbox series

[1/3] powf: Fix the hi+lo decomposition for 2/(3ln2) (FreeBSD)

Message ID 20240911140646.2143-2-fabian.schriever@gtd-gmbh.de
State New
Headers
Series Fix powf inaccuracies up to ~169 ULP reported by Paul Zimmermann |

Commit Message

Fabian Schriever Sept. 11, 2024, 2:06 p.m. UTC 
  The decomposition needs to be into 12+24 bits of precision for extra-
precision multiplication, but was into 13+24 bits. On i386 with -O1 the
bug was hidden by accidental extra precision, but on amd64, in 2^32
trials the bug caused about 200000 errors of more than 1 ulp, with a
maximum error of about 80 ulps. Now the maximum error in 2^32 trials
on amd64 is 0.8573 ulps. It is still 0.8316 ulps on i386 with -O1.

The nearby decomposition of 1/ln2 and the decomposition of 2/(3ln2) in
the double precision version seem to be sub-optimal but not broken.

Reference: https://github.com/freebsd/freebsd-src/commit/b4437c3d322a0f6d23d12b6f76d2fc72d2ff0ec2
Original Author: Bruce Evans
---
 newlib/libm/math/ef_pow.c | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)
diff mbox series

Patch

diff --git a/newlib/libm/math/ef_pow.c b/newlib/libm/math/ef_pow.c
index 07b225f8c..4ceba3d24 100644
--- a/newlib/libm/math/ef_pow.c
+++ b/newlib/libm/math/ef_pow.c
@@ -50,8 +50,8 @@  lg2_h  =  6.93145752e-01, /* 0x3f317200 */
 lg2_l  =  1.42860654e-06, /* 0x35bfbe8c */
 ovt =  4.2995665694e-08, /* -(128-log2(ovfl+.5ulp)) */
 cp    =  9.6179670095e-01, /* 0x3f76384f =2/(3ln2) */
-cp_h  =  9.6179199219e-01, /* 0x3f763800 =head of cp */
-cp_l  =  4.7017383622e-06, /* 0x369dc3a0 =tail of cp_h */
+cp_h  =  9.6191406250e-01f, /* 0x3f764000 =12b of cp */
+cp_l  = -1.1736857402e-04f, /* 0xb8f623c6 =tail of cp_h */
 ivln2    =  1.4426950216e+00, /* 0x3fb8aa3b =1/ln2 */
 ivln2_h  =  1.4426879883e+00, /* 0x3fb8aa00 =16b 1/ln2*/
 ivln2_l  =  7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/