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

On 5/7/26 22:55, Paul Zimmermann wrote:
>         Hi Steve,
> 
> I'll try to exercise your code when I return from vacation.

Thanks, I would appreciate someone else testing the code.

>
>> 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?

I haven't completely identified the cause of the 127
values that yield a not correctly rounded value.  My
testing framework is not set up to record each and
every value that does not yield 0.5 ULP.  If I
reduce the interval for exhaustive testing and change
the reported output I, for example, see

% ./tlibm rsqrt -fPE -x 8 -X 17
Interval tested for rsqrtf: [8,17]
       ulp <= 0.5: 100.000% 8912895 |  100.000% 8912895
0.5 < ulp <= 0.6:  0.000%       1 | 100.000% 8912896
Max ulp: 5.00000060e-01 at 1.59999981e+01, 0x1.fffffcp+3

Here, the rogue value is nearly 16, i.e., a power of 2.
The reference result is from mpfr_rec_sqrt() with the
precision of MPFR set to 4*24=96 bits.  I can increase
the precision, but the result does not change.

% ./tlibm rsqrt -fa 0x1.fffffcp3 -p 96
    x =  1.59999981e+01f, /* 0x417ffffe */
libm =  2.50000000e-01f, /* 0x3e800000 */
mpfr =  2.50000030e-01f, /* 0x3e800001 */
  ULP = 0.50000006
% ./tlibm rsqrt -fa 0x1.fffffcp3 -p 256
    x =  1.59999981e+01f, /* 0x417ffffe */
libm =  2.50000000e-01f, /* 0x3e800000 */
mpfr =  2.50000030e-01f, /* 0x3e800001 */
  ULP = 0.50000006

Note, I use strtof() to convert 0x1.fffffcp3
to the float value, and the reported mpfr
value is from use of mpfr_get_flt().  As x
is slightly less than 16, I'm incline to
accept MPFR as correct.

>> 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?

If I knew how, yes. :-).

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

The testing framework was setup to look at functions with
ULP's that exceed 0.5 by tenths.  For me, the 5th decimal
place is cosmetic.  For example,

% ./tlibm j0 -fPED -x 1 -X 2
Interval tested for j0f: [1,2]
       ulp <= 0.5: 64.738% 5430646 |  64.738% 5430646
0.5 < ulp <= 0.6:  8.870%  744080 |  73.608% 6174726
0.6 < ulp <= 0.7:  7.148%  599657 |  80.757% 6774383
0.7 < ulp <= 0.8:  5.760%  483165 |  86.517% 7257548
0.8 < ulp <= 0.9:  4.488%  376474 |  91.005% 7634022
0.9 < ulp <= 1.0:  3.110%  260863 |  94.114% 7894885
1.0 < ulp <= 1.5:  5.271%  442161 |  99.385% 8337046
1.5 < ulp <= 2.0:  0.481%   40354 |  99.866% 8377400
2.0 < ulp <= 3.0:  0.132%   11081 |  99.998% 8388481
3.0 < ulp       :  0.002%     127 | 100.000% 8388608
Max ulp: 3.66519642 at 1.99708652e+00, 0x1.ff4110p+0

-- 
steve