[libc] Add implementation for hypotf

[lntue]

Summary

Truncating the sum of squares, and then use shift-and-add algorithm to compute its square root.
Required MPFR testing infra is updated in https://reviews.llvm.org/D87514

[lntue]

Fix comments alignment and remove commented tests.

[hungptit]

libc/src/math/hypotf.cpp

88

Shall we use the return early pattern for this if block?

libc/test/src/math/hypotf_test.cpp

21

Does const/constexpr work for these global variables?

[lntue]

libc/test/src/math/hypotf_test.cpp

21

The template class FPBits<T> does not support constexpr yet.

[lntue]

Change to early return pattern.

[sivachandra]

libc/src/math/hypotf.cpp

2

With other implementations in libc/src, we have followed snake_case. I think we should do that here for consistency.

100

I understand why the use of a 64-bit integers is required here. But, I would ideally prefer to restrict a function which deals with 32-bit numbers not have to deal with wider numbers. It surely is not a problem on most modern architectures. So, we can revisit if it ever becomes a problem. FWIW, the AOR implementations also switch to double_t for higher precision with intermediate computations.

[lntue]

Change to snake_case.

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.hypotf
     libc.src.math.logb
     libc.src.math.logbf
     libc.src.math.logbl
     libc.src.math.modf
     libc.src.math.modff
     libc.src.math.modfl
     libc.src.math.round
     libc.src.math.roundf

libc/config/linux/api.td

def MathAPI : PublicAPI<"math.h"> {
    "fmaxf",
    "fmaxl",
    "fmin",
    "fminf",
    "fminl",
    "frexp",
    "frexpf",
    "frexpl",
+   "hypotf",
    "logb",
    "logbf",
    "logbl",
    "modf",
    "modff",
    "modfl",
    "expf",
    "exp2f",

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.hypotf
     libc.src.math.logb
     libc.src.math.logbf
     libc.src.math.logbl
     libc.src.math.modf
     libc.src.math.modff
     libc.src.math.modfl
     libc.src.math.remainderf
     libc.src.math.remainder

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<"hypotf", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<FloatType>]>,
+
           FunctionSpec<"logb", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
           FunctionSpec<"logbf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
           FunctionSpec<"logbl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>]>,
 
           FunctionSpec<"modf", RetValSpec<DoubleType>, [ArgSpec<DoubleType>, ArgSpec<DoublePtr>]>,
           FunctionSpec<"modff", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<FloatPtr>]>,
           FunctionSpec<"modfl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>, ArgSpec<LongDoublePtr>]>,
 

libc/src/math/CMakeLists.txt

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

libc/src/math/hypotf.h

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

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"
+
+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<float> xBits(x), yBits(y);
+
+  if (xBits.isInf() || yBits.isInf()) {
+    return FPBits<float>::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<float>::value + 2) ||
+        (y == 0)) {
+      return abs(x);
+    } else if ((yBits.exponent >=
+                xBits.exponent + MantissaWidth<float>::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<float>::value + 1);
+
+    aMant <<= 1;
+    bMant <<= 1;
+
+    uint32_t leadingOne;
+    int yMantWidth;
+    if (aExp != 0) {
+      leadingOne = One;
+      aMant |= One;
+      yMantWidth = MantissaWidth<float>::value + 1;
+    } else {
+      leadingOne = findLeadingOne(aMant, yMantWidth);
+    }
+
+    if (bExp != 0) {
+      bMant |= One;
+    }
+
+    aMantSq = static_cast<uint64_t>(aMant) * aMant;
+    bMantSq = static_cast<uint64_t>(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<float>::maxExponent) {
+          return FPBits<float>::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<uint32_t>(sum >> yMantWidth) - leadingOne;
+    uint32_t tailBits = static_cast<uint32_t>(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<float>::maxExponent) {
+        return FPBits<float>::inf();
+      }
+    }
+
+    Y |= static_cast<uint32_t>(outExp) << MantissaWidth<float>::value;
+    return *reinterpret_cast<float *>(&Y);
+  }
+}
+
+} // namespace __llvm_libc

libc/test/src/math/CMakeLists.txt

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

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 "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<float>;
+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);
+}
+
+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);
+  }
+}