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

On 2024-09-11 16:06, Fabian Schriever wrote:
> (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(-)
> 
> diff --git a/newlib/libm/math/ef_pow.c b/newlib/libm/math/ef_pow.c
> index 4ceba3d24..c7b1975c3 100644
> --- a/newlib/libm/math/ef_pow.c
> +++ b/newlib/libm/math/ef_pow.c
> @@ -173,7 +173,8 @@ ivln2_l  =  7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/
>   	    GET_FLOAT_WORD(is,s_h);
>   	    SET_FLOAT_WORD(s_h,is&0xfffff000);
>   	/* t_h=ax+bp[k] High */
> -	    SET_FLOAT_WORD(t_h,((ix>>1)|0x20000000)+0x0040000+(k<<21));
> +        is = ((ix >> 1) & 0xfffff000U) | 0x20000000;
> +        SET_FLOAT_WORD(t_h, is + 0x00400000 + (k << 21));

The indentation of the above 2 added lines appears to be incorrect.

Kind regards,
Torbjörn

>   	    t_l = ax - (t_h-bp[k]);
>   	    s_l = v*((u-s_h*t_h)-s_h*t_l);
>   	/* compute log(ax) */