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/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/architectures.h"
#include "src/__support/macros/properties/cpu_features.h"

#ifdef LIBC_TARGET_CPU_HAS_SSE2
#include <immintrin.h>
#elif defined(LIBC_TARGET_ARCH_IS_AARCH64) && defined(__ARM_NEON)
#include <arm_neon.h>
#endif // LIBC_TARGET_CPU_HAS_SSE2

#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) {
  using FPBits = typename fputil::FPBits<double>;

  // Scaling factor = 1/(min normal number) = 2^1022
  int64_t exp_hi =
      static_cast<int64_t>(hi + 1022) * (1LL << FPBits::FRACTION_LEN);
  double mid_hi = cpp::bit_cast<double>(exp_hi + cpp::bit_cast<int64_t>(mid));
  double lo_scaled =
      (lo != 0.0) ? cpp::bit_cast<double>(exp_hi + cpp::bit_cast<int64_t>(lo))
                  : 0.0;

  double extra_factor = 0.0;
  uint64_t scale_down = 0x3FE0'0000'0000'0000; // 1022 in the exponent field.

  // 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 = 0x3FF0'0000'0000'0000; // 1023 in the exponent field.
  }

  // 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 cpp::bit_cast<double>(cpp::bit_cast<uint64_t>(r) - scale_down);
  } else {
    double err_scaled =
        cpp::bit_cast<double>(exp_hi + cpp::bit_cast<int64_t>(err));

    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 cpp::bit_cast<double>(cpp::bit_cast<uint64_t>(upper) - scale_down);
    }

    return cpp::nullopt;
  }
}

namespace math {
namespace double_eval {

using fputil::DoubleDouble;

// 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};

// Directly add an integer `hi` to the exponent field of a double precision `r`.
// It is equivalent to scaling `r` by 2^hi when both `r` and the resulting bits
// are normal.
// The default implementation is simply bit cast `r` to int64`, shift `hi` to
// the exponent field, then add `r` bits and shifted `hi`, before bit cast the
// result back to double.
// The addition is done in integer, and most likely generated as addition in the
// general purpose registers (GPR), causing a full round trip from
// floating-point/SIMD register -> GPR -> FP/SIMD.
// When SSE2 or NEON are available, we can use the compiler intrinsic to keep
// `r` in the FP/SIMD register, and perform the integer addition in the FP/SIMD
// registers instead.  This in general results in a reduction in latency since
// `r` is normally on the critical path.
LIBC_ALWAYS_INLINE double scale_exp2(double r, int hi) {
  using FPBits = typename fputil::FPBits<double>;
#ifdef LIBC_TARGET_CPU_HAS_SSE2
  // Move `hi` to SIMD register and shift to exponent field position.
  __m128i v_exp = _mm_slli_epi64(_mm_cvtsi32_si128(hi), FPBits::FRACTION_LEN);
  // Load `r` and "bit_cast" it to integer SIMD type.
  __m128i v_r = _mm_castpd_si128(_mm_set_sd(r));
  // Perform integer addition in SIMD registers and bit cast the result back to
  // double.
  return _mm_cvtsd_f64(_mm_castsi128_pd(_mm_add_epi64(v_exp, v_r)));
#elif defined(LIBC_TARGET_ARCH_IS_AARCH64) && defined(__ARM_NEON)
  // Shift `hi` to exponent field position and move it to SIMD register.
  int64x1_t v_exp =
      vdup_n_s64(static_cast<int64_t>(hi) << FPBits::FRACTION_LEN);
  // Bit cast `r` to integer SIMD type.
  int64x1_t v_r = vreinterpret_s64_f64(vdup_n_f64(r));
  // Perform integer addition in SIMD registers and bit cast the result back to
  // double.
  return vget_lane_f64(vreinterpret_f64_s64(vadd_s64(v_exp, v_r)), 0);
#else
  int64_t exp_hi = static_cast<int64_t>(hi) * (1LL << FPBits::FRACTION_LEN);
  return cpp::bit_cast<double>(exp_hi + cpp::bit_cast<int64_t>(r));
#endif // LIBC_TARGET_CPU_HAS_SSE2
}

// 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 = typename 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.66.
//
// 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.66.
//
// 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.36.
//
// 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.36,
// reducing the maximum error compared to P0 by (2 + sqrt(3))^3 / 16 ~ 3.25x.
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);
}

#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^-22, 2^-18) and pre-split
// 12-bit + 53-bit logarithm constants for exact non-FMA multiplication without
// Veltkamp splitting:
LIBC_INLINE_VAR constexpr double SPLIT_22 = 0x1.8p30; // 1.5 * 2^(-22 + 52)
LIBC_INLINE_VAR constexpr double SPLIT_18 = 0x1.8p34; // 1.5 * 2^(-18 + 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 bound for double precision fast path (N = 7, Ziv failure rate <
// 0.05%):
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},
};

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 = typename 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);
}

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 cross 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^-23:
  double m_hi = (exp2_mid.hi + SPLIT_22) - SPLIT_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 [0.997, 2.0), 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