[libc-commits] [libc] [libc][stdfix] Add sqrt for fixed point types. (PR #83042)
via libc-commits
libc-commits at lists.llvm.org
Mon Feb 26 15:32:19 PST 2024
https://github.com/lntue updated https://github.com/llvm/llvm-project/pull/83042
>From 71a4e5221a676bfa709cee40260e5815ab35fee6 Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Mon, 26 Feb 2024 17:56:39 +0000
Subject: [PATCH 1/4] [libc][stdfix] Add sqrt for fixed point types.
---
libc/config/linux/x86_64/entrypoints.txt | 6 +
libc/docs/math/stdfix.rst | 2 +-
libc/spec/stdc_ext.td | 8 ++
libc/src/__support/fixed_point/CMakeLists.txt | 13 ++
libc/src/__support/fixed_point/sqrt.h | 131 ++++++++++++++++++
libc/src/stdfix/CMakeLists.txt | 15 ++
libc/src/stdfix/sqrtuhk.cpp | 19 +++
libc/src/stdfix/sqrtuhk.h | 20 +++
libc/src/stdfix/sqrtuhr.cpp | 19 +++
libc/src/stdfix/sqrtuhr.h | 20 +++
libc/src/stdfix/sqrtuk.cpp | 19 +++
libc/src/stdfix/sqrtuk.h | 20 +++
libc/src/stdfix/sqrtulr.cpp | 19 +++
libc/src/stdfix/sqrtulr.h | 20 +++
libc/src/stdfix/sqrtur.cpp | 19 +++
libc/src/stdfix/sqrtur.h | 20 +++
libc/test/src/stdfix/CMakeLists.txt | 22 +++
libc/test/src/stdfix/SqrtTest.h | 60 ++++++++
libc/test/src/stdfix/sqrtuhk_test.cpp | 13 ++
libc/test/src/stdfix/sqrtuhr_test.cpp | 13 ++
libc/test/src/stdfix/sqrtuk_test.cpp | 13 ++
libc/test/src/stdfix/sqrtulr_test.cpp | 13 ++
libc/test/src/stdfix/sqrtur_test.cpp | 13 ++
23 files changed, 516 insertions(+), 1 deletion(-)
create mode 100644 libc/src/__support/fixed_point/sqrt.h
create mode 100644 libc/src/stdfix/sqrtuhk.cpp
create mode 100644 libc/src/stdfix/sqrtuhk.h
create mode 100644 libc/src/stdfix/sqrtuhr.cpp
create mode 100644 libc/src/stdfix/sqrtuhr.h
create mode 100644 libc/src/stdfix/sqrtuk.cpp
create mode 100644 libc/src/stdfix/sqrtuk.h
create mode 100644 libc/src/stdfix/sqrtulr.cpp
create mode 100644 libc/src/stdfix/sqrtulr.h
create mode 100644 libc/src/stdfix/sqrtur.cpp
create mode 100644 libc/src/stdfix/sqrtur.h
create mode 100644 libc/test/src/stdfix/SqrtTest.h
create mode 100644 libc/test/src/stdfix/sqrtuhk_test.cpp
create mode 100644 libc/test/src/stdfix/sqrtuhr_test.cpp
create mode 100644 libc/test/src/stdfix/sqrtuk_test.cpp
create mode 100644 libc/test/src/stdfix/sqrtulr_test.cpp
create mode 100644 libc/test/src/stdfix/sqrtur_test.cpp
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index f2a224d45bbae7..cd44b5c52b58cd 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -473,6 +473,12 @@ if(LIBC_COMPILER_HAS_FIXED_POINT)
libc.src.stdfix.roundur
libc.src.stdfix.roundulk
libc.src.stdfix.roundulr
+ libc.src.stdfix.sqrtuhk
+ libc.src.stdfix.sqrtuhr
+ libc.src.stdfix.sqrtuk
+ libc.src.stdfix.sqrtur
+ # libc.src.stdfix.sqrtulk
+ libc.src.stdfix.sqrtulr
)
endif()
diff --git a/libc/docs/math/stdfix.rst b/libc/docs/math/stdfix.rst
index 080066e53bd2f1..79f499e61f121c 100644
--- a/libc/docs/math/stdfix.rst
+++ b/libc/docs/math/stdfix.rst
@@ -78,7 +78,7 @@ Fixed-point Arithmetics
+---------------+----------------+-------------+---------------+------------+----------------+-------------+----------------+-------------+---------------+------------+----------------+-------------+
| round | |check| | |check| | |check| | |check| | |check| | |check| | |check| | |check| | |check| | |check| | |check| | |check| |
+---------------+----------------+-------------+---------------+------------+----------------+-------------+----------------+-------------+---------------+------------+----------------+-------------+
-| sqrt | | | | | | | | | | | | |
+| sqrt | |check| | | |check| | | |check| | | |check| | | |check| | | | |
+---------------+----------------+-------------+---------------+------------+----------------+-------------+----------------+-------------+---------------+------------+----------------+-------------+
================== =========
diff --git a/libc/spec/stdc_ext.td b/libc/spec/stdc_ext.td
index 6620142146c471..be1e6d4ba2fcdc 100644
--- a/libc/spec/stdc_ext.td
+++ b/libc/spec/stdc_ext.td
@@ -47,6 +47,14 @@ def StdcExt : StandardSpec<"stdc_ext"> {
GuardedFunctionSpec<"rounduhk", RetValSpec<UnsignedShortAccumType>, [ArgSpec<UnsignedShortAccumType>, ArgSpec<IntType>], "LIBC_COMPILER_HAS_FIXED_POINT">,
GuardedFunctionSpec<"rounduk", RetValSpec<UnsignedAccumType>, [ArgSpec<UnsignedAccumType>, ArgSpec<IntType>], "LIBC_COMPILER_HAS_FIXED_POINT">,
GuardedFunctionSpec<"roundulk", RetValSpec<UnsignedLongAccumType>, [ArgSpec<UnsignedLongAccumType>, ArgSpec<IntType>], "LIBC_COMPILER_HAS_FIXED_POINT">,
+
+ GuardedFunctionSpec<"sqrtuhr", RetValSpec<UnsignedShortFractType>, [ArgSpec<UnsignedShortFractType>], "LIBC_COMPILER_HAS_FIXED_POINT">,
+ GuardedFunctionSpec<"sqrtur", RetValSpec<UnsignedFractType>, [ArgSpec<UnsignedFractType>], "LIBC_COMPILER_HAS_FIXED_POINT">,
+ GuardedFunctionSpec<"sqrtulr", RetValSpec<UnsignedLongFractType>, [ArgSpec<UnsignedLongFractType>], "LIBC_COMPILER_HAS_FIXED_POINT">,
+
+ GuardedFunctionSpec<"sqrtuhk", RetValSpec<UnsignedShortAccumType>, [ArgSpec<UnsignedShortAccumType>], "LIBC_COMPILER_HAS_FIXED_POINT">,
+ GuardedFunctionSpec<"sqrtuk", RetValSpec<UnsignedAccumType>, [ArgSpec<UnsignedAccumType>], "LIBC_COMPILER_HAS_FIXED_POINT">,
+ GuardedFunctionSpec<"sqrtulk", RetValSpec<UnsignedLongAccumType>, [ArgSpec<UnsignedLongAccumType>], "LIBC_COMPILER_HAS_FIXED_POINT">,
]
>;
diff --git a/libc/src/__support/fixed_point/CMakeLists.txt b/libc/src/__support/fixed_point/CMakeLists.txt
index 64f9dacc7ba5f0..0ed118f2408849 100644
--- a/libc/src/__support/fixed_point/CMakeLists.txt
+++ b/libc/src/__support/fixed_point/CMakeLists.txt
@@ -21,3 +21,16 @@ add_header_library(
libc.src.__support.CPP.bit
libc.src.__support.math_extras
)
+
+add_header_library(
+ sqrt
+ HDRS
+ sqrt.h
+ DEPENDS
+ .fx_rep
+ libc.include.llvm-libc-macros.stdfix_macros
+ libc.src.__support.macros.attributes
+ libc.src.__support.macros.optimization
+ libc.src.__support.CPP.bit
+ libc.src.__support.CPP.type_traits
+)
diff --git a/libc/src/__support/fixed_point/sqrt.h b/libc/src/__support/fixed_point/sqrt.h
new file mode 100644
index 00000000000000..c53baf08446409
--- /dev/null
+++ b/libc/src/__support/fixed_point/sqrt.h
@@ -0,0 +1,131 @@
+//===-- Calculate square root of fixed point numbers. -----*- 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___SUPPORT_FIXEDPOINT_SQRT_H
+#define LLVM_LIBC_SRC___SUPPORT_FIXEDPOINT_SQRT_H
+
+#include "include/llvm-libc-macros/stdfix-macros.h"
+#include "src/__support/CPP/bit.h"
+#include "src/__support/CPP/type_traits.h"
+#include "src/__support/macros/attributes.h" // LIBC_INLINE
+#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
+
+#include "fx_rep.h"
+
+#ifdef LIBC_COMPILER_HAS_FIXED_POINT
+
+namespace LIBC_NAMESPACE::fixed_point {
+
+namespace internal {
+
+template <typename T> struct SqrtConfig;
+
+template <> struct SqrtConfig<unsigned short fract> {
+ using Type = unsigned short fract;
+ static constexpr int EXTRA_STEPS = 0;
+};
+
+template <> struct SqrtConfig<unsigned fract> {
+ using Type = unsigned fract;
+ static constexpr int EXTRA_STEPS = 1;
+};
+
+template <> struct SqrtConfig<unsigned long fract> {
+ using Type = unsigned long fract;
+ static constexpr int EXTRA_STEPS = 2;
+};
+
+template <>
+struct SqrtConfig<unsigned short accum> : SqrtConfig<unsigned fract> {};
+
+template <>
+struct SqrtConfig<unsigned accum> : SqrtConfig<unsigned long fract> {};
+
+// TODO: unsigned long accum type is 64-bit, and will need 64-bit fract type.
+// Probably we will use DyadicFloat<64> for intermediate computations instead.
+
+// Linear approximation for the initial values, with errors bounded by:
+// max(1.5 * 2^-11, eps)
+// Generated with Sollya:
+// > for i from 4 to 15 do {
+// P = fpminimax(sqrt(x), 1, [|8, 8|], [i * 2^-4, (i + 1)*2^-4],
+// fixed, absolute);
+// print("{", coeff(P, 1), "uhr,", coeff(P, 0), "uhr},");
+// };
+static constexpr unsigned short fract SQRT_FIRST_APPROX[12][2] = {
+ {0x1.e8p-1uhr, 0x1.0cp-2uhr}, {0x1.bap-1uhr, 0x1.28p-2uhr},
+ {0x1.94p-1uhr, 0x1.44p-2uhr}, {0x1.74p-1uhr, 0x1.6p-2uhr},
+ {0x1.6p-1uhr, 0x1.74p-2uhr}, {0x1.4ep-1uhr, 0x1.88p-2uhr},
+ {0x1.3ep-1uhr, 0x1.9cp-2uhr}, {0x1.32p-1uhr, 0x1.acp-2uhr},
+ {0x1.22p-1uhr, 0x1.c4p-2uhr}, {0x1.18p-1uhr, 0x1.d4p-2uhr},
+ {0x1.08p-1uhr, 0x1.fp-2uhr}, {0x1.04p-1uhr, 0x1.f8p-2uhr},
+};
+
+} // namespace internal
+
+template <typename T>
+LIBC_INLINE constexpr cpp::enable_if_t<
+ cpp::is_fixed_point_v<T> && cpp::is_unsigned_v<T>, T>
+sqrt(T x) {
+ using BitType = typename FXRep<T>::StorageType;
+ BitType x_bit = cpp::bit_cast<BitType>(x);
+
+ if (LIBC_UNLIKELY(x_bit == 0))
+ return FXRep<T>::ZERO();
+
+ int leading_zeros = cpp::countl_zero(x_bit);
+ constexpr int STORAGE_LENGTH = sizeof(BitType) * CHAR_BIT;
+ constexpr int EXP_ADJUSTMENT = STORAGE_LENGTH - FXRep<T>::FRACTION_LEN - 1;
+ // x_exp is the real exponent of the leading bit of x.
+ int x_exp = EXP_ADJUSTMENT - leading_zeros;
+ int shift = EXP_ADJUSTMENT - 1 - (x_exp & (~1));
+ // Normalize.
+ x_bit <<= shift;
+ using FracType = typename internal::SqrtConfig<T>::Type;
+ FracType x_frac = cpp::bit_cast<FracType>(x_bit);
+
+ // Use use Newton method to approximate sqrt(a):
+ // x_{n + 1} = 1/2 (x_n + a / x_n)
+ // For the initial values, we choose x_0
+
+ // Use the leading 4 bits to do look up for sqrt(x).
+ int index = (x_bit >> (STORAGE_LENGTH - 4)) - 4;
+ FracType a = static_cast<FracType>(internal::SQRT_FIRST_APPROX[index][0]);
+ FracType b = static_cast<FracType>(internal::SQRT_FIRST_APPROX[index][1]);
+
+ // Initial approximation step.
+ // Estimated error bounds: | r - sqrt(x_frac) | < max(1.5 * 2^-11, eps).
+ FracType r = a * x_frac + b;
+
+ // Further Newton-method iterations for square-root:
+ // x_{n + 1} = 0.5 * (x_n + a / x_n)
+ // We distribute and do the multiplication by 0.5 first to avoid overflow.
+ // TODO: Investigate the performance and accuracy of using division-free
+ // iterations from:
+ // Blanchard, J. D. and Chamberland, M., "Newton's Method Without Division",
+ // The American Mathematical Monthly (2023).
+ // https://chamberland.math.grinnell.edu/papers/newton.pdf
+ for (int i = 0; i < internal::SqrtConfig<T>::EXTRA_STEPS; ++i) {
+ r = (r >> 1) + (x_frac >> 1) / r;
+ }
+
+ int r_exp = (x_exp >> 1);
+ int shift_back = EXP_ADJUSTMENT - r_exp;
+
+ // Re-scaling
+ r >>= shift_back;
+
+ // Return result.
+ return cpp::bit_cast<T>(r);
+}
+
+} // namespace LIBC_NAMESPACE::fixed_point
+
+#endif // LIBC_COMPILER_HAS_FIXED_POINT
+
+#endif // LLVM_LIBC_SRC___SUPPORT_FIXEDPOINT_SQRT_H
diff --git a/libc/src/stdfix/CMakeLists.txt b/libc/src/stdfix/CMakeLists.txt
index 6e2ed1bdfeafe7..cb2134fe33cf9c 100644
--- a/libc/src/stdfix/CMakeLists.txt
+++ b/libc/src/stdfix/CMakeLists.txt
@@ -17,6 +17,21 @@ foreach(suffix IN ITEMS hr r lr hk k lk)
)
endforeach()
+foreach(suffix IN ITEMS uhr ur ulr uhk uk)
+ add_entrypoint_object(
+ sqrt${suffix}
+ HDRS
+ sqrt${suffix}.h
+ SRCS
+ sqrt${suffix}.cpp
+ COMPILE_OPTIONS
+ -O3
+ -ffixed-point
+ DEPENDS
+ libc.src.__support.fixed_point.sqrt
+ )
+endforeach()
+
foreach(suffix IN ITEMS hr r lr hk k lk uhr ur ulr uhk uk ulk)
add_entrypoint_object(
round${suffix}
diff --git a/libc/src/stdfix/sqrtuhk.cpp b/libc/src/stdfix/sqrtuhk.cpp
new file mode 100644
index 00000000000000..e8dc842c8a9980
--- /dev/null
+++ b/libc/src/stdfix/sqrtuhk.cpp
@@ -0,0 +1,19 @@
+//===-- Implementation of sqrtuhk 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 "sqrtuhk.h"
+#include "src/__support/common.h"
+#include "src/__support/fixed_point/sqrt.h"
+
+namespace LIBC_NAMESPACE {
+
+LLVM_LIBC_FUNCTION(unsigned short accum, sqrtuhk, (unsigned short accum x)) {
+ return fixed_point::sqrt(x);
+}
+
+} // namespace LIBC_NAMESPACE
diff --git a/libc/src/stdfix/sqrtuhk.h b/libc/src/stdfix/sqrtuhk.h
new file mode 100644
index 00000000000000..80000a0079696d
--- /dev/null
+++ b/libc/src/stdfix/sqrtuhk.h
@@ -0,0 +1,20 @@
+//===-- Implementation header for sqrtuhk -----------------------*- 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_STDFIX_SQRTUHK_H
+#define LLVM_LIBC_SRC_STDFIX_SQRTUHK_H
+
+#include "include/llvm-libc-macros/stdfix-macros.h"
+
+namespace LIBC_NAMESPACE {
+
+unsigned short accum sqrtuhk(unsigned short accum x);
+
+} // namespace LIBC_NAMESPACE
+
+#endif // LLVM_LIBC_SRC_STDFIX_SQRTUHK_H
diff --git a/libc/src/stdfix/sqrtuhr.cpp b/libc/src/stdfix/sqrtuhr.cpp
new file mode 100644
index 00000000000000..6bba07aa20d595
--- /dev/null
+++ b/libc/src/stdfix/sqrtuhr.cpp
@@ -0,0 +1,19 @@
+//===-- Implementation of sqrtuhr 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 "sqrtuhr.h"
+#include "src/__support/common.h"
+#include "src/__support/fixed_point/sqrt.h"
+
+namespace LIBC_NAMESPACE {
+
+LLVM_LIBC_FUNCTION(unsigned short fract, sqrtuhr, (unsigned short fract x)) {
+ return fixed_point::sqrt(x);
+}
+
+} // namespace LIBC_NAMESPACE
diff --git a/libc/src/stdfix/sqrtuhr.h b/libc/src/stdfix/sqrtuhr.h
new file mode 100644
index 00000000000000..fd95f0924e8d48
--- /dev/null
+++ b/libc/src/stdfix/sqrtuhr.h
@@ -0,0 +1,20 @@
+//===-- Implementation header for sqrtuhr -----------------------*- 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_STDFIX_SQRTUHR_H
+#define LLVM_LIBC_SRC_STDFIX_SQRTUHR_H
+
+#include "include/llvm-libc-macros/stdfix-macros.h"
+
+namespace LIBC_NAMESPACE {
+
+unsigned short fract sqrtuhr(unsigned short fract x);
+
+} // namespace LIBC_NAMESPACE
+
+#endif // LLVM_LIBC_SRC_STDFIX_SQRTUHR_H
diff --git a/libc/src/stdfix/sqrtuk.cpp b/libc/src/stdfix/sqrtuk.cpp
new file mode 100644
index 00000000000000..6e5d8118c83b73
--- /dev/null
+++ b/libc/src/stdfix/sqrtuk.cpp
@@ -0,0 +1,19 @@
+//===-- Implementation of sqrtuk 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 "sqrtuk.h"
+#include "src/__support/common.h"
+#include "src/__support/fixed_point/sqrt.h"
+
+namespace LIBC_NAMESPACE {
+
+LLVM_LIBC_FUNCTION(unsigned accum, sqrtuk, (unsigned accum x)) {
+ return fixed_point::sqrt(x);
+}
+
+} // namespace LIBC_NAMESPACE
diff --git a/libc/src/stdfix/sqrtuk.h b/libc/src/stdfix/sqrtuk.h
new file mode 100644
index 00000000000000..04d0adadde9ad2
--- /dev/null
+++ b/libc/src/stdfix/sqrtuk.h
@@ -0,0 +1,20 @@
+//===-- Implementation header for sqrtuk ------------------------*- 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_STDFIX_SQRTUK_H
+#define LLVM_LIBC_SRC_STDFIX_SQRTUK_H
+
+#include "include/llvm-libc-macros/stdfix-macros.h"
+
+namespace LIBC_NAMESPACE {
+
+unsigned accum sqrtuk(unsigned accum x);
+
+} // namespace LIBC_NAMESPACE
+
+#endif // LLVM_LIBC_SRC_STDFIX_SQRTUK_H
diff --git a/libc/src/stdfix/sqrtulr.cpp b/libc/src/stdfix/sqrtulr.cpp
new file mode 100644
index 00000000000000..c9e5cd51f66bc5
--- /dev/null
+++ b/libc/src/stdfix/sqrtulr.cpp
@@ -0,0 +1,19 @@
+//===-- Implementation of sqrtulr 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 "sqrtulr.h"
+#include "src/__support/common.h"
+#include "src/__support/fixed_point/sqrt.h"
+
+namespace LIBC_NAMESPACE {
+
+LLVM_LIBC_FUNCTION(unsigned long fract, sqrtulr, (unsigned long fract x)) {
+ return fixed_point::sqrt(x);
+}
+
+} // namespace LIBC_NAMESPACE
diff --git a/libc/src/stdfix/sqrtulr.h b/libc/src/stdfix/sqrtulr.h
new file mode 100644
index 00000000000000..284adaaf35bf59
--- /dev/null
+++ b/libc/src/stdfix/sqrtulr.h
@@ -0,0 +1,20 @@
+//===-- Implementation header for sqrtulr -----------------------*- 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_STDFIX_SQRTULR_H
+#define LLVM_LIBC_SRC_STDFIX_SQRTULR_H
+
+#include "include/llvm-libc-macros/stdfix-macros.h"
+
+namespace LIBC_NAMESPACE {
+
+unsigned long fract sqrtulr(unsigned long fract x);
+
+} // namespace LIBC_NAMESPACE
+
+#endif // LLVM_LIBC_SRC_STDFIX_SQRTULR_H
diff --git a/libc/src/stdfix/sqrtur.cpp b/libc/src/stdfix/sqrtur.cpp
new file mode 100644
index 00000000000000..ac5be84910849f
--- /dev/null
+++ b/libc/src/stdfix/sqrtur.cpp
@@ -0,0 +1,19 @@
+//===-- Implementation of sqrtur 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 "sqrtur.h"
+#include "src/__support/common.h"
+#include "src/__support/fixed_point/sqrt.h"
+
+namespace LIBC_NAMESPACE {
+
+LLVM_LIBC_FUNCTION(unsigned fract, sqrtur, (unsigned fract x)) {
+ return fixed_point::sqrt(x);
+}
+
+} // namespace LIBC_NAMESPACE
diff --git a/libc/src/stdfix/sqrtur.h b/libc/src/stdfix/sqrtur.h
new file mode 100644
index 00000000000000..df9dfe5a0bf39e
--- /dev/null
+++ b/libc/src/stdfix/sqrtur.h
@@ -0,0 +1,20 @@
+//===-- Implementation header for sqrtur ------------------------*- 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_STDFIX_SQRTUR_H
+#define LLVM_LIBC_SRC_STDFIX_SQRTUR_H
+
+#include "include/llvm-libc-macros/stdfix-macros.h"
+
+namespace LIBC_NAMESPACE {
+
+unsigned fract sqrtur(unsigned fract x);
+
+} // namespace LIBC_NAMESPACE
+
+#endif // LLVM_LIBC_SRC_STDFIX_SQRTUR_H
diff --git a/libc/test/src/stdfix/CMakeLists.txt b/libc/test/src/stdfix/CMakeLists.txt
index b6e0256bb68800..4140b5b29f3b3c 100644
--- a/libc/test/src/stdfix/CMakeLists.txt
+++ b/libc/test/src/stdfix/CMakeLists.txt
@@ -22,6 +22,28 @@ foreach(suffix IN ITEMS hr r lr hk k lk)
)
endforeach()
+foreach(suffix IN ITEMS uhr ur ulr uhk uk)
+ add_libc_test(
+ sqrt${suffix}_test
+ SUITE
+ libc-stdfix-tests
+ HDRS
+ SqrtTest.h
+ SRCS
+ sqrt${suffix}_test.cpp
+ COMPILE_OPTIONS
+ -O3
+ -ffixed-point
+ DEPENDS
+ libc.src.stdfix.sqrt${suffix}
+ libc.src.__support.CPP.bit
+ libc.src.__support.fixed_point.fx_rep
+ libc.src.__support.fixed_point.sqrt
+ libc.src.__support.FPUtil.basic_operations
+ libc.src.__support.FPUtil.sqrt
+ )
+endforeach()
+
foreach(suffix IN ITEMS hr r lr hk k lk uhr ur ulr uhk uk ulk)
add_libc_test(
round${suffix}_test
diff --git a/libc/test/src/stdfix/SqrtTest.h b/libc/test/src/stdfix/SqrtTest.h
new file mode 100644
index 00000000000000..ed4a9dee017e29
--- /dev/null
+++ b/libc/test/src/stdfix/SqrtTest.h
@@ -0,0 +1,60 @@
+//===-- Utility class to test fixed-point sqrt ------------------*- 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
+//
+//===----------------------------------------------------------------------===//
+
+#include "test/UnitTest/Test.h"
+
+#include "src/__support/CPP/bit.h"
+#include "src/__support/FPUtil/BasicOperations.h"
+#include "src/__support/FPUtil/sqrt.h"
+#include "src/__support/fixed_point/fx_rep.h"
+#include "src/__support/fixed_point/sqrt.h"
+
+template <typename T> class SqrtTest : public LIBC_NAMESPACE::testing::Test {
+
+ using FXRep = LIBC_NAMESPACE::fixed_point::FXRep<T>;
+ static constexpr T zero = FXRep::ZERO();
+ static constexpr T min = FXRep::MIN();
+ static constexpr T max = FXRep::MAX();
+ static constexpr T half = static_cast<T>(0.5);
+ static constexpr T quarter = static_cast<T>(0.25);
+ static constexpr T one =
+ (FXRep::INTEGRAL_LEN > 0) ? static_cast<T>(1) : FXRep::MAX();
+ static constexpr T eps = FXRep::EPS();
+
+public:
+ typedef T (*SqrtFunc)(T);
+
+ void testSpecialNumbers(SqrtFunc func) {
+ EXPECT_EQ(zero, func(zero));
+ EXPECT_EQ(half, func(quarter));
+
+ if constexpr (FXRep::INTEGRAL_LEN) {
+ EXPECT_EQ(one, func(one));
+ EXPECT_EQ(static_cast<T>(2.0), func(static_cast<T>(4.0)));
+ }
+
+ using StorageType = FXRep<T>::StorageType;
+
+ constexpr size_t COUNT = 255;
+ constexpr StorageType STEP =
+ ~StorageType(0) / static_cast<StorageType>(COUNT);
+ constexpr double ERR = static_cast<double>(eps);
+ for (size_t i = 0, StorageType x = 0; i < COUNT; ++i, x += STEP) {
+ T v = cpp::bit_cast<T>(x);
+ double v_d = static_cast<double>(v);
+ double errors = LIBC_NAMESPACE::fputil::abs(
+ static_cast<double>(func(v)) - LIBC_NAMESPACE::fputil::sqrt(v_d));
+ ASSERT_LE(errors, ERR);
+ }
+ }
+};
+
+#define LIST_SQRT_TESTS(T, func) \
+ using LlvmLibcSqrtTest = SqrtTest<T>; \
+ TEST_F(LlvmLibcSqrtTest, SpecialNumbers) { testSpecialNumbers(&func); } \
+ static_assert(true, "Require semicolon.")
diff --git a/libc/test/src/stdfix/sqrtuhk_test.cpp b/libc/test/src/stdfix/sqrtuhk_test.cpp
new file mode 100644
index 00000000000000..d6ff5385cab863
--- /dev/null
+++ b/libc/test/src/stdfix/sqrtuhk_test.cpp
@@ -0,0 +1,13 @@
+//===-- Unittests for sqrtuhk ---------------------------------------------===//
+//
+// 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 "SqrtTest.h"
+
+#include "src/stdfix/sqrtuhk.h"
+
+LIST_SQRT_TESTS(unsigned short accum, LIBC_NAMESPACE::sqrtuhk);
diff --git a/libc/test/src/stdfix/sqrtuhr_test.cpp b/libc/test/src/stdfix/sqrtuhr_test.cpp
new file mode 100644
index 00000000000000..22f00a4231b312
--- /dev/null
+++ b/libc/test/src/stdfix/sqrtuhr_test.cpp
@@ -0,0 +1,13 @@
+//===-- Unittests for sqrtuhr ---------------------------------------------===//
+//
+// 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 "SqrtTest.h"
+
+#include "src/stdfix/sqrtuhr.h"
+
+LIST_SQRT_TESTS(unsigned short fract, LIBC_NAMESPACE::sqrtuhr);
diff --git a/libc/test/src/stdfix/sqrtuk_test.cpp b/libc/test/src/stdfix/sqrtuk_test.cpp
new file mode 100644
index 00000000000000..5a3105de1e0bff
--- /dev/null
+++ b/libc/test/src/stdfix/sqrtuk_test.cpp
@@ -0,0 +1,13 @@
+//===-- Unittests for sqrtuk ----------------------------------------------===//
+//
+// 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 "SqrtTest.h"
+
+#include "src/stdfix/sqrtuk.h"
+
+LIST_SQRT_TESTS(unsigned accum, LIBC_NAMESPACE::sqrtuk);
diff --git a/libc/test/src/stdfix/sqrtulr_test.cpp b/libc/test/src/stdfix/sqrtulr_test.cpp
new file mode 100644
index 00000000000000..1be4e2b5e0a6f9
--- /dev/null
+++ b/libc/test/src/stdfix/sqrtulr_test.cpp
@@ -0,0 +1,13 @@
+//===-- Unittests for sqrtulr ---------------------------------------------===//
+//
+// 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 "SqrtTest.h"
+
+#include "src/stdfix/sqrtulr.h"
+
+LIST_SQRT_TESTS(unsigned long fract, LIBC_NAMESPACE::sqrtulr);
diff --git a/libc/test/src/stdfix/sqrtur_test.cpp b/libc/test/src/stdfix/sqrtur_test.cpp
new file mode 100644
index 00000000000000..12b1c2211db054
--- /dev/null
+++ b/libc/test/src/stdfix/sqrtur_test.cpp
@@ -0,0 +1,13 @@
+//===-- Unittests for sqrtur ----------------------------------------------===//
+//
+// 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 "SqrtTest.h"
+
+#include "src/stdfix/sqrtur.h"
+
+LIST_SQRT_TESTS(unsigned fract, LIBC_NAMESPACE::sqrtur);
>From 17182015e9b6a43479ccf572e613713fde686a25 Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Mon, 26 Feb 2024 19:20:01 +0000
Subject: [PATCH 2/4] Fix issues with fixed point arrays.
---
libc/src/__support/fixed_point/sqrt.h | 30 ++++++++++++++++-----------
libc/test/src/stdfix/SqrtTest.h | 16 +++++++++-----
2 files changed, 29 insertions(+), 17 deletions(-)
diff --git a/libc/src/__support/fixed_point/sqrt.h b/libc/src/__support/fixed_point/sqrt.h
index c53baf08446409..b7e2153a0827c9 100644
--- a/libc/src/__support/fixed_point/sqrt.h
+++ b/libc/src/__support/fixed_point/sqrt.h
@@ -57,21 +57,25 @@ struct SqrtConfig<unsigned accum> : SqrtConfig<unsigned long fract> {};
// fixed, absolute);
// print("{", coeff(P, 1), "uhr,", coeff(P, 0), "uhr},");
// };
-static constexpr unsigned short fract SQRT_FIRST_APPROX[12][2] = {
- {0x1.e8p-1uhr, 0x1.0cp-2uhr}, {0x1.bap-1uhr, 0x1.28p-2uhr},
- {0x1.94p-1uhr, 0x1.44p-2uhr}, {0x1.74p-1uhr, 0x1.6p-2uhr},
- {0x1.6p-1uhr, 0x1.74p-2uhr}, {0x1.4ep-1uhr, 0x1.88p-2uhr},
- {0x1.3ep-1uhr, 0x1.9cp-2uhr}, {0x1.32p-1uhr, 0x1.acp-2uhr},
- {0x1.22p-1uhr, 0x1.c4p-2uhr}, {0x1.18p-1uhr, 0x1.d4p-2uhr},
- {0x1.08p-1uhr, 0x1.fp-2uhr}, {0x1.04p-1uhr, 0x1.f8p-2uhr},
+// Clang is having issues with this array:
+// static constexpr unsigned short fract SQRT_FIRST_APPROX[12][2] = {
+// {0x1.e8p-1uhr, 0x1.0cp-2uhr}, {0x1.bap-1uhr, 0x1.28p-2uhr},
+// {0x1.94p-1uhr, 0x1.44p-2uhr}, {0x1.74p-1uhr, 0x1.6p-2uhr},
+// {0x1.6p-1uhr, 0x1.74p-2uhr}, {0x1.4ep-1uhr, 0x1.88p-2uhr},
+// {0x1.3ep-1uhr, 0x1.9cp-2uhr}, {0x1.32p-1uhr, 0x1.acp-2uhr},
+// {0x1.22p-1uhr, 0x1.c4p-2uhr}, {0x1.18p-1uhr, 0x1.d4p-2uhr},
+// {0x1.08p-1uhr, 0x1.fp-2uhr}, {0x1.04p-1uhr, 0x1.f8p-2uhr},
+// };
+// We are using their storage type instead.
+static constexpr uint8_t SQRT_FIRST_APPROX[12][2] = {
+ {244, 67}, {221, 74}, {202, 81}, {186, 88}, {176, 93}, {167, 98},
+ {159, 103}, {153, 107}, {145, 113}, {140, 117}, {132, 124}, {130, 126},
};
} // namespace internal
template <typename T>
-LIBC_INLINE constexpr cpp::enable_if_t<
- cpp::is_fixed_point_v<T> && cpp::is_unsigned_v<T>, T>
-sqrt(T x) {
+LIBC_INLINE constexpr cpp::enable_if_t<cpp::is_fixed_point_v<T>, T> sqrt(T x) {
using BitType = typename FXRep<T>::StorageType;
BitType x_bit = cpp::bit_cast<BitType>(x);
@@ -95,8 +99,10 @@ sqrt(T x) {
// Use the leading 4 bits to do look up for sqrt(x).
int index = (x_bit >> (STORAGE_LENGTH - 4)) - 4;
- FracType a = static_cast<FracType>(internal::SQRT_FIRST_APPROX[index][0]);
- FracType b = static_cast<FracType>(internal::SQRT_FIRST_APPROX[index][1]);
+ FracType a = static_cast<FracType>(cpp::bit_cast<unsigned short fract>(
+ internal::SQRT_FIRST_APPROX[index][0]));
+ FracType b = static_cast<FracType>(cpp::bit_cast<unsigned short fract>(
+ internal::SQRT_FIRST_APPROX[index][1]));
// Initial approximation step.
// Estimated error bounds: | r - sqrt(x_frac) | < max(1.5 * 2^-11, eps).
diff --git a/libc/test/src/stdfix/SqrtTest.h b/libc/test/src/stdfix/SqrtTest.h
index ed4a9dee017e29..628be0deb770fb 100644
--- a/libc/test/src/stdfix/SqrtTest.h
+++ b/libc/test/src/stdfix/SqrtTest.h
@@ -38,18 +38,24 @@ template <typename T> class SqrtTest : public LIBC_NAMESPACE::testing::Test {
EXPECT_EQ(static_cast<T>(2.0), func(static_cast<T>(4.0)));
}
- using StorageType = FXRep<T>::StorageType;
+ using StorageType = typename FXRep::StorageType;
constexpr size_t COUNT = 255;
constexpr StorageType STEP =
~StorageType(0) / static_cast<StorageType>(COUNT);
- constexpr double ERR = static_cast<double>(eps);
- for (size_t i = 0, StorageType x = 0; i < COUNT; ++i, x += STEP) {
- T v = cpp::bit_cast<T>(x);
+ constexpr double ERR = 3.0 * static_cast<double>(eps);
+ StorageType x = 0;
+ for (size_t i = 0; i < COUNT; ++i, x += STEP) {
+ T v = LIBC_NAMESPACE::cpp::bit_cast<T>(x);
double v_d = static_cast<double>(v);
double errors = LIBC_NAMESPACE::fputil::abs(
static_cast<double>(func(v)) - LIBC_NAMESPACE::fputil::sqrt(v_d));
- ASSERT_LE(errors, ERR);
+ if (errors > ERR) {
+ // Print out the failure input and output.
+ EXPECT_EQ(v, zero);
+ EXPECT_EQ(func(v), zero);
+ }
+ ASSERT_TRUE(errors <= ERR);
}
}
};
>From 7f7759e19780fa90183f5e3460106bd780924824 Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Mon, 26 Feb 2024 19:26:49 +0000
Subject: [PATCH 3/4] Add link to the issue.
---
libc/src/__support/fixed_point/sqrt.h | 2 ++
1 file changed, 2 insertions(+)
diff --git a/libc/src/__support/fixed_point/sqrt.h b/libc/src/__support/fixed_point/sqrt.h
index b7e2153a0827c9..dcf2c715497130 100644
--- a/libc/src/__support/fixed_point/sqrt.h
+++ b/libc/src/__support/fixed_point/sqrt.h
@@ -67,6 +67,8 @@ struct SqrtConfig<unsigned accum> : SqrtConfig<unsigned long fract> {};
// {0x1.08p-1uhr, 0x1.fp-2uhr}, {0x1.04p-1uhr, 0x1.f8p-2uhr},
// };
// We are using their storage type instead.
+// TODO(https://github.com/llvm/llvm-project/issues/83050): Use fixed point
+// array when the bug is fixed.
static constexpr uint8_t SQRT_FIRST_APPROX[12][2] = {
{244, 67}, {221, 74}, {202, 81}, {186, 88}, {176, 93}, {167, 98},
{159, 103}, {153, 107}, {145, 113}, {140, 117}, {132, 124}, {130, 126},
>From 5a604c1966598782a0f5036899f9aa3845956817 Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Mon, 26 Feb 2024 23:31:49 +0000
Subject: [PATCH 4/4] Address comments.
---
libc/src/__support/fixed_point/sqrt.h | 40 ++++++++++-----------------
1 file changed, 15 insertions(+), 25 deletions(-)
diff --git a/libc/src/__support/fixed_point/sqrt.h b/libc/src/__support/fixed_point/sqrt.h
index dcf2c715497130..d8df294b18a1a8 100644
--- a/libc/src/__support/fixed_point/sqrt.h
+++ b/libc/src/__support/fixed_point/sqrt.h
@@ -57,21 +57,13 @@ struct SqrtConfig<unsigned accum> : SqrtConfig<unsigned long fract> {};
// fixed, absolute);
// print("{", coeff(P, 1), "uhr,", coeff(P, 0), "uhr},");
// };
-// Clang is having issues with this array:
-// static constexpr unsigned short fract SQRT_FIRST_APPROX[12][2] = {
-// {0x1.e8p-1uhr, 0x1.0cp-2uhr}, {0x1.bap-1uhr, 0x1.28p-2uhr},
-// {0x1.94p-1uhr, 0x1.44p-2uhr}, {0x1.74p-1uhr, 0x1.6p-2uhr},
-// {0x1.6p-1uhr, 0x1.74p-2uhr}, {0x1.4ep-1uhr, 0x1.88p-2uhr},
-// {0x1.3ep-1uhr, 0x1.9cp-2uhr}, {0x1.32p-1uhr, 0x1.acp-2uhr},
-// {0x1.22p-1uhr, 0x1.c4p-2uhr}, {0x1.18p-1uhr, 0x1.d4p-2uhr},
-// {0x1.08p-1uhr, 0x1.fp-2uhr}, {0x1.04p-1uhr, 0x1.f8p-2uhr},
-// };
-// We are using their storage type instead.
-// TODO(https://github.com/llvm/llvm-project/issues/83050): Use fixed point
-// array when the bug is fixed.
-static constexpr uint8_t SQRT_FIRST_APPROX[12][2] = {
- {244, 67}, {221, 74}, {202, 81}, {186, 88}, {176, 93}, {167, 98},
- {159, 103}, {153, 107}, {145, 113}, {140, 117}, {132, 124}, {130, 126},
+static constexpr unsigned short fract SQRT_FIRST_APPROX[12][2] = {
+ {0x1.e8p-1uhr, 0x1.0cp-2uhr}, {0x1.bap-1uhr, 0x1.28p-2uhr},
+ {0x1.94p-1uhr, 0x1.44p-2uhr}, {0x1.74p-1uhr, 0x1.6p-2uhr},
+ {0x1.6p-1uhr, 0x1.74p-2uhr}, {0x1.4ep-1uhr, 0x1.88p-2uhr},
+ {0x1.3ep-1uhr, 0x1.9cp-2uhr}, {0x1.32p-1uhr, 0x1.acp-2uhr},
+ {0x1.22p-1uhr, 0x1.c4p-2uhr}, {0x1.18p-1uhr, 0x1.d4p-2uhr},
+ {0x1.08p-1uhr, 0x1.fp-2uhr}, {0x1.04p-1uhr, 0x1.f8p-2uhr},
};
} // namespace internal
@@ -100,11 +92,13 @@ LIBC_INLINE constexpr cpp::enable_if_t<cpp::is_fixed_point_v<T>, T> sqrt(T x) {
// For the initial values, we choose x_0
// Use the leading 4 bits to do look up for sqrt(x).
+ // After normalization, 0.25 <= x_frac < 1, so the leading 4 bits of x_frac
+ // are between 0b0100 and 0b1111. Hence the lookup table only needs 12
+ // entries, and we can get the index by subtracting the leading 4 bits of
+ // x_frac by 4 = 0b0100.
int index = (x_bit >> (STORAGE_LENGTH - 4)) - 4;
- FracType a = static_cast<FracType>(cpp::bit_cast<unsigned short fract>(
- internal::SQRT_FIRST_APPROX[index][0]));
- FracType b = static_cast<FracType>(cpp::bit_cast<unsigned short fract>(
- internal::SQRT_FIRST_APPROX[index][1]));
+ FracType a = static_cast<FracType>(internal::SQRT_FIRST_APPROX[index][0]);
+ FracType b = static_cast<FracType>(internal::SQRT_FIRST_APPROX[index][1]);
// Initial approximation step.
// Estimated error bounds: | r - sqrt(x_frac) | < max(1.5 * 2^-11, eps).
@@ -118,15 +112,11 @@ LIBC_INLINE constexpr cpp::enable_if_t<cpp::is_fixed_point_v<T>, T> sqrt(T x) {
// Blanchard, J. D. and Chamberland, M., "Newton's Method Without Division",
// The American Mathematical Monthly (2023).
// https://chamberland.math.grinnell.edu/papers/newton.pdf
- for (int i = 0; i < internal::SqrtConfig<T>::EXTRA_STEPS; ++i) {
+ for (int i = 0; i < internal::SqrtConfig<T>::EXTRA_STEPS; ++i)
r = (r >> 1) + (x_frac >> 1) / r;
- }
-
- int r_exp = (x_exp >> 1);
- int shift_back = EXP_ADJUSTMENT - r_exp;
// Re-scaling
- r >>= shift_back;
+ r >>= EXP_ADJUSTMENT - (x_exp >> 1);
// Return result.
return cpp::bit_cast<T>(r);
More information about the libc-commits
mailing list