libc/src/__support/math/exp_utils.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
/// Shared utilities and tables for double-precision exponential functions.
///
//===----------------------------------------------------------------------===//

#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP_UTILS_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP_UTILS_H

#include "hdr/errno_macros.h"
#include "hdr/fenv_macros.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/fp_bitwise.h"
#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/FPUtil/rounding_mode.h"
#include "src/__support/common.h"
#include "src/__support/macros/attributes.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h"
#include "src/__support/macros/properties/cpu_features.h"

#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
#include "src/__support/math/exp2_accurate_128.h"
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS

namespace LIBC_NAMESPACE_DECL {

// Rounding tests for 2^hi * (mid + lo) when the output might be denormal. We
// assume further that 1 <= mid < 2, mid + lo < 2, and |lo| << mid.
// Notice that, if 0 < x < 2^-1022,
//   double(2^-1022 + x) - 2^-1022 = double(x).
// So if we scale x up by 2^1022, we can use
//   double(1.0 + 2^1022 * x) - 1.0 to test how x is rounded in denormal range.
template <bool SKIP_ZIV_TEST = false>
LIBC_INLINE constexpr cpp::optional<double>
ziv_test_denorm(int hi, double mid, double lo, double err) {
  // Scaling factor = 1/(min normal number) = 2^1022
  int scale_up = hi + 1022;
  double mid_hi = fputil::scale_exp2(mid, scale_up);
  double lo_scaled = (lo != 0.0) ? fputil::scale_exp2(lo, scale_up) : 0.0;

  double extra_factor = 0.0;
  int scale_down = -1022;

  // Result is denormal if (mid_hi + lo_scale < 1.0).
  if ((1.0 - mid_hi) > lo_scaled) {
    // Extra rounding step is needed, which adds more rounding errors.
    err += 0x1.0p-52;
    extra_factor = 1.0;
    scale_down = -1023;
  }

  // By adding 1.0, the results will have similar rounding points as denormal
  // outputs.
  if constexpr (SKIP_ZIV_TEST) {
    double r = extra_factor + (mid_hi + lo_scaled);
    return fputil::scale_exp2(r, scale_down);
  } else {
    double err_scaled = fputil::scale_exp2(err, scale_up);

    double lo_u = lo_scaled + err_scaled;
    double lo_l = lo_scaled - err_scaled;

    double upper = extra_factor + (mid_hi + lo_u);
    double lower = extra_factor + (mid_hi + lo_l);

    if (LIBC_LIKELY(upper == lower))
      return fputil::scale_exp2(upper, scale_down);

    return cpp::nullopt;
  }
}

namespace math {
namespace double_eval {

using fputil::DoubleDouble;
using fputil::scale_exp2;

// log(2) in double-double precision, generated by Sollya with:
// > display = hexadecimal!;
// > prec = 500;
// > hi = round(log(2), D, RN);
// > lo = round(log(2) - hi, D, RN);
// > print("{", lo, ",", hi, "};");
LIBC_INLINE_VAR constexpr DoubleDouble LOG_2 = {0x1.abc9e3b39803fp-56,
                                                0x1.62e42fefa39efp-1};

// log2(e) in double-double precision, generated by Sollya with:
// > display = hexadecimal!;
// > prec = 500;
// > hi = round(log2(exp(1)), D, RN);
// > lo = round(log2(exp(1)) - hi, D, RN);
// > print("{", lo, ",", hi, "};");
LIBC_INLINE_VAR constexpr DoubleDouble LOG2_E = {0x1.777d0ffda0d24p-56,
                                                 0x1.71547652b82fep+0};

// Check for exceptional cases (tiny, overflow, underflow, inf, or nan) for
// exp, exp2, and exp10. Callers have already verified that x is in the exact
// tiny, underflow, or overflow/inf/nan range for the specific function before
// calling this helper.
LIBC_NOINLINE LIBC_INLINE double set_exceptional(double x) {
  using FPBits = fputil::FPBits<double>;
  FPBits xbits(x);

  // |x| <= 0x1.0p-53 (caller has already checked the exact tiny bound):
  constexpr uint64_t TINY_BOUND = 0x3ca0'0000'0000'0000ULL;
  if (xbits.abs().uintval() <= TINY_BOUND) {
    // b^x ~ 1 + x
    return 1.0 + x;
  }

  // Negative underflow or -inf/nan (caller has already excluded
  // non-underflowing negative inputs):
  if (xbits.is_neg()) {
    // b^(-Inf) = 0
    if (xbits.is_inf())
      return 0.0;

    // b^(nan) = nan
    // Adding 0.0 to silence signaling nan.
    if (xbits.is_nan())
      return x + 0.0;

    fputil::set_errno_if_required(ERANGE);
    fputil::raise_underflow_except_if_required<double>();
#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
    if (fputil::quick_get_round() == FE_UPWARD)
      return FPBits::min_subnormal().get_val();
#endif
    return 0.0;
  }

  // Positive overflow or +inf/nan (caller has already checked the exact
  // overflow bound)

  // x is finite and positive
  if (xbits.uintval() < FPBits::inf().uintval()) {
    fputil::set_errno_if_required(ERANGE);
    fputil::raise_overflow_except_if_required<double>();
#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
    int rounding = fputil::quick_get_round();
    if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
      return FPBits::max_normal().get_val();
#endif
    return FPBits::inf().get_val();
  }
  // x is +inf or nan
  // Adding `inf` to silence signaling nan.
  return x + FPBits::inf().get_val();
}

#if !defined(LIBC_MATH_HAS_SKIP_ACCURATE_PASS) ||                              \
    !defined(LIBC_MATH_HAS_SMALL_TABLES)

LIBC_INLINE_VAR constexpr int EXP_N_BITS = 7;
LIBC_INLINE_VAR constexpr double EXP_SCALE = 128.0; // 2^EXP_N_BITS = 2^7
LIBC_INLINE_VAR constexpr double EXP_SCALE_INV = 0x1.0p-7;
LIBC_INLINE_VAR constexpr int EXP_MASK = 0x7f;

// Polynomial approximation for (exp(x) - 1 - x) / x^2 on
// [-B, B] with B = 2^-8 * log(2):
//
// 1) A degree-3 minimax polynomial P0 for (exp(x) - 1 - x) / x^2 on [-B, B]
// (which minimizes the error of exp(x) - (1 + x + x^2 * P0(x)) weighted by
// 1 / x^2, i.e., relative error with respect to x^2):
// > display = hexadecimal!;
// > prec = 500;
// > B = 2^(-8) * log(2);
// > P0 = fpminimax((exp(x) - 1 - x)/x^2, 3, [|D...|], [-B, B]);
// > dirtyinfnorm(exp(x) - (1 + x + x^2 * P0), [-B, B]);
//   0x1.4362...p-64 < 2^-63.
//
// Some explanation for the errors (from Chebyshev polynomials):
// Because P0 minimizes |(exp(x) - 1 - x) / x^2 - P0(x)|, its error is governed
// by the degree-4 Chebyshev polynomial:
//   T_4(x / B) = 8*(x / B)^4 - 8*(x / B)^2 + 1
// approximating the leading omitted Taylor term x^4 / 720, with magnitude:
//   B^4 / (720 * 8).
// Multiplying by x^2 gives an absolute error of:
//   x^2 * B^4 / (720 * 8) * T_4(x / B),
// peaking at x = +-B with:
//   B^6 / (720 * 8) ~ 0x1.4309...p-64 < 2^-63.
//
// 2) If we instead minimize the absolute error of exp(x) - (1 + x + P(x))
// directly on [-B, B] with monomials [|2, 3, 4, 5|]:
// > P = fpminimax(exp(x) - 1 - x, [|2, 3, 4, 5|], [|D...|], [-B, B],
//                 absolute);
// > dirtyinfnorm(exp(x) - (1 + x + P), [-B, B]);
//   0x1.8e05...p-66 < 2^-65.
//
// Some explanation for the better errors from Chebyshev polynomials:
// The free x^2 and x^4 coefficients in P(x) approximate the omitted x^6 / 720
// term on [-B, B] subject to vanishing at x = 0. Setting t = (x / B)^2 in
// [0, 1], this is equivalent to approximating t^3 on [0, 1] by:
//   a_2 * t^2 + a_1 * t
// vanishing at t = 0.
// The degree-3 Chebyshev polynomial:
//   T_3(z) = 4*z^3 - 3*z
// has a root at z = -sqrt(3)/2 and three alternating extrema +-1 on
// [-sqrt(3)/2, 1] at z = -1/2, 1/2, 1.
// Mapping t in [0, 1] affinely onto z in [-sqrt(3)/2, 1] via:
//   z = ((2 + sqrt(3)) / 2) * t - sqrt(3) / 2
// gives an error curve T_3(z) that vanishes at t = 0 and equioscillates three
// times on (0, 1], with leading coefficient:
//   4 * ((2 + sqrt(3)) / 2)^3 = (2 + sqrt(3))^3 / 2.
// Dividing by this leading coefficient gives a maximum error on [-B, B] of:
//   (2 / (2 + sqrt(3))^3) * B^6 / 720 = 2 * (2 - sqrt(3))^3 * B^6 / 720
//                                     ~ 0x1.8dba...p-66 < 2^-65,
// reducing the maximum error compared to P0 by:
//   (2 + sqrt(3))^3 / 16 > 3x.
LIBC_INLINE_VAR constexpr double EXP_P_COEFFS[4] = {
    0x1.ffffffffffdbdp-2, // p_0
    0x1.555555555543dp-3, // p_1
    0x1.55555cf16d6eep-5, // p_2
    0x1.111116782220ap-7, // p_3
};

// Evaluate P(u) ~ (e^u - 1 - u) / u^2 using Estrin's scheme:
LIBC_ALWAYS_INLINE double exp_poly(double u, double u2) {
  double c0 = fputil::multiply_add(u, EXP_P_COEFFS[1], EXP_P_COEFFS[0]);
  double c1 = fputil::multiply_add(u, EXP_P_COEFFS[3], EXP_P_COEFFS[2]);
  return fputil::multiply_add(u2, c1, c0);
}

// Splitting -2^-7 * log(2) ~ hi + lo with 35-bit-precision hi part so that
// kd * MLOG_2_EXP2_HI is exact for |kd| < 2^18:
// > hi = round(-2^-7 * log(2), 35, RN);
// > lo = round(-2^-7 * log(2) - hi, D, RN);
LIBC_INLINE_VAR constexpr double MLOG_2_EXP2_HI = -0x1.62e42fefcp-8;
LIBC_INLINE_VAR constexpr double MLOG_2_EXP2_MID = 0x1.c610ca86c3899p-44;

#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS || !LIBC_MATH_HAS_SMALL_TABLES

#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS

// Fixed-point splitting constants (ulp = 2^-18, 2^-30) and pre-split
// 12-bit + 53-bit logarithm constants for exact non-FMA multiplication without
// Veltkamp splitting:
LIBC_INLINE_VAR constexpr double SPLIT_18 = 0x1.8p34; // 1.5 * 2^(-18 + 52)
LIBC_INLINE_VAR constexpr double SPLIT_30 = 0x1.8p22; // 1.5 * 2^(-30 + 52)

// Splitting log(2) ~ hi + lo with 12-bit-precision hi part.
// > hi = round(log(2), 12, RN);
// > lo = round(log(2) - hi, D, RN);
LIBC_INLINE_VAR constexpr double LOG_2_HI_12 = 0x1.62ep-1;
LIBC_INLINE_VAR constexpr double LOG_2_LO_12 = 0x1.0bfbe8e7bcd5ep-15;

// Error bounds for double precision fast path (N = 7, Ziv failure rate <
// 0.05%):
LIBC_INLINE_VAR constexpr double EXP_ERR_D = 0x1.4p-62;
LIBC_INLINE_VAR constexpr double EXP2_ERR_D = 0x1.0p-64;

// Lookup table for 2^(i * 2^-7) with i = 0..127.
// Generated by Sollya with:
// > for i from 0 to 127 do {
//     a = 2^(i * 2^-7);
//     b = round(a, D, RN);
//     c = round(a - b, D, RN);
//     print("{", c, ",", b, "},");
//   };
LIBC_INLINE_VAR constexpr DoubleDouble EXP2_MID[128] = {
    {0, 0x1p0},
    {0x1.b61299ab8cdb7p-54, 0x1.0163da9fb3335p0},
    {-0x1.19083535b085dp-56, 0x1.02c9a3e778061p0},
    {-0x1.0a31c1977c96ep-54, 0x1.04315e86e7f85p0},
    {0x1.d73e2a475b465p-55, 0x1.059b0d3158574p0},
    {-0x1.c91dfe2b13c27p-55, 0x1.0706b29ddf6dep0},
    {0x1.186be4bb284ffp-57, 0x1.0874518759bc8p0},
    {0x1.1487818316136p-54, 0x1.09e3ecac6f383p0},
    {0x1.8a62e4adc610bp-54, 0x1.0b5586cf9890fp0},
    {0x1.01edc16e24f71p-54, 0x1.0cc922b7247f7p0},
    {0x1.03a1727c57b53p-59, 0x1.0e3ec32d3d1a2p0},
    {-0x1.b9bedc44ebd7bp-57, 0x1.0fb66affed31bp0},
    {-0x1.6c51039449b3ap-54, 0x1.11301d0125b51p0},
    {-0x1.1b514b36ca5c7p-58, 0x1.12abdc06c31ccp0},
    {-0x1.32fbf9af1369ep-54, 0x1.1429aaea92dep0},
    {0x1.2406ab9eeab0ap-55, 0x1.15a98c8a58e51p0},
    {-0x1.19041b9d78a76p-55, 0x1.172b83c7d517bp0},
    {-0x1.11023d1970f6cp-54, 0x1.18af9388c8deap0},
    {0x1.e5b4c7b4968e4p-55, 0x1.1a35beb6fcb75p0},
    {-0x1.95386352ef607p-54, 0x1.1bbe084045cd4p0},
    {0x1.e016e00a2643cp-54, 0x1.1d4873168b9aap0},
    {-0x1.1df98027bb78cp-54, 0x1.1ed5022fcd91dp0},
    {0x1.dc775814a8495p-55, 0x1.2063b88628cd6p0},
    {0x1.2a97e9494a5eep-55, 0x1.21f49917ddc96p0},
    {0x1.9b07eb6c70573p-54, 0x1.2387a6e756238p0},
    {0x1.ac155bef4f4a4p-55, 0x1.251ce4fb2a63fp0},
    {0x1.2bd339940e9d9p-55, 0x1.26b4565e27cddp0},
    {-0x1.a4c3a8c3f0d7ep-54, 0x1.284dfe1f56381p0},
    {0x1.612e8afad1255p-55, 0x1.29e9df51fdee1p0},
    {-0x1.10adcd6381aa4p-59, 0x1.2b87fd0dad99p0},
    {0x1.0024754db41d5p-54, 0x1.2d285a6e4030bp0},
    {0x1.1ca0f45d52383p-56, 0x1.2ecafa93e2f56p0},
    {0x1.6f46ad23182e4p-55, 0x1.306fe0a31b715p0},
    {0x1.a9ce78e18047cp-55, 0x1.32170fc4cd831p0},
    {0x1.32721843659a6p-54, 0x1.33c08b26416ffp0},
    {-0x1.b5cee5c4e4628p-55, 0x1.356c55f929ff1p0},
    {-0x1.63aeabf42eae2p-54, 0x1.371a7373aa9cbp0},
    {-0x1.e958d3c9904bdp-54, 0x1.38cae6d05d866p0},
    {-0x1.5e436d661f5e3p-56, 0x1.3a7db34e59ff7p0},
    {-0x1.efff8375d29c3p-54, 0x1.3c32dc313a8e5p0},
    {0x1.ada0911f09ebcp-55, 0x1.3dea64c123422p0},
    {-0x1.7d023f956f9f3p-54, 0x1.3fa4504ac801cp0},
    {-0x1.ef3691c309278p-58, 0x1.4160a21f72e2ap0},
    {-0x1.1c7dde35f7999p-55, 0x1.431f5d950a897p0},
    {0x1.89b7a04ef80dp-59, 0x1.44e086061892dp0},
    {0x1.c944bd1648a76p-54, 0x1.46a41ed1d0057p0},
    {0x1.3c1a3b69062fp-56, 0x1.486a2b5c13cdp0},
    {0x1.9cb62f3d1be56p-54, 0x1.4a32af0d7d3dep0},
    {0x1.d4397afec42e2p-56, 0x1.4bfdad5362a27p0},
    {0x1.8ecdbbc6a7833p-54, 0x1.4dcb299fddd0dp0},
    {-0x1.4b309d25957e3p-54, 0x1.4f9b2769d2ca7p0},
    {-0x1.f768569bd93efp-55, 0x1.516daa2cf6642p0},
    {-0x1.07abe1db13cadp-55, 0x1.5342b569d4f82p0},
    {-0x1.d689cefede59bp-55, 0x1.551a4ca5d920fp0},
    {0x1.9bb2c011d93adp-54, 0x1.56f4736b527dap0},
    {0x1.295e15b9a1de8p-55, 0x1.58d12d497c7fdp0},
    {0x1.6324c054647adp-54, 0x1.5ab07dd485429p0},
    {0x1.c4b1b816986a2p-60, 0x1.5c9268a5946b7p0},
    {0x1.ba6f93080e65ep-54, 0x1.5e76f15ad2148p0},
    {-0x1.3e2429b56de47p-54, 0x1.605e1b976dc09p0},
    {-0x1.383c17e40b497p-54, 0x1.6247eb03a5585p0},
    {-0x1.c483c759d8933p-55, 0x1.6434634ccc32p0},
    {-0x1.bb60987591c34p-54, 0x1.6623882552225p0},
    {0x1.038ae44f73e65p-57, 0x1.68155d44ca973p0},
    {-0x1.bdd3413b26456p-54, 0x1.6a09e667f3bcdp0},
    {-0x1.2895667ff0b0dp-56, 0x1.6c012750bdabfp0},
    {-0x1.bbe3a683c88abp-57, 0x1.6dfb23c651a2fp0},
    {-0x1.83c0f25860ef6p-55, 0x1.6ff7df9519484p0},
    {-0x1.16e4786887a99p-55, 0x1.71f75e8ec5f74p0},
    {-0x1.0a8d96c65d53cp-54, 0x1.73f9a48a58174p0},
    {-0x1.0245957316dd3p-54, 0x1.75feb564267c9p0},
    {0x1.866b80a02162dp-54, 0x1.780694fde5d3fp0},
    {-0x1.41577ee04992fp-55, 0x1.7a11473eb0187p0},
    {0x1.f124cd1164dd6p-54, 0x1.7c1ed0130c132p0},
    {0x1.05d02ba15797ep-56, 0x1.7e2f336cf4e62p0},
    {-0x1.27c86626d972bp-54, 0x1.80427543e1a12p0},
    {-0x1.d4c1dd41532d8p-54, 0x1.82589994cce13p0},
    {-0x1.8d684a341cdfbp-55, 0x1.8471a4623c7adp0},
    {-0x1.fc6f89bd4f6bap-54, 0x1.868d99b4492edp0},
    {0x1.994c2f37cb53ap-54, 0x1.88ac7d98a6699p0},
    {0x1.6e9f156864b27p-54, 0x1.8ace5422aa0dbp0},
    {-0x1.0d55e32e9e3aap-56, 0x1.8cf3216b5448cp0},
    {0x1.5cc13a2e3976cp-55, 0x1.8f1ae99157736p0},
    {-0x1.dd6792e582524p-54, 0x1.9145b0b91ffc6p0},
    {-0x1.75fc781b57ebcp-57, 0x1.93737b0cdc5e5p0},
    {-0x1.64b7c96a5f039p-56, 0x1.95a44cbc8520fp0},
    {-0x1.d185b7c1b85d1p-54, 0x1.97d829fde4e5p0},
    {-0x1.173bd91cee632p-54, 0x1.9a0f170ca07bap0},
    {0x1.c7c46b071f2bep-56, 0x1.9c49182a3f09p0},
    {0x1.824ca78e64c6ep-56, 0x1.9e86319e32323p0},
    {-0x1.359495d1cd533p-54, 0x1.a0c667b5de565p0},
    {0x1.6305c7ddc36abp-54, 0x1.a309bec4a2d33p0},
    {-0x1.d2f6edb8d41e1p-54, 0x1.a5503b23e255dp0},
    {0x1.bcb7ecac563c7p-54, 0x1.a799e1330b358p0},
    {0x1.0fac90ef7fd31p-54, 0x1.a9e6b5579fdbfp0},
    {-0x1.f9234cae76cdp-55, 0x1.ac36bbfd3f37ap0},
    {0x1.7a1cd345dcc81p-54, 0x1.ae89f995ad3adp0},
    {-0x1.bdef54c80e425p-54, 0x1.b0e07298db666p0},
    {-0x1.2805e3084d708p-57, 0x1.b33a2b84f15fbp0},
    {-0x1.c71dfbbba6de3p-54, 0x1.b59728de5593ap0},
    {-0x1.5584f7e54ac3bp-56, 0x1.b7f76f2fb5e47p0},
    {-0x1.efcd30e54292ep-54, 0x1.ba5b030a1064ap0},
    {0x1.23dd07a2d9e84p-55, 0x1.bcc1e904bc1d2p0},
    {-0x1.efdca3f6b9c73p-54, 0x1.bf2c25bd71e09p0},
    {0x1.11065895048ddp-55, 0x1.c199bdd85529cp0},
    {0x1.b4537e083c60ap-54, 0x1.c40ab5fffd07ap0},
    {0x1.2884dff483cadp-54, 0x1.c67f12e57d14bp0},
    {0x1.1acbc48805c44p-56, 0x1.c8f6d9406e7b5p0},
    {0x1.503cbd1e949dbp-56, 0x1.cb720dcef9069p0},
    {-0x1.dd83b53829d72p-55, 0x1.cdf0b555dc3fap0},
    {-0x1.cbc3743797a9cp-54, 0x1.d072d4a07897cp0},
    {-0x1.d487b719d8578p-54, 0x1.d2f87080d89f2p0},
    {0x1.2ed02d75b3707p-55, 0x1.d5818dcfba487p0},
    {-0x1.11ec18beddfe8p-54, 0x1.d80e316c98398p0},
    {0x1.c2300696db532p-54, 0x1.da9e603db3285p0},
    {0x1.2da5778f018c3p-54, 0x1.dd321f301b46p0},
    {-0x1.1a5cd4f184b5cp-54, 0x1.dfc97337b9b5fp0},
    {-0x1.7b627817a1496p-54, 0x1.e264614f5a129p0},
    {0x1.39e8980a9cc8fp-55, 0x1.e502ee78b3ff6p0},
    {0x1.2d522ca0c8de2p-54, 0x1.e7a51fbc74c83p0},
    {-0x1.e9c23179c2893p-54, 0x1.ea4afa2a490dap0},
    {-0x1.c93f3b411ad8cp-54, 0x1.ecf482d8e67f1p0},
    {0x1.dc7f486a4b6bp-54, 0x1.efa1bee615a27p0},
    {0x1.3a1a5bf0d8e43p-54, 0x1.f252b376bba97p0},
    {0x1.9d3e12dd8a18bp-54, 0x1.f50765b6e454p0},
    {-0x1.dbb12d006350ap-54, 0x1.f7bfdad9cbe14p0},
    {0x1.74853f3a5931ep-55, 0x1.fa7c1819e90d8p0},
    {0x1.2eb74966579e7p-57, 0x1.fd3c22b8f71f1p0},
};

// Cold path for scaling and rounding:
//   2^hi * (r_hi + r_lo) ~ b^x
// when the result may be subnormal (hi <= -1022). Falls back to 128-bit
// evaluation exp_accurate_128(x, log2_base) if Ziv's denormal rounding test
// fails.
//
// Input expectations:
//   - x: original input in the non-overflowing/non-underflowing range of b^x.
//   - hi: integer exponent k >> 7 in [-1075, -1022].
//   - r_hi, lo_u: high part and upper-bounded low part (r_lo + err_d) of the
//     unscaled mantissa:
//       2^(idx / 2^7) * exp(u) in [1 - 2^-8, 2),
//     with |lo_u - err_d| << r_hi.
//   - err_d: Ziv error bound on r_hi + r_lo.
//   - log2_base: pointer to 192-bit log2(b) table (or nullptr when b = 2).
LIBC_NOINLINE LIBC_INLINE double
exp_denorm(double x, int hi, double r_hi, double lo_u, double err_d,
           const ExpBase192 *log2_base = nullptr) {
  using FPBits = fputil::FPBits<double>;
  // Exact subnormal power of 2 when x is an integer in [-1074, -1022].
  if (log2_base == nullptr && LIBC_UNLIKELY(static_cast<double>(hi) == x))
    return FPBits(uint64_t(1) << (hi + 1074)).get_val();

  if (auto r_denorm = ziv_test_denorm(hi, r_hi, lo_u - err_d, err_d);
      LIBC_LIKELY(r_denorm.has_value())) {
    double res = r_denorm.value();
    if (LIBC_UNLIKELY(FPBits(res).is_normal()))
      return res;

    fputil::set_errno_if_required(ERANGE);
    fputil::raise_underflow_except_if_required<double>();
    return res;
  }
  return exp_accurate_128(x, log2_base);
}

// Evaluate:
//   b^x = 2^hi * 2^(idx / 2^7) * exp(u)
// in double-double precision and apply Ziv's rounding test, falling back to
// exp_denorm (when hi <= -1022) or exp_accurate_128 (when Ziv's test fails).
//
// Input expectations:
//   - x: original input after filtering out tiny, overflow, underflow, and
//     Inf/NaN cases.
//   - u: reduced argument:
//       u.hi + u.lo ~ x * log(b) - k * 2^-7 * log(2)
//     with |u.hi + u.lo| <= 2^-8 * log(2). On targets without FMA,
//       lsb(u.hi) >= 2^-30 and |u.lo| < 2^-18
//     so that m_hi * u.hi is exact in double.
//   - hi: integer exponent k >> 7 in [-1075, 1024].
//     When hi == 1024, k = 1024 * 2^7 so
//       idx = 0 (exp2_mid = 1.0) and u < -2^-46 < 0,
//     since x is strictly below the overflow threshold, which keeps upper in
//     [1 - 2^-8, 1) with unbiased exponent -1 (biased exponent 1022). Adding
//     hi = 1024 in scale_exp2(upper, hi) then produces biased exponent 2046
//     (unbiased exponent 1023) without overflowing.
//   - idx: table index k & 0x7f in [0, 127].
//   - err_d: Ziv error bound on the double-double approximation r_hi + r_lo.
//   - log2_base: pointer to 192-bit log2(b) table (or nullptr when b = 2).
LIBC_ALWAYS_INLINE double exp2_eval(double x, DoubleDouble u, int hi,
                                    unsigned idx, double err_d = EXP2_ERR_D,
                                    const ExpBase192 *log2_base = nullptr) {
  DoubleDouble exp2_mid = EXP2_MID[idx];

  // Reconstruction of 2^(idx / 2^7) * exp(u) ~ exp2_mid * (1 + u + u^2 * P(u))
  // where:
  //   exp2_mid = 2^(idx / 2^7) = exp2_mid.hi + exp2_mid.lo
  //   u        = dx * log(b)   = u.hi + u.lo
#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
  // Decompose exp2_mid * (1 + u + u^2 * P(u)) ~ r_hi + (r_lo + tail):
  //   r_hi + r_lo = exp2_mid.hi * (1 + u.hi)
  //   tail        = exp2_mid.hi * (u.lo + u^2 * P(u))
  //                 + exp2_mid.lo * (1 + u.hi)
  // where the omitted term:
  //   exp2_mid.lo * (u.lo + u^2 * P(u))
  // is bounded by 2^(-54 - 18) = 2^-72 and absorbed into err_d.

  // p ~ (exp(u) - 1 - u) / u^2
  double u2 = u.hi * u.hi;
  double p = exp_poly(u.hi, u2);

  // r_hi = round(exp2_mid.hi * (1 + u.hi))
  double r_hi = fputil::multiply_add(exp2_mid.hi, u.hi, exp2_mid.hi);
  // mid_lo ~ exp2_mid.lo * (1 + u.hi)
  double mid_lo = fputil::multiply_add(exp2_mid.lo, u.hi, exp2_mid.lo);
  // poly_lo ~ u.lo + u^2 * P(u)
  double poly_lo = fputil::multiply_add(u2, p, u.lo);

  // Fold +/-err_d into mid_lo so tail_u, tail_l, and r_lo can be evaluated in
  // parallel:
  //   mid_u = mid_lo + err_d
  //   mid_l = mid_lo - err_d
  double mid_u = mid_lo + err_d;
  double mid_l = mid_lo - err_d;
  // tail_{u,l} ~ exp2_mid.hi * (u.lo + u^2 * P(u)) + exp2_mid.lo * (1 + u.hi)
  //              +/- err_d
  double tail_u = fputil::multiply_add(exp2_mid.hi, poly_lo, mid_u);
  double tail_l = fputil::multiply_add(exp2_mid.hi, poly_lo, mid_l);
  // r_lo = exp2_mid.hi * (1 + u.hi) - r_hi (exact to within 2^-105)
  double r_lo = fputil::multiply_add(exp2_mid.hi, u.hi, exp2_mid.hi - r_hi);

  // lo_{u,l} ~ r_lo + tail +/- err_d
  double lo_u = r_lo + tail_u;
  double lo_l = r_lo + tail_l;
#else
  // p ~ (exp(u) - 1 - u) / u^2
  double u_sum = u.hi + u.lo;
  double u2 = u_sum * u_sum;
  double p = exp_poly(u_sum, u2);

  // Split exp2_mid = m_hi + m_lo with lsb(m_hi) >= 2^-22 and |m_lo| < 2^-22:
  LIBC_BIT_CAST_CONSTEXPR_VAR double HI_MASK_22 =
      cpp::bit_cast<double>(0xffff'ffff'c000'0000ULL);
  double m_hi = fputil::fp_and(exp2_mid.hi, HI_MASK_22); // lsb = 2^-22
  double m_lo = (exp2_mid.hi - m_hi) + exp2_mid.lo;

  // Since lsb(u.hi) >= 2^-30 and r_hi is in [1 - 2^-8, 2), both:
  //   m_hi * u.hi  and  r_hi = m_hi * (1 + u.hi)
  // are exact:
  double prod_hi = m_hi * u.hi; // exact, lsb = 2^-52
  double r_hi = m_hi + prod_hi; // exact

  // poly_lo ~ u.lo + u^2 * P(u)
  double u2_p = u2 * p;
  double poly_lo = u.lo + u2_p;

  // r_lo ~ exp2_mid.hi * (u.lo + u^2 * P(u)) + m_lo * (1 + u.hi)
  double l0 = exp2_mid.hi * poly_lo;
  double l1 = m_lo * (1.0 + u.hi);
  double r_lo = l0 + l1;

  // lo_{u,l} = r_lo +/- err_d
  double lo_u = r_lo + err_d;
  double lo_l = r_lo - err_d;
#endif
  // Ziv's rounding test in the current rounding mode:
  double upper = r_hi + lo_u;
  double lower = r_hi + lo_l;

  if (LIBC_LIKELY(upper == lower && hi > -1022))
    return scale_exp2(upper, hi);

  if (LIBC_UNLIKELY(hi <= -1022))
    return exp_denorm(x, hi, r_hi, lo_u, err_d, log2_base);

  return exp_accurate_128(x, log2_base);
}

#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS

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

#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP_UTILS_H