src/__support/mathvec/logf.h

//===----------------------------------------------------------------------===//
//
// 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
//
//===----------------------------------------------------------------------===//
///
/// \file
/// This file contains the implementation for single-precision SIMD log.
///
//===----------------------------------------------------------------------===//

#ifndef LLVM_LIBC_SRC___SUPPORT_MATHVEC_LOGF_H
#define LLVM_LIBC_SRC___SUPPORT_MATHVEC_LOGF_H

#include "src/__support/CPP/bit.h"
#include "src/__support/CPP/simd.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/macros/properties/cpu_features.h"
#include "src/__support/math/common_constants.h"
#include "src/__support/mathvec/common_constants.h"

namespace LIBC_NAMESPACE_DECL {

namespace mathvec {

template <size_t N>
LIBC_INLINE cpp::simd<float, N>
inline_logf(cpp::simd<uint32_t, N> ux, cpp::simd<int32_t, N> exponent_bias) {
  using namespace common_constants_internal;

  // Extract the 7 leading bits of the mantissa
  cpp::simd<uint32_t, N> index = (ux >> 16) & 0x7fU;

  // Compute unbiased exponent
  cpp::simd<int32_t, N> m = cpp::simd_cast<int32_t>(ux >> 23);
  m += exponent_bias;

  // Sets the exponent such that 1 <= u < 2
  cpp::simd<float, N> u =
      cpp::bit_cast<cpp::simd<float, N>>((ux & 0x007f'ffffU) | 0x3f80'0000U);

  cpp::simd<double, N> log_r =
      cpp::gather<cpp::simd<double, N>>(true, index, LOG_TABLE);

#ifdef LIBC_TARGET_CPU_HAS_FMA_FLOAT
  cpp::simd<float, N> r = cpp::gather<cpp::simd<float, N>>(true, index, R);
  cpp::simd<double, N> v = cpp::simd_cast<double>(
      cpp::multiply_add(u, r, cpp::simd<float, N>(-1.0f)));
#else
  cpp::simd<double, N> r = cpp::gather<cpp::simd<double, N>>(true, index, RD);
  cpp::simd<double, N> v = cpp::simd_cast<double>(u) * r - 1.0;
#endif // LIBC_TARGET_CPU_HAS_FMA_FLOAT

  // Degree-5 polynomial approximation of log generated by Sollya with:
  // > P = fpminimax(log(1 + x)/x, 4, [|1, D...|], [-2^-8, 2^-7]);
  constexpr cpp::simd<double, N> c0 = -0x1.000000000fe63p-1;
  constexpr cpp::simd<double, N> c1 = 0x1.555556e963c16p-2;
  constexpr cpp::simd<double, N> c2 = -0x1.000028dedf986p-2;
  constexpr cpp::simd<double, N> c3 = 0x1.966681bfda7f7p-3;
  constexpr cpp::simd<double, N> log2 = 0x1.62e42fefa39efp-1;

  // Evaluate log(1 + v) in double precision.
  cpp::simd<double, N> v_squared = v * v;
  cpp::simd<double, N> p2 = cpp::multiply_add(v, c3, c2);
  cpp::simd<double, N> p1 = cpp::multiply_add(v, c1, c0);

  // Accumulate the exponent independently of the polynomial's FMA chain.
  cpp::simd<double, N> exponent = cpp::simd_cast<double>(m);
  cpp::simd<double, N> p0 = cpp::multiply_add(exponent, log2, log_r) + v;
  cpp::simd<double, N> y_d =
      cpp::multiply_add(v_squared, cpp::multiply_add(v_squared, p2, p1), p0);
  cpp::simd<float, N> y = cpp::simd_cast<float>(y_d);

#if !(defined(LIBC_TARGET_CPU_HAS_FMA_FLOAT) ||                                \
      defined(LIBC_MATH_HAS_SKIP_ACCURATE_PASS))
  y = (ux == 0x6f31a8ec) ? 0x1.08b512p+6 : y;
  y = (ux == 0x79e7ec37) ? 0x1.43ff6ep+6 : y;
#endif // !(LIBC_TARGET_CPU_HAS_FMA_FLOAT || LIBC_MATH_HAS_SKIP_ACCURATE_PASS)

  return y;
}

template <size_t N>
LIBC_INLINE cpp::simd<float, N> special_case(cpp::simd<float, N> x) {
  using FPBits = fputil::FPBits<float>;
  cpp::simd<uint32_t, N> ux = cpp::bit_cast<cpp::simd<uint32_t, N>>(x);

  // If x is subnormal, multiply by 23 to produce a normalised x,
  // and subtract 23 from the exponent_bias to account for this scaling.
  cpp::simd<bool, N> is_subnormal = (ux < 0x0080'0000U);
  cpp::simd<uint32_t, N> normalized =
      cpp::bit_cast<cpp::simd<uint32_t, N>>(x * 0x1.0p23f);
  cpp::simd<int32_t, N> exponent_bias = is_subnormal ? -150 : -127;
  ux = is_subnormal ? normalized : ux;

  // Compute logf with the normalised input.
  cpp::simd<float, N> y = inline_logf(ux, exponent_bias);

  // Handle remaining specialc cases.
  cpp::simd<bool, N> is_negative = x < 0;
  cpp::simd<float, N> special = is_negative ? FPBits::quiet_nan().get_val() : x;
  special = x == 0 ? -FPBits::inf().get_val() : special;

  return ((ux >= 0x7f80'0000U) || (x == 0)) ? special : y;
}

template <size_t N>
LIBC_INLINE cpp::simd<float, N> logf(cpp::simd<float, N> x) {
  cpp::simd<uint32_t, N> ux = cpp::bit_cast<cpp::simd<uint32_t, N>>(x);

  // Handle zero, inf, nan, and subnormals
  cpp::simd<bool, N> is_special = (ux < 0x0080'0000U) || (ux >= 0x7f80'0000U);
  if (LIBC_UNLIKELY(cpp::any_of(is_special)))
    return special_case(x);

  return inline_logf(ux, cpp::simd<int32_t, N>(-127));
}
} // namespace mathvec
} // namespace LIBC_NAMESPACE_DECL

#endif // LLVM_LIBC_SRC___SUPPORT_MATHVEC_LOGF_H