libc/src/__support/math/powf_float_eval.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
/// Float-only evaluation implementation for powf(x, y).
///
//===----------------------------------------------------------------------===//

#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_POWF_FLOAT_EVAL_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_POWF_FLOAT_EVAL_H

#include "src/__support/CPP/bit.h"
#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/double_double.h"
#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/FPUtil/nearest_integer.h"
#include "src/__support/common.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h"
#include "src/__support/macros/properties/cpu_features.h"
#include "src/__support/math/exp2f_float_utils.h"
#include "src/__support/math/powf_utils.h"

namespace LIBC_NAMESPACE_DECL {
namespace math {
namespace float_eval {

namespace powf_internal = LIBC_NAMESPACE::math::powf_internal;

// Overview of float-precision pow(x, y) = x^y computation:
//
// Let x = 2^(e_x) * m_x > 0.  Then:
//   x^y = 2^( y * log2(x) )
//       = 2^( y * ( e_x + log2(m_x) ) )
//       = 2^( k + rem )
//       = 2^k * 2^rem,
// where:
//   k = round(y * log2(x)),
//   rem = y * log2(x) - k.
//
// In particular, k is an integer, and |rem| <= 0.5.
//
// For the final result to fit in single precision, the exponent field is
// bounded by:
//   -150 <= y * log2(x) <= 128,
// and anything outside that range overflows or underflows.
//
// Since |rem| <= 0.5:
//   2^-0.5 <= 2^rem <= 2^0.5,
// the relative error of x^y = 2^k * 2^rem is approximately:
//   ~ absolute error of rem * ln(2)
//   ~ absolute error of y * log2(x) * ln(2) + rel_err(exp2f_eval)
//   ~ relative error(log2(x)) * |y * log2(x)| * ln(2) + 2^-25.
//
// With exp2f_eval providing error ~ 2^-25, the budget for
// relative error(log2(x)) * |y * log2(x)| * ln(2) is <= 2^-25.
// For single-precision results, the 8-bit exponent field bounds:
//   |y * log2(x)| < 2^8 = 256.
// Thus:
//   relative error(log2(x)) <= 2^-25 / (2^8 * ln(2)) < 2^-33.
//
// To compute log2(x), we perform range reduction for log2(m_x) using a compact
// 5-bit lookup table (32 entries):
//   dx = r * m_x - 1 in [-2^-6, 2^-5].
// Then m_x = (1 + dx) / r, and:
//   log2(m_x) = log2(1 + dx) - log2(r),
// where -log2(r) is obtained from the lookup table LOG2_R_FF_32 with
// FloatFloat precision (~48 bits). The computation of dx = r * m_x - 1 is
// exact.
//
// The expansion of log2(1 + dx) is:
//   log2(1 + dx) = dx * log2(e) + dx^2 * Q(dx).
//
// We compute only the leading term dx * log2(e) in FloatFloat using
// exact_mult, and evaluate the tail Q(dx) using a degree-4 minimax polynomial
// purely in single precision:
//   Q(dx) ~ (log2(1 + dx) - dx * log2(e)) / dx^2.
// Generated by Sollya with:
// > log2_e = 0x1.715476p0 + 0x1.4ae0cp-26;
// > f = (log2(1 + x) - x * log2_e) / x^2;
// > Q = fpminimax(f, 4, [|SG...|], [-2^-6, 2^-5]);
// > dirtyinfnorm(((log2(1 + x) - x * log2_e) - x^2 * Q) / x, [-2^-6, 2^-5]);
//   0x1.0af1...p-34
//
// Finally, we compute the product y * log2(x) in FloatFloat:
//   y * log2(x) = y * ((e_x - log2(r)) + log2(1 + dx)).
// When x in [1 - 2^-5, 1) (where e_x = -1 and R[31] = 0.5 so
// -log2(R[31]) = 1.0f), e_x - log2(r) cancels to 0, avoiding cancellation.
// We decompose into integer k and fractional remainder rem:
//   k = round(total_hi),
//   rem = (total_hi - k) + total_lo,
// and evaluate 2^k * 2^rem via exp2f_eval(rem, k) from exp2f_float_utils.h,
// which provides a compact, code-size-friendly implementation without large
// exponent lookup tables.

LIBC_INLINE float powf(float x, float y) {
  using FPBits = fputil::FPBits<float>;
  using namespace powf_internal;

  FPBits xbits(x), ybits(y);
  uint32_t x_u = xbits.uintval();
  uint32_t y_u = ybits.uintval();
  uint32_t y_a = ybits.abs().uintval();

  if (LIBC_UNLIKELY((x_u & 0x007F'FFFF) == 0 || (y_u & 0x007F'FFFF) == 0)) {
    if (auto r = powf_internal::check_special_inputs(x, y);
        LIBC_UNLIKELY(r.has_value()))
      return r.value();
  }

  int e_x = 0;
  uint32_t sign = 0;
  if (LIBC_UNLIKELY(y_a <= powf_internal::Y_LOWER_BOUND ||
                    y_a >= powf_internal::Y_UPPER_BOUND ||
                    x_u >= FPBits::inf().uintval() ||
                    x_u < FPBits::min_normal().uintval())) {
    if (auto r = powf_internal::check_exceptional_cases(x, y, e_x, sign);
        LIBC_UNLIKELY(r.has_value()))
      return r.value();
  }

  Sign out_sign = (sign == 0) ? Sign::POS : Sign::NEG;

  // Extract exponent field and mantissa of x.
  xbits = FPBits(x);
  e_x += xbits.get_biased_exponent() - FPBits::EXP_BIAS;
  uint32_t x_mant = xbits.get_mantissa();
  unsigned idx_x = static_cast<unsigned>(x_mant >> (FPBits::FRACTION_LEN - 5));
  FPBits m_x = FPBits(x_mant | 0x3f80'0000U);

  // log2(e) in FloatFloat:
  constexpr fputil::FloatFloat LOG2_E = {0x1.4ae0cp-26f, 0x1.715476p0f};

  // Degree-4 Sollya minimax polynomial for:
  //   Q(dx) ~ (log2(1 + dx) - dx * log2(e)) / dx^2
  // Generated by Sollya with:
  // > log2_e = 0x1.715476p0 + 0x1.4ae0cp-26;
  // > f = (log2(1 + x) - x * log2_e) / x^2;
  // > Q = fpminimax(f, 4, [|SG...|], [-2^-6, 2^-5]);
  // > dirtyinfnorm(((log2(1 + x) - x * log2_e) - x^2 * Q) / x, [-2^-6, 2^-5]);
  //   0x1.0af1...p-34
  constexpr float COEFFS[] = {-0x1.715476p-1f, 0x1.ec709ep-2f, -0x1.71618ap-2f,
                              0x1.2810e4p-2f, -0x1.ae42b2p-3f};

  // Range reduction for log2(m_x):
  //   dx = r * m_x - 1 in [-2^-6, 2^-5).
  // The computation is exact: m_x = (1 + dx) / r.
  float r = powf_internal::R_32[idx_x];
#if defined(LIBC_TARGET_CPU_HAS_FMA_FLOAT)
  float dx = fputil::multiply_add(r, m_x.get_val(), -1.0f); // Exact
#else  // !LIBC_TARGET_CPU_HAS_FMA_FLOAT
  // Without FMA, split m_x = c + cd where c retains the high 14 bits
  // to perform exact multiplication: r * c - 1.0f and r * cd.
  float c = FPBits(m_x.uintval() & 0x3fff'e000U).get_val();
  float cd = m_x.get_val() - c;
  float dx_hi = fputil::multiply_add(r, c, -1.0f);
  float dx = fputil::multiply_add(r, cd, dx_hi); // Exact
#endif // LIBC_TARGET_CPU_HAS_FMA_FLOAT

  fputil::FloatFloat log2_r = powf_internal::LOG2_R_FF_32[idx_x];
  float e_xf = static_cast<float>(e_x);
  // (e_xf + log2_r.hi) is exact since log2_r.hi has 16 bits of precision and
  // |e_x| <= 150 < 2^8.
  float e_x_hi = e_xf + log2_r.hi;

  // Early overflow and underflow check:
  //   z = y * log2(x) ~ y * (e_x_hi + dx * log2(e))
  // With |dx| < 2^-5, |log2(1 + dx) - dx * log2(e)| <= (dx^2 / 2) * log2(e) <
  // 2^-10. For |y * log2(x)| <= 128, the error |z - y * log2(x)| < 4. Hence:
  //   z > 132  => overflow
  //   z < -155 => underflow
  float z = y * fputil::multiply_add(dx, LOG2_E.hi, e_x_hi);

  if (LIBC_UNLIKELY(z > 132.0f))
    return set_overflow(out_sign);

  if (LIBC_UNLIKELY(z < -155.0f))
    return set_underflow(out_sign);

  // Approximate the tail Q(dx) using Estrin's scheme:
  float dx2 = dx * dx;
  float dx4 = dx2 * dx2;

  float p0 = fputil::multiply_add(dx, COEFFS[1], COEFFS[0]);
  float p1 = fputil::multiply_add(dx, COEFFS[3], COEFFS[2]);

  float q0 = fputil::multiply_add(dx2, p1, p0);
  float p = fputil::multiply_add(dx4, COEFFS[4], q0);
  float p_tail = dx2 * p;

  // Compute dx * log2(e) in FloatFloat.
  fputil::FloatFloat dx_log2_e = fputil::exact_mult<float>(dx, LOG2_E.hi);
  float dx_log2_e_lo = fputil::multiply_add(dx, LOG2_E.lo, dx_log2_e.lo);

  // Combine with log2(1 + dx) in FloatFloat:
  fputil::FloatFloat log2_x_hi =
      fputil::exact_add<false, float>(e_x_hi, dx_log2_e.hi);
  float log2_x_lo = (dx_log2_e_lo + p_tail) + log2_r.lo + log2_x_hi.lo;
  fputil::FloatFloat log2_x =
      fputil::exact_add<false, float>(log2_x_hi.hi, log2_x_lo);

  fputil::FloatFloat s = fputil::exact_mult<float>(y, log2_x.hi);
  float total_hi = s.hi;
  float total_lo = fputil::multiply_add(y, log2_x.lo, s.lo);

  // Exponent reduction for 2^(total_hi + total_lo):
  //   k = round(total_hi)
  //   rem = (total_hi - k) + total_lo, with |rem| <= 0.5.
  float kf = fputil::nearest_integer(total_hi);
  int64_t k = static_cast<int64_t>(kf);

  float rem_hi = total_hi - kf;
  float rem = rem_hi + total_lo;

  // Normal range: -125 <= k <= 127.
  uint64_t k_u = static_cast<uint64_t>(k + 125);
  if (LIBC_LIKELY(k_u <= 252)) {
    float result = exp2f_eval(rem, static_cast<int>(k));
    return cpp::bit_cast<float>(cpp::bit_cast<uint32_t>(result) | sign);
  }

  // Exceptional path: k > 128 (overflow), k < -150 (underflow), or subnormal.
  if (k > 128)
    return set_overflow(out_sign);

  if (k < -150)
    return set_underflow(out_sign);

  // Subnormal range: -150 <= k <= -126.
  if (k <= -126) {
    float result = exp2f_eval(rem, static_cast<int>(k));
    float res = cpp::bit_cast<float>(cpp::bit_cast<uint32_t>(result) | sign);

    if (LIBC_UNLIKELY(FPBits(res).is_normal()))
      return res;

    fputil::set_errno_if_required(ERANGE);
    fputil::raise_underflow_except_if_required<float>();
    return res;
  }

  // Boundary case: k = 128.
  float result = exp2f_eval(rem, static_cast<int>(k));
  if (LIBC_UNLIKELY(FPBits(result).is_inf()))
    return set_overflow(out_sign);

  return cpp::bit_cast<float>(cpp::bit_cast<uint32_t>(result) | sign);
}

} // namespace float_eval
} // namespace math
} // namespace LIBC_NAMESPACE_DECL

#endif // LLVM_LIBC_SRC___SUPPORT_MATH_POWF_FLOAT_EVAL_H