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

#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_POWF_DOUBLE_EVAL_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_POWF_DOUBLE_EVAL_H

#include "src/__support/CPP/bit.h"
#include "src/__support/CPP/optional.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/common_constants.h"
#include "src/__support/math/exp_constants.h"
#include "src/__support/math/powf_utils.h"

namespace LIBC_NAMESPACE_DECL {
namespace math {
namespace double_eval {

namespace powf_internal = LIBC_NAMESPACE::math::powf_internal;

#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// Accurate evaluation in DoubleDouble precision when Ziv's test fails.
LIBC_INLINE double powf_accurate(float x, float y, double e_x, uint64_t sign) {
  using FloatBits = fputil::FPBits<float>;
  using DoubleBits = fputil::FPBits<double>;
  using fputil::DoubleDouble;

  uint32_t x_u = FloatBits(x).uintval();
  uint32_t x_mant = x_u & FloatBits::FRACTION_MASK;
  int idx_x = static_cast<int>(x_mant >> (FloatBits::FRACTION_LEN - 7));
  // Append hidden bit.
  float m_x = cpp::bit_cast<float>(x_mant | 0x3f80'0000);
  double dx = fputil::multiply_add(
      static_cast<double>(m_x),
      static_cast<double>(common_constants_internal::R[idx_x]), -1.0);

  // Degree-11 minimax polynomial for log2(1 + x)/x on [-2^-8, 2^-7] in
  // double-double, generated by Sollya with:
  // > display = hexadecimal;
  // > prec = 256;
  // > P = fpminimax(log2(1 + x)/x, 11, [|DD...|], [-2^(-8), 2^(-7)]);
  // > for i from 0 to 11 do {
  //     c = coeff(P, i);
  //     hi = round(c, D, RN);
  //     lo = round(c - hi, D, RN);
  //     print("{", lo, ", ", hi, "},");
  //   };
  // > dirtyinfnorm((log2(1 + x) - x * P) / x, [-2^(-8), 2^(-7)]);
  //   0x1.c598...p-104 < 2^-103.
  constexpr DoubleDouble LOG2_COEFFS[] = {
      {0x1.777d0ffda0d34p-56, 0x1.71547652b82fep0},
      {-0x1.777d0ffd88ca8p-57, -0x1.71547652b82fep-1},
      {0x1.d27f052d9d41dp-56, 0x1.ec709dc3a03fdp-2},
      {-0x1.777f259c6c0ap-58, -0x1.71547652b82fep-2},
      {0x1.e5c055599afdbp-56, 0x1.2776c50ef9bfep-2},
      {0x1.e923056eef156p-58, -0x1.ec709dc3a03fdp-3},
      {0x1.1f24e26b5ad0dp-57, 0x1.a61762a7add76p-3},
      {0x1.16d9b5013b55ep-58, -0x1.71547652bf304p-3},
      {-0x1.1551148c3d1d8p-57, 0x1.484b13f00425ep-3},
      {-0x1.a1cf327c6703p-60, -0x1.2776ca937f335p-3},
      {-0x1.d24f6254c2e6ap-57, 0x1.0c9211efe6d7ep-3},
      {0x1.a81cff2946d7cp-59, -0x1.e1f45baaf7f1dp-4},
  };

  // Evaluate P(dx) ~ log2(1 + dx) / dx using Horner's scheme.
  DoubleDouble p = LOG2_COEFFS[11];
  for (int i = 10; i >= 0; --i)
    p = fputil::add(LOG2_COEFFS[i], fputil::quick_mult(dx, p));

  // log2(1 + dx) = dx * P(dx) in DoubleDouble.
  DoubleDouble log2_1p = fputil::quick_mult(dx, p);

  // -log2(r) is represented in TripleDouble as {lo, mid, hi}.
  // We decompose log2(x) = e_x - log2(r) + log2(1 + dx) into:
  //   log2(x) = (e_x + hi) + (log2_1p + mid + lo)
  //           = (e_x + hi) + log2_tail.
  // The lower part log2_tail accumulates log2_1p, mid, and lo in DoubleDouble:
  DoubleDouble log2_tail =
      fputil::add<false>(log2_1p, powf_internal::LOG2_R_TD[idx_x].mid);
  log2_tail = fputil::add<false>(log2_tail, powf_internal::LOG2_R_TD[idx_x].lo);

  // Range reduction for 2^(y * log2(x)):
  // Scale by 2^6 = 64 to use the 64-entry EXP2_MID1 table:
  //   64 * y * log2(x) = y6 * (e_x + hi) + y6 * log2_tail.
  //
  // Since e_x + hi has at most 8 + 20 = 28 bits of precision and y has 24 bits,
  // the product y6 * (e_x + hi) has at most 28 + 24 = 52 bits, making y6_hi
  // exact in double precision.
  constexpr double SCALE = 0x1.0p6;
  double y6 = static_cast<double>(y) * SCALE;
  DoubleDouble prod = fputil::quick_mult(y6, log2_tail);

  double e_x_hi = e_x + powf_internal::LOG2_R_TD[idx_x].hi;
  DoubleDouble y6_hi = fputil::exact_mult(y6, e_x_hi);

  // hm = round(64 * y * log2(x)) is the nearest integer index.
  double hm_d = y6_hi.hi + prod.hi;
  double hm = fputil::nearest_integer(hm_d);

  // lo6 = 64 * y * log2(x) - hm = (y6_hi - hm) + prod:
  DoubleDouble lo6_hi = fputil::exact_add<false>(y6_hi.hi, -hm);
  lo6_hi.lo += y6_hi.lo;

  DoubleDouble lo6 = fputil::add<false>(prod, lo6_hi);

  // lo = lo6 / 64, with |lo| <= 2^-7.
  DoubleDouble lo = fputil::quick_mult(0x1.0p-6, lo6);

  int64_t hm_i = static_cast<int64_t>(hm);
  int idx_y = static_cast<int>(hm_i & 0x3f);
  int hi = static_cast<int>(hm_i >> 6);
  int64_t exp_hi_i = static_cast<int64_t>(static_cast<uint64_t>(hi)
                                          << DoubleBits::FRACTION_LEN);

  // Degree-8 minimax polynomial for (2^x - 1)/x on [-2^-7, 2^-7] in
  // double-double, generated by Sollya with:
  // > display = hexadecimal;
  // > prec = 256;
  // > P = fpminimax((2^x - 1)/x, 8, [|DD...|], [-2^(-7), 2^(-7)]);
  // > for i from 0 to 8 do {
  //     c = coeff(P, i);
  //     hi = round(c, D, RN);
  //     lo = round(c - hi, D, RN);
  //     print("{", lo, ", ", hi, "},");
  //   };
  // > dirtyinfnorm((2^x - (1 + x * P)) / (2^x), [-2^(-7), 2^(-7)]);
  //   0x1.e62d...p-106 < 2^-105.
  constexpr DoubleDouble EXP2_COEFFS_DD[] = {
      {0x1.abc9e3b39803ep-56, 0x1.62e42fefa39efp-1},
      {-0x1.5e43a54066432p-57, 0x1.ebfbdff82c58fp-3},
      {-0x1.d331626e9fe04p-59, 0x1.c6b08d704a0cp-5},
      {0x1.4f492022eeb74p-62, 0x1.3b2ab6fba4e77p-7},
      {0x1.035733b00abdp-66, 0x1.5d87fe78a6731p-10},
      {0x1.9e4f81ec951dep-68, 0x1.430912f86c123p-13},
      {0x1.7be627c299948p-74, 0x1.ffcbfc588b997p-17},
      {-0x1.c8a006d5b39acp-74, 0x1.62c03345a64f6p-20},
      {-0x1.fd14495e4b96p-81, 0x1.b5253eef8c0a9p-24},
  };

  DoubleDouble exp2_p = EXP2_COEFFS_DD[8];
  for (int i = 7; i >= 0; --i)
    exp2_p = fputil::add(EXP2_COEFFS_DD[i], fputil::quick_mult(lo, exp2_p));

  DoubleDouble p_lo = fputil::quick_mult(lo, exp2_p);
  DoubleDouble exp2_lo = fputil::exact_add(1.0, p_lo.hi);
  exp2_lo.lo += p_lo.lo;

  DoubleDouble exp2_hi_mid_dd;
  exp2_hi_mid_dd.hi = cpp::bit_cast<double>(
      exp_hi_i + cpp::bit_cast<int64_t>(EXP2_MID1[idx_y].hi) + sign);
  exp2_hi_mid_dd.lo =
      (idx_y != 0)
          ? cpp::bit_cast<double>(
                exp_hi_i + cpp::bit_cast<int64_t>(EXP2_MID1[idx_y].mid) + sign)
          : 0.0;

  DoubleDouble rr = fputil::quick_mult(exp2_hi_mid_dd, exp2_lo);

  DoubleDouble r = fputil::exact_add(rr.hi, rr.lo);

  // Round to odd to avoid double rounding when converting to float.
  uint64_t r_bits = cpp::bit_cast<uint64_t>(r.hi);
  if (LIBC_UNLIKELY(((r_bits & 0x0fff'ffffULL) == 0) && (r.lo != 0.0))) {
    if (DoubleBits(r.hi).sign() == DoubleBits(r.lo).sign()) {
      ++r_bits;
    } else if ((r_bits & DoubleBits::FRACTION_MASK) > 0) {
      --r_bits;
    }
  }

  return cpp::bit_cast<double>(r_bits);
}
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS

#if !defined(LIBC_MATH_HAS_SKIP_ACCURATE_PASS) &&                              \
    !defined(LIBC_MATH_HAS_SMALL_TABLES)
// Check if x^y is an exact rounding boundary.
// When x^y = exact_m * 2^exact_exp with exact_m <= 2^25:
// - exact_m has at most 26 bits, so it fits in double precision without
// rounding.
// - For single precision, exact_exp is in [-200, 127], well within the normal
//   range [-1022, 1023] of double precision.
// We compute exact_d = exact_m * 2^exact_exp in double precision and cast to
// float for the final rounded result.
LIBC_INLINE cpp::optional<float>
check_exact_boundary(float x, float y, uint64_t sign,
                     Sign out_sign = Sign::POS) {
  uint32_t exact_m = 0;
  int exact_exp = 0;
  if (LIBC_UNLIKELY(powf_internal::is_exact_rounding_boundary(x, y, exact_m,
                                                              exact_exp))) {
    // Number of bits in exact_m: 2^(l - 1) <= exact_m < 2^l.
    int l = 32 - cpp::countl_zero(exact_m);
    // Unbiased exponent of exact_m * 2^exact_exp.
    int unbiased_exp = exact_exp + l - 1;

    if (LIBC_UNLIKELY(unbiased_exp > 127))
      return powf_internal::set_overflow(out_sign);

    // Below the minimum subnormal float (2^-149).
    if (LIBC_UNLIKELY(unbiased_exp < -149)) {
      fputil::set_errno_if_required(ERANGE);
      fputil::raise_underflow_except_if_required<float>();
    }

    // Scale factor = sign * 2^exact_exp.
    using DoubleBits = fputil::FPBits<double>;
    double scale = cpp::bit_cast<double>(
        (static_cast<uint64_t>(exact_exp + DoubleBits::EXP_BIAS)
         << DoubleBits::FRACTION_LEN) |
        sign);
    double exact_d = static_cast<double>(exact_m) * scale;
    return static_cast<float>(exact_d);
  }
  return cpp::nullopt;
}
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS && !LIBC_MATH_HAS_SMALL_TABLES

// Overview of powf(x, y) = x^y computation in double precision:
//
// Let x = 2^(e_x) * m_x > 0.  Then:
//   x^y = 2^( y * log2(x) )
//       = 2^( y * ( e_x + log2(m_x) ) )
//       = 2^( k + f )
//       = 2^k * 2^f,
// where:
//   k = round(y * log2(x)),
//   f = y * log2(x) - k.
//
// In particular, k is an integer, and |f| <= 0.5.
// For the final result to fit in single precision:
//   -150 <= y * log2(x) <= 128,
// and anything outside that range overflows or underflows.
//
// Fast pass:
// 1. Range reduction for log2(m_x) using a 32-entry lookup table:
//      dx = r * m_x - 1 in [-2^-6, 2^-5],
//      log2(m_x) = log2(1 + dx) - log2(r),
//    where -log2(r) is retrieved from LOG2_R_32[idx_x], and log2(1 + dx)/dx is
//    approximated by a degree-6 polynomial in double precision with Estrin's
//    scheme.
// 2. Exponent reduction:
//      z = y * log2(x),
//      k = round(z),
//      f = z - k in [-0.5, 0.5].
// 3. 2^f is evaluated using a degree-7 polynomial.
//    Absorbing +-ERR into the constant term 1.0 +- ERR of the final FMA allows
//    computing the upper and lower bounds for Ziv's rounding test in parallel.
// 4. If upper == lower in single precision, return upper.
// 5. If Ziv's test fails, fall back to check_exact_boundary and powf_accurate
//    (DoubleDouble / 128-entry table).
//
// Accurate pass:
// 1. Range reduction using a 128-entry table: -log2(r) = {lo, mid, hi}.
//    log2(1 + dx) is evaluated using a degree-11 polynomial in DoubleDouble.
// 2. Decompose: log2(x) = (e_x + hi) + log2_tail, where:
//      log2_tail = log2(1 + dx) + mid + lo.
//    Because e_x + hi has <= 28 bits of precision and y has 24 bits,
//    y * (e_x + hi) is exact in double precision.
// 3. Exponent reduction with 64-entry EXP2_MID1 table and degree-8 polynomial
//    in DoubleDouble.
LIBC_INLINE float powf(float x, float y) {
  using FloatBits = fputil::FPBits<float>;
  using DoubleBits = fputil::FPBits<double>;
  using fputil::DoubleDouble;
  using namespace common_constants_internal;

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

  // Quick filter for special inputs:
  // - x: 0, +-1, 2^k, +-Inf, NaN.
  // - y: 0, +-1, +-2, +-0.5, +-Inf, NaN.
  if (LIBC_UNLIKELY((x_u & 0x001F'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();
  }

#if !defined(LIBC_MATH_HAS_SKIP_ACCURATE_PASS) &&                              \
    !defined(LIBC_MATH_HAS_SMALL_TABLES)
  float orig_x_abs = xbits.abs().get_val();
  float orig_y = y;
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS && !LIBC_MATH_HAS_SMALL_TABLES

  int ex = -FloatBits::EXP_BIAS;
  uint64_t sign = 0;

  // Check for exceptional cases:
  // - |y| <= 2^-40 or |y| >= 2^40.
  // - x < 0, x is subnormal, zero, Inf, or NaN.
  if (LIBC_UNLIKELY(y_a <= powf_internal::Y_LOWER_BOUND ||
                    y_a >= powf_internal::Y_UPPER_BOUND ||
                    x_u >= FloatBits::inf().uintval() ||
                    x_u < FloatBits::min_normal().uintval())) {
    if (auto r = powf_internal::check_exceptional_cases(x, y, ex, sign);
        LIBC_UNLIKELY(r.has_value()))
      return r.value();
  }

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

  // x^y = 2^( y * log2(x) ) = 2^( y * ( e_x + log2(m_x) ) )
  // Compute log2(x) = e_x + log2(m_x)
  x_u = FloatBits(x).uintval();

  // Extract exponent field of x.
  ex += (x_u >> FloatBits::FRACTION_LEN);
  double e_x = static_cast<double>(ex);
  uint32_t x_mant = x_u & FloatBits::FRACTION_MASK;
  // Top 5 bits of mantissa for 32-entry table lookup.
  int idx_x = static_cast<int>(x_mant >> (FloatBits::FRACTION_LEN - 5));
  // Embed mantissa into double format: m_x in [1.0, 2.0).
  uint64_t m_x_u =
      (static_cast<uint64_t>(x_mant) << 29) | 0x3ff0'0000'0000'0000ULL;
  double m_x = cpp::bit_cast<double>(m_x_u);

  // Range reduction: dx = m_x * r - 1.0 in [-2^-6, 2^-5].
  // Since r has <= 7 bits and m_x has 24 bits, m_x * r - 1.0 is exact in
  // double.
  double dx =
      fputil::multiply_add(m_x, powf_internal::R_32_D[idx_x], -1.0); // Exact

  // Degree-6 minimax polynomial approximation for log2(1 + dx)/dx on
  // [-2^-6, 2^-5]. Generated by Sollya with:
  // > display = hexadecimal;
  // > p6 = fpminimax(log2(1 + x)/x, 6, [|D...|], [-2^-6, 2^-5]);
  // > for i from 0 to 6 do print(round(coeff(p6, i), D, RN), ",");
  // > dirtyinfnorm((log2(1 + x) - x * p6) / log2(1 + x), [-2^-6, 2^-5]);
  //   0x1.050a...p-47
  // > dirtyinfnorm(log2(1 + x) - x * p6, [-2^-6, 2^-5]);
  //   0x1.712b...p-52
  constexpr double COEFFS[] = {
      0x1.71547652b831fp0,   -0x1.71547652b2ff5p-1, 0x1.ec709dbd02463p-2,
      -0x1.7154787f73adep-2, 0x1.2777b7f57b29fp-2,  -0x1.ec55ea50eaba7p-3,
      0x1.92c78e7d5a931p-3,
  };

  // Evaluate P(dx) ~ log2(1 + dx)/dx with Estrin's scheme:
  double dx2 = dx * dx;
  double lp0 = fputil::multiply_add(dx, COEFFS[1], COEFFS[0]);
  double lp1 = fputil::multiply_add(dx, COEFFS[3], COEFFS[2]);
  double lp2 = fputil::multiply_add(dx, COEFFS[5], COEFFS[4]);

  double dx4 = dx2 * dx2;
  double lq0 = fputil::multiply_add(dx2, lp1, lp0);
  double lq1 = fputil::multiply_add(dx2, COEFFS[6], lp2);

  double p = fputil::multiply_add(dx4, lq1, lq0);

  // log2(x) = e_x - log2(r) + dx * P(dx).
  double log2_r_ex = powf_internal::LOG2_R_32[idx_x] + e_x;
  double s = fputil::multiply_add(dx, p, log2_r_ex);

  double y_d = static_cast<double>(y);
  double z = y_d * s;

  // y * log2(x) = k + f, with k integer and |f| <= 0.5.
  double kd = fputil::nearest_integer(z);
  int64_t k = static_cast<int64_t>(kd);
  double f = fputil::multiply_add(y_d, s, -kd);

  // Degree-7 polynomial approximation P(f) ~ (2^f - 1)/f on [-0.5, 0.5]
  // Generated by Sollya with:
  // > display = hexadecimal;
  // > p = fpminimax((2^x - 1)/x, 7, [|D...|], [-0.5, 0.5]);
  // > for i from 0 to 7 do print(round(coeff(p, i), D, RN), ",");
  // > dirtyinfnorm(2^x - (1 + x * p), [-0.5, 0.5]);
  //   0x1.0496...p-39
  // > dirtyinfnorm((2^x - (1 + x * p)) / 2^x, [-0.5, 0.5]);
  //   0x1.0496...p-39
  constexpr double EXP2_COEFFS[] = {
      0x1.62e42fef9cdf7p-1,  0x1.ebfbdff85f49p-3,   0x1.c6b08da69c963p-5,
      0x1.3b2ab6a0e1172p-7,  0x1.5d87762dfc134p-10, 0x1.43099b4e829a9p-13,
      0x1.00c080f7699bep-16, 0x1.62c0108a065d8p-20,
  };

  // Evaluate P(f) ~ (2^f - 1)/f on [-0.5, 0.5] with Estrin's scheme:
  double f2 = f * f;
  double f4 = f2 * f2;

  double p0 = fputil::multiply_add(f, EXP2_COEFFS[1], EXP2_COEFFS[0]);
  double p1 = fputil::multiply_add(f, EXP2_COEFFS[3], EXP2_COEFFS[2]);
  double p2 = fputil::multiply_add(f, EXP2_COEFFS[5], EXP2_COEFFS[4]);
  double p3 = fputil::multiply_add(f, EXP2_COEFFS[7], EXP2_COEFFS[6]);

  double q0 = fputil::multiply_add(f2, p1, p0);
  double q1 = fputil::multiply_add(f2, p3, p2);

  double poly = fputil::multiply_add(f4, q1, q0);

  // Normal range: -125 <= k <= 128.
  // (k + 125) as unsigned <= 253 covers all normal outputs.
  uint64_t k_u = static_cast<uint64_t>(k + 125);

  if (LIBC_LIKELY(k_u <= 253)) {
    // Scale by 2^k.
    uint64_t exp2_k_i =
        (static_cast<uint64_t>(static_cast<int>(k) + DoubleBits::EXP_BIAS)
         << DoubleBits::FRACTION_LEN) |
        sign;
    double exp2_k = cpp::bit_cast<double>(exp2_k_i);

#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
    double pp = fputil::multiply_add(f, poly, 1.0);
    double r_d = pp * exp2_k;
    return static_cast<float>(r_d);
#else // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
#ifdef LIBC_TARGET_CPU_HAS_FMA
    // Absorbing +-ERR into 1.0 computes both bounds using two parallel FMAs.
    constexpr double ERR = 0x1.5p-39;
    double pp_hi = fputil::multiply_add(f, poly, 1.0 + ERR);
    double pp_lo = fputil::multiply_add(f, poly, 1.0 - ERR);
    double r_d = pp_hi * exp2_k;
#else  // !LIBC_TARGET_CPU_HAS_FMA
    // Without FMA, compute unbiased pp = f * poly + 1.0 directly.
    double pp = fputil::multiply_add(f, poly, 1.0);
    double r_d = pp * exp2_k;
#endif // LIBC_TARGET_CPU_HAS_FMA

    float res = static_cast<float>(r_d);

    // Since 2^k is an exact power of 2, testing if the bounds round to the same
    // float is equivalent to testing the fully scaled bounds.
#ifdef LIBC_TARGET_CPU_HAS_FMA
    float upper = static_cast<float>(pp_hi);
    float lower = static_cast<float>(pp_lo);
#else  // !LIBC_TARGET_CPU_HAS_FMA
    constexpr double ERR = 0x1.5p-39;
    float upper = static_cast<float>(pp + ERR);
    float lower = static_cast<float>(pp - ERR);
#endif // LIBC_TARGET_CPU_HAS_FMA

    if (LIBC_LIKELY(upper == lower))
      return res;

    // Accurate fallback when Ziv's test fails:
#ifndef LIBC_MATH_HAS_SMALL_TABLES
    if (auto r = check_exact_boundary(orig_x_abs, orig_y, sign, out_sign);
        LIBC_UNLIKELY(r.has_value()))
      return r.value();
#endif // !LIBC_MATH_HAS_SMALL_TABLES

    double r_dd = powf_accurate(x, y, e_x, sign);
    res = static_cast<float>(r_dd);
    if (LIBC_UNLIKELY(FloatBits(res).is_inf()))
      return powf_internal::set_overflow(out_sign);
    return res;
#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
  }

#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
  // Exceptional path: k > 128 (overflow), k < -155 (underflow), or denormal.
  if (k > 128)
    return powf_internal::set_overflow(out_sign);

  if (k < -155)
    return powf_internal::set_underflow(out_sign);

  // Denormal and underflow path for -155 <= k <= -126:
#ifndef LIBC_MATH_HAS_SMALL_TABLES
  if (auto r = check_exact_boundary(orig_x_abs, orig_y, sign, out_sign);
      LIBC_UNLIKELY(r.has_value()))
    return r.value();
#endif // !LIBC_MATH_HAS_SMALL_TABLES

  int64_t exp2_k_i = (static_cast<uint64_t>(k + DoubleBits::EXP_BIAS)
                      << DoubleBits::FRACTION_LEN) |
                     sign;
  double exp2_k = cpp::bit_cast<double>(exp2_k_i);

  // Add +-2^(-126 + (53 - 24)) = +-2^-97 to mimic single precision denormal
  // rounding in double precision.
  // For -155 <= k <= -126, the least significant bit of the double precision
  // sum (r_d + denorm_bias) is aligned at 2^(-97 - 52) = 2^-149, matching the
  // least significant bit of single precision denormals 2^(-126 - 23) = 2^-149.
  // This keeps the values in the double precision normal range while performing
  // rounding according to the current rounding mode.
  double denorm_bias = cpp::bit_cast<double>(0x39e0'0000'0000'0000ULL | sign);

  constexpr double ERR = 0x1.5p-39;

#ifdef LIBC_TARGET_CPU_HAS_FMA
  // Absorbing +-ERR into 1.0 computes both bounds using two parallel FMAs.
  double pp_hi = fputil::multiply_add(f, poly, 1.0 + ERR);
  double pp_lo = fputil::multiply_add(f, poly, 1.0 - ERR);
  double u_hi = fputil::multiply_add(pp_hi, exp2_k, denorm_bias);
  double u_lo = fputil::multiply_add(pp_lo, exp2_k, denorm_bias);
#else  // !LIBC_TARGET_CPU_HAS_FMA
  // Without FMA, compute unbiased pp = f * poly + 1.0 directly.
  double pp = fputil::multiply_add(f, poly, 1.0);
  double r_d = pp * exp2_k;
  double err = ERR * r_d;
  double u_hi = (r_d + err) + denorm_bias;
  double u_lo = (r_d - err) + denorm_bias;
#endif // LIBC_TARGET_CPU_HAS_FMA

  if (LIBC_LIKELY(u_hi == u_lo)) {
    if (LIBC_UNLIKELY(u_hi == denorm_bias))
      return powf_internal::set_underflow(out_sign);

    float res = static_cast<float>(u_hi - denorm_bias);
    if (LIBC_UNLIKELY(FloatBits(res).is_normal())) {
#ifdef LIBC_TARGET_CPU_HAS_FMA
      double pp = fputil::multiply_add(f, poly, 1.0);
#endif // LIBC_TARGET_CPU_HAS_FMA
      if (static_cast<float>(pp) < 1.0f)
        fputil::raise_underflow_except_if_required<float>();
      return res;
    }

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

  // Ziv's test failed for denormal input, fall back to accurate pass.
  double r_dd = powf_accurate(x, y, e_x, sign);
  float res = static_cast<float>(r_dd);
  if (LIBC_UNLIKELY(FloatBits(res).is_normal())) {
    if (static_cast<float>(r_dd * 0x1.0p126) < 1.0f)
      fputil::raise_underflow_except_if_required<float>();
    return res;
  }
  fputil::set_errno_if_required(ERANGE);
  fputil::raise_underflow_except_if_required<float>();
  return res;
#else  // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
  if (k > 128)
    return powf_internal::set_overflow(out_sign);
  return powf_internal::set_underflow(out_sign);
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
}

} // namespace double_eval
} // namespace math
} // namespace LIBC_NAMESPACE_DECL

#endif // LLVM_LIBC_SRC___SUPPORT_MATH_POWF_DOUBLE_EVAL_H