[libc-commits] [libc] [llvm] [libc][math][c23] Improve rsqrtf16 function for targets without fp32 FPUs. (PR #160639)
Anton Shepelev via libc-commits
libc-commits at lists.llvm.org
Wed Jun 10 23:25:45 PDT 2026
https://github.com/amemov updated https://github.com/llvm/llvm-project/pull/160639
>From 5bc74a6af2720225a3152d9630b5b1ec44b90a6f Mon Sep 17 00:00:00 2001
From: Anton Shepelev <shepelev777 at gmail.com>
Date: Wed, 24 Sep 2025 20:59:06 -0700
Subject: [PATCH 01/10] - Draft of math approximation for targets that have
only int hardware - Fixed typo in +inf case. Should return +0 according to
F.10.4.9
---
libc/src/__support/math/rsqrtf16.h | 109 +++++++++++++++++++++++++++--
1 file changed, 105 insertions(+), 4 deletions(-)
diff --git a/libc/src/__support/math/rsqrtf16.h b/libc/src/__support/math/rsqrtf16.h
index ab7529682950d..551aa9710ea70 100644
--- a/libc/src/__support/math/rsqrtf16.h
+++ b/libc/src/__support/math/rsqrtf16.h
@@ -58,12 +58,11 @@ LIBC_INLINE constexpr float16 rsqrtf16(float16 x) {
return FPBits::quiet_nan().get_val();
}
- // x = +inf => rsqrt(x) = 0
- return FPBits::zero().get_val();
+ // x = +inf => rsqrt(x) = +0
+ return FPBits::zero(xbits.sign()).get_val();
}
- // TODO: add integer based implementation when LIBC_TARGET_CPU_HAS_FPU_FLOAT
- // is not defined
+#ifdef LIBC_TARGET_CPU_HAS_FPU_FLOAT
float result = 1.0f / fputil::sqrt<float>(fputil::cast<float>(x));
// Targeted post-corrections to ensure correct rounding in half for specific
@@ -76,6 +75,108 @@ LIBC_INLINE constexpr float16 rsqrtf16(float16 x) {
}
return fputil::cast<float16>(result);
+
+#else
+ // Range reduction:
+ // x can be expressed as m*2^e, where e - int exponent and m - mantissa
+ // rsqrtf16(x) = rsqrtf16(m*2^e)
+ // rsqrtf16(m*2^e) = 1/sqrt(m) * 1/sqrt(2^e) = 1/sqrt(m) * 1/2^(e/2)
+ // 1/sqrt(m) * 1/2^(e/2) = 1/sqrt(m) * 2^(-e/2)
+
+ // Compute reduction directly from half bits to avoid frexp/ldexp overhead.
+ int exponent = 0;
+ int signifcand = 0; // same as mantissa, but int
+ uint16_t eh = static_cast<uint16_t>((x_abs >> 10) & 0x1F);
+ uint16_t frac = static_cast<uint16_t>(x_abs & 0x3FF);
+
+ int result;
+ if (eh != 0) {
+ // ((2^-1 + frac/2^11) * 2) * 2^(eh-15)
+
+ // Normal: x = (1 + frac/2^10) * 2^(eh-15) = ((0.5 + frac/2^11) * 2) *
+ // 2^(eh-15)
+ // => mantissa in [0.5,1): m = 0.5 + frac/2^11, exponent = (eh - 15) + 1 =
+ // eh - 14
+ exponent = static_cast<int>(eh) - 14;
+ mantissa = 0.5f + static_cast<float>(frac) * 0x1.0p-11f;
+ } else {
+ // Subnormal: x = (frac/2^10) * 2^(1-15) = frac * 2^-24.
+ // Normalize frac so that bit 9 becomes 1; then mantissa m = (frac <<
+ // t)/2^10 ∈ [0.5,1) and exponent E = -14 - t so that x = m * 2^E.
+ if (LIBC_UNLIKELY(frac == 0)) {
+ // Should have been handled by zero check above, but keep safe.
+ return FPBits::inf(Sign::POS).get_val();
+ }
+ int shifts = 0;
+ while ((frac & 0x200u) == 0u) { // bring into [0x200, 0x3FF]
+ frac <<= 1;
+ ++shifts;
+ }
+ exponent = -14 - shifts;
+ mantissa = static_cast<float>(frac) * 0x1.0p-10f;
+ }
+
+ float result = 0.0f;
+ int exp_floored = -(exponent >> 1);
+
+ if (mantissa == 0.5f) {
+ // When mantissa is 0.5f, x was a power of 2 (or subnormal that normalizes
+ // this way). 1/sqrt(0.5f) = sqrt(2.0f).
+ // If exponent is odd (exponent = 2k + 1):
+ // rsqrt(x) = (1/sqrt(0.5)) * 2^(-(2k+1)/2) = sqrt(2) * 2^(-k-0.5)
+ // = sqrt(2) * 2^(-k) * (1/sqrt(2)) = 2^(-k)
+ // exp_floored = -((2k+1)>>1) = -(k) = -k
+ // So result = ldexp(1.0f, exp_floored)
+ // If exponent is even (exponent = 2k):
+ // rsqrt(x) = (1/sqrt(0.5)) * 2^(-2k/2) = sqrt(2) * 2^(-k)
+ // exp_floored = -((2k)>>1) = -(k) = -k
+ // So result = ldexp(sqrt(2.0f), exp_floored)
+ if (exponent & 1) {
+ result = fputil::ldexp(1.0f, exp_floored);
+ } else {
+ constexpr float SQRT_2_F = 0x1.6a09e6p0f; // sqrt(2.0f)
+ result = fputil::ldexp(SQRT_2_F, exp_floored);
+ }
+ } else {
+ // 4 Degree minimax polynomial (single-precision coefficients) generated
+ // with Sollya: P = fpminimax(1/sqrt(x), 4,
+ // [|single,single,single,single,single|], [0.5;1])
+ float y = fputil::polyeval(mantissa,
+ 0x1.771256p1f, // c0
+ -0x1.5e7c4ap2f, // c1
+ 0x1.b3851cp2f, // c2
+ -0x1.1a27ep2f, // c3
+ 0x1.265c66p0f); // c4
+
+ // Newton-Raphson iteration in float (use multiply_add to leverage FMA when
+ // available):
+ float y2 = y * y;
+ float factor = fputil::multiply_add(-0.5f * mantissa, y2, 1.5f);
+ y = y * factor;
+
+ result = fputil::ldexp(y, exp_floored);
+ if (exponent & 1) {
+ constexpr float ONE_OVER_SQRT2 = 0x1.6a09e6p-1f; // 1/sqrt(2)
+ result *= ONE_OVER_SQRT2;
+ }
+
+ // Targeted post-correction: for the specific half-precision mantissa
+ // pattern M == 0x011F we observe a consistent -1 ULP bias across exponents.
+ // Apply a tiny upward nudge to cross the rounding boundary in all modes.
+ const uint16_t half_mantissa = static_cast<uint16_t>(x_abs & 0x3ff);
+ if (half_mantissa == 0x011F) {
+ // Nudge up to fix consistent -1 ULP at that mantissa boundary
+ result = fputil::multiply_add(result, 0x1.0p-21f,
+ result); // result *= (1 + 2^-21)
+ } else if (half_mantissa == 0x0313) {
+ // Nudge down to fix +1 ULP under upward rounding at this mantissa
+ // boundary
+ result = fputil::multiply_add(result, -0x1.0p-21f,
+ result); // result *= (1 - 2^-21)
+ }
+ }
+ return fputil::cast<float16>(result);
+#endif
}
} // namespace math
>From d80b51ed0dcc5305526c57dbb61f47740a9c45ff Mon Sep 17 00:00:00 2001
From: amemov <shepelev777 at gmail.com>
Date: Mon, 17 Nov 2025 22:00:28 -0800
Subject: [PATCH 02/10] Added approximation from previous PR
---
libc/src/__support/math/rsqrtf16.h | 45 +++++-------------------------
1 file changed, 7 insertions(+), 38 deletions(-)
diff --git a/libc/src/__support/math/rsqrtf16.h b/libc/src/__support/math/rsqrtf16.h
index 551aa9710ea70..4b237950b0b17 100644
--- a/libc/src/__support/math/rsqrtf16.h
+++ b/libc/src/__support/math/rsqrtf16.h
@@ -16,6 +16,7 @@
#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/ManipulationFunctions.h"
+#include "src/__support/FPUtil/PolyEval.h"
#include "src/__support/FPUtil/cast.h"
#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/FPUtil/sqrt.h"
@@ -77,44 +78,10 @@ LIBC_INLINE constexpr float16 rsqrtf16(float16 x) {
return fputil::cast<float16>(result);
#else
- // Range reduction:
- // x can be expressed as m*2^e, where e - int exponent and m - mantissa
- // rsqrtf16(x) = rsqrtf16(m*2^e)
- // rsqrtf16(m*2^e) = 1/sqrt(m) * 1/sqrt(2^e) = 1/sqrt(m) * 1/2^(e/2)
- // 1/sqrt(m) * 1/2^(e/2) = 1/sqrt(m) * 2^(-e/2)
+ float xf = fputil::cast<float>(x);
- // Compute reduction directly from half bits to avoid frexp/ldexp overhead.
int exponent = 0;
- int signifcand = 0; // same as mantissa, but int
- uint16_t eh = static_cast<uint16_t>((x_abs >> 10) & 0x1F);
- uint16_t frac = static_cast<uint16_t>(x_abs & 0x3FF);
-
- int result;
- if (eh != 0) {
- // ((2^-1 + frac/2^11) * 2) * 2^(eh-15)
-
- // Normal: x = (1 + frac/2^10) * 2^(eh-15) = ((0.5 + frac/2^11) * 2) *
- // 2^(eh-15)
- // => mantissa in [0.5,1): m = 0.5 + frac/2^11, exponent = (eh - 15) + 1 =
- // eh - 14
- exponent = static_cast<int>(eh) - 14;
- mantissa = 0.5f + static_cast<float>(frac) * 0x1.0p-11f;
- } else {
- // Subnormal: x = (frac/2^10) * 2^(1-15) = frac * 2^-24.
- // Normalize frac so that bit 9 becomes 1; then mantissa m = (frac <<
- // t)/2^10 ∈ [0.5,1) and exponent E = -14 - t so that x = m * 2^E.
- if (LIBC_UNLIKELY(frac == 0)) {
- // Should have been handled by zero check above, but keep safe.
- return FPBits::inf(Sign::POS).get_val();
- }
- int shifts = 0;
- while ((frac & 0x200u) == 0u) { // bring into [0x200, 0x3FF]
- frac <<= 1;
- ++shifts;
- }
- exponent = -14 - shifts;
- mantissa = static_cast<float>(frac) * 0x1.0p-10f;
- }
+ float mantissa = fputil::frexp(xf, exponent);
float result = 0.0f;
int exp_floored = -(exponent >> 1);
@@ -139,8 +106,9 @@ LIBC_INLINE constexpr float16 rsqrtf16(float16 x) {
}
} else {
// 4 Degree minimax polynomial (single-precision coefficients) generated
- // with Sollya: P = fpminimax(1/sqrt(x), 4,
- // [|single,single,single,single,single|], [0.5;1])
+ // with Sollya:
+ // P = fpminimax(1/sqrt(x), 4,
+ // [|single,single,single,single,single|], [0.5;1])
float y = fputil::polyeval(mantissa,
0x1.771256p1f, // c0
-0x1.5e7c4ap2f, // c1
@@ -175,6 +143,7 @@ LIBC_INLINE constexpr float16 rsqrtf16(float16 x) {
result); // result *= (1 - 2^-21)
}
}
+
return fputil::cast<float16>(result);
#endif
}
>From 2a365b3db82877adbf8d72ca3796f1827c45ad7e Mon Sep 17 00:00:00 2001
From: amemov <shepelev777 at gmail.com>
Date: Sat, 25 Apr 2026 13:37:01 -0700
Subject: [PATCH 03/10] new approximation draft
---
libc/src/__support/math/CMakeLists.txt | 14 ++
libc/src/__support/math/rsqrtf16.h | 237 +++++++++++++-----
.../llvm-project-overlay/libc/BUILD.bazel | 1 -
3 files changed, 182 insertions(+), 70 deletions(-)
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index b6b397c88a256..87c925b4e0f8a 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -126,6 +126,20 @@ add_header_library(
libc.src.__support.macros.properties.types
)
+add_header_library(
+ rsqrtf16
+ HDRS
+ rsqrtf16.h
+ DEPENDS
+ libc.src.__support.FPUtil.cast
+ libc.src.__support.FPUtil.fenv_impl
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.multiply_add
+ libc.src.__support.FPUtil.sqrt
+ libc.src.__support.macros.optimization
+ libc.src.__support.macros.properties.types
+)
+
add_header_library(
asin_utils
HDRS
diff --git a/libc/src/__support/math/rsqrtf16.h b/libc/src/__support/math/rsqrtf16.h
index 4b237950b0b17..23c74a8515c15 100644
--- a/libc/src/__support/math/rsqrtf16.h
+++ b/libc/src/__support/math/rsqrtf16.h
@@ -15,8 +15,6 @@
#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FPBits.h"
-#include "src/__support/FPUtil/ManipulationFunctions.h"
-#include "src/__support/FPUtil/PolyEval.h"
#include "src/__support/FPUtil/cast.h"
#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/FPUtil/sqrt.h"
@@ -25,6 +23,173 @@
namespace LIBC_NAMESPACE_DECL {
namespace math {
+namespace rsqrtf16_internal {
+
+LIBC_INLINE_VAR constexpr int RSQRT_FRACTION_BITS = 29;
+LIBC_INLINE_VAR constexpr int64_t ONE = int64_t(1) << RSQRT_FRACTION_BITS;
+LIBC_INLINE_VAR constexpr int64_t THREE_HALVES = 3 * (ONE >> 1);
+
+// Degree-4 minimax polynomial generated with Sollya:
+// P = fpminimax(1/sqrt(x), 4,
+// [|single,single,single,single,single|], [0.5;1])
+// Coefficients are stored in Q29 fixed-point format.
+LIBC_INLINE_VAR constexpr int64_t COEFFS[5] = {
+ 1'573'164'416, -2'940'085'504, 3'653'406'208, -2'366'894'080, 617'319'616,
+};
+LIBC_INLINE_VAR constexpr int64_t ONE_OVER_SQRT2 = 0x16a09e60;
+
+LIBC_INLINE constexpr int floor_log2(uint64_t x) {
+ int result = -1;
+ while (x) {
+ x >>= 1;
+ ++result;
+ }
+ return result;
+}
+
+LIBC_INLINE constexpr int64_t eval_polynomial(uint32_t m) {
+ int64_t y = COEFFS[4];
+ y = COEFFS[3] + ((y * m) >> RSQRT_FRACTION_BITS);
+ y = COEFFS[2] + ((y * m) >> RSQRT_FRACTION_BITS);
+ y = COEFFS[1] + ((y * m) >> RSQRT_FRACTION_BITS);
+ y = COEFFS[0] + ((y * m) >> RSQRT_FRACTION_BITS);
+ return y;
+}
+
+LIBC_INLINE constexpr int64_t newton_raphson(uint32_t m, int64_t y) {
+ int64_t y2 = (y * y) >> RSQRT_FRACTION_BITS;
+ int64_t my2 = (static_cast<int64_t>(m) * y2) >> RSQRT_FRACTION_BITS;
+ int64_t factor = THREE_HALVES - (my2 >> 1);
+ return (y * factor) >> RSQRT_FRACTION_BITS;
+}
+
+LIBC_INLINE constexpr uint16_t fixed_to_half_bits(uint64_t y, int scale_exp) {
+ int y_log2 = floor_log2(y);
+ int out_exp = scale_exp + y_log2 - RSQRT_FRACTION_BITS;
+ int biased_exp = out_exp + 15;
+
+ uint32_t out_sig = y_log2 >= 10 ? static_cast<uint32_t>(y >> (y_log2 - 10))
+ : static_cast<uint32_t>(y << (10 - y_log2));
+
+ if (biased_exp <= 0)
+ return 0x0400;
+ if (biased_exp >= 31)
+ return 0x7bff;
+
+ return static_cast<uint16_t>((biased_exp << 10) | (out_sig & 0x3ff));
+}
+
+LIBC_INLINE constexpr uint16_t approximate_rsqrt(uint16_t x_abs) {
+ uint32_t x_mant = x_abs & 0x03ff;
+ int exponent = 0;
+
+ if (x_abs >= 0x0400) {
+ x_mant |= 0x0400;
+ exponent = static_cast<int>(x_abs >> 10) - 14;
+ } else {
+ exponent = -13;
+ while ((x_mant & 0x0400) == 0) {
+ x_mant <<= 1;
+ --exponent;
+ }
+ }
+
+ uint32_t m = x_mant << (RSQRT_FRACTION_BITS - 11);
+ int64_t y = newton_raphson(m, eval_polynomial(m));
+
+ int scale_exp = 0;
+ if (exponent & 1) {
+ y = (y * ONE_OVER_SQRT2) >> RSQRT_FRACTION_BITS;
+ scale_exp = -((exponent - 1) / 2);
+ } else {
+ scale_exp = -(exponent / 2);
+ }
+
+ return fixed_to_half_bits(static_cast<uint64_t>(y), scale_exp);
+}
+
+// Compare y = sig * 2^exp with 1 / sqrt(x_sig * 2^x_exp).
+// Return -1 if y is below the exact value, 0 if exact, and 1 if above.
+LIBC_INLINE constexpr int compare_with_rsqrt(uint32_t sig, int exp,
+ uint32_t x_sig, int x_exp) {
+ uint64_t lhs = static_cast<uint64_t>(sig) * sig * x_sig;
+ int scale = 2 * exp + x_exp;
+
+ if (scale >= 0)
+ return (scale == 0 && lhs == 1) ? 0 : 1;
+
+ int rshift = -scale;
+ if (rshift >= 64)
+ return -1;
+
+ uint64_t rhs = uint64_t(1) << rshift;
+ if (lhs < rhs)
+ return -1;
+ if (lhs > rhs)
+ return 1;
+ return 0;
+}
+
+LIBC_INLINE constexpr int compare_half_with_rsqrt(uint16_t y, uint32_t x_sig,
+ int x_exp) {
+ uint32_t y_sig = 0x0400 | (y & 0x03ff);
+ int y_exp = static_cast<int>(y >> 10) - 25;
+ return compare_with_rsqrt(y_sig, y_exp, x_sig, x_exp);
+}
+
+LIBC_INLINE constexpr uint16_t floor_rsqrt(uint16_t approx, uint32_t x_sig,
+ int x_exp) {
+ uint16_t y = approx < 0x0400 ? 0x0400 : approx;
+ while (compare_half_with_rsqrt(y, x_sig, x_exp) > 0)
+ --y;
+ while (y < 0x7bff && compare_half_with_rsqrt(y + 1, x_sig, x_exp) <= 0)
+ ++y;
+ return y;
+}
+
+LIBC_INLINE constexpr uint16_t round_result(uint16_t y, uint32_t x_sig,
+ int x_exp) {
+ if (compare_half_with_rsqrt(y, x_sig, x_exp) == 0)
+ return y;
+
+ int rounding_mode = FE_TONEAREST;
+ if (!cpp::is_constant_evaluated())
+ rounding_mode = fputil::get_round();
+ if (rounding_mode == FE_UPWARD)
+ return y + 1;
+ if (rounding_mode != FE_TONEAREST)
+ return y;
+
+ uint32_t y_sig = 0x0400 | (y & 0x03ff);
+ int y_exp = static_cast<int>(y >> 10) - 25;
+ uint32_t midpoint_sig = (y_sig << 1) | 1;
+ int midpoint_cmp = compare_with_rsqrt(midpoint_sig, y_exp - 1, x_sig, x_exp);
+
+ if (midpoint_cmp < 0)
+ return y + 1;
+ if (midpoint_cmp > 0)
+ return y;
+ return (y & 1) ? static_cast<uint16_t>(y + 1) : y;
+}
+
+LIBC_INLINE constexpr float16 rsqrtf16_no_float(uint16_t x_abs) {
+ uint32_t x_sig = 0;
+ int x_exp = 0;
+ if (x_abs >= 0x0400) {
+ x_sig = 0x0400 | (x_abs & 0x03ff);
+ x_exp = static_cast<int>(x_abs >> 10) - 25;
+ } else {
+ x_sig = x_abs;
+ x_exp = -24;
+ }
+
+ uint16_t approx = approximate_rsqrt(x_abs);
+ uint16_t y = floor_rsqrt(approx, x_sig, x_exp);
+ return fputil::FPBits<float16>(round_result(y, x_sig, x_exp)).get_val();
+}
+
+} // namespace rsqrtf16_internal
+
LIBC_INLINE constexpr float16 rsqrtf16(float16 x) {
using FPBits = fputil::FPBits<float16>;
FPBits xbits(x);
@@ -78,73 +243,7 @@ LIBC_INLINE constexpr float16 rsqrtf16(float16 x) {
return fputil::cast<float16>(result);
#else
- float xf = fputil::cast<float>(x);
-
- int exponent = 0;
- float mantissa = fputil::frexp(xf, exponent);
-
- float result = 0.0f;
- int exp_floored = -(exponent >> 1);
-
- if (mantissa == 0.5f) {
- // When mantissa is 0.5f, x was a power of 2 (or subnormal that normalizes
- // this way). 1/sqrt(0.5f) = sqrt(2.0f).
- // If exponent is odd (exponent = 2k + 1):
- // rsqrt(x) = (1/sqrt(0.5)) * 2^(-(2k+1)/2) = sqrt(2) * 2^(-k-0.5)
- // = sqrt(2) * 2^(-k) * (1/sqrt(2)) = 2^(-k)
- // exp_floored = -((2k+1)>>1) = -(k) = -k
- // So result = ldexp(1.0f, exp_floored)
- // If exponent is even (exponent = 2k):
- // rsqrt(x) = (1/sqrt(0.5)) * 2^(-2k/2) = sqrt(2) * 2^(-k)
- // exp_floored = -((2k)>>1) = -(k) = -k
- // So result = ldexp(sqrt(2.0f), exp_floored)
- if (exponent & 1) {
- result = fputil::ldexp(1.0f, exp_floored);
- } else {
- constexpr float SQRT_2_F = 0x1.6a09e6p0f; // sqrt(2.0f)
- result = fputil::ldexp(SQRT_2_F, exp_floored);
- }
- } else {
- // 4 Degree minimax polynomial (single-precision coefficients) generated
- // with Sollya:
- // P = fpminimax(1/sqrt(x), 4,
- // [|single,single,single,single,single|], [0.5;1])
- float y = fputil::polyeval(mantissa,
- 0x1.771256p1f, // c0
- -0x1.5e7c4ap2f, // c1
- 0x1.b3851cp2f, // c2
- -0x1.1a27ep2f, // c3
- 0x1.265c66p0f); // c4
-
- // Newton-Raphson iteration in float (use multiply_add to leverage FMA when
- // available):
- float y2 = y * y;
- float factor = fputil::multiply_add(-0.5f * mantissa, y2, 1.5f);
- y = y * factor;
-
- result = fputil::ldexp(y, exp_floored);
- if (exponent & 1) {
- constexpr float ONE_OVER_SQRT2 = 0x1.6a09e6p-1f; // 1/sqrt(2)
- result *= ONE_OVER_SQRT2;
- }
-
- // Targeted post-correction: for the specific half-precision mantissa
- // pattern M == 0x011F we observe a consistent -1 ULP bias across exponents.
- // Apply a tiny upward nudge to cross the rounding boundary in all modes.
- const uint16_t half_mantissa = static_cast<uint16_t>(x_abs & 0x3ff);
- if (half_mantissa == 0x011F) {
- // Nudge up to fix consistent -1 ULP at that mantissa boundary
- result = fputil::multiply_add(result, 0x1.0p-21f,
- result); // result *= (1 + 2^-21)
- } else if (half_mantissa == 0x0313) {
- // Nudge down to fix +1 ULP under upward rounding at this mantissa
- // boundary
- result = fputil::multiply_add(result, -0x1.0p-21f,
- result); // result *= (1 - 2^-21)
- }
- }
-
- return fputil::cast<float16>(result);
+ return rsqrtf16_internal::rsqrtf16_no_float(x_abs);
#endif
}
diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
index ed27d4c6e7bf0..e3efddb726acd 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -3716,7 +3716,6 @@ libc_support_library(
":__support_fputil_cast",
":__support_fputil_fenv_impl",
":__support_fputil_fp_bits",
- ":__support_fputil_manipulation_functions",
":__support_fputil_multiply_add",
":__support_fputil_sqrt",
":__support_macros_optimization",
>From fe03cd5a60cceb843ab241109296dc820cbe5910 Mon Sep 17 00:00:00 2001
From: amemov <shepelev777 at gmail.com>
Date: Sat, 25 Apr 2026 14:01:24 -0700
Subject: [PATCH 04/10] fix: remove duplicate target from cmake
---
libc/src/__support/math/CMakeLists.txt | 16 ++--------------
1 file changed, 2 insertions(+), 14 deletions(-)
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 87c925b4e0f8a..7c580742b04b0 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -126,20 +126,6 @@ add_header_library(
libc.src.__support.macros.properties.types
)
-add_header_library(
- rsqrtf16
- HDRS
- rsqrtf16.h
- DEPENDS
- libc.src.__support.FPUtil.cast
- libc.src.__support.FPUtil.fenv_impl
- libc.src.__support.FPUtil.fp_bits
- libc.src.__support.FPUtil.multiply_add
- libc.src.__support.FPUtil.sqrt
- libc.src.__support.macros.optimization
- libc.src.__support.macros.properties.types
-)
-
add_header_library(
asin_utils
HDRS
@@ -4596,8 +4582,10 @@ add_header_library(
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.multiply_add
libc.src.__support.FPUtil.sqrt
libc.src.__support.macros.optimization
+ libc.src.__support.macros.properties.types
)
add_header_library(
>From 918e6ca9eb479b9362286d57b8ac2cddb8996d9d Mon Sep 17 00:00:00 2001
From: amemov <shepelev777 at gmail.com>
Date: Sun, 26 Apr 2026 19:39:05 -0700
Subject: [PATCH 05/10] reduced branching
---
libc/src/__support/math/CMakeLists.txt | 4 ++++
libc/src/__support/math/rsqrtf16.h | 21 +++++++------------
.../llvm-project-overlay/libc/BUILD.bazel | 1 +
3 files changed, 13 insertions(+), 13 deletions(-)
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 7c580742b04b0..05e50f45caf83 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -4595,11 +4595,15 @@ add_header_library(
DEPENDS
libc.hdr.errno_macros
libc.hdr.fenv_macros
+ libc.include.llvm-libc-macros.float16_macros
+ libc.src.__support.CPP.bit
libc.src.__support.FPUtil.cast
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.multiply_add
libc.src.__support.FPUtil.sqrt
libc.src.__support.macros.optimization
+ libc.src.__support.macros.properties.types
)
add_header_library(
diff --git a/libc/src/__support/math/rsqrtf16.h b/libc/src/__support/math/rsqrtf16.h
index 23c74a8515c15..67eaab1e0f652 100644
--- a/libc/src/__support/math/rsqrtf16.h
+++ b/libc/src/__support/math/rsqrtf16.h
@@ -13,6 +13,7 @@
#ifdef LIBC_TYPES_HAS_FLOAT16
+#include "src/__support/CPP/bit.h"
#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/cast.h"
@@ -39,12 +40,7 @@ LIBC_INLINE_VAR constexpr int64_t COEFFS[5] = {
LIBC_INLINE_VAR constexpr int64_t ONE_OVER_SQRT2 = 0x16a09e60;
LIBC_INLINE constexpr int floor_log2(uint64_t x) {
- int result = -1;
- while (x) {
- x >>= 1;
- ++result;
- }
- return result;
+ return 63 - cpp::countl_zero(x);
}
LIBC_INLINE constexpr int64_t eval_polynomial(uint32_t m) {
@@ -87,11 +83,9 @@ LIBC_INLINE constexpr uint16_t approximate_rsqrt(uint16_t x_abs) {
x_mant |= 0x0400;
exponent = static_cast<int>(x_abs >> 10) - 14;
} else {
- exponent = -13;
- while ((x_mant & 0x0400) == 0) {
- x_mant <<= 1;
- --exponent;
- }
+ int shift = cpp::countl_zero(x_mant) - (32 - 11);
+ x_mant <<= shift;
+ exponent = -13 - shift;
}
uint32_t m = x_mant << (RSQRT_FRACTION_BITS - 11);
@@ -140,9 +134,10 @@ LIBC_INLINE constexpr int compare_half_with_rsqrt(uint16_t y, uint32_t x_sig,
LIBC_INLINE constexpr uint16_t floor_rsqrt(uint16_t approx, uint32_t x_sig,
int x_exp) {
uint16_t y = approx < 0x0400 ? 0x0400 : approx;
- while (compare_half_with_rsqrt(y, x_sig, x_exp) > 0)
+ if (LIBC_UNLIKELY(compare_half_with_rsqrt(y, x_sig, x_exp) > 0))
--y;
- while (y < 0x7bff && compare_half_with_rsqrt(y + 1, x_sig, x_exp) <= 0)
+ if (LIBC_UNLIKELY(y < 0x7bff &&
+ compare_half_with_rsqrt(y + 1, x_sig, x_exp) <= 0))
++y;
return y;
}
diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
index e3efddb726acd..ddc62c2b26232 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -3713,6 +3713,7 @@ libc_support_library(
name = "__support_math_rsqrtf16",
hdrs = ["src/__support/math/rsqrtf16.h"],
deps = [
+ ":__support_cpp_bit",
":__support_fputil_cast",
":__support_fputil_fenv_impl",
":__support_fputil_fp_bits",
>From efd70509251e305e29cab8c89957b40b01e09af6 Mon Sep 17 00:00:00 2001
From: amemov <shepelev777 at gmail.com>
Date: Fri, 5 Jun 2026 20:19:54 -0700
Subject: [PATCH 06/10] chore: added g-benchmark for rsqrtf16 / improved
approximation
---
libc/benchmarks/CMakeLists.txt | 26 ++++-
.../LibcRsqrtf16GoogleBenchmarkMain.cpp | 63 +++++++++++
libc/src/__support/math/rsqrtf16.h | 103 ++++++++++--------
3 files changed, 143 insertions(+), 49 deletions(-)
create mode 100644 libc/benchmarks/LibcRsqrtf16GoogleBenchmarkMain.cpp
diff --git a/libc/benchmarks/CMakeLists.txt b/libc/benchmarks/CMakeLists.txt
index 65c7cd76fad29..3e1619a940b26 100644
--- a/libc/benchmarks/CMakeLists.txt
+++ b/libc/benchmarks/CMakeLists.txt
@@ -19,9 +19,15 @@ set(LLVM_LINK_COMPONENTS
# Add Unit Testing Support
#==============================================================================
-make_gtest_target()
+if(COMMAND make_gtest_target)
+ make_gtest_target()
+endif()
function(add_libc_benchmark_unittest target_name)
+ if(NOT COMMAND make_gtest_target)
+ return()
+ endif()
+
if(NOT LLVM_INCLUDE_TESTS)
return()
endif()
@@ -214,3 +220,21 @@ target_link_libraries(libc.benchmarks.memory_functions.opt_host
benchmark_main
)
llvm_update_compile_flags(libc.benchmarks.memory_functions.opt_host)
+
+add_executable(libc.benchmarks.rsqrtf16.opt_host
+ EXCLUDE_FROM_ALL
+ LibcRsqrtf16GoogleBenchmarkMain.cpp
+)
+target_include_directories(libc.benchmarks.rsqrtf16.opt_host
+ PRIVATE
+ ${LIBC_SOURCE_DIR}
+ ${LIBC_INCLUDE_DIR}
+)
+target_link_libraries(libc.benchmarks.rsqrtf16.opt_host
+ PRIVATE
+ benchmark_main
+ libc-benchmark
+ libc.src.errno.errno.__internal__
+ libc.src.__support.math.rsqrtf16
+)
+llvm_update_compile_flags(libc.benchmarks.rsqrtf16.opt_host)
diff --git a/libc/benchmarks/LibcRsqrtf16GoogleBenchmarkMain.cpp b/libc/benchmarks/LibcRsqrtf16GoogleBenchmarkMain.cpp
new file mode 100644
index 0000000000000..ef93c7142177e
--- /dev/null
+++ b/libc/benchmarks/LibcRsqrtf16GoogleBenchmarkMain.cpp
@@ -0,0 +1,63 @@
+#include "benchmark/benchmark.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/math/rsqrtf16.h"
+
+#include <stddef.h>
+#include <stdint.h>
+
+namespace {
+
+using FPBits = LIBC_NAMESPACE::fputil::FPBits<float16>;
+
+constexpr uint16_t INPUTS[] = {
+ // Subnormals.
+ 0x0001, 0x0002, 0x0004, 0x0008, 0x0010, 0x0020, 0x0040, 0x0080,
+ 0x0100, 0x0200, 0x03ff,
+ // Normals spread across all exponent ranges.
+ 0x0400, 0x0401, 0x047f, 0x0555, 0x07ff, 0x0800, 0x0c00, 0x1000,
+ 0x1400, 0x1800, 0x1c00, 0x2000, 0x2400, 0x2800, 0x2c00, 0x3000,
+ 0x3400, 0x3800, 0x3c00, 0x3c01, 0x3d55, 0x3fff, 0x4000, 0x4400,
+ 0x4800, 0x4c00, 0x5000, 0x5400, 0x5800, 0x5c00, 0x6000, 0x6400,
+ 0x6800, 0x6c00, 0x7000, 0x7400, 0x7800, 0x7bff,
+};
+
+constexpr size_t INPUT_COUNT = sizeof(INPUTS) / sizeof(INPUTS[0]);
+
+float16 get_input(uint16_t bits) { return FPBits(bits).get_val(); }
+
+float16 rsqrtf16_integer_finite(float16 x) {
+ FPBits xbits(x);
+ return LIBC_NAMESPACE::math::rsqrtf16_internal::rsqrtf16_no_float(
+ xbits.uintval() & 0x7fff);
+}
+
+void BM_Rsqrtf16HostFpu(benchmark::State &state) {
+ for (auto _ : state) {
+ for (uint16_t bits : INPUTS)
+ benchmark::DoNotOptimize(LIBC_NAMESPACE::math::rsqrtf16(get_input(bits)));
+ }
+ state.SetItemsProcessed(state.iterations() * INPUT_COUNT);
+}
+
+void BM_Rsqrtf16IntegerFallbackFiniteWrapper(benchmark::State &state) {
+ for (auto _ : state) {
+ for (uint16_t bits : INPUTS)
+ benchmark::DoNotOptimize(rsqrtf16_integer_finite(get_input(bits)));
+ }
+ state.SetItemsProcessed(state.iterations() * INPUT_COUNT);
+}
+
+void BM_Rsqrtf16IntegerFallback(benchmark::State &state) {
+ for (auto _ : state) {
+ for (uint16_t bits : INPUTS)
+ benchmark::DoNotOptimize(
+ LIBC_NAMESPACE::math::rsqrtf16_internal::rsqrtf16_no_float(bits));
+ }
+ state.SetItemsProcessed(state.iterations() * INPUT_COUNT);
+}
+
+} // namespace
+
+BENCHMARK(BM_Rsqrtf16HostFpu);
+BENCHMARK(BM_Rsqrtf16IntegerFallbackFiniteWrapper);
+BENCHMARK(BM_Rsqrtf16IntegerFallback);
diff --git a/libc/src/__support/math/rsqrtf16.h b/libc/src/__support/math/rsqrtf16.h
index 67eaab1e0f652..7e7375036b150 100644
--- a/libc/src/__support/math/rsqrtf16.h
+++ b/libc/src/__support/math/rsqrtf16.h
@@ -30,12 +30,14 @@ LIBC_INLINE_VAR constexpr int RSQRT_FRACTION_BITS = 29;
LIBC_INLINE_VAR constexpr int64_t ONE = int64_t(1) << RSQRT_FRACTION_BITS;
LIBC_INLINE_VAR constexpr int64_t THREE_HALVES = 3 * (ONE >> 1);
-// Degree-4 minimax polynomial generated with Sollya:
-// P = fpminimax(1/sqrt(x), 4,
-// [|single,single,single,single,single|], [0.5;1])
-// Coefficients are stored in Q29 fixed-point format.
-LIBC_INLINE_VAR constexpr int64_t COEFFS[5] = {
- 1'573'164'416, -2'940'085'504, 3'653'406'208, -2'366'894'080, 617'319'616,
+// Midpoint lookup table for 1/sqrt(x) on 16 sub-intervals of [0.5;1).
+// Values are stored in Q29 fixed-point format. The Newton step and exact
+// rounding below correct the seed before producing the final half result.
+LIBC_INLINE_VAR constexpr uint32_t RSQRT_APPROX[16] = {
+ 747'657'839, 725'981'977, 706'088'274, 687'745'184,
+ 670'761'200, 654'976'372, 640'255'922, 626'485'368,
+ 613'566'757, 601'415'717, 589'959'130, 579'133'272,
+ 568'882'316, 559'157'115, 549'914'212, 541'115'017,
};
LIBC_INLINE_VAR constexpr int64_t ONE_OVER_SQRT2 = 0x16a09e60;
@@ -43,13 +45,8 @@ LIBC_INLINE constexpr int floor_log2(uint64_t x) {
return 63 - cpp::countl_zero(x);
}
-LIBC_INLINE constexpr int64_t eval_polynomial(uint32_t m) {
- int64_t y = COEFFS[4];
- y = COEFFS[3] + ((y * m) >> RSQRT_FRACTION_BITS);
- y = COEFFS[2] + ((y * m) >> RSQRT_FRACTION_BITS);
- y = COEFFS[1] + ((y * m) >> RSQRT_FRACTION_BITS);
- y = COEFFS[0] + ((y * m) >> RSQRT_FRACTION_BITS);
- return y;
+LIBC_INLINE constexpr int64_t initial_approximation(uint32_t x_mant) {
+ return RSQRT_APPROX[(x_mant - 0x0400) >> 6];
}
LIBC_INLINE constexpr int64_t newton_raphson(uint32_t m, int64_t y) {
@@ -75,13 +72,24 @@ LIBC_INLINE constexpr uint16_t fixed_to_half_bits(uint64_t y, int scale_exp) {
return static_cast<uint16_t>((biased_exp << 10) | (out_sig & 0x3ff));
}
-LIBC_INLINE constexpr uint16_t approximate_rsqrt(uint16_t x_abs) {
+struct ApproxResult {
+ uint16_t value;
+ uint32_t x_sig;
+ int x_exp;
+};
+
+LIBC_INLINE constexpr ApproxResult approximate_rsqrt(uint16_t x_abs) {
uint32_t x_mant = x_abs & 0x03ff;
+ uint32_t x_sig = x_mant;
+ int x_exp = -24;
int exponent = 0;
if (x_abs >= 0x0400) {
+ int biased_exp = static_cast<int>(x_abs >> 10);
x_mant |= 0x0400;
- exponent = static_cast<int>(x_abs >> 10) - 14;
+ x_sig = x_mant;
+ x_exp = biased_exp - 25;
+ exponent = biased_exp - 14;
} else {
int shift = cpp::countl_zero(x_mant) - (32 - 11);
x_mant <<= shift;
@@ -89,7 +97,7 @@ LIBC_INLINE constexpr uint16_t approximate_rsqrt(uint16_t x_abs) {
}
uint32_t m = x_mant << (RSQRT_FRACTION_BITS - 11);
- int64_t y = newton_raphson(m, eval_polynomial(m));
+ int64_t y = newton_raphson(m, initial_approximation(x_mant));
int scale_exp = 0;
if (exponent & 1) {
@@ -99,7 +107,8 @@ LIBC_INLINE constexpr uint16_t approximate_rsqrt(uint16_t x_abs) {
scale_exp = -(exponent / 2);
}
- return fixed_to_half_bits(static_cast<uint64_t>(y), scale_exp);
+ return {fixed_to_half_bits(static_cast<uint64_t>(y), scale_exp), x_sig,
+ x_exp};
}
// Compare y = sig * 2^exp with 1 / sqrt(x_sig * 2^x_exp).
@@ -107,15 +116,9 @@ LIBC_INLINE constexpr uint16_t approximate_rsqrt(uint16_t x_abs) {
LIBC_INLINE constexpr int compare_with_rsqrt(uint32_t sig, int exp,
uint32_t x_sig, int x_exp) {
uint64_t lhs = static_cast<uint64_t>(sig) * sig * x_sig;
- int scale = 2 * exp + x_exp;
-
- if (scale >= 0)
- return (scale == 0 && lhs == 1) ? 0 : 1;
-
- int rshift = -scale;
- if (rshift >= 64)
- return -1;
-
+ // For all finite positive half inputs and candidates produced by this
+ // algorithm, 2 * exp + x_exp is in [-34, -20].
+ int rshift = -(2 * exp + x_exp);
uint64_t rhs = uint64_t(1) << rshift;
if (lhs < rhs)
return -1;
@@ -131,20 +134,32 @@ LIBC_INLINE constexpr int compare_half_with_rsqrt(uint16_t y, uint32_t x_sig,
return compare_with_rsqrt(y_sig, y_exp, x_sig, x_exp);
}
-LIBC_INLINE constexpr uint16_t floor_rsqrt(uint16_t approx, uint32_t x_sig,
- int x_exp) {
+struct FloorResult {
+ uint16_t value;
+ int cmp;
+};
+
+LIBC_INLINE constexpr FloorResult floor_rsqrt(uint16_t approx, uint32_t x_sig,
+ int x_exp) {
uint16_t y = approx < 0x0400 ? 0x0400 : approx;
- if (LIBC_UNLIKELY(compare_half_with_rsqrt(y, x_sig, x_exp) > 0))
+ int cmp = compare_half_with_rsqrt(y, x_sig, x_exp);
+ if (LIBC_UNLIKELY(cmp > 0)) {
--y;
- if (LIBC_UNLIKELY(y < 0x7bff &&
- compare_half_with_rsqrt(y + 1, x_sig, x_exp) <= 0))
- ++y;
- return y;
+ cmp = compare_half_with_rsqrt(y, x_sig, x_exp);
+ } else if (LIBC_UNLIKELY(y < 0x7bff)) {
+ int next_cmp = compare_half_with_rsqrt(y + 1, x_sig, x_exp);
+ if (LIBC_UNLIKELY(next_cmp <= 0)) {
+ ++y;
+ cmp = next_cmp;
+ }
+ }
+ return {y, cmp};
}
-LIBC_INLINE constexpr uint16_t round_result(uint16_t y, uint32_t x_sig,
+LIBC_INLINE constexpr uint16_t round_result(FloorResult floor, uint32_t x_sig,
int x_exp) {
- if (compare_half_with_rsqrt(y, x_sig, x_exp) == 0)
+ uint16_t y = floor.value;
+ if (floor.cmp == 0)
return y;
int rounding_mode = FE_TONEAREST;
@@ -168,19 +183,11 @@ LIBC_INLINE constexpr uint16_t round_result(uint16_t y, uint32_t x_sig,
}
LIBC_INLINE constexpr float16 rsqrtf16_no_float(uint16_t x_abs) {
- uint32_t x_sig = 0;
- int x_exp = 0;
- if (x_abs >= 0x0400) {
- x_sig = 0x0400 | (x_abs & 0x03ff);
- x_exp = static_cast<int>(x_abs >> 10) - 25;
- } else {
- x_sig = x_abs;
- x_exp = -24;
- }
-
- uint16_t approx = approximate_rsqrt(x_abs);
- uint16_t y = floor_rsqrt(approx, x_sig, x_exp);
- return fputil::FPBits<float16>(round_result(y, x_sig, x_exp)).get_val();
+ ApproxResult approx = approximate_rsqrt(x_abs);
+ FloorResult floor = floor_rsqrt(approx.value, approx.x_sig, approx.x_exp);
+ return fputil::FPBits<float16>(
+ round_result(floor, approx.x_sig, approx.x_exp))
+ .get_val();
}
} // namespace rsqrtf16_internal
>From f67434782f061d8a2082c6e3dffe6c2d2e02d60a Mon Sep 17 00:00:00 2001
From: amemov <shepelev777 at gmail.com>
Date: Fri, 5 Jun 2026 20:23:58 -0700
Subject: [PATCH 07/10] fix: reformatted benchmark .cpp file
---
.../LibcRsqrtf16GoogleBenchmarkMain.cpp | 56 ++++++++++++++++---
1 file changed, 49 insertions(+), 7 deletions(-)
diff --git a/libc/benchmarks/LibcRsqrtf16GoogleBenchmarkMain.cpp b/libc/benchmarks/LibcRsqrtf16GoogleBenchmarkMain.cpp
index ef93c7142177e..2c6e4b02018a4 100644
--- a/libc/benchmarks/LibcRsqrtf16GoogleBenchmarkMain.cpp
+++ b/libc/benchmarks/LibcRsqrtf16GoogleBenchmarkMain.cpp
@@ -11,14 +11,56 @@ using FPBits = LIBC_NAMESPACE::fputil::FPBits<float16>;
constexpr uint16_t INPUTS[] = {
// Subnormals.
- 0x0001, 0x0002, 0x0004, 0x0008, 0x0010, 0x0020, 0x0040, 0x0080,
- 0x0100, 0x0200, 0x03ff,
+ 0x0001,
+ 0x0002,
+ 0x0004,
+ 0x0008,
+ 0x0010,
+ 0x0020,
+ 0x0040,
+ 0x0080,
+ 0x0100,
+ 0x0200,
+ 0x03ff,
// Normals spread across all exponent ranges.
- 0x0400, 0x0401, 0x047f, 0x0555, 0x07ff, 0x0800, 0x0c00, 0x1000,
- 0x1400, 0x1800, 0x1c00, 0x2000, 0x2400, 0x2800, 0x2c00, 0x3000,
- 0x3400, 0x3800, 0x3c00, 0x3c01, 0x3d55, 0x3fff, 0x4000, 0x4400,
- 0x4800, 0x4c00, 0x5000, 0x5400, 0x5800, 0x5c00, 0x6000, 0x6400,
- 0x6800, 0x6c00, 0x7000, 0x7400, 0x7800, 0x7bff,
+ 0x0400,
+ 0x0401,
+ 0x047f,
+ 0x0555,
+ 0x07ff,
+ 0x0800,
+ 0x0c00,
+ 0x1000,
+ 0x1400,
+ 0x1800,
+ 0x1c00,
+ 0x2000,
+ 0x2400,
+ 0x2800,
+ 0x2c00,
+ 0x3000,
+ 0x3400,
+ 0x3800,
+ 0x3c00,
+ 0x3c01,
+ 0x3d55,
+ 0x3fff,
+ 0x4000,
+ 0x4400,
+ 0x4800,
+ 0x4c00,
+ 0x5000,
+ 0x5400,
+ 0x5800,
+ 0x5c00,
+ 0x6000,
+ 0x6400,
+ 0x6800,
+ 0x6c00,
+ 0x7000,
+ 0x7400,
+ 0x7800,
+ 0x7bff,
};
constexpr size_t INPUT_COUNT = sizeof(INPUTS) / sizeof(INPUTS[0]);
>From e96d8711566e5054b248cd78a3b2e67178880a92 Mon Sep 17 00:00:00 2001
From: amemov <shepelev777 at gmail.com>
Date: Fri, 5 Jun 2026 20:43:57 -0700
Subject: [PATCH 08/10] chore: added more annotations to the code
---
libc/src/__support/math/rsqrtf16.h | 19 +++++++++++++++++++
1 file changed, 19 insertions(+)
diff --git a/libc/src/__support/math/rsqrtf16.h b/libc/src/__support/math/rsqrtf16.h
index 7e7375036b150..cebdf0b4e7067 100644
--- a/libc/src/__support/math/rsqrtf16.h
+++ b/libc/src/__support/math/rsqrtf16.h
@@ -26,6 +26,9 @@ namespace math {
namespace rsqrtf16_internal {
+// Fixed-point computations below use Q29: the integer N represents
+// N * 2^-29. Multiplying two Q29 values produces a Q58 value, so products are
+// shifted right by RSQRT_FRACTION_BITS to return to Q29.
LIBC_INLINE_VAR constexpr int RSQRT_FRACTION_BITS = 29;
LIBC_INLINE_VAR constexpr int64_t ONE = int64_t(1) << RSQRT_FRACTION_BITS;
LIBC_INLINE_VAR constexpr int64_t THREE_HALVES = 3 * (ONE >> 1);
@@ -50,6 +53,9 @@ LIBC_INLINE constexpr int64_t initial_approximation(uint32_t x_mant) {
}
LIBC_INLINE constexpr int64_t newton_raphson(uint32_t m, int64_t y) {
+ // Refine y ~= 1/sqrt(m) with:
+ // y_{n+1} = y_n * (1.5 - 0.5 * m * y_n^2)
+ // where both m and y are stored in Q29.
int64_t y2 = (y * y) >> RSQRT_FRACTION_BITS;
int64_t my2 = (static_cast<int64_t>(m) * y2) >> RSQRT_FRACTION_BITS;
int64_t factor = THREE_HALVES - (my2 >> 1);
@@ -84,6 +90,10 @@ LIBC_INLINE constexpr ApproxResult approximate_rsqrt(uint16_t x_abs) {
int x_exp = -24;
int exponent = 0;
+ // Decompose the finite positive input as:
+ // x = m * 2^exponent, with 0.5 <= m < 1.
+ // `x_sig` and `x_exp` keep the exact input as x_sig * 2^x_exp for the integer
+ // rounding test below.
if (x_abs >= 0x0400) {
int biased_exp = static_cast<int>(x_abs >> 10);
x_mant |= 0x0400;
@@ -99,6 +109,8 @@ LIBC_INLINE constexpr ApproxResult approximate_rsqrt(uint16_t x_abs) {
uint32_t m = x_mant << (RSQRT_FRACTION_BITS - 11);
int64_t y = newton_raphson(m, initial_approximation(x_mant));
+ // Since rsqrt(m * 2^e) = rsqrt(m) * 2^(-e/2), odd exponents need one
+ // extra factor of 1/sqrt(2) before applying the integral power of two.
int scale_exp = 0;
if (exponent & 1) {
y = (y * ONE_OVER_SQRT2) >> RSQRT_FRACTION_BITS;
@@ -113,6 +125,8 @@ LIBC_INLINE constexpr ApproxResult approximate_rsqrt(uint16_t x_abs) {
// Compare y = sig * 2^exp with 1 / sqrt(x_sig * 2^x_exp).
// Return -1 if y is below the exact value, 0 if exact, and 1 if above.
+// Instead of computing a reciprocal square root, square both sides:
+// y <= 1/sqrt(x) <=> y^2 * x <= 1.
LIBC_INLINE constexpr int compare_with_rsqrt(uint32_t sig, int exp,
uint32_t x_sig, int x_exp) {
uint64_t lhs = static_cast<uint64_t>(sig) * sig * x_sig;
@@ -141,6 +155,8 @@ struct FloorResult {
LIBC_INLINE constexpr FloorResult floor_rsqrt(uint16_t approx, uint32_t x_sig,
int x_exp) {
+ // The table seed and Newton step have been validated exhaustively to produce
+ // a candidate at most one half-precision step below the exact floor.
uint16_t y = approx < 0x0400 ? 0x0400 : approx;
int cmp = compare_half_with_rsqrt(y, x_sig, x_exp);
if (LIBC_UNLIKELY(cmp > 0)) {
@@ -162,6 +178,9 @@ LIBC_INLINE constexpr uint16_t round_result(FloorResult floor, uint32_t x_sig,
if (floor.cmp == 0)
return y;
+ // Once `y` is the greatest half value below the exact result, directed
+ // rounding is immediate. Round-to-nearest compares against the midpoint
+ // between `y` and the next half value, then applies ties-to-even.
int rounding_mode = FE_TONEAREST;
if (!cpp::is_constant_evaluated())
rounding_mode = fputil::get_round();
>From f82cad9ce85f7440c9548443ac900ba7d36b1de3 Mon Sep 17 00:00:00 2001
From: amemov <shepelev777 at gmail.com>
Date: Fri, 5 Jun 2026 20:53:00 -0700
Subject: [PATCH 09/10] fix: remove tmp stub from cmake file for benchmarks
---
libc/benchmarks/CMakeLists.txt | 8 +-------
1 file changed, 1 insertion(+), 7 deletions(-)
diff --git a/libc/benchmarks/CMakeLists.txt b/libc/benchmarks/CMakeLists.txt
index 3e1619a940b26..13f1f00de9b4a 100644
--- a/libc/benchmarks/CMakeLists.txt
+++ b/libc/benchmarks/CMakeLists.txt
@@ -19,15 +19,9 @@ set(LLVM_LINK_COMPONENTS
# Add Unit Testing Support
#==============================================================================
-if(COMMAND make_gtest_target)
- make_gtest_target()
-endif()
+make_gtest_target()
function(add_libc_benchmark_unittest target_name)
- if(NOT COMMAND make_gtest_target)
- return()
- endif()
-
if(NOT LLVM_INCLUDE_TESTS)
return()
endif()
>From 2e3b72c32604db56becf3e8bb99bcb97f7ce7021 Mon Sep 17 00:00:00 2001
From: amemov <shepelev777 at gmail.com>
Date: Wed, 10 Jun 2026 23:25:27 -0700
Subject: [PATCH 10/10] chore: upd names for constants + add comments
---
libc/src/__support/math/rsqrtf16.h | 101 +++++++++++++++++++++--------
1 file changed, 73 insertions(+), 28 deletions(-)
diff --git a/libc/src/__support/math/rsqrtf16.h b/libc/src/__support/math/rsqrtf16.h
index cebdf0b4e7067..d3356d4dbc9e5 100644
--- a/libc/src/__support/math/rsqrtf16.h
+++ b/libc/src/__support/math/rsqrtf16.h
@@ -26,6 +26,8 @@ namespace math {
namespace rsqrtf16_internal {
+using FPBits = fputil::FPBits<float16>;
+
// Fixed-point computations below use Q29: the integer N represents
// N * 2^-29. Multiplying two Q29 values produces a Q58 value, so products are
// shifted right by RSQRT_FRACTION_BITS to return to Q29.
@@ -33,14 +35,44 @@ LIBC_INLINE_VAR constexpr int RSQRT_FRACTION_BITS = 29;
LIBC_INLINE_VAR constexpr int64_t ONE = int64_t(1) << RSQRT_FRACTION_BITS;
LIBC_INLINE_VAR constexpr int64_t THREE_HALVES = 3 * (ONE >> 1);
+LIBC_INLINE_VAR constexpr int HALF_FRACTION_LEN = FPBits::FRACTION_LEN;
+LIBC_INLINE_VAR constexpr int HALF_SIGNIFICAND_LEN = HALF_FRACTION_LEN + 1;
+LIBC_INLINE_VAR constexpr int HALF_EXP_BIAS = FPBits::EXP_BIAS;
+LIBC_INLINE_VAR constexpr uint16_t HALF_FRACTION_MASK = FPBits::FRACTION_MASK;
+LIBC_INLINE_VAR constexpr uint16_t HALF_MIN_NORMAL =
+ FPBits::min_normal().uintval();
+LIBC_INLINE_VAR constexpr uint16_t HALF_MAX_NORMAL =
+ FPBits::max_normal().uintval();
+LIBC_INLINE_VAR constexpr uint32_t HALF_HIDDEN_BIT =
+ uint32_t(1) << HALF_FRACTION_LEN;
+LIBC_INLINE_VAR constexpr int UINT32_BITS = 8 * static_cast<int>(sizeof(uint32_t));
+
+// Exact representation exponents for values stored as:
+// x = significand * 2^exponent,
+// where normal significands include the hidden bit.
+LIBC_INLINE_VAR constexpr int EXACT_NORMAL_EXP_OFFSET =
+ -HALF_EXP_BIAS - HALF_FRACTION_LEN;
+LIBC_INLINE_VAR constexpr int EXACT_SUBNORMAL_EXP =
+ 1 - HALF_EXP_BIAS - HALF_FRACTION_LEN;
+
+// Exponents for the reduced form used by the approximation:
+// x = m * 2^exponent, with 0.5 <= m < 1.
+LIBC_INLINE_VAR constexpr int REDUCED_NORMAL_EXP_OFFSET = 1 - HALF_EXP_BIAS;
+LIBC_INLINE_VAR constexpr int REDUCED_SUBNORMAL_EXP =
+ EXACT_SUBNORMAL_EXP + HALF_SIGNIFICAND_LEN;
+
+LIBC_INLINE_VAR constexpr int RSQRT_APPROX_BITS = 4;
+LIBC_INLINE_VAR constexpr int RSQRT_APPROX_INDEX_SHIFT =
+ HALF_FRACTION_LEN - RSQRT_APPROX_BITS;
+
// Midpoint lookup table for 1/sqrt(x) on 16 sub-intervals of [0.5;1).
// Values are stored in Q29 fixed-point format. The Newton step and exact
// rounding below correct the seed before producing the final half result.
LIBC_INLINE_VAR constexpr uint32_t RSQRT_APPROX[16] = {
- 747'657'839, 725'981'977, 706'088'274, 687'745'184,
- 670'761'200, 654'976'372, 640'255'922, 626'485'368,
- 613'566'757, 601'415'717, 589'959'130, 579'133'272,
- 568'882'316, 559'157'115, 549'914'212, 541'115'017,
+ 0x2c905a6f, 0x2b459b19, 0x2a160d52, 0x28fe28a0,
+ 0x27fb00f0, 0x270a2574, 0x262987b2, 0x25576878,
+ 0x24924925, 0x23d8e025, 0x232a0fda, 0x2284df58,
+ 0x21e8748c, 0x21540f7b, 0x20c70664, 0x2040c289,
};
LIBC_INLINE_VAR constexpr int64_t ONE_OVER_SQRT2 = 0x16a09e60;
@@ -49,7 +81,8 @@ LIBC_INLINE constexpr int floor_log2(uint64_t x) {
}
LIBC_INLINE constexpr int64_t initial_approximation(uint32_t x_mant) {
- return RSQRT_APPROX[(x_mant - 0x0400) >> 6];
+ return RSQRT_APPROX[(x_mant - HALF_HIDDEN_BIT) >>
+ RSQRT_APPROX_INDEX_SHIFT];
}
LIBC_INLINE constexpr int64_t newton_raphson(uint32_t m, int64_t y) {
@@ -63,21 +96,31 @@ LIBC_INLINE constexpr int64_t newton_raphson(uint32_t m, int64_t y) {
}
LIBC_INLINE constexpr uint16_t fixed_to_half_bits(uint64_t y, int scale_exp) {
+ // Convert y * 2^scale_exp, with y in Q29, to an approximate positive normal
+ // half bit pattern. This only creates a nearby candidate; exact rounding is
+ // handled by floor_rsqrt and round_result.
int y_log2 = floor_log2(y);
int out_exp = scale_exp + y_log2 - RSQRT_FRACTION_BITS;
- int biased_exp = out_exp + 15;
+ int biased_exp = out_exp + HALF_EXP_BIAS;
- uint32_t out_sig = y_log2 >= 10 ? static_cast<uint32_t>(y >> (y_log2 - 10))
- : static_cast<uint32_t>(y << (10 - y_log2));
+ uint32_t out_sig =
+ y_log2 >= HALF_FRACTION_LEN
+ ? static_cast<uint32_t>(y >> (y_log2 - HALF_FRACTION_LEN))
+ : static_cast<uint32_t>(y << (HALF_FRACTION_LEN - y_log2));
if (biased_exp <= 0)
- return 0x0400;
- if (biased_exp >= 31)
- return 0x7bff;
+ return HALF_MIN_NORMAL;
+ if (biased_exp >= FPBits::MAX_BIASED_EXPONENT)
+ return HALF_MAX_NORMAL;
- return static_cast<uint16_t>((biased_exp << 10) | (out_sig & 0x3ff));
+ return static_cast<uint16_t>((biased_exp << HALF_FRACTION_LEN) |
+ (out_sig & HALF_FRACTION_MASK));
}
+// `value` is the approximate positive half bit pattern produced by the table
+// seed, Newton step, and exponent scaling. `x_sig` and `x_exp` keep the exact
+// input as x_sig * 2^x_exp, which is needed to compare candidates against the
+// mathematical result without floating-point operations.
struct ApproxResult {
uint16_t value;
uint32_t x_sig;
@@ -85,28 +128,28 @@ struct ApproxResult {
};
LIBC_INLINE constexpr ApproxResult approximate_rsqrt(uint16_t x_abs) {
- uint32_t x_mant = x_abs & 0x03ff;
+ uint32_t x_mant = x_abs & HALF_FRACTION_MASK;
uint32_t x_sig = x_mant;
- int x_exp = -24;
+ int x_exp = EXACT_SUBNORMAL_EXP;
int exponent = 0;
// Decompose the finite positive input as:
// x = m * 2^exponent, with 0.5 <= m < 1.
// `x_sig` and `x_exp` keep the exact input as x_sig * 2^x_exp for the integer
// rounding test below.
- if (x_abs >= 0x0400) {
- int biased_exp = static_cast<int>(x_abs >> 10);
- x_mant |= 0x0400;
+ if (x_abs >= HALF_MIN_NORMAL) {
+ int biased_exp = static_cast<int>(x_abs >> HALF_FRACTION_LEN);
+ x_mant |= HALF_HIDDEN_BIT;
x_sig = x_mant;
- x_exp = biased_exp - 25;
- exponent = biased_exp - 14;
+ x_exp = biased_exp + EXACT_NORMAL_EXP_OFFSET;
+ exponent = biased_exp + REDUCED_NORMAL_EXP_OFFSET;
} else {
- int shift = cpp::countl_zero(x_mant) - (32 - 11);
+ int shift = cpp::countl_zero(x_mant) - (UINT32_BITS - HALF_SIGNIFICAND_LEN);
x_mant <<= shift;
- exponent = -13 - shift;
+ exponent = REDUCED_SUBNORMAL_EXP - shift;
}
- uint32_t m = x_mant << (RSQRT_FRACTION_BITS - 11);
+ uint32_t m = x_mant << (RSQRT_FRACTION_BITS - HALF_SIGNIFICAND_LEN);
int64_t y = newton_raphson(m, initial_approximation(x_mant));
// Since rsqrt(m * 2^e) = rsqrt(m) * 2^(-e/2), odd exponents need one
@@ -143,8 +186,9 @@ LIBC_INLINE constexpr int compare_with_rsqrt(uint32_t sig, int exp,
LIBC_INLINE constexpr int compare_half_with_rsqrt(uint16_t y, uint32_t x_sig,
int x_exp) {
- uint32_t y_sig = 0x0400 | (y & 0x03ff);
- int y_exp = static_cast<int>(y >> 10) - 25;
+ uint32_t y_sig = HALF_HIDDEN_BIT | (y & HALF_FRACTION_MASK);
+ int y_exp =
+ static_cast<int>(y >> HALF_FRACTION_LEN) + EXACT_NORMAL_EXP_OFFSET;
return compare_with_rsqrt(y_sig, y_exp, x_sig, x_exp);
}
@@ -157,12 +201,12 @@ LIBC_INLINE constexpr FloorResult floor_rsqrt(uint16_t approx, uint32_t x_sig,
int x_exp) {
// The table seed and Newton step have been validated exhaustively to produce
// a candidate at most one half-precision step below the exact floor.
- uint16_t y = approx < 0x0400 ? 0x0400 : approx;
+ uint16_t y = approx < HALF_MIN_NORMAL ? HALF_MIN_NORMAL : approx;
int cmp = compare_half_with_rsqrt(y, x_sig, x_exp);
if (LIBC_UNLIKELY(cmp > 0)) {
--y;
cmp = compare_half_with_rsqrt(y, x_sig, x_exp);
- } else if (LIBC_UNLIKELY(y < 0x7bff)) {
+ } else if (LIBC_UNLIKELY(y < HALF_MAX_NORMAL)) {
int next_cmp = compare_half_with_rsqrt(y + 1, x_sig, x_exp);
if (LIBC_UNLIKELY(next_cmp <= 0)) {
++y;
@@ -189,8 +233,9 @@ LIBC_INLINE constexpr uint16_t round_result(FloorResult floor, uint32_t x_sig,
if (rounding_mode != FE_TONEAREST)
return y;
- uint32_t y_sig = 0x0400 | (y & 0x03ff);
- int y_exp = static_cast<int>(y >> 10) - 25;
+ uint32_t y_sig = HALF_HIDDEN_BIT | (y & HALF_FRACTION_MASK);
+ int y_exp =
+ static_cast<int>(y >> HALF_FRACTION_LEN) + EXACT_NORMAL_EXP_OFFSET;
uint32_t midpoint_sig = (y_sig << 1) | 1;
int midpoint_cmp = compare_with_rsqrt(midpoint_sig, y_exp - 1, x_sig, x_exp);
More information about the libc-commits
mailing list