[libc-commits] [libc] [libc][math][c23] Add rsqrtf16() function (PR #137545)
via libc-commits
libc-commits at lists.llvm.org
Tue May 20 13:13:26 PDT 2025
https://github.com/amemov updated https://github.com/llvm/llvm-project/pull/137545
>From 43fefedc6b8c644fb71aa5a12b17e9c76aac6cef Mon Sep 17 00:00:00 2001
From: amemov <shepelev777 at gmail.com>
Date: Sun, 27 Apr 2025 19:54:02 +0000
Subject: [PATCH 1/5] Initial implementation of rsqrtf16 - Defined and declared
rsqrtf16 - Added entrypoints to support rsqrt in float16 format - Added tests
- both unit and exhaustive - Added MPFR support for Rsqrt to verify the test
results
TODO:
- Write computation as a polynomial in the end of rsqrtf16.cpp
- Check if errors are generated correctly (Reference from C23: The rsqrt functions compute the reciprocal of the nonnegative square root of the argument. A domain error occurs if the argument is less than zero. A pole error may occur if the argument equals zero.)
---
libc/config/linux/x86_64/entrypoints.txt | 1 +
libc/docs/headers/math/index.rst | 2 +-
libc/include/math.yaml | 7 +++
libc/src/math/CMakeLists.txt | 2 +
libc/src/math/generic/CMakeLists.txt | 19 ++++++
libc/src/math/generic/rsqrtf16.cpp | 67 ++++++++++++++++++++++
libc/src/math/rsqrtf16.h | 21 +++++++
libc/test/src/math/CMakeLists.txt | 11 ++++
libc/test/src/math/rsqrtf16_test.cpp | 42 ++++++++++++++
libc/test/src/math/smoke/CMakeLists.txt | 11 ++++
libc/test/src/math/smoke/rsqrtf16_test.cpp | 38 ++++++++++++
libc/utils/MPFRWrapper/MPCommon.cpp | 6 ++
libc/utils/MPFRWrapper/MPCommon.h | 1 +
libc/utils/MPFRWrapper/MPFRUtils.cpp | 2 +
libc/utils/MPFRWrapper/MPFRUtils.h | 1 +
15 files changed, 230 insertions(+), 1 deletion(-)
create mode 100644 libc/src/math/generic/rsqrtf16.cpp
create mode 100644 libc/src/math/rsqrtf16.h
create mode 100644 libc/test/src/math/rsqrtf16_test.cpp
create mode 100644 libc/test/src/math/smoke/rsqrtf16_test.cpp
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 0e4fbdd2f2f57..dc3120ca913a3 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -740,6 +740,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
libc.src.math.rintf16
libc.src.math.roundevenf16
libc.src.math.roundf16
+ libc.src.math.rsqrtf16
libc.src.math.scalblnf16
libc.src.math.scalbnf16
libc.src.math.setpayloadf16
diff --git a/libc/docs/headers/math/index.rst b/libc/docs/headers/math/index.rst
index cf88e6237d1e3..a29e70a7b625b 100644
--- a/libc/docs/headers/math/index.rst
+++ b/libc/docs/headers/math/index.rst
@@ -333,7 +333,7 @@ Higher Math Functions
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
| rootn | | | | | | 7.12.7.8 | F.10.4.8 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
-| rsqrt | | | | | | 7.12.7.9 | F.10.4.9 |
+| rsqrt | | | | |check| | | 7.12.7.9 | F.10.4.9 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
| sin | |check| | |check| | | |check| | | 7.12.4.6 | F.10.1.6 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
diff --git a/libc/include/math.yaml b/libc/include/math.yaml
index b725c33c0bb06..35aee5f54d8cd 100644
--- a/libc/include/math.yaml
+++ b/libc/include/math.yaml
@@ -2237,6 +2237,13 @@ functions:
return_type: long double
arguments:
- type: long double
+ - name: rsqrtf16
+ standards:
+ - stdc
+ return_type: _Float16
+ arguments:
+ - type: _Float16
+ guard: LIBC_TYPES_HAS_FLOAT16
- name: scalbln
standards:
- stdc
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index 5161b2b61aa94..689327b2f1249 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -467,6 +467,8 @@ add_math_entrypoint_object(roundevenl)
add_math_entrypoint_object(roundevenf16)
add_math_entrypoint_object(roundevenf128)
+add_math_entrypoint_object(rsqrtf16)
+
add_math_entrypoint_object(scalbln)
add_math_entrypoint_object(scalblnf)
add_math_entrypoint_object(scalblnl)
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 25f583035fbc2..ae0977ce255c2 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -955,6 +955,25 @@ add_entrypoint_object(
libc.src.__support.FPUtil.nearest_integer_operations
)
+add_entrypoint_object(
+ rsqrtf16
+ SRCS
+ rsqrtf16.cpp
+ HDRS
+ ../rsqrtf16.h
+ DEPENDS
+ libc.hdr.errno_macros
+ libc.hdr.fenv_macros
+ 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.polyeval
+ libc.src.__support.FPUtil.sqrt
+ libc.src.__support.macros.optimization
+ libc.src.__support.macros.properties.types
+)
+
add_entrypoint_object(
lround
SRCS
diff --git a/libc/src/math/generic/rsqrtf16.cpp b/libc/src/math/generic/rsqrtf16.cpp
new file mode 100644
index 0000000000000..930c8f6fc7b4e
--- /dev/null
+++ b/libc/src/math/generic/rsqrtf16.cpp
@@ -0,0 +1,67 @@
+//===-- Half-precision rsqrt function -------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception.
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/math/rsqrtf16.h"
+#include "hdr/errno_macros.h"
+#include "hdr/fenv_macros.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.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"
+#include "src/__support/macros/optimization.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+LLVM_LIBC_FUNCTION(float16, rsqrtf16, (float16 x)) {
+ using FPBits = fputil::FPBits<float16>;
+ FPBits xbits(x);
+
+ uint16_t x_u = xbits.uintval();
+ uint16_t x_abs = x_u & 0x7fff;
+ uint16_t x_sign = x_u >> 15;
+
+ // x is NaN
+ if (LIBC_UNLIKELY(xbits.is_nan())) {
+ if (xbits.is_signaling_nan()) {
+ fputil::raise_except_if_required(FE_INVALID);
+ return FPBits::quiet_nan().get_val();
+ }
+ return x;
+ }
+
+ // |x| = 0
+ if (LIBC_UNLIKELY(x_abs == 0x0)) {
+ fputil::raise_except_if_required(FE_DIVBYZERO);
+ fputil::set_errno_if_required(ERANGE);
+ return FPBits::quiet_nan().get_val();
+ }
+
+ // -inf <= x < 0
+ if (LIBC_UNLIKELY(x_sign == 1)) {
+ fputil::raise_except_if_required(FE_INVALID);
+ fputil::set_errno_if_required(EDOM);
+ return FPBits::quiet_nan().get_val();
+ }
+
+ // x = +inf => rsqrt(x) = 0
+ if (LIBC_UNLIKELY(xbits.is_inf())) {
+ return fputil::cast<float16>(0.0f);
+ }
+
+ // x = 1 => rsqrt(x) = 1
+ if (LIBC_UNLIKELY(x_u == 0x1)) {
+ return fputil::cast<float16>(1.0f);
+ }
+
+ // x is valid, estimate the result - below is temporary solution for just testing
+ float xf = x;
+ return fputil::cast<float16>(1.0f / xf);
+}
+} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/rsqrtf16.h b/libc/src/math/rsqrtf16.h
new file mode 100644
index 0000000000000..c88ab5256ce88
--- /dev/null
+++ b/libc/src/math/rsqrtf16.h
@@ -0,0 +1,21 @@
+//===-- Implementation header for rsqrtf16 ----------------------*- C++ -*-===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception.
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC_MATH_RSQRTF16_H
+#define LLVM_LIBC_SRC_MATH_RSQRTF16_H
+
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/properties/types.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+float16 rsqrtf16(float16 x);
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC_MATH_RSQRTF16_H
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index ebf9f1c86cf15..8daf59032622c 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -1560,6 +1560,17 @@ add_fp_unittest(
libc.src.math.sqrtl
)
+add_fp_unittest(
+ rsqrtf16_test
+ NEED_MPFR
+ SUITE
+ libc-math-unittests
+ SRCS
+ rsqrtf16_test.cpp
+ DEPENDS
+ libc.src.math.rsqrtf16
+)
+
add_fp_unittest(
sqrtf16_test
NEED_MPFR
diff --git a/libc/test/src/math/rsqrtf16_test.cpp b/libc/test/src/math/rsqrtf16_test.cpp
new file mode 100644
index 0000000000000..ec22df5bd31c7
--- /dev/null
+++ b/libc/test/src/math/rsqrtf16_test.cpp
@@ -0,0 +1,42 @@
+//===-- Exhaustive test for rsqrtf16 --------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/math/rsqrtf16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+using LlvmLibcRsqrtf16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+// Range: [0, Inf]
+static constexpr uint16_t POS_START = 0x0000U;
+static constexpr uint16_t POS_STOP = 0x7c00U;
+
+// Range: [-Inf, 0]
+static constexpr uint16_t NEG_START = 0x8000U;
+static constexpr uint16_t NEG_STOP = 0xfc00U;
+
+TEST_F(LlvmLibcRsqrtf16Test, PositiveRange) {
+ for (uint16_t v = POS_START; v <= POS_STOP; ++v) {
+ float16 x = FPBits(v).get_val();
+
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Rsqrt, x,
+ LIBC_NAMESPACE::rsqrtf16(x), 0.5);
+ }
+}
+
+TEST_F(LlvmLibcRsqrtf16Test, NegativeRange) {
+ for (uint16_t v = NEG_START; v <= NEG_STOP; ++v) {
+ float16 x = FPBits(v).get_val();
+
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Rsqrt, x,
+ LIBC_NAMESPACE::rsqrtf16(x), 0.5);
+ }
+}
\ No newline at end of file
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 79b189159e9d8..fe5d8a6b3b864 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -2958,6 +2958,17 @@ add_fp_unittest(
libc.src.math.sqrtl
)
+add_fp_unittest(
+ rsqrtf16_test
+ SUITE
+ libc-math-smoke-tests
+ SRCS
+ rsqrtf16_test.cpp
+ DEPENDS
+ libc.src.errno.errno
+ libc.src.math.rsqrtf16
+)
+
add_fp_unittest(
sqrtf16_test
SUITE
diff --git a/libc/test/src/math/smoke/rsqrtf16_test.cpp b/libc/test/src/math/smoke/rsqrtf16_test.cpp
new file mode 100644
index 0000000000000..eb429d3dbf6c8
--- /dev/null
+++ b/libc/test/src/math/smoke/rsqrtf16_test.cpp
@@ -0,0 +1,38 @@
+//===-- Unittests for rsqrtf16 --------------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception.
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/errno/libc_errno.h"
+#include "src/math/rsqrtf16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+using LlvmLibcRsqrtf16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+TEST_F(LlvmLibcRsqrtf16Test, SpecialNumbers) {
+ LIBC_NAMESPACE::libc_errno = 0;
+ EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::rsqrtf16(aNaN));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::rsqrtf16(sNaN),
+ FE_INVALID);
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::rsqrtf16(0.0f));
+ EXPECT_MATH_ERRNO(ERANGE);
+
+ EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::rsqrtf16(1.0f));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ(0.0f, LIBC_NAMESPACE::rsqrtf16(inf));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::rsqrtf16(neg_inf));
+ EXPECT_MATH_ERRNO(EDOM);
+
+ EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::rsqrtf16(-2.0f));
+ EXPECT_MATH_ERRNO(EDOM);
+}
\ No newline at end of file
diff --git a/libc/utils/MPFRWrapper/MPCommon.cpp b/libc/utils/MPFRWrapper/MPCommon.cpp
index ccd4d2d01a4e2..1a78ca5c24ba2 100644
--- a/libc/utils/MPFRWrapper/MPCommon.cpp
+++ b/libc/utils/MPFRWrapper/MPCommon.cpp
@@ -366,6 +366,12 @@ MPFRNumber MPFRNumber::rint(mpfr_rnd_t rnd) const {
return result;
}
+MPFRNumber MPFRNumber::rsqrt() const {
+ MPFRNumber result(*this);
+ mpfr_rec_sqrt(result.value, value, mpfr_rounding);
+ return result;
+}
+
MPFRNumber MPFRNumber::mod_2pi() const {
MPFRNumber result(0.0, 1280);
MPFRNumber _2pi(0.0, 1280);
diff --git a/libc/utils/MPFRWrapper/MPCommon.h b/libc/utils/MPFRWrapper/MPCommon.h
index 99cb7ec66a2ca..43218ee7662db 100644
--- a/libc/utils/MPFRWrapper/MPCommon.h
+++ b/libc/utils/MPFRWrapper/MPCommon.h
@@ -216,6 +216,7 @@ class MPFRNumber {
bool round_to_long(long &result) const;
bool round_to_long(mpfr_rnd_t rnd, long &result) const;
MPFRNumber rint(mpfr_rnd_t rnd) const;
+ MPFRNumber rsqrt() const;
MPFRNumber mod_2pi() const;
MPFRNumber mod_pi_over_2() const;
MPFRNumber mod_pi_over_4() const;
diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp
index 8853f96ef8f92..a68d21650f004 100644
--- a/libc/utils/MPFRWrapper/MPFRUtils.cpp
+++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp
@@ -86,6 +86,8 @@ unary_operation(Operation op, InputType input, unsigned int precision,
return mpfrInput.round();
case Operation::RoundEven:
return mpfrInput.roundeven();
+ case Operation::Rsqrt:
+ return mpfrInput.rsqrt();
case Operation::Sin:
return mpfrInput.sin();
case Operation::Sinpi:
diff --git a/libc/utils/MPFRWrapper/MPFRUtils.h b/libc/utils/MPFRWrapper/MPFRUtils.h
index c77a6aa3adeae..532b2c8d9d819 100644
--- a/libc/utils/MPFRWrapper/MPFRUtils.h
+++ b/libc/utils/MPFRWrapper/MPFRUtils.h
@@ -55,6 +55,7 @@ enum class Operation : int {
ModPIOver4,
Round,
RoundEven,
+ Rsqrt,
Sin,
Sinpi,
Sinh,
>From 10f24106d6947ce0c8408c527e81d9201589be10 Mon Sep 17 00:00:00 2001
From: amemov <shepelev777 at gmail.com>
Date: Sun, 27 Apr 2025 20:07:27 +0000
Subject: [PATCH 2/5] Clang-formated the files
---
libc/src/math/generic/rsqrtf16.cpp | 9 +++++----
libc/test/src/math/smoke/rsqrtf16_test.cpp | 3 +--
2 files changed, 6 insertions(+), 6 deletions(-)
diff --git a/libc/src/math/generic/rsqrtf16.cpp b/libc/src/math/generic/rsqrtf16.cpp
index 930c8f6fc7b4e..8e97ec896aa0c 100644
--- a/libc/src/math/generic/rsqrtf16.cpp
+++ b/libc/src/math/generic/rsqrtf16.cpp
@@ -35,7 +35,7 @@ LLVM_LIBC_FUNCTION(float16, rsqrtf16, (float16 x)) {
}
return x;
}
-
+
// |x| = 0
if (LIBC_UNLIKELY(x_abs == 0x0)) {
fputil::raise_except_if_required(FE_DIVBYZERO);
@@ -44,7 +44,7 @@ LLVM_LIBC_FUNCTION(float16, rsqrtf16, (float16 x)) {
}
// -inf <= x < 0
- if (LIBC_UNLIKELY(x_sign == 1)) {
+ if (LIBC_UNLIKELY(x_sign == 1)) {
fputil::raise_except_if_required(FE_INVALID);
fputil::set_errno_if_required(EDOM);
return FPBits::quiet_nan().get_val();
@@ -54,13 +54,14 @@ LLVM_LIBC_FUNCTION(float16, rsqrtf16, (float16 x)) {
if (LIBC_UNLIKELY(xbits.is_inf())) {
return fputil::cast<float16>(0.0f);
}
-
+
// x = 1 => rsqrt(x) = 1
if (LIBC_UNLIKELY(x_u == 0x1)) {
return fputil::cast<float16>(1.0f);
}
- // x is valid, estimate the result - below is temporary solution for just testing
+ // x is valid, estimate the result - below is temporary solution for just
+ // testing
float xf = x;
return fputil::cast<float16>(1.0f / xf);
}
diff --git a/libc/test/src/math/smoke/rsqrtf16_test.cpp b/libc/test/src/math/smoke/rsqrtf16_test.cpp
index eb429d3dbf6c8..d50d8ee0ad2ff 100644
--- a/libc/test/src/math/smoke/rsqrtf16_test.cpp
+++ b/libc/test/src/math/smoke/rsqrtf16_test.cpp
@@ -17,8 +17,7 @@ TEST_F(LlvmLibcRsqrtf16Test, SpecialNumbers) {
EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::rsqrtf16(aNaN));
EXPECT_MATH_ERRNO(0);
- EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::rsqrtf16(sNaN),
- FE_INVALID);
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::rsqrtf16(sNaN), FE_INVALID);
EXPECT_MATH_ERRNO(0);
EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::rsqrtf16(0.0f));
>From 8d15cfac9a39c497c951b9fdd6b73583d63020f6 Mon Sep 17 00:00:00 2001
From: amemov <shepelev777 at gmail.com>
Date: Sun, 27 Apr 2025 22:09:15 +0000
Subject: [PATCH 3/5] Replaced the computation for valid X with polynomial
approximation
---
libc/src/math/generic/rsqrtf16.cpp | 9 ++++++---
1 file changed, 6 insertions(+), 3 deletions(-)
diff --git a/libc/src/math/generic/rsqrtf16.cpp b/libc/src/math/generic/rsqrtf16.cpp
index 8e97ec896aa0c..2ae026499b89c 100644
--- a/libc/src/math/generic/rsqrtf16.cpp
+++ b/libc/src/math/generic/rsqrtf16.cpp
@@ -60,9 +60,12 @@ LLVM_LIBC_FUNCTION(float16, rsqrtf16, (float16 x)) {
return fputil::cast<float16>(1.0f);
}
- // x is valid, estimate the result - below is temporary solution for just
- // testing
+ // x is valid, estimate the result
+ // 3-degree polynomial generated using Sollya
+ // P = fpminimax(1/sqrt(x), [|1, 2, 3|], [|SG...|], [0.5, 1]);
float xf = x;
- return fputil::cast<float16>(1.0f / xf);
+ float result =
+ fputil::polyeval(xf, 0x1.d42408p2f, -0x1.7cc4fep3f, 0x1.66cb6ap2f);
+ return fputil::cast<float16>(result);
}
} // namespace LIBC_NAMESPACE_DECL
>From e49620105adb5eaab40a9c4a2cdade8469121bf3 Mon Sep 17 00:00:00 2001
From: amemov <shepelev777 at gmail.com>
Date: Tue, 29 Apr 2025 21:37:45 +0000
Subject: [PATCH 4/5] Added range reduction to the approximation
---
libc/src/math/generic/CMakeLists.txt | 2 +-
libc/src/math/generic/rsqrtf16.cpp | 39 ++++++++++++++++++----
libc/test/src/math/rsqrtf16_test.cpp | 2 +-
libc/test/src/math/smoke/rsqrtf16_test.cpp | 4 +--
4 files changed, 36 insertions(+), 11 deletions(-)
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index ae0977ce255c2..9335e7cbcb270 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -968,8 +968,8 @@ add_entrypoint_object(
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
libc.src.__support.FPUtil.multiply_add
+ libc.src.__support.FPUtil.manipulation_functions
libc.src.__support.FPUtil.polyeval
- libc.src.__support.FPUtil.sqrt
libc.src.__support.macros.optimization
libc.src.__support.macros.properties.types
)
diff --git a/libc/src/math/generic/rsqrtf16.cpp b/libc/src/math/generic/rsqrtf16.cpp
index 2ae026499b89c..5bf4ab3e46fda 100644
--- a/libc/src/math/generic/rsqrtf16.cpp
+++ b/libc/src/math/generic/rsqrtf16.cpp
@@ -11,10 +11,10 @@
#include "hdr/fenv_macros.h"
#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"
+#include "src/__support/FPUtil/multiply_add.h" // to remove
#include "src/__support/macros/optimization.h"
namespace LIBC_NAMESPACE_DECL {
@@ -40,7 +40,7 @@ LLVM_LIBC_FUNCTION(float16, rsqrtf16, (float16 x)) {
if (LIBC_UNLIKELY(x_abs == 0x0)) {
fputil::raise_except_if_required(FE_DIVBYZERO);
fputil::set_errno_if_required(ERANGE);
- return FPBits::quiet_nan().get_val();
+ return FPBits::inf(Sign::POS).get_val();
}
// -inf <= x < 0
@@ -61,11 +61,36 @@ LLVM_LIBC_FUNCTION(float16, rsqrtf16, (float16 x)) {
}
// x is valid, estimate the result
- // 3-degree polynomial generated using Sollya
- // P = fpminimax(1/sqrt(x), [|1, 2, 3|], [|SG...|], [0.5, 1]);
+ // 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 = x;
- float result =
- fputil::polyeval(xf, 0x1.d42408p2f, -0x1.7cc4fep3f, 0x1.66cb6ap2f);
+ int exponent;
+ float mantissa = fputil::frexp(xf, exponent);
+
+ // 6-degree polynomial generated using Sollya
+ // P = fpminimax(1/sqrt(x), [|0,1,2,3,4,5|], [|SG...|], [0.5, 1]);
+ float interm =
+ fputil::polyeval(mantissa, 0x1.9c81c4p1f, -0x1.e2c57cp2f, 0x1.91e8bp3f,
+ -0x1.899954p3f, 0x1.9edcp2f, -0x1.6bd93cp0f);
+
+ // Round (-e/2)
+ int exp_floored = -(exponent >> 1);
+
+ // rsqrt(x) = 1/sqrt(mantissa) * 2^(-e/2)
+ // rsqrt(x) = P(mantissa) * 2*(exp_floored)
+ float result = fputil::ldexp(interm, exp_floored);
+
+ // Handle the case where exponent is odd
+ if (exponent & 1) {
+ const float ONE_OVER_SQRT2 =
+ 0x1.6a09e667f3bcc908b2fb1366ea957d3e3adec1751p-1f;
+ result *= ONE_OVER_SQRT2;
+ }
+
return fputil::cast<float16>(result);
}
} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/test/src/math/rsqrtf16_test.cpp b/libc/test/src/math/rsqrtf16_test.cpp
index ec22df5bd31c7..d2f3fe8f49b92 100644
--- a/libc/test/src/math/rsqrtf16_test.cpp
+++ b/libc/test/src/math/rsqrtf16_test.cpp
@@ -39,4 +39,4 @@ TEST_F(LlvmLibcRsqrtf16Test, NegativeRange) {
EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Rsqrt, x,
LIBC_NAMESPACE::rsqrtf16(x), 0.5);
}
-}
\ No newline at end of file
+}
diff --git a/libc/test/src/math/smoke/rsqrtf16_test.cpp b/libc/test/src/math/smoke/rsqrtf16_test.cpp
index d50d8ee0ad2ff..8e69027e67e13 100644
--- a/libc/test/src/math/smoke/rsqrtf16_test.cpp
+++ b/libc/test/src/math/smoke/rsqrtf16_test.cpp
@@ -20,7 +20,7 @@ TEST_F(LlvmLibcRsqrtf16Test, SpecialNumbers) {
EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::rsqrtf16(sNaN), FE_INVALID);
EXPECT_MATH_ERRNO(0);
- EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::rsqrtf16(0.0f));
+ EXPECT_FP_EQ(inf, LIBC_NAMESPACE::rsqrtf16(0.0f));
EXPECT_MATH_ERRNO(ERANGE);
EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::rsqrtf16(1.0f));
@@ -34,4 +34,4 @@ TEST_F(LlvmLibcRsqrtf16Test, SpecialNumbers) {
EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::rsqrtf16(-2.0f));
EXPECT_MATH_ERRNO(EDOM);
-}
\ No newline at end of file
+}
>From ba9bdd2681f48dd353398c22501c8fa5390e83af Mon Sep 17 00:00:00 2001
From: amemov <shepelev777 at gmail.com>
Date: Tue, 20 May 2025 19:56:14 +0000
Subject: [PATCH 5/5] Added Newton-Raphson iterations -The accuracy improved
drastically, but it still fails
---
libc/src/math/generic/CMakeLists.txt | 2 +-
libc/src/math/generic/rsqrtf16.cpp | 35 +++++++++++++++++++---------
2 files changed, 25 insertions(+), 12 deletions(-)
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 9335e7cbcb270..11ffd8d69829f 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -967,7 +967,7 @@ add_entrypoint_object(
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.fma
libc.src.__support.FPUtil.manipulation_functions
libc.src.__support.FPUtil.polyeval
libc.src.__support.macros.optimization
diff --git a/libc/src/math/generic/rsqrtf16.cpp b/libc/src/math/generic/rsqrtf16.cpp
index 5bf4ab3e46fda..806328b60aca2 100644
--- a/libc/src/math/generic/rsqrtf16.cpp
+++ b/libc/src/math/generic/rsqrtf16.cpp
@@ -10,11 +10,11 @@
#include "hdr/errno_macros.h"
#include "hdr/fenv_macros.h"
#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FMA.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" // to remove
#include "src/__support/macros/optimization.h"
namespace LIBC_NAMESPACE_DECL {
@@ -55,11 +55,6 @@ LLVM_LIBC_FUNCTION(float16, rsqrtf16, (float16 x)) {
return fputil::cast<float16>(0.0f);
}
- // x = 1 => rsqrt(x) = 1
- if (LIBC_UNLIKELY(x_u == 0x1)) {
- return fputil::cast<float16>(1.0f);
- }
-
// x is valid, estimate the result
// Range reduction:
// x can be expressed as m*2^e, where e - int exponent and m - mantissa
@@ -72,23 +67,41 @@ LLVM_LIBC_FUNCTION(float16, rsqrtf16, (float16 x)) {
float mantissa = fputil::frexp(xf, exponent);
// 6-degree polynomial generated using Sollya
- // P = fpminimax(1/sqrt(x), [|0,1,2,3,4,5|], [|SG...|], [0.5, 1]);
+ // bigger polynomial doesn't generate better results-> the current one
+ // produces the least number of errors but still errors are presents P =
+ // fpminimax(1/(sqrt(x)), [|0,1,2,3,4,5|], [|SG...|], [0.5, 1]);
float interm =
fputil::polyeval(mantissa, 0x1.9c81c4p1f, -0x1.e2c57cp2f, 0x1.91e8bp3f,
-0x1.899954p3f, 0x1.9edcp2f, -0x1.6bd93cp0f);
+ // Apply one Newton-Raphson iteration to refine the approximation of
+ // 1/sqrt(mantissa) y_new = y_old * (1.5 - 0.5 * mantissa * y_old^2) Using
+ // fputil::fma for potential precision benefits in the factor calculation
+ float interm_sq = interm * interm;
+ float factor = fputil::fma<float>(-0.5f * mantissa, interm_sq, 1.5f);
+ float interm_refined = interm * factor; // Final multiplication
+
+ // Apply a second Newton-Raphson iteration
+ // y_new = y_old * (1.5 - 0.5 * mantissa * y_old^2)
+ // y_old is now interm_refined
+ float interm_refined_sq = interm_refined * interm_refined;
+ float factor2 = fputil::fma<float>(-0.5f * mantissa, interm_refined_sq, 1.5f);
+ float interm_refined2 = interm_refined * factor2;
+
// Round (-e/2)
int exp_floored = -(exponent >> 1);
// rsqrt(x) = 1/sqrt(mantissa) * 2^(-e/2)
// rsqrt(x) = P(mantissa) * 2*(exp_floored)
- float result = fputil::ldexp(interm, exp_floored);
+ // float result = fputil::ldexp(interm, exp_floored);
+ float result = fputil::ldexp(interm_refined2, exp_floored);
// Handle the case where exponent is odd
if (exponent & 1) {
- const float ONE_OVER_SQRT2 =
- 0x1.6a09e667f3bcc908b2fb1366ea957d3e3adec1751p-1f;
- result *= ONE_OVER_SQRT2;
+ const float ONE_OVER_SQRT2 = 0x1.6a09e6p-1f;
+ // result *= ONE_OVER_SQRT2;
+ result = fputil::fma<float>(result, ONE_OVER_SQRT2,
+ 0.0f); // Use FMA for multiplication
}
return fputil::cast<float16>(result);
More information about the libc-commits
mailing list