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

  (1) The bit for the 1.0 part of bp[k] was right shifted by 4.  This
    seems to have been caused by a typo in converting e_pow.c to
    e_powf.c.
(2) The lower 12 bits of ax+bp[k] were not discarded, so t_h was
    actually plain ax+bp[k].  This seems to have been caused by a logic
    error in the conversion.

These bugs gave wrong results like:

    powf(-1.1, 101.0) = -15158.703 (should be -15158.707)
      hex values: BF8CCCCD 42CA0000 C66CDAD0 C66CDAD4

Fixing (1) gives a result wrong in the opposite direction
(hex C66CDAD8), and fixing (2) gives the correct result.

ucbtest has been reporting this particular wrong result on i386 systems
with unpatched libraries for 9 years.  I finally figured out the extent
of the bugs.  On i386's they are normally hidden by extra precision.
We use the trick of representing floats as a sum of 2 floats (one much
smaller) to get extra precision in intermediate calculations without
explicitly using more than float precision.  This trick is just a
pessimization when extra precision is available naturally (as it always
is when dealing with IEEE single precision, so the float precision part
of the library is mostly misimplemented).  (1) and (2) break the trick
in different ways, except on i386's it turns out that the intermediate
calculations are done in enough precision to mask both the bugs and
the limited precision of the float variables (as far as ucbtest can
check).

ucbtest detects the bugs because it forces float precision, but this
is not a normal mode of operation so the bug normally has little effect
on i386's.

On systems that do float arithmetic in float precision, e.g., amd64's,
there is no accidental extra precision and the bugs just give wrong
results.

Reference: https://github.com/freebsd/freebsd-src/commit/12be4e0d5a54a6750913aee2564d164baa71f0dc
Original Author: Bruce Evans
---
 newlib/libm/math/ef_pow.c | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)
  

(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 <keithp@keithp.com>
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 <keithp@keithp.com>
---
 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*/
  

On 26.09.24 04:11, Keith Packard wrote:
> From 3dd0fce84efd8a45e3501175c05086f8a1cbd832 Mon Sep 17 00:00:00 2001
> From: Keith Packard <keithp@keithp.com>
> 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 <keithp@keithp.com>
> ---
>   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