Re: [PATCH 2/3 v2] powf: Fixed 2 bugs in the computation /* t_h=ax+bp[k] High */. (FreeBSD)

Keith Packard <[email protected]> Wed, 25 Sep 2024 19:11:17 -0700

Newsgroups	gmane.comp.lib.newlib
Message-ID	<[email protected]>

(Sorry for the previous accidental message)

I was running the glibc test cases and found this one:

        powf(0x1.059c76p+0, 0x1.ff80bep+11)

The desired result is 0x1.fffcaap+127, but newlib yields
0x1.fffbf8p+127.

I tracked this down to an imprecise constant 'cp' for 2/(3*ln(2)) used in the
computation.

From 3dd0fce84efd8a45e3501175c05086f8a1cbd832 Mon Sep 17 00:00:00 2001
From: Keith Packard <[email protected]>
Date: Wed, 25 Sep 2024 18:34:28 -0700
Subject: [PATCH] libm: Improve powf accuracy

One of the constants, 'cp', was not correctly converted to extended
precision (using two floats) leading to reduced precision in the
computation of log2(x).

The new values use a larger value for cp_h which is closer to the
actual value.  That leaves a smaller magnitude for cp_l resulting in
increased precision for the result. These values were computed
according to the Veltkamp/Dekker method.

2/(3*ln2): 0x1.ec709dc3a03fd748p-1
Old cp_h:  0x1.ec7p-1
New cp_h:  0x1.ec8p-1

The desired value is much closer to the higher value than the lower.

Signed-off-by: Keith Packard <[email protected]>
---
 newlib/libm/math/ef_pow.c | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/newlib/libm/math/ef_pow.c b/newlib/libm/math/ef_pow.c
index a2aba1054..d69f1ae48 100644
--- a/newlib/libm/math/ef_pow.c
+++ b/newlib/libm/math/ef_pow.c
@@ -47,8 +47,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.61914062e-01,   /* 0x3f764000 =head of cp */
+cp_l  = -1.17368574e-04,   /* 0xb8f623c6 =tail of cp */
 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*/
-- 
2.45.2

-- 
-keith

signature.asc (application/pgp-signature, 832 B)