llvm-project/libc/utils/FPUtil/Sqrt.h

//===-- Square root of IEEE 754 floating point numbers ----------*- 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_UTILS_FPUTIL_SQRT_H
#define LLVM_LIBC_UTILS_FPUTIL_SQRT_H

#include "FPBits.h"

#include "utils/CPP/TypeTraits.h"

namespace __llvm_libc {
namespace fputil {

namespace internal {

template <typename T>
static inline void normalize(int &exponent,
                             typename FPBits<T>::UIntType &mantissa);

template <> inline void normalize<float>(int &exponent, uint32_t &mantissa) {
  // Use binary search to shift the leading 1 bit.
  // With MantissaWidth<float> = 23, it will take
  // ceil(log2(23)) = 5 steps checking the mantissa bits as followed:
  // Step 1: 0000 0000 0000 XXXX XXXX XXXX
  // Step 2: 0000 00XX XXXX XXXX XXXX XXXX
  // Step 3: 000X XXXX XXXX XXXX XXXX XXXX
  // Step 4: 00XX XXXX XXXX XXXX XXXX XXXX
  // Step 5: 0XXX XXXX XXXX XXXX XXXX XXXX
  constexpr int nsteps = 5; // = ceil(log2(MantissaWidth))
  constexpr uint32_t bounds[nsteps] = {1 << 12, 1 << 18, 1 << 21, 1 << 22,
                                       1 << 23};
  constexpr int shifts[nsteps] = {12, 6, 3, 2, 1};

  for (int i = 0; i < nsteps; ++i) {
    if (mantissa < bounds[i]) {
      exponent -= shifts[i];
      mantissa <<= shifts[i];
    }
  }
}

template <> inline void normalize<double>(int &exponent, uint64_t &mantissa) {
  // Use binary search to shift the leading 1 bit similar to float.
  // With MantissaWidth<double> = 52, it will take
  // ceil(log2(52)) = 6 steps checking the mantissa bits.
  constexpr int nsteps = 6; // = ceil(log2(MantissaWidth))
  constexpr uint64_t bounds[nsteps] = {1ULL << 26, 1ULL << 39, 1ULL << 46,
                                       1ULL << 49, 1ULL << 51, 1ULL << 52};
  constexpr int shifts[nsteps] = {27, 14, 7, 4, 2, 1};

  for (int i = 0; i < nsteps; ++i) {
    if (mantissa < bounds[i]) {
      exponent -= shifts[i];
      mantissa <<= shifts[i];
    }
  }
}

#if !(defined(__x86_64__) || defined(__i386__))
template <>
inline void normalize<long double>(int &exponent, __uint128_t &mantissa) {
  // Use binary search to shift the leading 1 bit similar to float.
  // With MantissaWidth<long double> = 112, it will take
  // ceil(log2(112)) = 7 steps checking the mantissa bits.
  constexpr int nsteps = 7; // = ceil(log2(MantissaWidth))
  constexpr __uint128_t bounds[nsteps] = {
      __uint128_t(1) << 56,  __uint128_t(1) << 84,  __uint128_t(1) << 98,
      __uint128_t(1) << 105, __uint128_t(1) << 109, __uint128_t(1) << 111,
      __uint128_t(1) << 112};
  constexpr int shifts[nsteps] = {57, 29, 15, 8, 4, 2, 1};

  for (int i = 0; i < nsteps; ++i) {
    if (mantissa < bounds[i]) {
      exponent -= shifts[i];
      mantissa <<= shifts[i];
    }
  }
}
#endif

} // namespace internal

// Correctly rounded IEEE 754 SQRT with round to nearest, ties to even.
// Shift-and-add algorithm.
template <typename T,
          cpp::EnableIfType<cpp::IsFloatingPointType<T>::Value, int> = 0>
static inline T sqrt(T x) {
  using UIntType = typename FPBits<T>::UIntType;
  constexpr UIntType One = UIntType(1) << MantissaWidth<T>::value;

  FPBits<T> bits(x);

  if (bits.isInfOrNaN()) {
    if (bits.sign && (bits.mantissa == 0)) {
      // sqrt(-Inf) = NaN
      return FPBits<T>::buildNaN(One >> 1);
    } else {
      // sqrt(NaN) = NaN
      // sqrt(+Inf) = +Inf
      return x;
    }
  } else if (bits.isZero()) {
    // sqrt(+0) = +0
    // sqrt(-0) = -0
    return x;
  } else if (bits.sign) {
    // sqrt( negative numbers ) = NaN
    return FPBits<T>::buildNaN(One >> 1);
  } else {
    int xExp = bits.getExponent();
    UIntType xMant = bits.mantissa;

    // Step 1a: Normalize denormal input and append hiddent bit to the mantissa
    if (bits.exponent == 0) {
      ++xExp; // let xExp be the correct exponent of One bit.
      internal::normalize<T>(xExp, xMant);
    } else {
      xMant |= One;
    }

    // Step 1b: Make sure the exponent is even.
    if (xExp & 1) {
      --xExp;
      xMant <<= 1;
    }

    // After step 1b, x = 2^(xExp) * xMant, where xExp is even, and
    // 1 <= xMant < 4.  So sqrt(x) = 2^(xExp / 2) * y, with 1 <= y < 2.
    // Notice that the output of sqrt is always in the normal range.
    // To perform shift-and-add algorithm to find y, let denote:
    //   y(n) = 1.y_1 y_2 ... y_n, we can define the nth residue to be:
    //   r(n) = 2^n ( xMant - y(n)^2 ).
    // That leads to the following recurrence formula:
    //   r(n) = 2*r(n-1) - y_n*[ 2*y(n-1) + 2^(-n-1) ]
    // with the initial conditions: y(0) = 1, and r(0) = x - 1.
    // So the nth digit y_n of the mantissa of sqrt(x) can be found by:
    //   y_n = 1 if 2*r(n-1) >= 2*y(n - 1) + 2^(-n-1)
    //         0 otherwise.
    UIntType y = One;
    UIntType r = xMant - One;

    for (UIntType current_bit = One >> 1; current_bit; current_bit >>= 1) {
      r <<= 1;
      UIntType tmp = (y << 1) + current_bit; // 2*y(n - 1) + 2^(-n-1)
      if (r >= tmp) {
        r -= tmp;
        y += current_bit;
      }
    }

    // We compute one more iteration in order to round correctly.
    bool lsb = y & 1; // Least significant bit
    bool rb = false;  // Round bit
    r <<= 2;
    UIntType tmp = (y << 2) + 1;
    if (r >= tmp) {
      r -= tmp;
      rb = true;
    }

    // Remove hidden bit and append the exponent field.
    xExp = ((xExp >> 1) + FPBits<T>::exponentBias);

    y = (y - One) | (static_cast<UIntType>(xExp) << MantissaWidth<T>::value);
    // Round to nearest, ties to even
    if (rb && (lsb || (r != 0))) {
      ++y;
    }

    return *reinterpret_cast<T *>(&y);
  }
}

} // namespace fputil
} // namespace __llvm_libc

#if (defined(__x86_64__) || defined(__i386__))
#include "SqrtLongDoubleX86.h"
#endif // defined(__x86_64__) || defined(__i386__)

#endif // LLVM_LIBC_UTILS_FPUTIL_SQRT_H