Re: [libm] Implementations for rsqrt[fl].
Steve Kargl <[email protected]>
| Newsgroups | gmane.os.freebsd.devel.hackers |
|---|---|
| Message-ID | <[email protected]> |
On 5/7/26 22:55, Paul Zimmermann wrote: > 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? > Don't know how to fix the issue, but the 127 are of the form: 0x1.fffffcp-127 ULP = 0.50000006 0x1.fffffcp-126 ULP = 0.08986171 0x1.fffffcp-125 ULP = 0.50000006 0x1.fffffcp-124 ULP = 0.08986171 0x1.fffffcp-123 ULP = 0.50000006 0x1.fffffcp-122 ULP = 0.08986171 0x1.fffffcp-121 ULP = 0.50000006 0x1.fffffcp-120 ULP = 0.08986171 After splitting x into significand and exponent: x = f*2^m, if I have an odd exponent I scale via f /= 2, m+=1. I played with various rounding modes but that did not help. -- steve