libc/src/__support/math/sinhcosh_utils.h - llvm-project.git - Git at Google


            
            //===----------------------------------------------------------------------===//

            
            //

            
            // 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 for double-precision sinh and cosh functions.

            
            ///

            
            //===----------------------------------------------------------------------===//


            #ifndef LLVM_LIBC_SRC___SUPPORT_MATH_SINHCOSH_UTILS_H

            
            #define LLVM_LIBC_SRC___SUPPORT_MATH_SINHCOSH_UTILS_H


            #include "src/__support/CPP/bit.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/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"

            
            #include "src/__support/math/exp_utils.h"


            #ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS

            
            #include "src/__support/math/sinhcosh_accurate_128.h"

            
            #endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS


            namespace LIBC_NAMESPACE_DECL {

            
            namespace math {

            
            namespace double_eval {


            #ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS


            // Splitting -log(2) ~ hi + lo with 35-bit-precision hi part so that

            
            // kd * MLOG_2_HI is exact for |kd| < 2^18:

            
            // > hi = round(-log(2), 35, RN);

            
            // > lo = round(-log(2) - hi, D, RN);

            
            LIBC_INLINE_VAR constexpr double MLOG_2_HI = -0x1.62e42fefcp-1;

            
            LIBC_INLINE_VAR constexpr double MLOG_2_LO = 0x1.c610ca86c3899p-37;


            // Degree-4 minimax polynomials in u^2 for (cosh(u) - 1) / u^2 and

            
            // (sinh(u) - u) / u^3 on [-log(2) / 2, log(2) / 2]:

            
            //   cosh(u) - 1 ~ u^2 * P_even(u^2)

            
            //   sinh(u) - u ~ u^3 * P_odd(u^2).

            
            //

            
            // Generated by Sollya with:

            
            // > display = hexadecimal!;

            
            // > prec = 500;

            
            // > B = (log(2) / 2)^2;

            
            // > P_even = fpminimax((cosh(sqrt(x)) - 1) / x, 4, [|D...|], [2^(-100), B]);

            
            // > for i from 0 to 4 do print(coeff(P_even, i), ",");

            
            // > dirtyinfnorm((cosh(sqrt(x)) - 1) - x * P_even, [0, B]);

            
            //   0x1.c244...p-57 < 2^-56.

            
            LIBC_INLINE_VAR constexpr double EXP_EVEN_N0_COEFFS[5] = {

            
                0x1.0000000000001p-1,  // ~ 1/2!

            
                0x1.5555555553ce4p-5,  // ~ 1/4!

            
                0x1.6c16c178a71cep-10, // ~ 1/6!

            
                0x1.a019b8c6ec68dp-16, // ~ 1/8!

            
                0x1.2891c21cf0c03p-22, // ~ 1/10!

            
            };


            // > P_odd = fpminimax((sinh(sqrt(x)) / sqrt(x) - 1) / x, 4, [|D...|],

            
            //                     [2^(-100), B]);

            
            // > for i from 0 to 4 do print(coeff(P_odd, i), ",");

            
            // > dirtyinfnorm((sinh(sqrt(x)) - sqrt(x)) - x * sqrt(x) * P_odd, [0, B]);

            
            //   0x1.4882...p-62 < 2^-61.

            
            LIBC_INLINE_VAR constexpr double EXP_ODD_N0_COEFFS[5] = {

            
                0x1.5555555555556p-3,  // ~ 1/3!

            
                0x1.111111111070cp-7,  // ~ 1/5!

            
                0x1.a01a01a901652p-13, // ~ 1/7!

            
                0x1.71de0a3fd7069p-19, // ~ 1/9!

            
                0x1.af3effc3463e9p-26, // ~ 1/11!

            
            };


            // Tableless (N = 0) evaluation of:

            
            //   sinh(x) = (e^x - e^-x) / 2   (when is_sinh is true)

            
            //   cosh(x) = (e^x + e^-x) / 2   (when is_sinh is false, with x >= 0)

            
            // for 2^-26 < |x| < 0x1.633ce8fb9f87ep+9 (sinh) or

            
            // 2^-27 < x < 0x1.633ce8fb9f87ep+9 (cosh).

            
            //

            
            // Range reduction:

            
            // Let |x| = k * log(2) + u, where:

            
            //   k   = round(|x| * log2(e))

            
            //   |u| <= log(2) / 2

            
            // then:

            
            //   e^|x|  = 2^k * e^u

            
            //   e^-|x| = 2^-k * e^-u = 2^k * 2^(-2 * k) * e^-u.

            
            //

            
            // Factoring out 2^(k - 1) and separating e^(+/- u) = cosh(u) +/- sinh(u):

            
            //   (e^|x| -/+ e^-|x|) / 2 = 2^(k - 1) * (e^u + m_m * e^-u)

            
            //                          = 2^(k - 1) * (s * cosh(u) + d * sinh(u))

            
            // where:

            
            //   m_m = (-/+ 1) * 2^(-2 * k)

            
            //   s   = 1 + m_m

            
            //   d   = 1 - m_m.

            
            //

            
            // Key properties across the domain:

            
            // - For |x| < 0x1.633ce8fb9f87ep+9, k is in [0, 1025] so k - 1 is in

            
            //   [-1, 1024]. When k = 1025, u < 0 keeps:

            
            //     r = e^u + m_m * e^-u < 1.0,

            
            //   so 2^1024 * r stays in the normal range and scale_exp2(r, k - 1) handles

            
            //   the entire domain without a separate top-octave branch.

            
            //

            
            // - For |x| < log(2) / 2, k = 0 gives u = |x| (exact, t2 = 0) and m_m = -/+ 1:

            
            //   * For sinh (m_m = -1): s = 0 and d = 2, reducing the evaluation to:

            
            //       2^-1 * 2 * (u + u^3 * P_odd(u^2)) = x + x^3 * P_odd(x^2)

            
            //     with no cancellation down to |x| = 2^-26.

            
            //   * For cosh (m_m = +1): s = 2 and d = 0, reducing to:

            
            //       1 + u^2 * P_even(u^2).

            
            //

            
            // - For 1 <= k <= 26, 2^(-2 * k) is in [2^-52, 2^-2], so both:

            
            //     s = 1.0 + m_m

            
            //     d = 1.0 - m_m

            
            //   are exact in all rounding modes. For k >= 27:

            
            //     |m_m| <= 2^-54 < 0.5 * ulp(1.0),

            
            //   so clamping m_m to 0.0 avoids a 1-ULP rounding error in directed rounding

            
            //   modes.

            
            LIBC_MATH_SMALL_TABLES_INLINE double exp_pm_eval(double x, bool is_sinh) {

            
              using FPBits = fputil::FPBits<double>;

            
              FPBits xbits(x);

            
              double x_a = xbits.abs().get_val();


              // Using ROUND_BIAS, ROUND_SHIFT, and ROUND_MASK (see exp_utils.h), we compute

            
              // k_biased = fma(x_a, LOG2_E.hi, ROUND_BIAS) and extract in parallel:

            
              //   kd = (k_biased & ROUND_MASK) - ROUND_SHIFT

            
              //   k  = (k_bits >> 32) & 0x3'ffff.

            
              double k_biased = fputil::multiply_add(x_a, LOG2_E.hi, ROUND_BIAS);

            
              double kd = fputil::fp_and(k_biased, ROUND_MASK) - ROUND_SHIFT;

            
              int k = static_cast<int>(cpp::bit_cast<uint64_t>(k_biased) >> 32) & 0x3'ffff;


              // Secondary scale factor:

            
              //   m_m = (-/+ 1) * 2^(-2 * k)   for k < 27

            
              //   m_m = 0.0                    for k >= 27 (where 2^(-2 * k) <= 2^-54).

            
              //

            
              // Subtracting (1 << 11) = 0x800 (for sinh) sets bit 11 of the 12-bit

            
              // two's-complement exponent field, flipping the sign bit of m_m when shifted

            
              // left by 52 in scale_exp2:

            
              int sign_shift = static_cast<int>(is_sinh) << 11;

            
              int exp_diff = -2 * k - sign_shift;

            
              double m_m = (k < 27) ? scale_exp2(1.0, exp_diff) : 0.0;


              // Reduced argument:

            
              //   u = (x_a + kd * MLOG_2_HI) + kd * MLOG_2_LO.

            
              // Since MLOG_2_HI has 35 significant bits and |kd| <= 1025,

            
              //   t1 = fma(kd, MLOG_2_HI, x_a)

            
              // is exact in all rounding modes with or without FMA.

            
              double t1 = fputil::multiply_add(kd, MLOG_2_HI, x_a); // exact

            
              double t2 = kd * MLOG_2_LO;

            
              double u = t1 + t2;

            
              double u2 = u * u;

            
              double u4 = u2 * u2;


              // Evaluate degree-4 even and odd minimax polynomials P_even(u^2) and

            
              // P_odd(u^2) in parallel using Estrin's scheme:

            
              //   cosh(u) ~ 1 + u^2 * P_even(u^2)

            
              //   sinh(u) ~ t1 + (t2 + u^3 * P_odd(u^2))


              // First Estrin layer: c_0 + u^2 * c_1 and c_2 + u^2 * c_3

            
              double e01 =

            
                  fputil::multiply_add(u2, EXP_EVEN_N0_COEFFS[1], EXP_EVEN_N0_COEFFS[0]);

            
              double o01 =

            
                  fputil::multiply_add(u2, EXP_ODD_N0_COEFFS[1], EXP_ODD_N0_COEFFS[0]);

            
              double e23 =

            
                  fputil::multiply_add(u2, EXP_EVEN_N0_COEFFS[3], EXP_EVEN_N0_COEFFS[2]);

            
              double o23 =

            
                  fputil::multiply_add(u2, EXP_ODD_N0_COEFFS[3], EXP_ODD_N0_COEFFS[2]);


              // Second Estrin layer: (c_2 + u^2 * c_3) + u^4 * c_4

            
              double e24 = fputil::multiply_add(u4, EXP_EVEN_N0_COEFFS[4], e23);

            
              double o24 = fputil::multiply_add(u4, EXP_ODD_N0_COEFFS[4], o23);


              // Final polynomial values P_even(u^2) and P_odd(u^2):

            
              double p_even = fputil::multiply_add(u4, e24, e01);

            
              double p_odd = fputil::multiply_add(u4, o24, o01);


              // Since m_m is either 0 or +/- 2^(-2 * k) with 0 <= k <= 26, both s and d are

            
              // exact in all rounding modes:

            
              double s = 1.0 + m_m; // exact

            
              double d = 1.0 - m_m; // exact


              // Apply sign(x) to the exact weights s and d before accumulating the

            
              // polynomial so every subsequent rounding operation rounds in the active

            
              // rounding direction for negative x as well (for cosh, the caller passes

            
              // x >= 0 so sign_x == +0.0):

            
              double sign_x = fputil::fp_xor(x, x_a);

            
              s = fputil::fp_xor(s, sign_x);

            
              d = fputil::fp_xor(d, sign_x);


              // Accumulate from smallest to largest magnitude before adding the leading

            
              // constant s. Adding s last avoids a double rounding at the leading ULP

            
              // position of s when |d * t1| (< 2^-1) and |s * u^2 * P_even(u^2)| (< 2^-4)

            
              // are added:

            
              double d_t2 = d * t2;

            
              double s_u2 = s * u2;

            
              double d_u = d * u;

            
              double d_u3 = d_u * u2;


              // even_lo ~ s * u^2 * P_even(u^2) + d * t2

            
              double even_lo = fputil::multiply_add(s_u2, p_even, d_t2);

            
              // poly_lo ~ d * u^3 * P_odd(u^2) + even_lo

            
              double poly_lo = fputil::multiply_add(d_u3, p_odd, even_lo);

            
              // poly ~ d * t1 + poly_lo

            
              double poly = fputil::multiply_add(d, t1, poly_lo);


              // r ~ s * cosh(u) + d * sinh(u)

            
              double r = s + poly;

            
              return scale_exp2(r, k - 1);

            
            }


            #else // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS


            // Fixed-point splitting constants (ulp = 2^-22, 2^-29) for exact multiplication

            
            // without Veltkamp splitting:

            
            LIBC_INLINE_VAR constexpr double SPLIT_22 = 0x1.8p30; // 1.5 * 2^(-22 + 52)

            
            LIBC_INLINE_VAR constexpr double SPLIT_29 = 0x1.8p23; // 1.5 * 2^(-29 + 52)


            // Evaluate:

            
            //   sinh(x) = (e^x - e^-x) / 2   (when is_sinh is true)

            
            //   cosh(x) = (e^x + e^-x) / 2   (when is_sinh is false, with x >= 0)

            
            // in double-double precision for 2^-2 <= |x| < 0x1.633ce8fb9f87ep+9 (sinh) or

            
            // 2^-27 < x < 0x1.633ce8fb9f87ep+9 (cosh), falling back to

            
            // exp_pm_accurate_128(x, is_sinh) when Ziv's rounding test fails.

            
            //

            
            // Range reduction:

            
            // Let |x| = k_p * (log(2) / 2^7) + u

            
            // in which:

            
            //   k_p   = round(|x| * 2^7 * log2(e)),  k_m = -k_p

            
            //   hi_p  = k_p >> 7,                    idx_p = k_p & 0x7f

            
            //   hi_m  = k_m >> 7,                    idx_m = k_m & 0x7f

            
            //   |u|  <= 2^-8 * log(2)

            
            // then:

            
            //   e^|x|  = 2^hi_p * 2^(idx_p / 2^7) * e^u

            
            //   e^-|x| = 2^hi_m * 2^(idx_m / 2^7) * e^-u

            
            //          = 2^hi_p * (2^(hi_m - hi_p) * 2^(idx_m / 2^7)) * e^-u.

            
            //

            
            // Factoring out 2^(hi_p - 1):

            
            //   (e^|x| -/+ e^-|x|) / 2 = 2^(hi_p - 1) * (m_p * e^u + m_m * e^-u)

            
            // where:

            
            //   m_p = 2^(idx_p / 2^7)

            
            //   m_m = (-/+ 1) * 2^(hi_m - hi_p) * 2^(idx_m / 2^7).

            
            //

            
            // For |x| < 0x1.633ce8fb9f87ep+9, hi_p is in [0, 1025] so hi_p - 1 is in

            
            // [-1, 1024]. When hi_p == 1025, idx_p == 0 and u < 0, which keeps:

            
            //   r = m_p * e^u + m_m * e^-u < 1.0

            
            // (biased exponent 1022). Adding hi_p - 1 = 1024 in scale_exp2(r, hi_p - 1)

            
            // produces biased exponent 2046, handling the top octave

            
            // [1024 * log(2), 1025 * log(2)) without a separate branch.

            
            LIBC_MATH_SMALL_TABLES_INLINE double exp_pm_eval(double x, bool is_sinh) {

            
              using FPBits = fputil::FPBits<double>;

            
              FPBits xbits(x);

            
              double x_a = xbits.abs().get_val();


              constexpr double LOG2_E_SCALE = EXP_SCALE * LOG2_E.hi;

            
              // Using ROUND_BIAS, ROUND_SHIFT, and ROUND_MASK (see exp_utils.h), we compute

            
              // k_biased = fma(x_a, LOG2_E_SCALE, ROUND_BIAS) and extract in parallel:

            
              //   kd  = (k_biased & ROUND_MASK) - ROUND_SHIFT

            
              //   k_p = (k_bits >> 32) & 0x3'ffff.

            
              double k_biased = fputil::multiply_add(x_a, LOG2_E_SCALE, ROUND_BIAS);

            
              double kd = fputil::fp_and(k_biased, ROUND_MASK) - ROUND_SHIFT;

            
              int k_p =

            
                  static_cast<int>(cpp::bit_cast<uint64_t>(k_biased) >> 32) & 0x3'ffff;

            
              int k_m = -k_p;


              int hi_p = k_p >> EXP_N_BITS;

            
              unsigned idx_p = static_cast<unsigned>(k_p & EXP_MASK);


              int hi_m = k_m >> EXP_N_BITS;

            
              unsigned idx_m = static_cast<unsigned>(k_m & EXP_MASK);


              // Secondary exponent offset:

            
              //   exp_diff = hi_m - hi_p

            
              // clamped at -100 (where e^-|x| <= 2^-99 * e^|x| << SINH_ERR_D) to avoid

            
              // underflow in m_m:

            
              int exp_diff = hi_m - hi_p;

            
              int clamped_diff = (exp_diff > -100) ? exp_diff : -100;


              // Subtracting (1 << 11) = 0x800 (for sinh) sets bit 11 of the 12-bit

            
              // two's-complement exponent field, flipping the sign bit of m_m when shifted

            
              // left by 52 in scale_exp2:

            
              int sign_shift = static_cast<int>(is_sinh) << 11;

            
              int scale_m = clamped_diff - sign_shift;


              DoubleDouble m_p = EXP2_MID[idx_p];

            
              // Scale m_m = EXP2_MID[idx_m] by (-/+ 1) * 2^clamped_diff:

            
              // - Since m_m.hi is in [1, 2), we use scale_exp2(m_m.hi, scale_m) directly.

            
              // - Since EXP2_MID[0].lo == 0.0, we scale m_m.lo by multiplying by:

            
              //     s_m = scale_exp2(1.0, scale_m).

            
              double s_m = scale_exp2(1.0, scale_m);

            
              DoubleDouble m_m = EXP2_MID[idx_m];

            
              m_m.lo *= s_m;

            
              m_m.hi = scale_exp2(m_m.hi, scale_m);


              // Reduced argument:

            
              //   u = (x_a + kd * MLOG_2_EXP2_HI) + kd * MLOG_2_EXP2_MID.

            
              // Since MLOG_2_EXP2_HI has 35 significant bits and |kd| < 2^18,

            
              //   t1 = fma(kd, MLOG_2_EXP2_HI, x_a)

            
              // is exact in all rounding modes with or without FMA.

            
              double t1 = fputil::multiply_add(kd, MLOG_2_EXP2_HI, x_a); // exact

            
              double u = fputil::multiply_add(kd, MLOG_2_EXP2_MID, t1);

            
              double u2 = u * u;


              // Separate e^(+/-u) into even and odd parts:

            
              //   e^u  ~ (1 + u^2 * P_even(u^2)) + (u + u^3 * P_odd(u^2))

            
              //   e^-u ~ (1 + u^2 * P_even(u^2)) - (u + u^3 * P_odd(u^2))

            
              // where:

            
              //   P_even(u^2) = p_0 + p_2 * u^2

            
              //   P_odd(u^2)  = p_1 + p_3 * u^2

            
              // Then:

            
              //   m_p * e^u + m_m * e^-u

            
              //     ~ (m_p + m_m) * (1 + u^2 * P_even(u^2))

            
              //       + (m_p - m_m) * (u + u^3 * P_odd(u^2)).

            
              double p_even = fputil::multiply_add(u2, EXP_P_COEFFS[2], EXP_P_COEFFS[0]);

            
              double p_odd = fputil::multiply_add(u2, EXP_P_COEFFS[3], EXP_P_COEFFS[1]);


            #ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE

            
              // Decompose m_p * e^u + m_m * e^-u ~ r_hi + r_lo:

            
              //   s.hi + s.lo  = m_p.hi + m_m.hi   (exact FastTwoSum)

            
              //   d.hi + d.lo  = m_p.hi - m_m.hi   (exact FastTwoSum)

            
              //   r_hi + r_err = s.hi + d.hi * t1

            
              //   r_lo         ~ r_err + (s.lo + (m_p.lo - m_m.lo + d.lo) * t1)

            
              //                  + ((m_p.lo + m_m.lo) + s.hi * u^2 * P_even(u^2))

            
              //                  + d.hi * (kd * MLOG_2_EXP2_MID + u^3 * P_odd(u^2)).

            
              //

            
              // For |x| >= 2^-2, k_p >= 46 so hi_p >= 0 and hi_m <= -1, hence

            
              // m_p.hi >= 1.0 > |m_m.hi|. By Sterbenz's lemma, both (s.hi - m_p.hi) and

            
              // (m_p.hi - d.hi) are exact in all rounding modes.

            
              DoubleDouble s = fputil::exact_add(m_p.hi, m_m.hi);

            
              DoubleDouble d = fputil::exact_sub(m_p.hi, m_m.hi);


              // Leading term r_hi = round(s.hi + d.hi * t1) and its rounding error r_err:

            
              // Since |d.hi * t1| <= 2^-7 * log(2) < 2^-7 and s.hi > 2^-2, (s.hi - r_hi) is

            
              // exact by Sterbenz's lemma.

            
              double r_hi = fputil::multiply_add(d.hi, t1, s.hi);

            
              double r_t = s.hi - r_hi; // exact

            
              double r_err = fputil::multiply_add(d.hi, t1, r_t);


              // Odd polynomial and low reduced-argument correction:

            
              //   odd_lo ~ (u - t1) + u^3 * P_odd(u^2) = kd * MLOG_2_EXP2_MID + u^3 * p_odd

            
              double t2 = kd * MLOG_2_EXP2_MID;

            
              double u3 = u * u2;

            
              double odd_lo = fputil::multiply_add(u3, p_odd, t2);


              // Combine table low parts and polynomial terms in parallel with the

            
              // FastTwoSum residuals (s.lo, d.lo) and r_hi rounding error (r_err):

            
              //   mp_lo   ~ m_p.lo * (1 + t1)

            
              //   mm_lo   ~ m_m.lo * (1 - t1)

            
              //   even_lo ~ mp_lo + (s.hi * u^2) * P_even(u^2)

            
              //   poly_lo ~ even_lo + d.hi * odd_lo

            
              double mp_lo = fputil::multiply_add(m_p.lo, t1, m_p.lo);

            
              double mm_lo = fputil::multiply_add(-m_m.lo, t1, m_m.lo);

            
              double s_u2 = s.hi * u2;

            
              double even_lo = fputil::multiply_add(s_u2, p_even, mp_lo);

            
              double poly_lo = fputil::multiply_add(d.hi, odd_lo, even_lo);


              // Total linear error from FastTwoSum residuals and r_hi rounding:

            
              //   sd_lo   ~ s.lo + d.lo * t1

            
              //   err_tot ~ r_err + sd_lo

            
              double sd_lo = fputil::multiply_add(d.lo, t1, s.lo);

            
              double err_tot = r_err + sd_lo;


              // Apply sign(x) to r_hi, mm_lo, poly_lo, and err_tot before the final Ziv

            
              // test additions so that (r_hi + lo_u) and (r_hi + lo_l) round in the active

            
              // rounding direction for negative inputs (for cosh, the caller passes x >= 0

            
              // so sign_x == +0.0):

            
              double sign_x = fputil::fp_xor(x, x_a);

            
              r_hi = fputil::fp_xor(r_hi, sign_x);

            
              mm_lo = fputil::fp_xor(mm_lo, sign_x);

            
              poly_lo = fputil::fp_xor(poly_lo, sign_x);

            
              err_tot = fputil::fp_xor(err_tot, sign_x);


              // Fold +/-SINH_ERR_D into mm_lo off the critical path so tail_{u,l} and

            
              // err_tot finish in parallel:

            
              constexpr double SINH_ERR_D = 0x1.4p-64;

            
              double mm_lo_u = mm_lo + SINH_ERR_D;

            
              double mm_lo_l = mm_lo - SINH_ERR_D;

            
              double tail_u = poly_lo + mm_lo_u;

            
              double tail_l = poly_lo + mm_lo_l;


              double lo_u = err_tot + tail_u;

            
              double lo_l = err_tot + tail_l;

            
            #else  // LIBC_TARGET_CPU_HAS_FMA_DOUBLE

            
              // Fixed-point grid splitting for exact leading reconstruction in all rounding

            
              // modes without FMA:

            
              //   lsb(mp_hi), lsb(mm_hi) >= 2^-22  =>  lsb(s_hi), lsb(d_hi) >= 2^-22

            
              //   lsb(u_hi)              >= 2^-29

            
              // Since s_hi = mp_hi + mm_hi > 2^-2, both s_hi, d_hi, and

            
              //   r_hi = s_hi + d_hi * u_hi  (with lsb >= 2^-51)

            
              // are exact in double.


              // Split u = u_hi + u_lo with lsb(u_hi) >= 2^-29 and |u_lo| <= 2^-30:

            
              double u_hi_t = t1 + SPLIT_29;

            
              double u_hi = u_hi_t - SPLIT_29; // exact, lsb = 2^-29

            
              double u_lo_t = t1 - u_hi;       // exact

            
              double u_lo = fputil::multiply_add(kd, MLOG_2_EXP2_MID, u_lo_t);


              // Split m_p = mp_hi + mp_lo and m_m = mm_hi + mm_lo with lsb >= 2^-22:

            
              double mp_hi_t = m_p.hi + SPLIT_22;

            
              double mp_hi = mp_hi_t - SPLIT_22; // exact, lsb = 2^-22

            
              double mp_lo_t = m_p.hi - mp_hi;   // exact

            
              double mp_lo = mp_lo_t + m_p.lo;


              double mm_hi_t = m_m.hi + SPLIT_22;

            
              double mm_hi = mm_hi_t - SPLIT_22; // exact, lsb = 2^-22

            
              double mm_lo_t = m_m.hi - mm_hi;   // exact

            
              double mm_lo = mm_lo_t + m_m.lo;


              // Exact leading sum, difference, and linear term r_hi = s_hi + d_hi * u_hi:

            
              double s_hi = mp_hi + mm_hi;                          // exact, lsb = 2^-22

            
              double d_hi = mp_hi - mm_hi;                          // exact, lsb = 2^-22

            
              double r_hi = fputil::multiply_add(d_hi, u_hi, s_hi); // exact, lsb = 2^-51


              // Low parts of the split table values and their linear contribution:

            
              //   s_lo = mp_lo + mm_lo,  d_lo = mp_lo - mm_lo

            
              //   m_lo ~ s_lo + d_lo * u_hi

            
              double s_lo = mp_lo + mm_lo;

            
              double d_lo = mp_lo - mm_lo;

            
              double m_lo = fputil::multiply_add(d_lo, u_hi, s_lo);


              // Full-precision leading weights for the quadratic and cubic terms:

            
              double s_m_hi = m_p.hi + m_m.hi;

            
              double d_m_hi = m_p.hi - m_m.hi;


              // Apply sign(x) before the final polynomial accumulation so that (r_hi +

            
              // lo_u) and (r_hi + lo_l) round in the active rounding direction for negative

            
              // x:

            
              double sign_x = fputil::fp_xor(x, x_a);

            
              r_hi = fputil::fp_xor(r_hi, sign_x);

            
              s_m_hi = fputil::fp_xor(s_m_hi, sign_x);

            
              d_m_hi = fputil::fp_xor(d_m_hi, sign_x);

            
              m_lo = fputil::fp_xor(m_lo, sign_x);


              // Combine even and odd polynomial terms with m_lo:

            
              //   even_lo ~ m_lo + s_m_hi * u^2 * P_even(u^2)

            
              //   odd_lo  ~ u_lo + u^3 * P_odd(u^2)

            
              //   r_lo    ~ even_lo + d_m_hi * odd_lo

            
              double u2_peven = u2 * p_even;

            
              double even_lo = fputil::multiply_add(s_m_hi, u2_peven, m_lo);

            
              double u3 = u * u2;

            
              double odd_lo = fputil::multiply_add(u3, p_odd, u_lo);

            
              double r_lo = fputil::multiply_add(d_m_hi, odd_lo, even_lo);


              // Total approximation and rounding error is bounded by SINH_ERR_D =

            
              // 0x1.4p-64:

            
              constexpr double SINH_ERR_D = 0x1.4p-64;

            
              double lo_u = r_lo + SINH_ERR_D;

            
              double lo_l = r_lo - SINH_ERR_D;

            
            #endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE


              // 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))

            
                return scale_exp2(upper, hi_p - 1);


              return exp_pm_accurate_128(x, is_sinh);

            
            }


            #endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS


            } // namespace double_eval

            
            } // namespace math

            
            } // namespace LIBC_NAMESPACE_DECL


            #endif // LLVM_LIBC_SRC___SUPPORT_MATH_SINHCOSH_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 for double-precision sinh and cosh functions.
	///
	//===----------------------------------------------------------------------===//

	#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_SINHCOSH_UTILS_H
	#define LLVM_LIBC_SRC___SUPPORT_MATH_SINHCOSH_UTILS_H

	#include "src/__support/CPP/bit.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/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"
	#include "src/__support/math/exp_utils.h"

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

	namespace LIBC_NAMESPACE_DECL {
	namespace math {
	namespace double_eval {

	#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS

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

	// Degree-4 minimax polynomials in u^2 for (cosh(u) - 1) / u^2 and
	// (sinh(u) - u) / u^3 on [-log(2) / 2, log(2) / 2]:
	// cosh(u) - 1 ~ u^2 * P_even(u^2)
	// sinh(u) - u ~ u^3 * P_odd(u^2).
	//
	// Generated by Sollya with:
	// > display = hexadecimal!;
	// > prec = 500;
	// > B = (log(2) / 2)^2;
	// > P_even = fpminimax((cosh(sqrt(x)) - 1) / x, 4, [\|D...\|], [2^(-100), B]);
	// > for i from 0 to 4 do print(coeff(P_even, i), ",");
	// > dirtyinfnorm((cosh(sqrt(x)) - 1) - x * P_even, [0, B]);
	// 0x1.c244...p-57 < 2^-56.
	LIBC_INLINE_VAR constexpr double EXP_EVEN_N0_COEFFS[5] = {
	0x1.0000000000001p-1, // ~ 1/2!
	0x1.5555555553ce4p-5, // ~ 1/4!
	0x1.6c16c178a71cep-10, // ~ 1/6!
	0x1.a019b8c6ec68dp-16, // ~ 1/8!
	0x1.2891c21cf0c03p-22, // ~ 1/10!
	};

	// > P_odd = fpminimax((sinh(sqrt(x)) / sqrt(x) - 1) / x, 4, [\|D...\|],
	// [2^(-100), B]);
	// > for i from 0 to 4 do print(coeff(P_odd, i), ",");
	// > dirtyinfnorm((sinh(sqrt(x)) - sqrt(x)) - x * sqrt(x) * P_odd, [0, B]);
	// 0x1.4882...p-62 < 2^-61.
	LIBC_INLINE_VAR constexpr double EXP_ODD_N0_COEFFS[5] = {
	0x1.5555555555556p-3, // ~ 1/3!
	0x1.111111111070cp-7, // ~ 1/5!
	0x1.a01a01a901652p-13, // ~ 1/7!
	0x1.71de0a3fd7069p-19, // ~ 1/9!
	0x1.af3effc3463e9p-26, // ~ 1/11!
	};

	// Tableless (N = 0) evaluation of:
	// sinh(x) = (e^x - e^-x) / 2 (when is_sinh is true)
	// cosh(x) = (e^x + e^-x) / 2 (when is_sinh is false, with x >= 0)
	// for 2^-26 < \|x\| < 0x1.633ce8fb9f87ep+9 (sinh) or
	// 2^-27 < x < 0x1.633ce8fb9f87ep+9 (cosh).
	//
	// Range reduction:
	// Let \|x\| = k * log(2) + u, where:
	// k = round(\|x\| * log2(e))
	// \|u\| <= log(2) / 2
	// then:
	// e^\|x\| = 2^k * e^u
	// e^-\|x\| = 2^-k * e^-u = 2^k * 2^(-2 * k) * e^-u.
	//
	// Factoring out 2^(k - 1) and separating e^(+/- u) = cosh(u) +/- sinh(u):
	// (e^\|x\| -/+ e^-\|x\|) / 2 = 2^(k - 1) * (e^u + m_m * e^-u)
	// = 2^(k - 1) * (s * cosh(u) + d * sinh(u))
	// where:
	// m_m = (-/+ 1) * 2^(-2 * k)
	// s = 1 + m_m
	// d = 1 - m_m.
	//
	// Key properties across the domain:
	// - For \|x\| < 0x1.633ce8fb9f87ep+9, k is in [0, 1025] so k - 1 is in
	// [-1, 1024]. When k = 1025, u < 0 keeps:
	// r = e^u + m_m * e^-u < 1.0,
	// so 2^1024 * r stays in the normal range and scale_exp2(r, k - 1) handles
	// the entire domain without a separate top-octave branch.
	//
	// - For \|x\| < log(2) / 2, k = 0 gives u = \|x\| (exact, t2 = 0) and m_m = -/+ 1:
	// * For sinh (m_m = -1): s = 0 and d = 2, reducing the evaluation to:
	// 2^-1 * 2 * (u + u^3 * P_odd(u^2)) = x + x^3 * P_odd(x^2)
	// with no cancellation down to \|x\| = 2^-26.
	// * For cosh (m_m = +1): s = 2 and d = 0, reducing to:
	// 1 + u^2 * P_even(u^2).
	//
	// - For 1 <= k <= 26, 2^(-2 * k) is in [2^-52, 2^-2], so both:
	// s = 1.0 + m_m
	// d = 1.0 - m_m
	// are exact in all rounding modes. For k >= 27:
	// \|m_m\| <= 2^-54 < 0.5 * ulp(1.0),
	// so clamping m_m to 0.0 avoids a 1-ULP rounding error in directed rounding
	// modes.
	LIBC_MATH_SMALL_TABLES_INLINE double exp_pm_eval(double x, bool is_sinh) {
	using FPBits = fputil::FPBits<double>;
	FPBits xbits(x);
	double x_a = xbits.abs().get_val();

	// Using ROUND_BIAS, ROUND_SHIFT, and ROUND_MASK (see exp_utils.h), we compute
	// k_biased = fma(x_a, LOG2_E.hi, ROUND_BIAS) and extract in parallel:
	// kd = (k_biased & ROUND_MASK) - ROUND_SHIFT
	// k = (k_bits >> 32) & 0x3'ffff.
	double k_biased = fputil::multiply_add(x_a, LOG2_E.hi, ROUND_BIAS);
	double kd = fputil::fp_and(k_biased, ROUND_MASK) - ROUND_SHIFT;
	int k = static_cast<int>(cpp::bit_cast<uint64_t>(k_biased) >> 32) & 0x3'ffff;

	// Secondary scale factor:
	// m_m = (-/+ 1) * 2^(-2 * k) for k < 27
	// m_m = 0.0 for k >= 27 (where 2^(-2 * k) <= 2^-54).
	//
	// Subtracting (1 << 11) = 0x800 (for sinh) sets bit 11 of the 12-bit
	// two's-complement exponent field, flipping the sign bit of m_m when shifted
	// left by 52 in scale_exp2:
	int sign_shift = static_cast<int>(is_sinh) << 11;
	int exp_diff = -2 * k - sign_shift;
	double m_m = (k < 27) ? scale_exp2(1.0, exp_diff) : 0.0;

	// Reduced argument:
	// u = (x_a + kd * MLOG_2_HI) + kd * MLOG_2_LO.
	// Since MLOG_2_HI has 35 significant bits and \|kd\| <= 1025,
	// t1 = fma(kd, MLOG_2_HI, x_a)
	// is exact in all rounding modes with or without FMA.
	double t1 = fputil::multiply_add(kd, MLOG_2_HI, x_a); // exact
	double t2 = kd * MLOG_2_LO;
	double u = t1 + t2;
	double u2 = u * u;
	double u4 = u2 * u2;

	// Evaluate degree-4 even and odd minimax polynomials P_even(u^2) and
	// P_odd(u^2) in parallel using Estrin's scheme:
	// cosh(u) ~ 1 + u^2 * P_even(u^2)
	// sinh(u) ~ t1 + (t2 + u^3 * P_odd(u^2))

	// First Estrin layer: c_0 + u^2 * c_1 and c_2 + u^2 * c_3
	double e01 =
	fputil::multiply_add(u2, EXP_EVEN_N0_COEFFS[1], EXP_EVEN_N0_COEFFS[0]);
	double o01 =
	fputil::multiply_add(u2, EXP_ODD_N0_COEFFS[1], EXP_ODD_N0_COEFFS[0]);
	double e23 =
	fputil::multiply_add(u2, EXP_EVEN_N0_COEFFS[3], EXP_EVEN_N0_COEFFS[2]);
	double o23 =
	fputil::multiply_add(u2, EXP_ODD_N0_COEFFS[3], EXP_ODD_N0_COEFFS[2]);

	// Second Estrin layer: (c_2 + u^2 * c_3) + u^4 * c_4
	double e24 = fputil::multiply_add(u4, EXP_EVEN_N0_COEFFS[4], e23);
	double o24 = fputil::multiply_add(u4, EXP_ODD_N0_COEFFS[4], o23);

	// Final polynomial values P_even(u^2) and P_odd(u^2):
	double p_even = fputil::multiply_add(u4, e24, e01);
	double p_odd = fputil::multiply_add(u4, o24, o01);

	// Since m_m is either 0 or +/- 2^(-2 * k) with 0 <= k <= 26, both s and d are
	// exact in all rounding modes:
	double s = 1.0 + m_m; // exact
	double d = 1.0 - m_m; // exact

	// Apply sign(x) to the exact weights s and d before accumulating the
	// polynomial so every subsequent rounding operation rounds in the active
	// rounding direction for negative x as well (for cosh, the caller passes
	// x >= 0 so sign_x == +0.0):
	double sign_x = fputil::fp_xor(x, x_a);
	s = fputil::fp_xor(s, sign_x);
	d = fputil::fp_xor(d, sign_x);

	// Accumulate from smallest to largest magnitude before adding the leading
	// constant s. Adding s last avoids a double rounding at the leading ULP
	// position of s when \|d * t1\| (< 2^-1) and \|s * u^2 * P_even(u^2)\| (< 2^-4)
	// are added:
	double d_t2 = d * t2;
	double s_u2 = s * u2;
	double d_u = d * u;
	double d_u3 = d_u * u2;

	// even_lo ~ s * u^2 * P_even(u^2) + d * t2
	double even_lo = fputil::multiply_add(s_u2, p_even, d_t2);
	// poly_lo ~ d * u^3 * P_odd(u^2) + even_lo
	double poly_lo = fputil::multiply_add(d_u3, p_odd, even_lo);
	// poly ~ d * t1 + poly_lo
	double poly = fputil::multiply_add(d, t1, poly_lo);

	// r ~ s * cosh(u) + d * sinh(u)
	double r = s + poly;
	return scale_exp2(r, k - 1);
	}

	#else // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS

	// Fixed-point splitting constants (ulp = 2^-22, 2^-29) for exact multiplication
	// without Veltkamp splitting:
	LIBC_INLINE_VAR constexpr double SPLIT_22 = 0x1.8p30; // 1.5 * 2^(-22 + 52)
	LIBC_INLINE_VAR constexpr double SPLIT_29 = 0x1.8p23; // 1.5 * 2^(-29 + 52)

	// Evaluate:
	// sinh(x) = (e^x - e^-x) / 2 (when is_sinh is true)
	// cosh(x) = (e^x + e^-x) / 2 (when is_sinh is false, with x >= 0)
	// in double-double precision for 2^-2 <= \|x\| < 0x1.633ce8fb9f87ep+9 (sinh) or
	// 2^-27 < x < 0x1.633ce8fb9f87ep+9 (cosh), falling back to
	// exp_pm_accurate_128(x, is_sinh) when Ziv's rounding test fails.
	//
	// Range reduction:
	// Let \|x\| = k_p * (log(2) / 2^7) + u
	// in which:
	// k_p = round(\|x\| * 2^7 * log2(e)), k_m = -k_p
	// hi_p = k_p >> 7, idx_p = k_p & 0x7f
	// hi_m = k_m >> 7, idx_m = k_m & 0x7f
	// \|u\| <= 2^-8 * log(2)
	// then:
	// e^\|x\| = 2^hi_p * 2^(idx_p / 2^7) * e^u
	// e^-\|x\| = 2^hi_m * 2^(idx_m / 2^7) * e^-u
	// = 2^hi_p * (2^(hi_m - hi_p) * 2^(idx_m / 2^7)) * e^-u.
	//
	// Factoring out 2^(hi_p - 1):
	// (e^\|x\| -/+ e^-\|x\|) / 2 = 2^(hi_p - 1) * (m_p * e^u + m_m * e^-u)
	// where:
	// m_p = 2^(idx_p / 2^7)
	// m_m = (-/+ 1) * 2^(hi_m - hi_p) * 2^(idx_m / 2^7).
	//
	// For \|x\| < 0x1.633ce8fb9f87ep+9, hi_p is in [0, 1025] so hi_p - 1 is in
	// [-1, 1024]. When hi_p == 1025, idx_p == 0 and u < 0, which keeps:
	// r = m_p * e^u + m_m * e^-u < 1.0
	// (biased exponent 1022). Adding hi_p - 1 = 1024 in scale_exp2(r, hi_p - 1)
	// produces biased exponent 2046, handling the top octave
	// [1024 * log(2), 1025 * log(2)) without a separate branch.
	LIBC_MATH_SMALL_TABLES_INLINE double exp_pm_eval(double x, bool is_sinh) {
	using FPBits = fputil::FPBits<double>;
	FPBits xbits(x);
	double x_a = xbits.abs().get_val();

	constexpr double LOG2_E_SCALE = EXP_SCALE * LOG2_E.hi;
	// Using ROUND_BIAS, ROUND_SHIFT, and ROUND_MASK (see exp_utils.h), we compute
	// k_biased = fma(x_a, LOG2_E_SCALE, ROUND_BIAS) and extract in parallel:
	// kd = (k_biased & ROUND_MASK) - ROUND_SHIFT
	// k_p = (k_bits >> 32) & 0x3'ffff.
	double k_biased = fputil::multiply_add(x_a, LOG2_E_SCALE, ROUND_BIAS);
	double kd = fputil::fp_and(k_biased, ROUND_MASK) - ROUND_SHIFT;
	int k_p =
	static_cast<int>(cpp::bit_cast<uint64_t>(k_biased) >> 32) & 0x3'ffff;
	int k_m = -k_p;

	int hi_p = k_p >> EXP_N_BITS;
	unsigned idx_p = static_cast<unsigned>(k_p & EXP_MASK);

	int hi_m = k_m >> EXP_N_BITS;
	unsigned idx_m = static_cast<unsigned>(k_m & EXP_MASK);

	// Secondary exponent offset:
	// exp_diff = hi_m - hi_p
	// clamped at -100 (where e^-\|x\| <= 2^-99 * e^\|x\| << SINH_ERR_D) to avoid
	// underflow in m_m:
	int exp_diff = hi_m - hi_p;
	int clamped_diff = (exp_diff > -100) ? exp_diff : -100;

	// Subtracting (1 << 11) = 0x800 (for sinh) sets bit 11 of the 12-bit
	// two's-complement exponent field, flipping the sign bit of m_m when shifted
	// left by 52 in scale_exp2:
	int sign_shift = static_cast<int>(is_sinh) << 11;
	int scale_m = clamped_diff - sign_shift;

	DoubleDouble m_p = EXP2_MID[idx_p];
	// Scale m_m = EXP2_MID[idx_m] by (-/+ 1) * 2^clamped_diff:
	// - Since m_m.hi is in [1, 2), we use scale_exp2(m_m.hi, scale_m) directly.
	// - Since EXP2_MID[0].lo == 0.0, we scale m_m.lo by multiplying by:
	// s_m = scale_exp2(1.0, scale_m).
	double s_m = scale_exp2(1.0, scale_m);
	DoubleDouble m_m = EXP2_MID[idx_m];
	m_m.lo *= s_m;
	m_m.hi = scale_exp2(m_m.hi, scale_m);

	// Reduced argument:
	// u = (x_a + kd * MLOG_2_EXP2_HI) + kd * MLOG_2_EXP2_MID.
	// Since MLOG_2_EXP2_HI has 35 significant bits and \|kd\| < 2^18,
	// t1 = fma(kd, MLOG_2_EXP2_HI, x_a)
	// is exact in all rounding modes with or without FMA.
	double t1 = fputil::multiply_add(kd, MLOG_2_EXP2_HI, x_a); // exact
	double u = fputil::multiply_add(kd, MLOG_2_EXP2_MID, t1);
	double u2 = u * u;

	// Separate e^(+/-u) into even and odd parts:
	// e^u ~ (1 + u^2 * P_even(u^2)) + (u + u^3 * P_odd(u^2))
	// e^-u ~ (1 + u^2 * P_even(u^2)) - (u + u^3 * P_odd(u^2))
	// where:
	// P_even(u^2) = p_0 + p_2 * u^2
	// P_odd(u^2) = p_1 + p_3 * u^2
	// Then:
	// m_p * e^u + m_m * e^-u
	// ~ (m_p + m_m) * (1 + u^2 * P_even(u^2))
	// + (m_p - m_m) * (u + u^3 * P_odd(u^2)).
	double p_even = fputil::multiply_add(u2, EXP_P_COEFFS[2], EXP_P_COEFFS[0]);
	double p_odd = fputil::multiply_add(u2, EXP_P_COEFFS[3], EXP_P_COEFFS[1]);

	#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
	// Decompose m_p * e^u + m_m * e^-u ~ r_hi + r_lo:
	// s.hi + s.lo = m_p.hi + m_m.hi (exact FastTwoSum)
	// d.hi + d.lo = m_p.hi - m_m.hi (exact FastTwoSum)
	// r_hi + r_err = s.hi + d.hi * t1
	// r_lo ~ r_err + (s.lo + (m_p.lo - m_m.lo + d.lo) * t1)
	// + ((m_p.lo + m_m.lo) + s.hi * u^2 * P_even(u^2))
	// + d.hi * (kd * MLOG_2_EXP2_MID + u^3 * P_odd(u^2)).
	//
	// For \|x\| >= 2^-2, k_p >= 46 so hi_p >= 0 and hi_m <= -1, hence
	// m_p.hi >= 1.0 > \|m_m.hi\|. By Sterbenz's lemma, both (s.hi - m_p.hi) and
	// (m_p.hi - d.hi) are exact in all rounding modes.
	DoubleDouble s = fputil::exact_add(m_p.hi, m_m.hi);
	DoubleDouble d = fputil::exact_sub(m_p.hi, m_m.hi);

	// Leading term r_hi = round(s.hi + d.hi * t1) and its rounding error r_err:
	// Since \|d.hi * t1\| <= 2^-7 * log(2) < 2^-7 and s.hi > 2^-2, (s.hi - r_hi) is
	// exact by Sterbenz's lemma.
	double r_hi = fputil::multiply_add(d.hi, t1, s.hi);
	double r_t = s.hi - r_hi; // exact
	double r_err = fputil::multiply_add(d.hi, t1, r_t);

	// Odd polynomial and low reduced-argument correction:
	// odd_lo ~ (u - t1) + u^3 * P_odd(u^2) = kd * MLOG_2_EXP2_MID + u^3 * p_odd
	double t2 = kd * MLOG_2_EXP2_MID;
	double u3 = u * u2;
	double odd_lo = fputil::multiply_add(u3, p_odd, t2);

	// Combine table low parts and polynomial terms in parallel with the
	// FastTwoSum residuals (s.lo, d.lo) and r_hi rounding error (r_err):
	// mp_lo ~ m_p.lo * (1 + t1)
	// mm_lo ~ m_m.lo * (1 - t1)
	// even_lo ~ mp_lo + (s.hi * u^2) * P_even(u^2)
	// poly_lo ~ even_lo + d.hi * odd_lo
	double mp_lo = fputil::multiply_add(m_p.lo, t1, m_p.lo);
	double mm_lo = fputil::multiply_add(-m_m.lo, t1, m_m.lo);
	double s_u2 = s.hi * u2;
	double even_lo = fputil::multiply_add(s_u2, p_even, mp_lo);
	double poly_lo = fputil::multiply_add(d.hi, odd_lo, even_lo);

	// Total linear error from FastTwoSum residuals and r_hi rounding:
	// sd_lo ~ s.lo + d.lo * t1
	// err_tot ~ r_err + sd_lo
	double sd_lo = fputil::multiply_add(d.lo, t1, s.lo);
	double err_tot = r_err + sd_lo;

	// Apply sign(x) to r_hi, mm_lo, poly_lo, and err_tot before the final Ziv
	// test additions so that (r_hi + lo_u) and (r_hi + lo_l) round in the active
	// rounding direction for negative inputs (for cosh, the caller passes x >= 0
	// so sign_x == +0.0):
	double sign_x = fputil::fp_xor(x, x_a);
	r_hi = fputil::fp_xor(r_hi, sign_x);
	mm_lo = fputil::fp_xor(mm_lo, sign_x);
	poly_lo = fputil::fp_xor(poly_lo, sign_x);
	err_tot = fputil::fp_xor(err_tot, sign_x);

	// Fold +/-SINH_ERR_D into mm_lo off the critical path so tail_{u,l} and
	// err_tot finish in parallel:
	constexpr double SINH_ERR_D = 0x1.4p-64;
	double mm_lo_u = mm_lo + SINH_ERR_D;
	double mm_lo_l = mm_lo - SINH_ERR_D;
	double tail_u = poly_lo + mm_lo_u;
	double tail_l = poly_lo + mm_lo_l;

	double lo_u = err_tot + tail_u;
	double lo_l = err_tot + tail_l;
	#else // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
	// Fixed-point grid splitting for exact leading reconstruction in all rounding
	// modes without FMA:
	// lsb(mp_hi), lsb(mm_hi) >= 2^-22 => lsb(s_hi), lsb(d_hi) >= 2^-22
	// lsb(u_hi) >= 2^-29
	// Since s_hi = mp_hi + mm_hi > 2^-2, both s_hi, d_hi, and
	// r_hi = s_hi + d_hi * u_hi (with lsb >= 2^-51)
	// are exact in double.

	// Split u = u_hi + u_lo with lsb(u_hi) >= 2^-29 and \|u_lo\| <= 2^-30:
	double u_hi_t = t1 + SPLIT_29;
	double u_hi = u_hi_t - SPLIT_29; // exact, lsb = 2^-29
	double u_lo_t = t1 - u_hi; // exact
	double u_lo = fputil::multiply_add(kd, MLOG_2_EXP2_MID, u_lo_t);

	// Split m_p = mp_hi + mp_lo and m_m = mm_hi + mm_lo with lsb >= 2^-22:
	double mp_hi_t = m_p.hi + SPLIT_22;
	double mp_hi = mp_hi_t - SPLIT_22; // exact, lsb = 2^-22
	double mp_lo_t = m_p.hi - mp_hi; // exact
	double mp_lo = mp_lo_t + m_p.lo;

	double mm_hi_t = m_m.hi + SPLIT_22;
	double mm_hi = mm_hi_t - SPLIT_22; // exact, lsb = 2^-22
	double mm_lo_t = m_m.hi - mm_hi; // exact
	double mm_lo = mm_lo_t + m_m.lo;

	// Exact leading sum, difference, and linear term r_hi = s_hi + d_hi * u_hi:
	double s_hi = mp_hi + mm_hi; // exact, lsb = 2^-22
	double d_hi = mp_hi - mm_hi; // exact, lsb = 2^-22
	double r_hi = fputil::multiply_add(d_hi, u_hi, s_hi); // exact, lsb = 2^-51

	// Low parts of the split table values and their linear contribution:
	// s_lo = mp_lo + mm_lo, d_lo = mp_lo - mm_lo
	// m_lo ~ s_lo + d_lo * u_hi
	double s_lo = mp_lo + mm_lo;
	double d_lo = mp_lo - mm_lo;
	double m_lo = fputil::multiply_add(d_lo, u_hi, s_lo);

	// Full-precision leading weights for the quadratic and cubic terms:
	double s_m_hi = m_p.hi + m_m.hi;
	double d_m_hi = m_p.hi - m_m.hi;

	// Apply sign(x) before the final polynomial accumulation so that (r_hi +
	// lo_u) and (r_hi + lo_l) round in the active rounding direction for negative
	// x:
	double sign_x = fputil::fp_xor(x, x_a);
	r_hi = fputil::fp_xor(r_hi, sign_x);
	s_m_hi = fputil::fp_xor(s_m_hi, sign_x);
	d_m_hi = fputil::fp_xor(d_m_hi, sign_x);
	m_lo = fputil::fp_xor(m_lo, sign_x);

	// Combine even and odd polynomial terms with m_lo:
	// even_lo ~ m_lo + s_m_hi * u^2 * P_even(u^2)
	// odd_lo ~ u_lo + u^3 * P_odd(u^2)
	// r_lo ~ even_lo + d_m_hi * odd_lo
	double u2_peven = u2 * p_even;
	double even_lo = fputil::multiply_add(s_m_hi, u2_peven, m_lo);
	double u3 = u * u2;
	double odd_lo = fputil::multiply_add(u3, p_odd, u_lo);
	double r_lo = fputil::multiply_add(d_m_hi, odd_lo, even_lo);

	// Total approximation and rounding error is bounded by SINH_ERR_D =
	// 0x1.4p-64:
	constexpr double SINH_ERR_D = 0x1.4p-64;
	double lo_u = r_lo + SINH_ERR_D;
	double lo_l = r_lo - SINH_ERR_D;
	#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE

	// 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))
	return scale_exp2(upper, hi_p - 1);

	return exp_pm_accurate_128(x, is_sinh);
	}

	#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS

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

	#endif // LLVM_LIBC_SRC___SUPPORT_MATH_SINHCOSH_UTILS_H