diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt --- a/libc/config/linux/aarch64/entrypoints.txt +++ b/libc/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/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt --- a/libc/config/linux/x86_64/entrypoints.txt +++ b/libc/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/libc/spec/stdc.td b/libc/spec/stdc.td --- a/libc/spec/stdc.td +++ b/libc/spec/stdc.td @@ -280,6 +280,7 @@ FunctionSpec<"frexpf", RetValSpec, [ArgSpec, ArgSpec]>, FunctionSpec<"frexpl", RetValSpec, [ArgSpec, ArgSpec]>, + FunctionSpec<"hypot", RetValSpec, [ArgSpec, ArgSpec]>, FunctionSpec<"hypotf", RetValSpec, [ArgSpec, ArgSpec]>, FunctionSpec<"ilogb", RetValSpec, [ArgSpec]>, diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt --- a/libc/src/math/CMakeLists.txt +++ b/libc/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/libc/src/math/hypot.h b/libc/src/math/hypot.h new file mode 100644 --- /dev/null +++ b/libc/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/libc/src/math/hypot.cpp b/libc/src/math/hypot.cpp new file mode 100644 --- /dev/null +++ b/libc/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/libc/src/math/hypotf.cpp b/libc/src/math/hypotf.cpp --- a/libc/src/math/hypotf.cpp +++ b/libc/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 x_bits(x), y_bits(y); - - if (x_bits.isInf() || y_bits.isInf()) { - return FPBits::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::value + 2) || - (y == 0)) { - return abs(x); - } else if ((y_bits.exponent >= - x_bits.exponent + MantissaWidth::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::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::value + 1; - } else { - leading_one = findLeadingOne(a_mant, y_mant_width); - } - - if (b_exp != 0) { - b_mant |= one; - } - - a_mant_sq = static_cast(a_mant) * a_mant; - b_mant_sq = static_cast(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::maxExponent) { - return FPBits::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(sum >> y_mant_width) - leading_one; - uint32_t tailBits = static_cast(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::maxExponent) { - return FPBits::inf(); - } - } - - Y |= static_cast(out_exp) << MantissaWidth::value; - return *reinterpret_cast(&Y); + return __llvm_libc::fputil::hypot(x, y); } } // namespace __llvm_libc diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt --- a/libc/test/src/math/CMakeLists.txt +++ b/libc/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/libc/test/src/math/HypotTest.h b/libc/test/src/math/HypotTest.h new file mode 100644 --- /dev/null +++ b/libc/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 +class HypotTestTemplate : public __llvm_libc::testing::Test { +private: + using Func = T (*)(T, T); + using FPBits = __llvm_libc::fputil::FPBits; + using UIntType = typename FPBits::UIntType; + const T nan = __llvm_libc::fputil::FPBits::buildNaN(1); + const T inf = __llvm_libc::fputil::FPBits::inf(); + const T negInf = __llvm_libc::fputil::FPBits::negInf(); + const T zero = __llvm_libc::fputil::FPBits::zero(); + const T negZero = __llvm_libc::fputil::FPBits::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 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 input{x, y}; + ASSERT_MPFR_MATCH(mpfr::Operation::Hypot, input, result, 0.5); + } + } +}; + +#endif // LLVM_LIBC_TEST_SRC_MATH_HYPOTTEST_H diff --git a/libc/test/src/math/hypot_test.cpp b/libc/test/src/math/hypot_test.cpp new file mode 100644 --- /dev/null +++ b/libc/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; + +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/libc/test/src/math/hypotf_test.cpp b/libc/test/src/math/hypotf_test.cpp --- a/libc/test/src/math/hypotf_test.cpp +++ b/libc/test/src/math/hypotf_test.cpp @@ -6,56 +6,15 @@ // //===----------------------------------------------------------------------===// -#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 - -using FPBits = __llvm_libc::fputil::FPBits; -using UIntType = FPBits::UIntType; - -namespace mpfr = __llvm_libc::testing::mpfr; +#include "HypotTest.h" -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); +#include "include/math.h" +#include "src/math/hypotf.h" - 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); +using HypotfTest = HypotTestTemplate; - EXPECT_FP_EQ(__llvm_libc::hypotf(negZero, zero), zero); -} +TEST_F(HypotfTest, SpecialNumbers) { testSpecialNumbers(&__llvm_libc::hypotf); } -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 input{x, y}; - ASSERT_MPFR_MATCH(mpfr::Operation::Hypot, input, result, 0.5); - } -} +TEST_F(HypotfTest, SubnormalRange) { testSubnormalRange(&__llvm_libc::hypotf); } -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 input{x, y}; - ASSERT_MPFR_MATCH(mpfr::Operation::Hypot, input, result, 0.5); - } -} +TEST_F(HypotfTest, NormalRange) { testNormalRange(&__llvm_libc::hypotf); } diff --git a/libc/src/math/hypotf.cpp b/libc/utils/FPUtil/Hypot.h copy from libc/src/math/hypotf.cpp copy to libc/utils/FPUtil/Hypot.h --- a/libc/src/math/hypotf.cpp +++ b/libc/utils/FPUtil/Hypot.h @@ -5,15 +5,23 @@ // 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 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 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}; @@ -27,6 +35,32 @@ return 1U << shift_length; } +template <> +inline uint64_t findLeadingOne(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 struct DoubleLength; + +template <> struct DoubleLength { using Type = uint32_t; }; + +template <> struct DoubleLength { using Type = uint64_t; }; + +template <> struct DoubleLength { using Type = __uint128_t; }; + // Correctly rounded IEEE 754 HYPOT(x, y) with round to nearest, ties to even. // // Algorithm: @@ -81,11 +115,17 @@ // - 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 x_bits(x), y_bits(y); +template ::Value, int> = 0> +static inline T hypot(T x, T y) { + using FPBits_t = FPBits; + using UIntType = typename FPBits::UIntType; + using DUIntType = typename DoubleLength::Type; + + FPBits_t x_bits(x), y_bits(y); if (x_bits.isInf() || y_bits.isInf()) { - return FPBits::inf(); + return FPBits_t::inf(); } if (x_bits.isNaN()) { return x; @@ -95,15 +135,15 @@ } 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::value + 2) || + if ((x_bits.exponent >= y_bits.exponent + MantissaWidth::value + 2) || (y == 0)) { return abs(x); } else if ((y_bits.exponent >= - x_bits.exponent + MantissaWidth::value + 2) || + x_bits.exponent + MantissaWidth::value + 2) || (x == 0)) { y_bits.sign = 0; return abs(y); @@ -124,27 +164,27 @@ out_exp = a_exp; // Add an extra bit to simplify the final rounding bit computation. - constexpr uint32_t one = 1U << (MantissaWidth::value + 1); + constexpr UIntType one = UIntType(1) << (MantissaWidth::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::value + 1; + y_mant_width = MantissaWidth::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(a_mant) * a_mant; - b_mant_sq = static_cast(b_mant) * b_mant; + a_mant_sq = static_cast(a_mant) * a_mant; + b_mant_sq = static_cast(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. @@ -152,11 +192,13 @@ // 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 @@ -164,8 +206,8 @@ sticky_bits = sticky_bits || ((sum & 0x3U) != 0); sum >>= 2; ++out_exp; - if (out_exp >= FPBits::maxExponent) { - return FPBits::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 @@ -175,22 +217,22 @@ } } - uint32_t Y = leading_one; - uint32_t R = static_cast(sum >> y_mant_width) - leading_one; - uint32_t tailBits = static_cast(sum) & (leading_one - 1); + UIntType Y = leading_one; + UIntType R = static_cast(sum >> y_mant_width) - leading_one; + UIntType tailBits = static_cast(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; @@ -210,13 +252,16 @@ if (Y >= (one >> 1)) { Y -= one >> 1; ++out_exp; - if (out_exp >= FPBits::maxExponent) { - return FPBits::inf(); + if (out_exp >= FPBits_t::maxExponent) { + return FPBits_t::inf(); } } - Y |= static_cast(out_exp) << MantissaWidth::value; - return *reinterpret_cast(&Y); + Y |= static_cast(out_exp) << MantissaWidth::value; + return *reinterpret_cast(&Y); } +} // namespace fputil } // namespace __llvm_libc + +#endif // LLVM_LIBC_UTILS_FPUTIL_HYPOT_H