Re: [libm] Implementations for rsqrt[fl].

       Hi Steve,

I'll try to exercise your code when I return from vacation.

> Date: Thu, 7 May 2026 16:05:16 -0700
> From: Steve Kargl <[email protected]>
> 
> All,
> 
> I submitted https://bugs.freebsd.org/bugzilla/show_bug.cgi?id=295089
> which contains an implementation of rsqrt[fl](x)=1/sqrt[fl](x).  These
> appear in C23 7.12.7.9 and F.10.4.9 as well as IEEE-754 2008.  Under
> IEEE-754, the functions are required(?)/recommended to be correctly
> rounded in all rounding modes.  I have only checked round-to-nearest.
> Using my testing franework, I see
> 
> % ./tlibm rsqrt -fPED -x 0x1p-127 -X 0x1p126 -s 0
> Interval tested for rsqrtf: [5.87747e-39,8.50706e+37]
> 1000000 calls, 0.011062 secs, 0.01106 usecs/call
>         ulp <= 0.5: 100.000% 2118123393 | 100.000% 2118123393
> 0.5 <  ulp <  0.6:  0.000%     127 | 100.000% 2118123520

what does this 127 mean?

> Max ulp: 0.500000 at 1.17549421e-38
> %./tlibm rsqrt -fa  1.17549421e-38
>     x =  1.17549421e-38f, /* 0x007fffff */
> libm =  9.22337204e+18f, /* 0x5f000000 */
> mpfr =  9.22337314e+18f, /* 0x5f000001 */
>   ULP = 0.50000

it seems libm is not correctly rounded in this example.
Do you plan to fix that?

Side remark: the ULP value should be rounded up to 0.50001.

Paul

> %./tlibm rsqrt -dP -x 0 -X 0x1p1001 -N700 -s 0
> Interval tested for rsqrt: [0,2.14302e+301]
> 700000000 calls, 9.325613 secs, 0.01332 usecs/call
>         ulp <= 0.5: 100.000% 699999999 | 100.000% 699999999
> Max ulp: 0.500000 at 1.2611313594916888e+301
> 
> % ./tlibm rsqrt -lP -x 0 -X 0x1p16380 -N700 -s 0
> Interval tested for rsqrtl: [0,7.43582e+4930]
> 700000000 calls, 53.109852 secs, 0.07587 usecs/call
>         ulp <= 0.5: 100.000% 699999999 | 100.000% 699999999
> Max ulp: 0.500000 at 6.157706927493729414081e+4930
> 
> If a committer wants to pick up the diff and apply it
> to src/  I can answer any questions you may have.
> -- 
> steve
> 
>