diff mbox series

[3/3] powf: Fixed another precision bug in powf() (FreeBSD)

Message ID 20240911140646.2143-4-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 
  Fixed another precision bug in powf(). This one is in the computation
[t=p_l+p_h High]. We multiply t by lg2_h, and want the result to be
exact. For the bogus float case of the high-low decomposition trick, we
normally discard the lowest 12 bits of the fraction for the high part,
keeping 12 bits of precision. That was used for t here, but it doesnt't
work because for some reason we only discard the lowest 9 bits in the
fraction for lg2_h.  Discard another 3 bits of the fraction for t to
compensate.

This bug gave wrong results like:

powf(0.9999999, -2.9999995) = 1.0000002 (should be 1.0000001)
hex values: 3F7FFFFF C03FFFFE 3F800002 3F800001

As explained in the log for the previous commit, the bug is normally
masked by doing float calculations in extra precision on i386's, but is
easily detected by ucbtest on systems that don't have accidental extra
precision.

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

Patch

diff --git a/newlib/libm/math/ef_pow.c b/newlib/libm/math/ef_pow.c
index c7b1975c3..8cf58dbb5 100644
--- a/newlib/libm/math/ef_pow.c
+++ b/newlib/libm/math/ef_pow.c
@@ -242,7 +242,7 @@  ivln2_l  =  7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/
 	} 
 	t = p_l+p_h;
 	GET_FLOAT_WORD(is,t);
-	SET_FLOAT_WORD(t,is&0xfffff000);
+	SET_FLOAT_WORD(t,is&0xffff8000);
 	u = t*lg2_h;
 	v = (p_l-(t-p_h))*lg2+t*lg2_l;
 	z = u+v;