[libc] Add implementation of hypot.
Refactor src/math/hypotf.cpp and test/src/math/hypotf_test.cpp and reuse them for hypot and hypot_test
Differential Revision: https://reviews.llvm.org/D91831
GitOrigin-RevId: 3b487d51e2ec699c27387fc30374f0d035b2a482
diff --git a/config/linux/aarch64/entrypoints.txt b/config/linux/aarch64/entrypoints.txt
index 1d8e5dd..3a3b050 100644
--- a/config/linux/aarch64/entrypoints.txt
+++ b/config/linux/aarch64/entrypoints.txt
@@ -68,6 +68,7 @@
libc.src.math.frexp
libc.src.math.frexpf
libc.src.math.frexpl
+ libc.src.math.hypot
libc.src.math.hypotf
libc.src.math.ilogb
libc.src.math.ilogbf
diff --git a/config/linux/x86_64/entrypoints.txt b/config/linux/x86_64/entrypoints.txt
index d6d56f2..7401715 100644
--- a/config/linux/x86_64/entrypoints.txt
+++ b/config/linux/x86_64/entrypoints.txt
@@ -104,6 +104,7 @@
libc.src.math.frexp
libc.src.math.frexpf
libc.src.math.frexpl
+ libc.src.math.hypot
libc.src.math.hypotf
libc.src.math.ilogb
libc.src.math.ilogbf
diff --git a/spec/stdc.td b/spec/stdc.td
index 7275f1e..d40fe8d 100644
--- a/spec/stdc.td
+++ b/spec/stdc.td
@@ -280,6 +280,7 @@
FunctionSpec<"frexpf", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<IntPtr>]>,
FunctionSpec<"frexpl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>, ArgSpec<IntPtr>]>,
+ FunctionSpec<"hypot", RetValSpec<DoubleType>, [ArgSpec<DoubleType>, ArgSpec<DoubleType>]>,
FunctionSpec<"hypotf", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<FloatType>]>,
FunctionSpec<"ilogb", RetValSpec<IntType>, [ArgSpec<DoubleType>]>,
diff --git a/src/math/CMakeLists.txt b/src/math/CMakeLists.txt
index fd75a3b..8201d73 100644
--- a/src/math/CMakeLists.txt
+++ b/src/math/CMakeLists.txt
@@ -713,3 +713,15 @@
COMPILE_OPTIONS
-O2
)
+
+add_entrypoint_object(
+ hypot
+ SRCS
+ hypot.cpp
+ HDRS
+ hypot.h
+ DEPENDS
+ libc.utils.FPUtil.fputil
+ COMPILE_OPTIONS
+ -O2
+)
diff --git a/src/math/hypot.cpp b/src/math/hypot.cpp
new file mode 100644
index 0000000..9d59365
--- /dev/null
+++ b/src/math/hypot.cpp
@@ -0,0 +1,18 @@
+//===-- Implementation of hypot 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 "utils/FPUtil/Hypot.h"
+#include "src/__support/common.h"
+
+namespace __llvm_libc {
+
+double LLVM_LIBC_ENTRYPOINT(hypot)(double x, double y) {
+ return __llvm_libc::fputil::hypot(x, y);
+}
+
+} // namespace __llvm_libc
diff --git a/src/math/hypot.h b/src/math/hypot.h
new file mode 100644
index 0000000..6c901ee
--- /dev/null
+++ b/src/math/hypot.h
@@ -0,0 +1,18 @@
+//===-- Implementation header for hypot -------------------------*- 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_HYPOT_H
+#define LLVM_LIBC_SRC_MATH_HYPOT_H
+
+namespace __llvm_libc {
+
+double hypot(double x, double y);
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_HYPOT_H
diff --git a/src/math/hypotf.cpp b/src/math/hypotf.cpp
index 10ebbb1..ebe7e97 100644
--- a/src/math/hypotf.cpp
+++ b/src/math/hypotf.cpp
@@ -6,217 +6,12 @@
//
//===----------------------------------------------------------------------===//
#include "src/__support/common.h"
-#include "utils/FPUtil/BasicOperations.h"
-#include "utils/FPUtil/FPBits.h"
+#include "utils/FPUtil/Hypot.h"
namespace __llvm_libc {
-using namespace fputil;
-
-uint32_t findLeadingOne(uint32_t mant, int &shift_length) {
- shift_length = 0;
- constexpr int nsteps = 5;
- constexpr uint32_t bounds[nsteps] = {1 << 16, 1 << 8, 1 << 4, 1 << 2, 1 << 1};
- constexpr int shifts[nsteps] = {16, 8, 4, 2, 1};
- for (int i = 0; i < nsteps; ++i) {
- if (mant >= bounds[i]) {
- shift_length += shifts[i];
- mant >>= shifts[i];
- }
- }
- return 1U << shift_length;
-}
-
-// Correctly rounded IEEE 754 HYPOT(x, y) with round to nearest, ties to even.
-//
-// Algorithm:
-// - Let a = max(|x|, |y|), b = min(|x|, |y|), then we have that:
-// a <= sqrt(a^2 + b^2) <= min(a + b, a*sqrt(2))
-// 1. So if b < eps(a)/2, then HYPOT(x, y) = a.
-//
-// - Moreover, the exponent part of HYPOT(x, y) is either the same or 1 more
-// than the exponent part of a.
-//
-// 2. For the remaining cases, we will use the digit-by-digit (shift-and-add)
-// algorithm to compute SQRT(Z):
-//
-// - For Y = y0.y1...yn... = SQRT(Z),
-// let Y(n) = y0.y1...yn be the first n fractional digits of Y.
-//
-// - The nth scaled residual R(n) is defined to be:
-// R(n) = 2^n * (Z - Y(n)^2)
-//
-// - Since Y(n) = Y(n - 1) + yn * 2^(-n), the scaled residual
-// satisfies the following recurrence formula:
-// R(n) = 2*R(n - 1) - yn*(2*Y(n - 1) + 2^(-n)),
-// with the initial conditions:
-// Y(0) = y0, and R(0) = Z - y0.
-//
-// - So the nth fractional digit of Y = SQRT(Z) can be decided by:
-// yn = 1 if 2*R(n - 1) >= 2*Y(n - 1) + 2^(-n),
-// 0 otherwise.
-//
-// 3. Precision analysis:
-//
-// - Notice that in the decision function:
-// 2*R(n - 1) >= 2*Y(n - 1) + 2^(-n),
-// the right hand side only uses up to the 2^(-n)-bit, and both sides are
-// non-negative, so R(n - 1) can be truncated at the 2^(-(n + 1))-bit, so
-// that 2*R(n - 1) is corrected up to the 2^(-n)-bit.
-//
-// - Thus, in order to round SQRT(a^2 + b^2) correctly up to n-fractional
-// bits, we need to perform the summation (a^2 + b^2) correctly up to (2n +
-// 2)-fractional bits, and the remaining bits are sticky bits (i.e. we only
-// care if they are 0 or > 0), and the comparisons, additions/subtractions
-// can be done in n-fractional bits precision.
-//
-// - For single precision (float), we can use uint64_t to store the sum a^2 +
-// b^2 exact up to (2n + 2)-fractional bits.
-//
-// - Then we can feed this sum into the digit-by-digit algorithm for SQRT(Z)
-// described above.
-//
-//
-// Special cases:
-// - HYPOT(x, y) is +Inf if x or y is +Inf or -Inf; else
-// - HYPOT(x, y) is NaN if x or y is NaN.
-//
float LLVM_LIBC_ENTRYPOINT(hypotf)(float x, float y) {
- FPBits<float> x_bits(x), y_bits(y);
-
- if (x_bits.isInf() || y_bits.isInf()) {
- return FPBits<float>::inf();
- }
- if (x_bits.isNaN()) {
- return x;
- }
- if (y_bits.isNaN()) {
- return y;
- }
-
- uint16_t a_exp, b_exp, out_exp;
- uint32_t a_mant, b_mant;
- uint64_t a_mant_sq, b_mant_sq;
- bool sticky_bits;
-
- if ((x_bits.exponent >= y_bits.exponent + MantissaWidth<float>::value + 2) ||
- (y == 0)) {
- return abs(x);
- } else if ((y_bits.exponent >=
- x_bits.exponent + MantissaWidth<float>::value + 2) ||
- (x == 0)) {
- y_bits.sign = 0;
- return abs(y);
- }
-
- if (x >= y) {
- a_exp = x_bits.exponent;
- a_mant = x_bits.mantissa;
- b_exp = y_bits.exponent;
- b_mant = y_bits.mantissa;
- } else {
- a_exp = y_bits.exponent;
- a_mant = y_bits.mantissa;
- b_exp = x_bits.exponent;
- b_mant = x_bits.mantissa;
- }
-
- out_exp = a_exp;
-
- // Add an extra bit to simplify the final rounding bit computation.
- constexpr uint32_t one = 1U << (MantissaWidth<float>::value + 1);
-
- a_mant <<= 1;
- b_mant <<= 1;
-
- uint32_t leading_one;
- int y_mant_width;
- if (a_exp != 0) {
- leading_one = one;
- a_mant |= one;
- y_mant_width = MantissaWidth<float>::value + 1;
- } else {
- leading_one = findLeadingOne(a_mant, y_mant_width);
- }
-
- if (b_exp != 0) {
- b_mant |= one;
- }
-
- a_mant_sq = static_cast<uint64_t>(a_mant) * a_mant;
- b_mant_sq = static_cast<uint64_t>(b_mant) * b_mant;
-
- // At this point, a_exp >= b_exp > a_exp - 25, so in order to line up aSqMant
- // and bSqMant, we need to shift bSqMant to the right by (a_exp - b_exp) bits.
- // But before that, remember to store the losing bits to sticky.
- // The shift length is for a^2 and b^2, so it's double of the exponent
- // difference between a and b.
- uint16_t shift_length = 2 * (a_exp - b_exp);
- sticky_bits = ((b_mant_sq & ((1ULL << shift_length) - 1)) != 0);
- b_mant_sq >>= shift_length;
-
- uint64_t sum = a_mant_sq + b_mant_sq;
- if (sum >= (1ULL << (2 * y_mant_width + 2))) {
- // a^2 + b^2 >= 4* leading_one^2, so we will need an extra bit to the left.
- if (leading_one == one) {
- // For normal result, we discard the last 2 bits of the sum and increase
- // the exponent.
- sticky_bits = sticky_bits || ((sum & 0x3U) != 0);
- sum >>= 2;
- ++out_exp;
- if (out_exp >= FPBits<float>::maxExponent) {
- return FPBits<float>::inf();
- }
- } else {
- // For denormal result, we simply move the leading bit of the result to
- // the left by 1.
- leading_one <<= 1;
- ++y_mant_width;
- }
- }
-
- uint32_t Y = leading_one;
- uint32_t R = static_cast<uint32_t>(sum >> y_mant_width) - leading_one;
- uint32_t tailBits = static_cast<uint32_t>(sum) & (leading_one - 1);
-
- for (uint32_t current_bit = leading_one >> 1; current_bit;
- current_bit >>= 1) {
- R = (R << 1) + ((tailBits & current_bit) ? 1 : 0);
- uint32_t tmp = (Y << 1) + current_bit; // 2*y(n - 1) + 2^(-n)
- if (R >= tmp) {
- R -= tmp;
- Y += current_bit;
- }
- }
-
- bool round_bit = Y & 1U;
- bool lsb = Y & 2U;
-
- if (Y >= one) {
- Y -= one;
-
- if (out_exp == 0) {
- out_exp = 1;
- }
- }
-
- Y >>= 1;
-
- // Round to the nearest, tie to even.
- if (round_bit && (lsb || sticky_bits || (R != 0))) {
- ++Y;
- }
-
- if (Y >= (one >> 1)) {
- Y -= one >> 1;
- ++out_exp;
- if (out_exp >= FPBits<float>::maxExponent) {
- return FPBits<float>::inf();
- }
- }
-
- Y |= static_cast<uint32_t>(out_exp) << MantissaWidth<float>::value;
- return *reinterpret_cast<float *>(&Y);
+ return __llvm_libc::fputil::hypot(x, y);
}
} // namespace __llvm_libc
diff --git a/test/src/math/CMakeLists.txt b/test/src/math/CMakeLists.txt
index cdffe73..8635e7a 100644
--- a/test/src/math/CMakeLists.txt
+++ b/test/src/math/CMakeLists.txt
@@ -736,3 +736,16 @@
libc.src.math.hypotf
libc.utils.FPUtil.fputil
)
+
+add_fp_unittest(
+ hypot_test
+ NEED_MPFR
+ SUITE
+ libc_math_unittests
+ SRCS
+ hypot_test.cpp
+ DEPENDS
+ libc.include.math
+ libc.src.math.hypot
+ libc.utils.FPUtil.fputil
+)
diff --git a/test/src/math/HypotTest.h b/test/src/math/HypotTest.h
new file mode 100644
index 0000000..f90807b
--- /dev/null
+++ b/test/src/math/HypotTest.h
@@ -0,0 +1,75 @@
+//===-- Utility class to test different flavors of hypot ------------------===//
+//
+// 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_TEST_SRC_MATH_HYPOTTEST_H
+#define LLVM_LIBC_TEST_SRC_MATH_HYPOTTEST_H
+
+#include "include/math.h"
+#include "utils/FPUtil/FPBits.h"
+#include "utils/FPUtil/Hypot.h"
+#include "utils/FPUtil/TestHelpers.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include "utils/UnitTest/Test.h"
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+template <typename T>
+class HypotTestTemplate : public __llvm_libc::testing::Test {
+private:
+ using Func = T (*)(T, T);
+ using FPBits = __llvm_libc::fputil::FPBits<T>;
+ using UIntType = typename FPBits::UIntType;
+ const T nan = __llvm_libc::fputil::FPBits<T>::buildNaN(1);
+ const T inf = __llvm_libc::fputil::FPBits<T>::inf();
+ const T negInf = __llvm_libc::fputil::FPBits<T>::negInf();
+ const T zero = __llvm_libc::fputil::FPBits<T>::zero();
+ const T negZero = __llvm_libc::fputil::FPBits<T>::negZero();
+
+public:
+ void testSpecialNumbers(Func func) {
+ EXPECT_FP_EQ(func(inf, nan), inf);
+ EXPECT_FP_EQ(func(nan, negInf), inf);
+ EXPECT_FP_EQ(func(zero, inf), inf);
+ EXPECT_FP_EQ(func(negInf, negZero), inf);
+
+ EXPECT_FP_EQ(func(nan, nan), nan);
+ EXPECT_FP_EQ(func(nan, zero), nan);
+ EXPECT_FP_EQ(func(negZero, nan), nan);
+
+ EXPECT_FP_EQ(func(negZero, zero), zero);
+ }
+
+ void testSubnormalRange(Func func) {
+ constexpr UIntType count = 1000001;
+ constexpr UIntType step =
+ (FPBits::maxSubnormal - FPBits::minSubnormal) / count;
+ for (UIntType v = FPBits::minSubnormal, w = FPBits::maxSubnormal;
+ v <= FPBits::maxSubnormal && w >= FPBits::minSubnormal;
+ v += step, w -= step) {
+ T x = FPBits(v), y = FPBits(w);
+ T result = func(x, y);
+ mpfr::BinaryInput<T> input{x, y};
+ ASSERT_MPFR_MATCH(mpfr::Operation::Hypot, input, result, 0.5);
+ }
+ }
+
+ void testNormalRange(Func func) {
+ constexpr UIntType count = 1000001;
+ constexpr UIntType step = (FPBits::maxNormal - FPBits::minNormal) / count;
+ for (UIntType v = FPBits::minNormal, w = FPBits::maxNormal;
+ v <= FPBits::maxNormal && w >= FPBits::minNormal;
+ v += step, w -= step) {
+ T x = FPBits(v), y = FPBits(w);
+ T result = func(x, y);
+ mpfr::BinaryInput<T> input{x, y};
+ ASSERT_MPFR_MATCH(mpfr::Operation::Hypot, input, result, 0.5);
+ }
+ }
+};
+
+#endif // LLVM_LIBC_TEST_SRC_MATH_HYPOTTEST_H
diff --git a/test/src/math/hypot_test.cpp b/test/src/math/hypot_test.cpp
new file mode 100644
index 0000000..d723f52
--- /dev/null
+++ b/test/src/math/hypot_test.cpp
@@ -0,0 +1,20 @@
+//===-- Unittests for hypot -----------------------------------------------===//
+//
+// 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 "HypotTest.h"
+
+#include "include/math.h"
+#include "src/math/hypot.h"
+
+using HypotTest = HypotTestTemplate<double>;
+
+TEST_F(HypotTest, SpecialNumbers) { testSpecialNumbers(&__llvm_libc::hypot); }
+
+TEST_F(HypotTest, SubnormalRange) { testSubnormalRange(&__llvm_libc::hypot); }
+
+TEST_F(HypotTest, NormalRange) { testNormalRange(&__llvm_libc::hypot); }
diff --git a/test/src/math/hypotf_test.cpp b/test/src/math/hypotf_test.cpp
index 1769307..21d1bea 100644
--- a/test/src/math/hypotf_test.cpp
+++ b/test/src/math/hypotf_test.cpp
@@ -6,56 +6,15 @@
//
//===----------------------------------------------------------------------===//
+#include "HypotTest.h"
+
+#include "include/math.h"
#include "src/math/hypotf.h"
-#include "utils/FPUtil/FPBits.h"
-#include "utils/FPUtil/TestHelpers.h"
-#include "utils/MPFRWrapper/MPFRUtils.h"
-#include "utils/UnitTest/Test.h"
-#include <math.h>
-using FPBits = __llvm_libc::fputil::FPBits<float>;
-using UIntType = FPBits::UIntType;
+using HypotfTest = HypotTestTemplate<float>;
-namespace mpfr = __llvm_libc::testing::mpfr;
+TEST_F(HypotfTest, SpecialNumbers) { testSpecialNumbers(&__llvm_libc::hypotf); }
-DECLARE_SPECIAL_CONSTANTS(float)
+TEST_F(HypotfTest, SubnormalRange) { testSubnormalRange(&__llvm_libc::hypotf); }
-TEST(HypotfTest, SpecialNumbers) {
- EXPECT_FP_EQ(__llvm_libc::hypotf(inf, nan), inf);
- EXPECT_FP_EQ(__llvm_libc::hypotf(nan, negInf), inf);
- EXPECT_FP_EQ(__llvm_libc::hypotf(zero, inf), inf);
- EXPECT_FP_EQ(__llvm_libc::hypotf(negInf, negZero), inf);
-
- EXPECT_FP_EQ(__llvm_libc::hypotf(nan, nan), nan);
- EXPECT_FP_EQ(__llvm_libc::hypotf(nan, zero), nan);
- EXPECT_FP_EQ(__llvm_libc::hypotf(negZero, nan), nan);
-
- EXPECT_FP_EQ(__llvm_libc::hypotf(negZero, zero), zero);
-}
-
-TEST(HypotfTest, SubnormalRange) {
- constexpr UIntType count = 1000001;
- constexpr UIntType step =
- (FPBits::maxSubnormal - FPBits::minSubnormal) / count;
- for (UIntType v = FPBits::minSubnormal, w = FPBits::maxSubnormal;
- v <= FPBits::maxSubnormal && w >= FPBits::minSubnormal;
- v += step, w -= step) {
- float x = FPBits(v), y = FPBits(w);
- float result = __llvm_libc::hypotf(x, y);
- mpfr::BinaryInput<float> input{x, y};
- ASSERT_MPFR_MATCH(mpfr::Operation::Hypot, input, result, 0.5);
- }
-}
-
-TEST(HypotfTest, NormalRange) {
- constexpr UIntType count = 1000001;
- constexpr UIntType step = (FPBits::maxNormal - FPBits::minNormal) / count;
- for (UIntType v = FPBits::minNormal, w = FPBits::maxNormal;
- v <= FPBits::maxNormal && w >= FPBits::minNormal; v += step, w -= step) {
- float x = FPBits(v), y = FPBits(w);
- float result = __llvm_libc::hypotf(x, y);
- ;
- mpfr::BinaryInput<float> input{x, y};
- ASSERT_MPFR_MATCH(mpfr::Operation::Hypot, input, result, 0.5);
- }
-}
+TEST_F(HypotfTest, NormalRange) { testNormalRange(&__llvm_libc::hypotf); }
diff --git a/utils/FPUtil/Hypot.h b/utils/FPUtil/Hypot.h
new file mode 100644
index 0000000..6795f9d
--- /dev/null
+++ b/utils/FPUtil/Hypot.h
@@ -0,0 +1,267 @@
+//===-- Implementation of hypotf 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
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_UTILS_FPUTIL_HYPOT_H
+#define LLVM_LIBC_UTILS_FPUTIL_HYPOT_H
+
+#include "BasicOperations.h"
+#include "FPBits.h"
+#include "utils/CPP/TypeTraits.h"
+
+namespace __llvm_libc {
+namespace fputil {
+
+namespace internal {
+
+template <typename T> static inline T findLeadingOne(T mant, int &shift_length);
+
+template <>
+inline uint32_t findLeadingOne<uint32_t>(uint32_t mant, int &shift_length) {
+ shift_length = 0;
+ constexpr int nsteps = 5;
+ constexpr uint32_t bounds[nsteps] = {1 << 16, 1 << 8, 1 << 4, 1 << 2, 1 << 1};
+ constexpr int shifts[nsteps] = {16, 8, 4, 2, 1};
+ for (int i = 0; i < nsteps; ++i) {
+ if (mant >= bounds[i]) {
+ shift_length += shifts[i];
+ mant >>= shifts[i];
+ }
+ }
+ return 1U << shift_length;
+}
+
+template <>
+inline uint64_t findLeadingOne<uint64_t>(uint64_t mant, int &shift_length) {
+ shift_length = 0;
+ constexpr int nsteps = 6;
+ constexpr uint64_t bounds[nsteps] = {1ULL << 32, 1ULL << 16, 1ULL << 8,
+ 1ULL << 4, 1ULL << 2, 1ULL << 1};
+ constexpr int shifts[nsteps] = {32, 16, 8, 4, 2, 1};
+ for (int i = 0; i < nsteps; ++i) {
+ if (mant >= bounds[i]) {
+ shift_length += shifts[i];
+ mant >>= shifts[i];
+ }
+ }
+ return 1ULL << shift_length;
+}
+
+} // namespace internal
+
+template <typename T> struct DoubleLength;
+
+template <> struct DoubleLength<uint16_t> { using Type = uint32_t; };
+
+template <> struct DoubleLength<uint32_t> { using Type = uint64_t; };
+
+template <> struct DoubleLength<uint64_t> { using Type = __uint128_t; };
+
+// Correctly rounded IEEE 754 HYPOT(x, y) with round to nearest, ties to even.
+//
+// Algorithm:
+// - Let a = max(|x|, |y|), b = min(|x|, |y|), then we have that:
+// a <= sqrt(a^2 + b^2) <= min(a + b, a*sqrt(2))
+// 1. So if b < eps(a)/2, then HYPOT(x, y) = a.
+//
+// - Moreover, the exponent part of HYPOT(x, y) is either the same or 1 more
+// than the exponent part of a.
+//
+// 2. For the remaining cases, we will use the digit-by-digit (shift-and-add)
+// algorithm to compute SQRT(Z):
+//
+// - For Y = y0.y1...yn... = SQRT(Z),
+// let Y(n) = y0.y1...yn be the first n fractional digits of Y.
+//
+// - The nth scaled residual R(n) is defined to be:
+// R(n) = 2^n * (Z - Y(n)^2)
+//
+// - Since Y(n) = Y(n - 1) + yn * 2^(-n), the scaled residual
+// satisfies the following recurrence formula:
+// R(n) = 2*R(n - 1) - yn*(2*Y(n - 1) + 2^(-n)),
+// with the initial conditions:
+// Y(0) = y0, and R(0) = Z - y0.
+//
+// - So the nth fractional digit of Y = SQRT(Z) can be decided by:
+// yn = 1 if 2*R(n - 1) >= 2*Y(n - 1) + 2^(-n),
+// 0 otherwise.
+//
+// 3. Precision analysis:
+//
+// - Notice that in the decision function:
+// 2*R(n - 1) >= 2*Y(n - 1) + 2^(-n),
+// the right hand side only uses up to the 2^(-n)-bit, and both sides are
+// non-negative, so R(n - 1) can be truncated at the 2^(-(n + 1))-bit, so
+// that 2*R(n - 1) is corrected up to the 2^(-n)-bit.
+//
+// - Thus, in order to round SQRT(a^2 + b^2) correctly up to n-fractional
+// bits, we need to perform the summation (a^2 + b^2) correctly up to (2n +
+// 2)-fractional bits, and the remaining bits are sticky bits (i.e. we only
+// care if they are 0 or > 0), and the comparisons, additions/subtractions
+// can be done in n-fractional bits precision.
+//
+// - For single precision (float), we can use uint64_t to store the sum a^2 +
+// b^2 exact up to (2n + 2)-fractional bits.
+//
+// - Then we can feed this sum into the digit-by-digit algorithm for SQRT(Z)
+// described above.
+//
+//
+// Special cases:
+// - HYPOT(x, y) is +Inf if x or y is +Inf or -Inf; else
+// - HYPOT(x, y) is NaN if x or y is NaN.
+//
+template <typename T,
+ cpp::EnableIfType<cpp::IsFloatingPointType<T>::Value, int> = 0>
+static inline T hypot(T x, T y) {
+ using FPBits_t = FPBits<T>;
+ using UIntType = typename FPBits<T>::UIntType;
+ using DUIntType = typename DoubleLength<UIntType>::Type;
+
+ FPBits_t x_bits(x), y_bits(y);
+
+ if (x_bits.isInf() || y_bits.isInf()) {
+ return FPBits_t::inf();
+ }
+ if (x_bits.isNaN()) {
+ return x;
+ }
+ if (y_bits.isNaN()) {
+ return y;
+ }
+
+ uint16_t a_exp, b_exp, out_exp;
+ UIntType a_mant, b_mant;
+ DUIntType a_mant_sq, b_mant_sq;
+ bool sticky_bits;
+
+ if ((x_bits.exponent >= y_bits.exponent + MantissaWidth<T>::value + 2) ||
+ (y == 0)) {
+ return abs(x);
+ } else if ((y_bits.exponent >=
+ x_bits.exponent + MantissaWidth<T>::value + 2) ||
+ (x == 0)) {
+ y_bits.sign = 0;
+ return abs(y);
+ }
+
+ if (x >= y) {
+ a_exp = x_bits.exponent;
+ a_mant = x_bits.mantissa;
+ b_exp = y_bits.exponent;
+ b_mant = y_bits.mantissa;
+ } else {
+ a_exp = y_bits.exponent;
+ a_mant = y_bits.mantissa;
+ b_exp = x_bits.exponent;
+ b_mant = x_bits.mantissa;
+ }
+
+ out_exp = a_exp;
+
+ // Add an extra bit to simplify the final rounding bit computation.
+ constexpr UIntType one = UIntType(1) << (MantissaWidth<T>::value + 1);
+
+ a_mant <<= 1;
+ b_mant <<= 1;
+
+ UIntType leading_one;
+ int y_mant_width;
+ if (a_exp != 0) {
+ leading_one = one;
+ a_mant |= one;
+ y_mant_width = MantissaWidth<T>::value + 1;
+ } else {
+ leading_one = internal::findLeadingOne(a_mant, y_mant_width);
+ }
+
+ if (b_exp != 0) {
+ b_mant |= one;
+ }
+
+ a_mant_sq = static_cast<DUIntType>(a_mant) * a_mant;
+ b_mant_sq = static_cast<DUIntType>(b_mant) * b_mant;
+
+ // At this point, a_exp >= b_exp > a_exp - 25, so in order to line up aSqMant
+ // and bSqMant, we need to shift bSqMant to the right by (a_exp - b_exp) bits.
+ // But before that, remember to store the losing bits to sticky.
+ // The shift length is for a^2 and b^2, so it's double of the exponent
+ // difference between a and b.
+ uint16_t shift_length = 2 * (a_exp - b_exp);
+ sticky_bits =
+ ((b_mant_sq & ((DUIntType(1) << shift_length) - DUIntType(1))) !=
+ DUIntType(0));
+ b_mant_sq >>= shift_length;
+
+ DUIntType sum = a_mant_sq + b_mant_sq;
+ if (sum >= (DUIntType(1) << (2 * y_mant_width + 2))) {
+ // a^2 + b^2 >= 4* leading_one^2, so we will need an extra bit to the left.
+ if (leading_one == one) {
+ // For normal result, we discard the last 2 bits of the sum and increase
+ // the exponent.
+ sticky_bits = sticky_bits || ((sum & 0x3U) != 0);
+ sum >>= 2;
+ ++out_exp;
+ if (out_exp >= FPBits_t::maxExponent) {
+ return FPBits_t::inf();
+ }
+ } else {
+ // For denormal result, we simply move the leading bit of the result to
+ // the left by 1.
+ leading_one <<= 1;
+ ++y_mant_width;
+ }
+ }
+
+ UIntType Y = leading_one;
+ UIntType R = static_cast<UIntType>(sum >> y_mant_width) - leading_one;
+ UIntType tailBits = static_cast<UIntType>(sum) & (leading_one - 1);
+
+ for (UIntType current_bit = leading_one >> 1; current_bit;
+ current_bit >>= 1) {
+ R = (R << 1) + ((tailBits & current_bit) ? 1 : 0);
+ UIntType tmp = (Y << 1) + current_bit; // 2*y(n - 1) + 2^(-n)
+ if (R >= tmp) {
+ R -= tmp;
+ Y += current_bit;
+ }
+ }
+
+ bool round_bit = Y & UIntType(1);
+ bool lsb = Y & UIntType(2);
+
+ if (Y >= one) {
+ Y -= one;
+
+ if (out_exp == 0) {
+ out_exp = 1;
+ }
+ }
+
+ Y >>= 1;
+
+ // Round to the nearest, tie to even.
+ if (round_bit && (lsb || sticky_bits || (R != 0))) {
+ ++Y;
+ }
+
+ if (Y >= (one >> 1)) {
+ Y -= one >> 1;
+ ++out_exp;
+ if (out_exp >= FPBits_t::maxExponent) {
+ return FPBits_t::inf();
+ }
+ }
+
+ Y |= static_cast<UIntType>(out_exp) << MantissaWidth<T>::value;
+ return *reinterpret_cast<T *>(&Y);
+}
+
+} // namespace fputil
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_UTILS_FPUTIL_HYPOT_H
