maxjhandsome
/
llvm-project
forked from OSchip/llvm-project
1
0
Fork 0
Code Issues Pull Requests Packages Releases Wiki Activity

[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
This commit is contained in:
Tue Ly 2020-09-17 23:22:18 -04:00
parent c00516d520
commit 3b487d51e2
12 changed files with 435 additions and 255 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
libc/config/linux/aarch64/entrypoints.txt
View File

@ -68,6 +68,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.frexp libc.src.math.frexp
libc.src.math.frexpf libc.src.math.frexpf
libc.src.math.frexpl libc.src.math.frexpl
libc.src.math.hypot
libc.src.math.hypotf libc.src.math.hypotf
libc.src.math.ilogb libc.src.math.ilogb
libc.src.math.ilogbf libc.src.math.ilogbf

1
libc/config/linux/x86_64/entrypoints.txt
View File

@ -104,6 +104,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.frexp libc.src.math.frexp
libc.src.math.frexpf libc.src.math.frexpf
libc.src.math.frexpl libc.src.math.frexpl
libc.src.math.hypot
libc.src.math.hypotf libc.src.math.hypotf
libc.src.math.ilogb libc.src.math.ilogb
libc.src.math.ilogbf libc.src.math.ilogbf

1
libc/spec/stdc.td
View File

@ -280,6 +280,7 @@ def StdC : StandardSpec<"stdc"> {
FunctionSpec<"frexpf", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<IntPtr>]>, FunctionSpec<"frexpf", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<IntPtr>]>,
FunctionSpec<"frexpl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>, 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<"hypotf", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<FloatType>]>,
FunctionSpec<"ilogb", RetValSpec<IntType>, [ArgSpec<DoubleType>]>, FunctionSpec<"ilogb", RetValSpec<IntType>, [ArgSpec<DoubleType>]>,

12
libc/src/math/CMakeLists.txt
View File

@ -713,3 +713,15 @@ add_entrypoint_object(
COMPILE_OPTIONS COMPILE_OPTIONS
-O2 -O2
) )
add_entrypoint_object(
hypot
SRCS
hypot.cpp
HDRS
hypot.h
DEPENDS
libc.utils.FPUtil.fputil
COMPILE_OPTIONS
-O2
)

18
libc/src/math/hypot.cpp Normal file
View File

@ -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

18
libc/src/math/hypot.h Normal file
View File

@ -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

209
libc/src/math/hypotf.cpp
View File

@ -6,217 +6,12 @@
// //
//===----------------------------------------------------------------------===// //===----------------------------------------------------------------------===//
#include "src/__support/common.h" #include "src/__support/common.h"
#include "utils/FPUtil/BasicOperations.h" #include "utils/FPUtil/Hypot.h"
#include "utils/FPUtil/FPBits.h"
namespace __llvm_libc { 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) { float LLVM_LIBC_ENTRYPOINT(hypotf)(float x, float y) {
FPBits<float> x_bits(x), y_bits(y); return __llvm_libc::fputil::hypot(x, 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);
} }
} // namespace __llvm_libc } // namespace __llvm_libc

13
libc/test/src/math/CMakeLists.txt
View File

@ -736,3 +736,16 @@ add_fp_unittest(
libc.src.math.hypotf libc.src.math.hypotf
libc.utils.FPUtil.fputil 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
)

75
libc/test/src/math/HypotTest.h Normal file
View File

@ -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

20
libc/test/src/math/hypot_test.cpp Normal file
View File

@ -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); }

55
libc/test/src/math/hypotf_test.cpp
View File

@ -6,56 +6,15 @@
// //
//===----------------------------------------------------------------------===// //===----------------------------------------------------------------------===//
#include "HypotTest.h"
#include "include/math.h"
#include "src/math/hypotf.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 HypotfTest = HypotTestTemplate<float>;
using UIntType = FPBits::UIntType;
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) { TEST_F(HypotfTest, NormalRange) { testNormalRange(&__llvm_libc::hypotf); }
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);
}
}

267
libc/utils/FPUtil/Hypot.h Normal file
View File

@ -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