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

Thomas Wucher <[email protected]>
Newsgroups gmane.comp.lib.newlib
Message-ID <[email protected]>
On 26.09.24 04:11, Keith Packard wrote:
> 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*/

This looks correct, thanks Keith! A very similar change has also been 
done in FreeBSD:

https://github.com/freebsd/freebsd-src/commit/b4437c3d322a0f6d23d12b6f76d2fc72d2ff0ec2
lmpx.com only provides a reader for public news (NNTP) servers. It is not affiliated with the servers or forums shown here and is not responsible for the content of articles, which is written by their respective authors.