[libc] Add implementation of hypot.

[lntue]

Summary

Refactor src/math/hypotf.cpp and test/src/math/hypotf_test.cpp and reuse them for hypot and hypot_test

[sivachandra]

libc/utils/FPUtil/Hypot.h

121

In other functions, we have not used the _t suffixes. To be consistent, can you follow the same convention here? If you feel strongly about the _t style, you can do a follow up cleanup pass changing everything to follow the _t style.

122

Same here: other places do not have the _t but use UIntType with the template type.

123

Same.

[sivachandra]

Just realized I missed selecting accept.

[lntue]

libc/utils/FPUtil/Hypot.h

121

Changed to UIntType

[lntue]

Change UInt_t to UIntType.

libc/config/linux/aarch64/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.fmaxf
     libc.src.math.fmaxl
     libc.src.math.fmin
     libc.src.math.fminf
     libc.src.math.fminl
     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
     libc.src.math.ilogbl
     libc.src.math.ldexp
     libc.src.math.ldexpf
     libc.src.math.ldexpl
     libc.src.math.logb

libc/config/linux/x86_64/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.fminf
     libc.src.math.fminl
     libc.src.math.fmax
     libc.src.math.fmaxf
     libc.src.math.fmaxl
     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
     libc.src.math.ilogbl
     libc.src.math.ldexp
     libc.src.math.ldexpf
     libc.src.math.ldexpl
     libc.src.math.logb

libc/spec/stdc.td

HeaderSpec Math = HeaderSpec<
           FunctionSpec<"fmax", RetValSpec<DoubleType>, [ArgSpec<DoubleType>, ArgSpec<DoubleType>]>,
           FunctionSpec<"fmaxf", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<FloatType>]>,
           FunctionSpec<"fmaxl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>, ArgSpec<LongDoubleType>]>,
 
           FunctionSpec<"frexp", RetValSpec<DoubleType>, [ArgSpec<DoubleType>, ArgSpec<IntPtr>]>,
           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>]>,
           FunctionSpec<"ilogbf", RetValSpec<IntType>, [ArgSpec<FloatType>]>,
           FunctionSpec<"ilogbl", RetValSpec<IntType>, [ArgSpec<LongDoubleType>]>,
 
           FunctionSpec<"ldexp", RetValSpec<DoubleType>, [ArgSpec<DoubleType>, ArgSpec<IntType>]>,
           FunctionSpec<"ldexpf", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<IntType>]>,

libc/src/math/CMakeLists.txt

SRCS
     fdiml.cpp
   HDRS
     fdiml.h
   DEPENDS
     libc.utils.FPUtil.fputil
   COMPILE_OPTIONS
     -O2
 )
+
+add_entrypoint_object(
+  hypot
+  SRCS
+    hypot.cpp
+  HDRS
+    hypot.h
+  DEPENDS
+    libc.utils.FPUtil.fputil
+  COMPILE_OPTIONS
+    -O2
+)

libc/src/math/hypot.h

+//===-- 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

libc/src/math/hypot.cpp

+//===-- 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

libc/src/math/hypotf.cpp

 //===-- 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
 //
 //===----------------------------------------------------------------------===//
 #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

libc/test/src/math/CMakeLists.txt

SUITE
     libc_math_unittests
   SRCS
     hypotf_test.cpp
   DEPENDS
     libc.include.math
     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
+)

libc/test/src/math/HypotTest.h

+//===-- 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

libc/test/src/math/hypot_test.cpp

+//===-- 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); }

libc/test/src/math/hypotf_test.cpp

 //===-- Unittests for hypotf ----------------------------------------------===//
 //
 // 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/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;
-
-namespace mpfr = __llvm_libc::testing::mpfr;
-
-DECLARE_SPECIAL_CONSTANTS(float)
-
-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);
-  }
-}
+
+using HypotfTest = HypotTestTemplate<float>;
+
+TEST_F(HypotfTest, SpecialNumbers) { testSpecialNumbers(&__llvm_libc::hypotf); }
+
+TEST_F(HypotfTest, SubnormalRange) { testSubnormalRange(&__llvm_libc::hypotf); }
+
+TEST_F(HypotfTest, NormalRange) { testNormalRange(&__llvm_libc::hypotf); }

libc/utils/FPUtil/Hypot.h

 //===-- 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
 //
 //===----------------------------------------------------------------------===//
-#include "src/__support/common.h"
-#include "utils/FPUtil/BasicOperations.h"
-#include "utils/FPUtil/FPBits.h"
+
+#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 {
 
-using namespace fputil;
+template <typename T> static inline T findLeadingOne(T mant, int &shift_length);
 
-uint32_t findLeadingOne(uint32_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

 //   -  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);
+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<float>::inf();
+    return FPBits_t::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;
+  UIntType a_mant, b_mant;
+  DUIntType a_mant_sq, b_mant_sq;
   bool sticky_bits;
 
-  if ((x_bits.exponent >= y_bits.exponent + MantissaWidth<float>::value + 2) ||
+  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<float>::value + 2) ||
+              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 uint32_t one = 1U << (MantissaWidth<float>::value + 1);
+  constexpr UIntType one = UIntType(1) << (MantissaWidth<T>::value + 1);
 
   a_mant <<= 1;
   b_mant <<= 1;
 
-  uint32_t leading_one;
+  UIntType leading_one;
   int y_mant_width;
   if (a_exp != 0) {
     leading_one = one;
     a_mant |= one;
-    y_mant_width = MantissaWidth<float>::value + 1;
+    y_mant_width = MantissaWidth<T>::value + 1;
   } else {
-    leading_one = findLeadingOne(a_mant, y_mant_width);
+    leading_one = internal::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;
+  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 & ((1ULL << shift_length) - 1)) != 0);
+  sticky_bits =
+      ((b_mant_sq & ((DUIntType(1) << shift_length) - DUIntType(1))) !=
+       DUIntType(0));
   b_mant_sq >>= shift_length;
 
-  uint64_t sum = a_mant_sq + b_mant_sq;
-  if (sum >= (1ULL << (2 * y_mant_width + 2))) {
+  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<float>::maxExponent) {
-        return FPBits<float>::inf();
+      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;
     }
   }
 
-  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);
+  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 (uint32_t current_bit = leading_one >> 1; current_bit;
+  for (UIntType 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)
+    UIntType 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;
+  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<float>::maxExponent) {
-      return FPBits<float>::inf();
+    if (out_exp >= FPBits_t::maxExponent) {
+      return FPBits_t::inf();
     }
   }
 
-  Y |= static_cast<uint32_t>(out_exp) << MantissaWidth<float>::value;
-  return *reinterpret_cast<float *>(&Y);
+  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