Created attachment 270500 [details] Patch with libm changes The attached diff implements the inverse square root function, i.e, rsqrt(x) = 1 / sqrt(x). Exhaustive testing of the float version suggests that it is correctly rounded in round-to-nearest for all test values in the range [0x1p-127,0x1p126] except for 127 values. Exhaustive testing of rsqrt and rsqrtl cannot be done, but 1100M values of x for rsqrt and 400M values for rsqrtl were tested. All tested values were correctly rounded. I do not have access to LD128 (i.e., IEEE 128-bit floating point) hardware, so the implementation of rsqrtl() is untested. The following is a summary of changes to source code. * lib/msun/Makefile: . Add s_rsqrt.c and s_rsqrtf.c to COMMON_SRCS. . For non-53-bit long double targets, add s_rsqrtl.c to COMMON_SRCS. . Add MLINKS for rsqrt.3, rsqrtf.3, and rsqrtl.3 to sqrt.3. * lib/msun/Symbol.map: . Add rsqrt, rsqrtf, and rsqrtl to the Symbol map for shared libm.so. * lib/msun/man/sqrt.3: . Update the sqrt.3 manual page to include information for rsqrt[fl]. . Note, these function come from ISO C23 (and IEEE-754 2008). * lib/msun/src/math.h: . Add prototypes for new functions. * lib/msun/src/math_private.h: . Add _SPLIT, _FAST2SUM, _SLOW2SUM, _XADD, _MUL, and _XMUL macros to perform type-type arthimetic (i.e., float-float). * src/s_rsqrt.c: . New file with the implementation of 'double rsqrt(double)'. . For 53-bit long double targets, add a weak reference for rsqrtl. * src/s_rsqrtf.c: . New file with the implementation of 'float rsqrt(float)'. * src/s_rsqrtl.c . New file with the implementation of 'long double rsqrt(long double)'. Note, the LD80 version uses bit twiddling and LD128 version is a straight C language implementation. The LD128 is untested due to lack of hardware.
Created attachment 270556 [details] New diff The new diff corrects two (silly) issues found by kib@. It also includes fixes for correct rounding of rsqrt and rsqrtf (in at least round-to-nearest mode).
A commit in branch main references this bug: URL: https://cgit.FreeBSD.org/src/commit/?id=3085fc9d97bd83785ba3ba43e0378d7d67987d1f commit 3085fc9d97bd83785ba3ba43e0378d7d67987d1f Author: Steve Kargl <kargl@FreeBSD.org> AuthorDate: 2026-05-08 14:06:08 +0000 Commit: Konstantin Belousov <kib@FreeBSD.org> CommitDate: 2026-05-10 16:36:33 +0000 [libm] implementation of rsqrt, rsqrtf, and rsqrtl From the PR: The attached diff implements the inverse square root function, i.e, rsqrt(x) = 1 / sqrt(x). Exhaustive testing of the float version suggests that it is correctly rounded in round-to-nearest for all test values in the range [0x1p-127,0x1p126]. Exhaustive testing of rsqrt and rsqrtl cannot be done, but 1100M values of x for rsqrt and 400M values for rsqrtl were tested. All tested values were correctly rounded. I do not have access to LD128 (i.e., IEEE 128-bit floating point) hardware, so the implementation of rsqrtl() is untested. The following is a summary of changes to source code. * lib/msun/Makefile: . Add s_rsqrt.c and s_rsqrtf.c to COMMON_SRCS. . For non-53-bit long double targets, add s_rsqrtl.c to COMMON_SRCS. . Add MLINKS for rsqrt.3, rsqrtf.3, and rsqrtl.3 to sqrt.3. * lib/msun/Symbol.map: . Add rsqrt, rsqrtf, and rsqrtl to the Symbol map for shared libm.so. * lib/msun/man/sqrt.3: . Update the sqrt.3 manual page to include information for rsqrt[fl]. . Note, these function come from ISO C23 (and IEEE-754 2008). * lib/msun/src/math.h: . Add prototypes for new functions. * lib/msun/src/math_private.h: . Add _SPLIT, _FAST2SUM, _SLOW2SUM, _XADD, _MUL, and _XMUL macros to perform type-type arthimetic (i.e., float-float). * src/s_rsqrt.c: . New file with the implementation of 'double rsqrt(double)'. . For 53-bit long double targets, add a weak reference for rsqrtl. * src/s_rsqrtf.c: . New file with the implementation of 'float rsqrt(float)'. * src/s_rsqrtl.c . New file with the implementation of 'long double rsqrt(long double)'. Note, the LD80 version uses bit twiddling and LD128 version is a straight C language implementation. The LD128 is untested due to lack of hardware. PR: 295089 MFC after: 1 week lib/msun/Makefile | 6 +- lib/msun/Symbol.map | 3 + lib/msun/man/sqrt.3 | 53 ++++++++++- lib/msun/src/math.h | 3 + lib/msun/src/math_private.h | 83 +++++++++++++++++ lib/msun/src/s_rsqrt.c (new) | 153 +++++++++++++++++++++++++++++++ lib/msun/src/s_rsqrtf.c (new) | 155 ++++++++++++++++++++++++++++++++ lib/msun/src/s_rsqrtl.c (new) | 203 ++++++++++++++++++++++++++++++++++++++++++ 8 files changed, 654 insertions(+), 5 deletions(-)
Oh, didn't realise you had filed a PR! Great to see it merged!
A commit in branch stable/15 references this bug: URL: https://cgit.FreeBSD.org/src/commit/?id=0093bb670537f19e99fff2d46b4831d7b9d44b4c commit 0093bb670537f19e99fff2d46b4831d7b9d44b4c Author: Steve Kargl <kargl@FreeBSD.org> AuthorDate: 2026-05-08 14:06:08 +0000 Commit: Konstantin Belousov <kib@FreeBSD.org> CommitDate: 2026-05-17 00:27:48 +0000 [libm] implementation of rsqrt, rsqrtf, and rsqrtl PR: 295089 (cherry picked from commit 3085fc9d97bd83785ba3ba43e0378d7d67987d1f) lib/msun/Makefile | 6 +- lib/msun/Symbol.map | 3 + lib/msun/man/sqrt.3 | 53 ++++++++++- lib/msun/src/math.h | 3 + lib/msun/src/math_private.h | 83 +++++++++++++++++ lib/msun/src/s_rsqrt.c (new) | 153 +++++++++++++++++++++++++++++++ lib/msun/src/s_rsqrtf.c (new) | 155 ++++++++++++++++++++++++++++++++ lib/msun/src/s_rsqrtl.c (new) | 203 ++++++++++++++++++++++++++++++++++++++++++ 8 files changed, 654 insertions(+), 5 deletions(-)
with attachment 270556 [details], I still get incorrect rounding for rsqrtf: ``` libm wrong by up to 5.01e-01 ulp(s) [1] for x=0x1.13e07p-129 rsqrt gives 0x1.5cc0a8p+64 mpfr_rsqrt gives 0x1.5cc0aap+64 Total: errors=129 (0.00%) errors2=0 maxerr=5.01e-01 ulp(s) ``` (I copied the auxiliary functions from https://github.com/freebsd/freebsd-src.git.)
I also get errors with the double version: ``` zimmerma@thym:~/svn/core-math$ CORE_MATH_CHECK_STD=true LIBM=/tmp/libm.a ./check.sh rsqrt Running worst cases check in --rndn mode... FAIL x=0x1.e60431782ff31p-1 ref=0x1.06c10c599fca7p+0 z=0x0.000002b344e03p-1022 ``` but maybe I did a mistake in copying the code. How can I build a libm.a with just rsqrtf, rsqrt and rsqrtl?