[PATCH 4/4] math: Use rsqrt from CORE-MATH
Adhemerval Zanella
adhemerval.zanella@linaro.org
Tue Sep 22 14:28:38 GMT 2026
The generic implementation computes 1.0 / sqrt (x), which rounds twice
and shows up to 1 ulp for rounding to nearest (for about 26% of the
inputs) and up to 2 ulp for the directed roundings. The CORE-MATH
implementation is correctly rounded in all rounding modes.
The code is adapted to glibc style and to use the math_config.h
definitions, along with:
- The 128-bit integer arithmetic of the refinement uses
math_uint128.h (as s_fma.c)
- The rounding mode is read with get_rounding_mode instead of the
hand-written MXCSR and FPCR readers.
- The check that avoids a spurious underflow in 1/x for x > 2^1022
used 0x7fd000000000000 (2^-896) instead of 0x7fd0000000000000. Both
branches compute the same value in between, so this only saves a
multiplication for most inputs.
Benchtests on x86_64 (Ryzen 9 5900X), aarch64 (Neoverse-N1), powerpc
(POWER10, built for power8), and loongarch64 (Loongson-3C5000L-LL):
* latency
- input: [0.5,2)
master patched improvement
x86_64 30.9951 44.6305 -30.55%
x86_64-v2 30.2757 44.8485 -32.49%
x86_64-v3 29.2565 33.6567 -13.07%
aarch64 11.4766 13.3578 -14.08%
powerpc 9.1164 9.8520 -7.47%
loongarch64 3.2031 3.2529 -1.53%
- input: normal
master patched improvement
x86_64 30.9749 44.6203 -30.58%
x86_64-v2 30.2693 44.8381 -32.49%
x86_64-v3 29.1302 33.6053 -13.32%
aarch64 11.4736 13.3584 -14.11%
powerpc 9.1027 9.8705 -7.78%
loongarch64 3.2027 3.2530 -1.55%
- input: subnormal
master patched improvement
x86_64 42.2356 86.7159 -51.29%
x86_64-v2 40.9051 85.9510 -52.41%
x86_64-v3 39.7846 59.5851 -33.23%
aarch64 11.7711 17.3402 -32.12%
powerpc 9.1276 11.8302 -22.84%
loongarch64 3.2530 4.3537 -25.28%
* reciprocal-throughput
- input: [0.5,2)
master patched improvement
x86_64 10.4331 18.7492 -44.35%
x86_64-v2 10.6881 19.4536 -45.06%
x86_64-v3 10.7814 12.2191 -11.77%
aarch64 5.0839 7.4389 -31.66%
powerpc 1.8437 4.5477 -59.46%
loongarch64 1.8517 2.2027 -15.94%
- input: normal
master patched improvement
x86_64 10.4175 18.8101 -44.62%
x86_64-v2 10.6821 19.4530 -45.09%
x86_64-v3 10.6977 12.2145 -12.42%
aarch64 5.0832 7.4384 -31.66%
powerpc 1.8375 4.5313 -59.45%
loongarch64 1.8517 2.2026 -15.93%
- input: subnormal
master patched improvement
x86_64 21.1591 60.7068 -65.15%
x86_64-v2 21.4036 60.2986 -64.50%
x86_64-v3 19.3992 32.4437 -40.21%
aarch64 5.1698 7.9184 -34.71%
powerpc 1.8438 4.6799 -60.60%
loongarch64 1.9019 2.2028 -13.66%
Checked on aarch64-linux-gnu, x86_64-linux-gnu, powerpc64le-linux-gnu,
and loongarch64-linux-gnuf64.
---
NEWS | 3 +-
SHARED-FILES | 2 +
sysdeps/i386/Makefile | 1 +
sysdeps/ieee754/dbl-64/libm-test-ulps | 12 +++
sysdeps/ieee754/dbl-64/s_rsqrt.c | 144 ++++++++++++++++++++++++++
5 files changed, 161 insertions(+), 1 deletion(-)
create mode 100644 sysdeps/ieee754/dbl-64/s_rsqrt.c
diff --git a/NEWS b/NEWS
index 56c6b581df0..9297ee47e74 100644
--- a/NEWS
+++ b/NEWS
@@ -9,7 +9,8 @@ Version 2.45
Major new features:
- [Add new features here]
+* The rsqrt and rsqrtf functions are now correctly rounded in all rounding
+ modes. The rsqrt implementation is imported from the CORE-MATH project.
Deprecated and removed features, and other changes affecting compatibility:
diff --git a/SHARED-FILES b/SHARED-FILES
index 86ad2283728..611744c3073 100644
--- a/SHARED-FILES
+++ b/SHARED-FILES
@@ -281,6 +281,8 @@ core-math:
sysdeps/ieee754/dbl-64/s_erf.c
# src/binary64/erfc/erfc.c, revision 55e9869e
sysdeps/ieee754/dbl-64/s_erfc.c
+ # src/binary64/rsqrt/rsqrt.c, revision 8ea8ea35
+ sysdeps/ieee754/dbl-64/s_rsqrt.c
# src/binary64/tanh/tanh.c, revision dbd377a9
sysdeps/ieee754/dbl-64/s_tanh.c
# src/binary32/acos/acosf.c, revision 50864ddf
diff --git a/sysdeps/i386/Makefile b/sysdeps/i386/Makefile
index f84dda84018..f60295d2513 100644
--- a/sysdeps/i386/Makefile
+++ b/sysdeps/i386/Makefile
@@ -17,6 +17,7 @@ CFLAGS-s_erfc.c += -fexcess-precision=standard
CFLAGS-s_erf_common.c += -fexcess-precision=standard
CFLAGS-s_fma.c += -fexcess-precision=standard
CFLAGS-s_fmaf.c += -fexcess-precision=standard
+CFLAGS-s_rsqrt.c += -fexcess-precision=standard
CFLAGS-s_tanh.c += -fexcess-precision=standard
endif
diff --git a/sysdeps/ieee754/dbl-64/libm-test-ulps b/sysdeps/ieee754/dbl-64/libm-test-ulps
index c46cfe0030a..cc1de492bd6 100644
--- a/sysdeps/ieee754/dbl-64/libm-test-ulps
+++ b/sysdeps/ieee754/dbl-64/libm-test-ulps
@@ -83,6 +83,18 @@ double: 0
Function: "lgamma_upward":
double: 0
+Function: "rsqrt":
+double: 0
+
+Function: "rsqrt_downward":
+double: 0
+
+Function: "rsqrt_towardzero":
+double: 0
+
+Function: "rsqrt_upward":
+double: 0
+
Function: "sinh":
double: 0
diff --git a/sysdeps/ieee754/dbl-64/s_rsqrt.c b/sysdeps/ieee754/dbl-64/s_rsqrt.c
new file mode 100644
index 00000000000..6b3f683efb3
--- /dev/null
+++ b/sysdeps/ieee754/dbl-64/s_rsqrt.c
@@ -0,0 +1,144 @@
+/* Correctly-rounded reciprocal square root of binary64 value.
+
+Copyright (c) 2022-2025 Alexei Sibidanov.
+
+The original version of this file was copied from the CORE-MATH
+project (file src/binary64/rsqrt/rsqrt.c, revision 8ea8ea35).
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+*/
+
+#include <fenv.h>
+#include <get-rounding-mode.h>
+#include <libm-alias-double.h>
+#include <math.h>
+#include <stdint.h>
+#include "math_config.h"
+#include <math_uint128.h>
+
+static __attribute__ ((noinline)) double
+as_rsqrt_refine (double rf, double a)
+{
+ uint64_t ir = asuint64 (rf), ia = asuint64 (a);
+ if (ia < UINT64_C (1) << MANTISSA_WIDTH)
+ {
+ /* Normalize a subnormal A. Only the parity of the exponent is used
+ below. */
+ int nz = stdc_leading_zeros (ia);
+ ia <<= nz - EXPONENT_WIDTH;
+ ia &= MANTISSA_MASK;
+ ia |= (uint64_t) (nz - 12) << MANTISSA_WIDTH;
+ }
+ /* An even power of 2 gives an exact result. */
+ if ((ia << EXPONENT_WIDTH) == UINT64_C (1) << 63)
+ return rf;
+
+ int mode = get_rounding_mode ();
+ int e = (ia >> MANTISSA_WIDTH) & 1;
+ uint64_t rm = (ir << EXPONENT_WIDTH | UINT64_C (1) << 63) >> EXPONENT_WIDTH;
+ uint64_t am = (get_mantissa (ia) | UINT64_C (1) << MANTISSA_WIDTH)
+ << (5 - e);
+ u128 rt = u128_mul (u128_from_u64 (rm), u128_from_u64 (am));
+ uint64_t rth = u128_high (rt), rtl = u128_low (rt);
+ u128 rrt = u128_mul (u128_from_u64 (rtl), u128_from_u64 (rm));
+ uint64_t t0 = u128_low (rrt), t1 = u128_high (rrt) + rth * rm;
+ rrt = u128_from_hl (t1, t0);
+ int64_t s = u128_high (rrt) >> 63, dd = 1 - 2 * s;
+ /* rts = s ? -(rt << 1) : rt << 1. */
+ u128 rt2 = u128_lshift (rt, 1);
+ uint64_t ms = -(uint64_t) s;
+ u128 rts = u128_from_hl (u128_high (rt2) ^ ms, u128_low (rt2) ^ ms);
+ rts = u128_add (rts, u128_from_u64 (s));
+ u128 prrt;
+ uint64_t am2 = am << 1, am20 = -am;
+ do
+ {
+ ir -= dd;
+ prrt = rrt;
+ am20 += am2;
+ u128 tt = u128_sub (rts, u128_from_u64 (am20));
+ rrt = u128_sub (rrt, tt);
+ }
+ while (__glibc_unlikely (!((u128_high (prrt) ^ u128_high (rrt)) >> 63)));
+ if (!(u128_high (rrt) >> 63))
+ {
+ ir += dd;
+ rrt = prrt;
+ }
+ if (__glibc_likely (mode == FE_TONEAREST))
+ {
+ rm = (ir << EXPONENT_WIDTH | UINT64_C (1) << 63) >> EXPONENT_WIDTH;
+ rt = u128_mul (u128_from_u64 (rm), u128_from_u64 (am));
+ rrt = u128_add (rrt, u128_from_u64 (am >> 2));
+ rrt = u128_add (rrt, rt);
+ ir += u128_high (rrt) >> 63;
+ }
+ else
+ ir += mode == FE_UPWARD;
+ return asdouble (ir);
+}
+
+double
+__rsqrt (double x)
+{
+ uint64_t ix = asuint64 (x);
+ double r;
+ if (__glibc_unlikely (ix < UINT64_C (1) << MANTISSA_WIDTH))
+ {
+ /* 0 <= x < 0x1p-1022. */
+ if (__glibc_unlikely (ix == 0))
+ return __math_divzero (0);
+ r = sqrt (x) / x;
+ }
+ else if (__glibc_unlikely (ix >= EXPONENT_MASK))
+ {
+ /* NaN, Inf, x <= 0. */
+ if (!(ix << 1))
+ return __math_divzero (1); /* x = -0. */
+ if (ix > UINT64_C (0xfff0000000000000))
+ return x + x; /* -NaN. */
+ if (ix >> 63)
+ return __math_invalid (x); /* x < 0. */
+ if (!(ix << 12))
+ return 0.0; /* +Inf. */
+ return x + x; /* +NaN. */
+ }
+ else
+ {
+ /* 0x1p-1022 <= x < 2^1024. */
+ if (__glibc_unlikely (ix > UINT64_C (0x7fd0000000000000)))
+ /* x > 2^1022: avoid a spurious underflow in 1/x. */
+ r = (4.0 / x) * (0.25 * sqrt (x));
+ else
+ r = (1.0 / x) * sqrt (x);
+ }
+ double rx = r * x, drx = fma (r, x, -rx);
+ double h = fma (r, rx, -1.0) + r * drx, dr = (r * 0.5) * h;
+ double rf = r - dr;
+ dr -= r - rf;
+ uint64_t idr = asuint64 (dr), irf = asuint64 (rf);
+ uint64_t aidr = (idr & EXP_MANT_MASK) - (irf & EXPONENT_MASK)
+ + (UINT64_C (0x3fe) << MANTISSA_WIDTH);
+ uint64_t mid = (aidr - UINT64_C (0x3c90000000000000) + 16) >> 5;
+ if (__glibc_unlikely (mid == 0 || aidr < UINT64_C (0x39b0000000000000)
+ || aidr > UINT64_C (0x3c9fffffffffff80)))
+ rf = as_rsqrt_refine (rf, x);
+ return rf;
+}
+libm_alias_double (__rsqrt, rsqrt)
--
2.53.0
More information about the Libc-alpha
mailing list