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 @@ -64,6 +64,7 @@ libc.src.math.frexp libc.src.math.frexpf libc.src.math.frexpl + libc.src.math.hypotf libc.src.math.logb libc.src.math.logbf libc.src.math.logbl diff --git a/libc/config/linux/api.td b/libc/config/linux/api.td --- a/libc/config/linux/api.td +++ b/libc/config/linux/api.td @@ -191,6 +191,7 @@ "frexp", "frexpf", "frexpl", + "hypotf", "logb", "logbf", "logbl", 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 @@ -97,6 +97,7 @@ libc.src.math.frexp libc.src.math.frexpf libc.src.math.frexpl + libc.src.math.hypotf libc.src.math.logb libc.src.math.logbf libc.src.math.logbl diff --git a/libc/spec/stdc.td b/libc/spec/stdc.td --- a/libc/spec/stdc.td +++ b/libc/spec/stdc.td @@ -296,6 +296,8 @@ FunctionSpec<"frexpf", RetValSpec, [ArgSpec, ArgSpec]>, FunctionSpec<"frexpl", RetValSpec, [ArgSpec, ArgSpec]>, + FunctionSpec<"hypotf", RetValSpec, [ArgSpec, ArgSpec]>, + FunctionSpec<"logb", RetValSpec, [ArgSpec]>, FunctionSpec<"logbf", RetValSpec, [ArgSpec]>, FunctionSpec<"logbl", 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 @@ -593,3 +593,15 @@ COMPILE_OPTIONS -O2 ) + +add_entrypoint_object( + hypotf + SRCS + hypotf.cpp + HDRS + hypotf.h + DEPENDS + libc.utils.FPUtil.fputil + COMPILE_OPTIONS + -O2 +) diff --git a/libc/src/math/hypotf.h b/libc/src/math/hypotf.h new file mode 100644 --- /dev/null +++ b/libc/src/math/hypotf.h @@ -0,0 +1,18 @@ +//===-- Implementation header for hypotf ------------------------*- 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_HYPOTF_H +#define LLVM_LIBC_SRC_MATH_HYPOTF_H + +namespace __llvm_libc { + +float hypotf(float x, float y); + +} // namespace __llvm_libc + +#endif // LLVM_LIBC_SRC_MATH_HYPOTF_H diff --git a/libc/src/math/hypotf.cpp b/libc/src/math/hypotf.cpp new file mode 100644 --- /dev/null +++ b/libc/src/math/hypotf.cpp @@ -0,0 +1,222 @@ +//===-- 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" + +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 (n + 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 xBits(x), yBits(y); + + if (xBits.isInf() || yBits.isInf()) { + return FPBits::inf(); + } else if (xBits.isNaN()) { + return x; + } else if (yBits.isNaN()) { + return y; + } else { + uint16_t aExp, bExp, outExp; + uint32_t aMant, bMant; + uint64_t aMantSq, bMantSq; + bool stickyBits; + + if ((xBits.exponent >= yBits.exponent + MantissaWidth::value + 2) || + (y == 0)) { + return abs(x); + } else if ((yBits.exponent >= + xBits.exponent + MantissaWidth::value + 2) || + (x == 0)) { + yBits.sign = 0; + return abs(y); + } + + if (x >= y) { + aExp = xBits.exponent; + aMant = xBits.mantissa; + bExp = yBits.exponent; + bMant = yBits.mantissa; + } else { + aExp = yBits.exponent; + aMant = yBits.mantissa; + bExp = xBits.exponent; + bMant = xBits.mantissa; + } + + outExp = aExp; + + // Add an extra bit to simplify computing rounding bit of the final result. + constexpr uint32_t One = 1U << (MantissaWidth::value + 1); + + aMant <<= 1; + bMant <<= 1; + + uint32_t leadingOne; + int yMantWidth; + if (aExp != 0) { + leadingOne = One; + aMant |= One; + yMantWidth = MantissaWidth::value + 1; + } else { + leadingOne = findLeadingOne(aMant, yMantWidth); + } + + if (bExp != 0) { + bMant |= One; + } + + aMantSq = static_cast(aMant) * aMant; + bMantSq = static_cast(bMant) * bMant; + + // At this point, aExp >= bExp > aExp - 25, so in order to line up aSqMant + // and bSqMant, we need to shift bSqMant to the right by (aExp - bExp) 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 shiftLength = 2 * (aExp - bExp); + stickyBits = ((bMantSq & ((1ULL << shiftLength) - 1)) != 0); + bMantSq >>= shiftLength; + + uint64_t sum = aMantSq + bMantSq; + if (sum >= (1ULL << (2 * yMantWidth + 2))) { + // a^2 + b^2 >= 4* leadingOne^2, so we will need an extra bit to the left. + if (leadingOne == One) { + // For normal result, we discard the last 2 bits of the sum and increase + // the exponent. + stickyBits = stickyBits || ((sum & 0x3U) != 0); + sum >>= 2; + ++outExp; + if (outExp >= FPBits::maxExponent) { + return FPBits::inf(); + } + } else { + // For denormal result, we simply move the leading bit of the result to + // the left by 1. + leadingOne <<= 1; + ++yMantWidth; + } + } + + uint32_t Y = leadingOne; + uint32_t R = static_cast(sum >> yMantWidth) - leadingOne; + uint32_t tailBits = static_cast(sum) & (leadingOne - 1); + + for (uint32_t currentBit = leadingOne >> 1; currentBit; currentBit >>= 1) { + R = (R << 1) + ((tailBits & currentBit) ? 1 : 0); + uint32_t tmp = (Y << 1) + currentBit; // 2*y(n - 1) + 2^(-n-1) + if (R >= tmp) { + R -= tmp; + Y += currentBit; + } + } + + bool roundBit = Y & 1U; + bool lsb = Y & 2U; + + if (Y >= One) { + Y -= One; + + if (outExp == 0) { + outExp = 1; + } + } + + Y >>= 1; + + // Round to the nearest, tie to even. + if (roundBit && (lsb || stickyBits || (R != 0))) { + ++Y; + } + + if (Y >= (One >> 1)) { + Y -= One >> 1; + ++outExp; + if (outExp >= FPBits::maxExponent) { + return FPBits::inf(); + } + } + + Y |= static_cast(outExp) << MantissaWidth::value; + return *reinterpret_cast(&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 @@ -591,3 +591,16 @@ libc.src.math.remquol libc.utils.FPUtil.fputil ) + +add_fp_unittest( + hypotf_test + NEED_MPFR + SUITE + libc_math_unittests + SRCS + hypotf_test.cpp + DEPENDS + libc.include.math + libc.src.math.hypotf + libc.utils.FPUtil.fputil +) diff --git a/libc/test/src/math/hypotf_test.cpp b/libc/test/src/math/hypotf_test.cpp new file mode 100644 --- /dev/null +++ b/libc/test/src/math/hypotf_test.cpp @@ -0,0 +1,73 @@ +//===-- 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 "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" + +using FPBits = __llvm_libc::fputil::FPBits; +using UIntType = FPBits::UIntType; + +namespace mpfr = __llvm_libc::testing::mpfr; + +static const float zero = FPBits::zero(); +static const float negZero = FPBits::negZero(); +static const float nan = FPBits::buildNaN(1); +static const float inf = FPBits::inf(); +static const float negInf = FPBits::negInf(); + +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); + + // uint32_t a = 0x3c26c596U; + // uint32_t b = 0x45d93a69U; + // float x = *reinterpret_cast(&a); + // float y = *reinterpret_cast(&b); + // float result = __llvm_libc::hypotf(x, y); + // mpfr::BinaryInput input{x, y}; + // ASSERT_MPFR_MATCH(mpfr::Operation::Hypot, input, result, 0.5); +} + +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(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); + } +}