ported/excel_numeric

fn sqrt(x: &Ext80, cw: u16) -> Ext80 { $($m)::+::ext_sqrt(x, cw) }

/// `FPATAN(x, 1)`.
fn atan_with<F: Fpu>(x: f64) -> f64 {
    F::store(&F::fpatan(&load(x), &ext_one(), CW), CW)
}

/// For `|x| > 1`: `pi/2 - FPATAN(1, |x|)` in the registers, the sign of `x` after the store.
fn atan_reduced_with<F: Fpu>(x: f64) -> f64 {
    let magnitude = x.abs();
    if magnitude <= 1.0 {
        return atan_with::<F>(x);
    }
    let angle = F::fpatan(&ext_one(), &load(magnitude), CW);
    let half_pi = F::div(&ext_pi(), &load(2.0), CW);
    F::store(&F::sub(&half_pi, &angle, CW), CW).copysign(x)
}

/// One x87 operation on two binary64 operands at PC=64, then the binary64 store.
fn binary_with<F: Fpu>(op: fn(&Ext80, &Ext80, u16) -> Ext80, a: f64, b: f64) -> f64 {
    F::store(&op(&load(a), &load(b), CW), CW)
}

/// The SIN chain: `FPREM1(x, FLDPI)`, `FSIN`, the `(-1)^Q` sign.
fn sin_with<F: Fpu>(x: f64) -> f64 {
    let (r, q) = F::prem1_quo(&load(x), &ext_pi(), CW);
    let mut s = F::sin(&r, CW);
    if q & 1 == 1 {
        s = ext_chs(&s, CW);
    }
    F::store(&s, CW)
}

/// The COS chain (see [`excel_cos`]).
fn cos_with<F: Fpu>(x: f64) -> f64 {
    let minus_one = ext_chs(&ext_one(), CW);
    let pi_half = F::scale(&ext_pi(), &minus_one, CW); // exact halving
    let xa = ext_abs(&load(x), CW);
    let (r, q) = F::prem1_quo(&xa, &pi_half, CW);
    let odd_quadrant_sine = || {
        let tangent = F::tan(&r, CW);
        let square = F::mul(&tangent, &tangent, CW);
        let denominator = F::add(&ext_one(), &square, CW);
        let ratio = F::div(&square, &denominator, CW);
        let mut sine = F::sqrt(&ratio, CW);
        if tangent.0[9] & 0x80 != 0 {
            sine = ext_chs(&sine, CW);
        }
        sine
    };
    let v = match q & 3 {
        0 => F::cos(&r, CW),
        1 => ext_chs(&odd_quadrant_sine(), CW),
        2 => ext_chs(&F::cos(&r, CW), CW),
        _ => odd_quadrant_sine(),
    };
    F::store(&v, CW)
}

/// The TAN chain: `FPREM1(x, FLDPI/2)`, `FPTAN`, and `-1/tan(r)` on odd quadrants.
fn tan_with<F: Fpu>(x: f64) -> f64 {
    let minus_one = ext_chs(&ext_one(), CW);
    let pi_half = F::scale(&ext_pi(), &minus_one, CW);
    let (r, q) = F::prem1_quo(&load(x), &pi_half, CW);
    let t = F::tan(&r, CW);
    let v = if q & 1 == 0 { t } else { ext_chs(&F::div(&ext_one(), &t, CW), CW) };
    F::store(&v, CW)
}

#[inline]
fn asuint64(f: f64) -> u64 {
    f.to_bits()
}

#[inline]
fn asdouble(i: u64) -> f64 {
    f64::from_bits(i)
}

#[inline]
fn top12(x: f64) -> u32 {
    (asuint64(x) >> 52) as u32
}

#[inline]
fn top16(x: f64) -> u32 {
    (asuint64(x) >> 48) as u32
}

#[inline]
fn c2() -> f64 {
    EXP_POLY[0]
}

#[inline]
fn c3() -> f64 {
    EXP_POLY[1]
}

#[inline]
fn c4() -> f64 {
    EXP_POLY[2]
}

#[inline]
fn c5() -> f64 {
    EXP_POLY[3]
}

/// Handle results that may overflow/underflow (subnormal range): Arm's `specialcase`
/// (`math/exp.c`), with the errno-only `check_oflow/uflow` tails reduced to identity (they do
/// not change the returned value).
#[inline]
fn specialcase(tmp: f64, mut sbits: u64, ki: u64) -> f64 {
    let scale: f64;
    let mut y: f64;
    if (ki & 0x80000000) == 0 {
        // k > 0, exponent of scale might have overflowed by <= 460.
        sbits = sbits.wrapping_sub(1009u64 << 52);
        scale = asdouble(sbits);
        y = f64::from_bits(0x7f00000000000000) /* 0x1p1009 */ * (scale + scale * tmp);
        return y;
    }
    // k < 0, need special care in the subnormal range.
    sbits = sbits.wrapping_add(1022u64 << 52);
    scale = asdouble(sbits);
    y = scale + scale * tmp;
    if y < 1.0 {
        // Round y to the right precision before scaling into subnormal range
        // to avoid double rounding.
        let hi: f64;
        let mut lo: f64;
        lo = scale - y + scale * tmp;
        hi = 1.0 + y;
        lo = 1.0 - hi + y + lo;
        y = (hi + lo) - 1.0;
        // Avoid -0.0 with downward rounding (WANT_ROUNDING).
        if y == 0.0 {
            y = 0.0;
        }
    }
    y = f64::from_bits(0x0010000000000000) /* 0x1p-1022 */ * y;
    y
}

/// Excel's internal expm1 (W109 identification, 17,992/18,000 on the
/// EXPON/GAMMA.DIST-a=1 windows): msvcr100 exports no C99 `expm1`, so Excel
/// builds it from its own primitives via Kahan's cancellation-free correction
/// in DOUBLE arithmetic — `u = exp(t)` (the x87 chain, RN53); `t` when
/// `u == 1`; `(u-1)*t/ln(u)` for `|t| < 1` (ln = the fyl2x chain, == CR on
/// everything probed); plain `u - 1` otherwise. The numerator is one double
/// product, then the divide (alternate orderings refuted at 82.9%/84.9%).
pub(crate) fn excel_expm1_internal(t: f64) -> f64 {
    expm1_internal_graph(t, false)
}

/// The same cancellation correction with each arithmetic assignment stored from x87 PC64.
/// EXPON.DIST's recorded rows distinguish both the numerator product and the final division
/// from binary64 arithmetic; its density also stores the rate product through x87.
pub(crate) fn excel_expm1_spilled(t: f64) -> f64 {
    expm1_internal_graph(t, true)
}

fn expm1_internal_graph(t: f64, spilled: bool) -> f64 {
    let u = excel_exp(t);
    if u == 1.0 {
        return t;
    }
    let difference = if spilled { excel_x87_sub(u, 1.0) } else { u - 1.0 };
    if t.abs() < 1.0 {
        if spilled {
            excel_x87_div(excel_x87_mul(difference, t), excel_log(u))
        } else {
            difference * t / excel_log(u)
        }
    } else {
        difference
    }
}

/// Portable faithful (`~0.502` ULP) double-precision `e^x` — Arm optimized-routines' `exp`
/// without fused multiply-add: a correctly-rounded reference in tests and the deterministic substrate the
/// financial family composes against (W108 Bead C, pending its x87 migration).
pub(crate) fn exp_portable(x: f64) -> f64 {
    let mut abstop: u32 = top12(x) & 0x7ff;

    // 0x1p-54 top12, 512.0 top12
    let t_2m54 = top12(f64::from_bits(0x3c90000000000000)); // 0x1p-54
    let t_512 = top12(512.0);
    if abstop.wrapping_sub(t_2m54) >= t_512.wrapping_sub(t_2m54) {
        if abstop.wrapping_sub(t_2m54) >= 0x80000000 {
            // Avoid spurious underflow for tiny x. 0 is common input.
            return 1.0 + x;
        }
        if abstop >= top12(1024.0) {
            if asuint64(x) == asuint64(f64::NEG_INFINITY) {
                return 0.0;
            }
            if abstop >= top12(f64::INFINITY) {
                return 1.0 + x;
            }
            if asuint64(x) >> 63 != 0 {
                return 0.0; // __math_uflow
            } else {
                return f64::INFINITY; // __math_oflow
            }
        }
        // Large x is special cased below.
        abstop = 0;
    }

    // exp(x) = 2^(k/N) * exp(r)
    let z = EXP_INVLN2N * x;
    // !TOINT_INTRINSICS path
    let kd0 = z + EXP_SHIFT;
    let ki = asuint64(kd0);
    let kd = kd0 - EXP_SHIFT;

    let r = x + kd * EXP_NEGLN2HIN + kd * EXP_NEGLN2LON;
    let idx = (2 * (ki % N_EXP)) as usize;
    let top = ki << (52 - EXP_TABLE_BITS);
    let tail = asdouble(EXP_TAB[idx]);
    let sbits = EXP_TAB[idx + 1].wrapping_add(top);
    let r2 = r * r;
    let tmp = tail + r + r2 * (c2() + r * c3()) + r2 * r2 * (c4() + r * c5());
    if abstop == 0 {
        return specialcase(tmp, sbits, ki);
    }
    let scale = asdouble(sbits);
    scale + scale * tmp
}

/// The positive-power staging of Excel's `POWER(x, y)`: `p = exp(y · ln x)` via
/// the x87 `exp`/`ln`, with the product formed as the x87 **double-rounded**
/// `t = RN53(RN64(y · ln x))`. Bit-exact to 64-bit Excel on `x86_64`.
///
/// `base` must be finite and `> 0`, `exp` finite and `> 0`. The caller
/// (`power_fn`) owns the surrounding algorithm: sign of a negative base, the
/// `y < 0` reciprocal via [`excel_x87_recip`], the zero/integer dispatch, and
/// the worksheet error/subnormal publication rules. The returned `p` is the raw
/// positive power and may be `+Inf`, `+0`, or subnormal.
pub(crate) fn excel_pow_positive(base: f64, exp: f64) -> f64 {
    // Excel special-cases exponent 0.5 to the correctly-rounded hardware `sqrt`
    // (SSE2 SQRTSD), NOT `exp(0.5·ln x)` — the latter is 1 ULP off (e.g.
    // POWER(2, 0.5) = √2 exactly). `exp` here is `|y|`, so this also covers the
    // `y = -0.5` case, whose positive power is then reciprocated by the caller.
    // (Live-Excel confirmed 400/400; corrects the reciprocal-staging reference.)
    //
    // The 0.5 shortcut belongs to the POWER wrapper, NOT the underlying pow:
    // the distribution substrate's pow sites hit the raw chain even at
    // exponent 0.5 (W109 b24 s0p5: CR sqrt 412/800 vs chain 800/800).
    if exp == 0.5 {
        return base.sqrt();
    }
    excel_pow_chain(base, exp)
}

/// The raw CRT `pow` chain `exp(RN53(RN64(y · ln x)))` with NO special cases —
/// the distribution-substrate pow (W109 lane 1). Site-proven: BINOM `q^n`
/// (b24 11,045/11,045), WEIBULL cdf/pdf pow sites (b24+b27 all blocks 100%,
/// b27D product-staging discriminator 113/113 for the double-rounded product,
/// no `y == 0.5 → sqrt` and no integer-exponent shortcut at these sites).
/// `y` may be negative or zero (`exp(0·ln x) = 1` exactly, as Excel's chain).
/// `base` must be finite and `> 0`.
pub(crate) fn excel_pow_chain(base: f64, exp: f64) -> f64 {
    let t = excel_x87_mul(exp, excel_log(base)); // RN53(RN64(exp * ln base))
    excel_exp(t)
}

/// Excel `SIN` — bit-exact to 64-bit Excel on `x86_64` (W109 G4-01, validated
/// 1020/1020 live rows incl. held-out, build 20131): the SIN
/// chain — `FPREM1(x, FLDPI)` (the 64-bit ROM π, NOT the hardware's internal
/// 66-bit reduction), `FSIN` on the residue, `(-1)^Q` parity fixup from the
/// quotient low bit. This is the source of Excel's large-argument trig
/// "error": the reduction constant is the fldpi double-rounding of π.
pub(crate) fn excel_sin(x: f64) -> f64 {
    if !x.is_finite() {
        return f64::NAN;
    }
    sin_with::<Native>(x)
}

/// Excel `COS` — bit-exact on `x86_64` (W109 G4-01/G4-07, 2561/2561 across
/// the original corpus, a focused/randomized discovery battery, and a fresh
/// oracle-blind held-out): `|x| < 2^-26` publishes exactly `1.0` (live-probed;
/// `2^-26` itself takes the chain). Otherwise Excel reduces `|x|` with
/// `FPREM1(|x|, FLDPI/2)`. Even quadrants use `FCOS`; odd quadrants reconstruct
/// the sine magnitude as `sqrt(tan(r)^2 / (1 + tan(r)^2))` in continuous x87
/// PC64/RN arithmetic, then apply the residue and quadrant signs.
pub(crate) fn excel_cos(x: f64) -> f64 {
    if !x.is_finite() {
        return f64::NAN;
    }
    if x.abs() < f64::from_bits(0x3E50000000000000) {
        // 2^-26
        return 1.0;
    }
    cos_with::<Native>(x)
}

/// Excel `TAN` — bit-exact on `x86_64` (W109 G4-01, 1020/1020): the TAN
/// π/2-quadrant chain — `FPREM1(x, FLDPI/2)`, `FPTAN` on the residue, and on
/// odd quadrants `-1/tan(r)` with the reciprocal taken in EXTENDED before the
/// single binary64 store.
pub(crate) fn excel_tan(x: f64) -> f64 {
    if !x.is_finite() {
        return f64::NAN;
    }
    tan_with::<Native>(x)
}

/// `RN53(RN64(a + b))` — the x87 double-rounded sum of Excel's legacy
/// x87-compiled function bodies (W109: XNPV's `1 + rate` and running-total
/// accumulation; dc-6qki: the worksheet `+` and SUM's running total). Bit-exact to 64-bit Excel
/// on `x86_64`; elsewhere the same operation in the software x87.
pub(crate) fn excel_x87_add(a: f64, b: f64) -> f64 {
    #[cfg(all(target_arch = "x86_64", not(feature = "soft-x87")))]
    {
        x87::add(a, b)
    }
    #[cfg(not(all(target_arch = "x86_64", not(feature = "soft-x87"))))]
    {
        binary_with::<Soft>(Soft::add, a, b)
    }
}

/// `RN53(RN64(a - b))` — the x87 double-rounded difference (W109: NPER's
/// numerator). Bit-exact to 64-bit Excel on `x86_64`.
pub(crate) fn excel_x87_sub(a: f64, b: f64) -> f64 {
    #[cfg(all(target_arch = "x86_64", not(feature = "soft-x87")))]
    {
        x87::sub(a, b)
    }
    #[cfg(not(all(target_arch = "x86_64", not(feature = "soft-x87"))))]
    {
        binary_with::<Soft>(Soft::sub, a, b)
    }
}

/// `RN53(RN64(a * b))` — the x87 double-rounded product as a `pub(crate)`
/// spill-loop helper (the same op `POWER` uses internally via [`x87::mul`]).
pub(crate) fn excel_x87_mul(a: f64, b: f64) -> f64 {
    #[cfg(all(target_arch = "x86_64", not(feature = "soft-x87")))]
    {
        x87::mul(a, b)
    }
    #[cfg(not(all(target_arch = "x86_64", not(feature = "soft-x87"))))]
    {
        binary_with::<Soft>(Soft::mul, a, b)
    }
}

/// `RN53(RN64(RN64(a * b) + c))`: carry the PC64 product into the sum,
/// then store once to binary64. Both x87 operations use the same backend
/// on native and portable targets, with no binary64 store between them.
pub(crate) fn excel_x87_mul_add(a: f64, b: f64, c: f64) -> f64 {
    let product = Native::mul(&load(a), &load(b), CW);
    Native::store(&Native::add(&product, &load(c), CW), CW)
}

/// Temme's `c * (w +/- coefficient * t / rta)`, retaining each PC64 result
/// until the final binary64 store. The complement starts from that stored result.
pub(crate) fn excel_x87_temme_finish(
    c: f64, w: f64, coefficient: f64, t: f64, rta: f64, subtract: bool,
) -> (f64, f64) {
    let correction = Native::div(&Native::mul(&load(coefficient), &load(t), CW), &load(rta), CW);
    let inner = if subtract {
        Native::sub(&load(w), &correction, CW)
    } else {
        Native::add(&load(w), &correction, CW)
    };
    let answer = Native::store(&Native::mul(&load(c), &inner, CW), CW);
    let half = load(0.5);
    let complement = Native::store(&Native::add(&half, &Native::sub(&half, &load(answer), CW), CW), CW);
    (answer, complement)
}

/// `RN53(RN64(a / b))` — the x87 double-rounded quotient (W109: XNPV's
/// per-term `value / pow`). Bit-exact to 64-bit Excel on `x86_64`.
pub(crate) fn excel_x87_div(a: f64, b: f64) -> f64 {
    #[cfg(all(target_arch = "x86_64", not(feature = "soft-x87")))]
    {
        x87::div(a, b)
    }
    #[cfg(not(all(target_arch = "x86_64", not(feature = "soft-x87"))))]
    {
        binary_with::<Soft>(Soft::div, a, b)
    }
}

/// `RN53(RN64(1.0 / x))` — the x87 double-rounded reciprocal Excel uses to stage
/// `POWER(x, y)` for `y < 0` (positive power, stored to `binary64`, then one
/// reciprocal). Bit-exact to 64-bit Excel on `x86_64`.
pub(crate) fn excel_x87_recip(x: f64) -> f64 {
    #[cfg(all(target_arch = "x86_64", not(feature = "soft-x87")))]
    {
        x87::recip(x)
    }
    #[cfg(not(all(target_arch = "x86_64", not(feature = "soft-x87"))))]
    {
        binary_with::<Soft>(Soft::div, 1.0, x)
    }
}

/// Magnitude-reduced ATAN graph characterized by W111 live observations.
pub(crate) fn excel_atan_reduced(x: f64) -> f64 {
    #[cfg(all(target_arch = "x86_64", not(feature = "soft-x87")))]
    {
        x87::atan_reduced(x)
    }
    #[cfg(not(all(target_arch = "x86_64", not(feature = "soft-x87"))))]
    {
        atan_reduced_with::<Soft>(x)
    }
}

/// Internal: the main `__log` table path, returning the result as an unrounded
/// (hi, lo_rest) pair such that `log(x) == hi + lo_rest` (before final rounding).
/// `ix` is `asuint64(x)` for a NORMAL, positive x already outside the close-to-1
/// window (subnormals must be pre-normalized by the caller as in `excel_log`).
/// Exposing the split lets `excel_log1p` inject the `1+x` compensation term
/// into `lo_rest` and round only once.
#[inline]
fn log_general_hi_lo(ix: u64) -> (f64, f64) {
    let tmp = ix.wrapping_sub(OFF);
    let i = ((tmp >> (52 - LOG_TABLE_BITS)) % N_LOG) as usize;
    let k = (tmp as i64) >> 52; // arithmetic shift
    let iz = ix.wrapping_sub(tmp & (0xfffu64 << 52));
    let (invc, logc) = LOG_TAB[i];
    let z = asdouble(iz);

    // no-FMA path: r = (z - chi - clo) * invc
    let (chi, clo) = LOG_TAB2[i];
    let r = (z - chi - clo) * invc;
    let kd = k as f64;

    let w = kd * LOG_LN2HI + logc;
    let hi = w + r;
    let lo = w - hi + r + kd * LOG_LN2LO;

    let r2 = r * r;
    let a = &LOG_POLY;
    // lo_rest = lo + polynomial (everything except the dominant `hi` term)
    let lo_rest = lo + r2 * a[0] + r * r2 * (a[1] + r * a[2] + r2 * (a[3] + r * a[4]));
    (hi, lo_rest)
}

/// e^x - 1 — faithful (~0.5 ULP) deterministic double-precision expm1.
pub(crate) fn excel_expm1(x: f64) -> f64 {
    let ix = asuint64(x);
    let abstop = top12(x) & 0x7ff;

    // |x| < 0x1p-53: expm1(x) == x to nearest (including +-0 and subnormals),
    // because the next term x^2/2 is below 0.5 ULP of x.
    if abstop < top12(f64::from_bits(0x3ca0000000000000)) {
        // 0x1p-53
        return x;
    }

    // Special / out-of-range handling.
    // top12(0x1p6)=... use |x| >= 0x1.63108c75a1937p+9 ~ 708 overflow; and
    // x <= -0x1.2b708872320e2p+5 (~ -37.4) => expm1 == -1 to nearest.
    if abstop >= top12(708.0_f64) {
        // near the exp overflow / large-magnitude regime, or inf/nan
        if abstop >= 0x7ff {
            // inf or nan
            if ix == asuint64(f64::NEG_INFINITY) {
                return -1.0;
            }
            // +inf -> +inf ; nan -> nan (x propagates, +1 keeps nan)
            return x + x;
        }
        if x > 0.0 {
            // exp overflows -> expm1 overflows to +inf
            return f64::INFINITY;
        }
        // very negative: expm1(x) == -1 (exp(x) underflows below 0.5 ULP of 1)
        return -1.0;
    }
    // For x <= about -37 (still finite, abstop below 708 top12 only if |x|<512),
    // handle the "expm1 == -1" saturation explicitly: exp(x) < 2^-54 => e^x-1
    // rounds to -1.
    if x < -37.0 {
        return -1.0;
    }

    // ---- exp's argument reduction (no-FMA path, matches excel_exp) ----
    let z = EXP_INVLN2N * x;
    let kd0 = z + EXP_SHIFT;
    let ki = asuint64(kd0);
    let kd = kd0 - EXP_SHIFT;
    let r = x + kd * EXP_NEGLN2HIN + kd * EXP_NEGLN2LON;
    let idx = (2 * (ki % N_EXP)) as usize;
    let top = ki << (52 - EXP_TABLE_BITS);
    let tail = asdouble(EXP_TAB[idx]);
    let sbits = EXP_TAB[idx + 1].wrapping_add(top);
    let r2 = r * r;
    // tmp such that exp(x) = scale * (1 + tmp)
    let tmp = tail + r + r2 * (c2() + r * c3()) + r2 * r2 * (c4() + r * c5());
    let scale = asdouble(sbits);

    // expm1(x) = scale*(1+tmp) - 1 = (scale - 1) + scale*tmp
    // Compute (scale - 1) with a compensated split so we keep the low bits when
    // scale is close to 1 (small k). Use a two-sum of (scale - 1).
    let sm1 = scale - 1.0;
    // exact error of (scale - 1): scale - 1 is exact when |k| small... not
    // always. Recover the rounding remainder of scale-1 via Dekker two-sum with
    // the known operands (scale, -1.0).
    // two_sum(scale, -1.0): s = scale + (-1); bb = s - scale; err = (scale-(s-bb)) + (-1 - bb)
    let s = sm1;
    let bb = s - scale;
    let err = (scale - (s - bb)) + (-1.0 - bb);
    // result = s + (scale*tmp + err), added in increasing-magnitude order.
    let y = scale * tmp + err;
    s + y
}

/// log(1 + x) — faithful (~0.5 ULP) deterministic double-precision log1p.
pub(crate) fn excel_log1p(x: f64) -> f64 {
    let ix = asuint64(x);
    let top = top16(x);

    // ---- edge cases ----
    // log1p(0) = +0 ; log1p(-0) = -0 (preserve sign of zero).
    if x == 0.0 {
        return x;
    }
    if x.is_nan() {
        return x + x;
    }
    if x == f64::INFINITY {
        return f64::INFINITY;
    }
    // x < -1 -> NaN ; x == -1 -> -inf.
    if x < -1.0 {
        return (x - x) / (x - x); // NaN
    }
    if x == -1.0 {
        return f64::NEG_INFINITY;
    }

    // ---- small |x|: direct poly1 double-double kernel with f = x ----
    // The __log close-to-1 branch computes log(1+r) with r = (input) - 1. Here
    // the input IS 1+x, so r = x exactly (no rounding), giving a cancellation-
    // free log1p for |x| in [-0x1p-4, 0x1.09p-4].
    // Bounds mirror excel_log's close-to-1 window, expressed on x directly.
    let x_lo = -f64::from_bits(0x3fb0000000000000); // -0x1p-4
    let x_hi = f64::from_bits(0x3fb0900000000000); //  0x1.09p-4
    if x >= x_lo && x <= x_hi {
        let r = x;
        let r2 = r * r;
        let r3 = r * r2;
        let b = &LOG_POLY1;
        let mut y = r3
            * (b[1]
                + r * b[2]
                + r2 * b[3]
                + r3 * (b[4]
                    + r * b[5]
                    + r2 * b[6]
                    + r3 * (b[7] + r * b[8] + r2 * b[9] + r3 * b[10])));
        let w = r * f64::from_bits(0x41a0000000000000); // 0x1p27
        let rhi = r + w - w;
        let rlo = r - rhi;
        let w2 = rhi * rhi * b[0]; // b[0] == -0.5
        let hi = r + w2;
        let mut lo = r - hi + w2;
        lo += b[0] * rlo * (rhi + r);
        y += lo;
        y += hi;
        return y;
    }

    // ---- general path: form 1+x as compensated (uhi, ulo), then
    //      log1p(x) = log(uhi) + ulo/uhi. ----
    let _ = (ix, top);
    // two_sum(1.0, x) (Knuth, no ordering assumption) so uhi + ulo == 1 + x
    // EXACTLY. ulo captures the bits of x lost when rounding 1+x.
    let uhi = 1.0 + x;
    let bv = uhi - 1.0;
    let ulo = (1.0 - (uhi - bv)) + (x - bv);

    // uhi is in (0, 2) and normal here (|x| may be up to 1e6). If uhi happens to
    // land in the close-to-1 window we still get a correct result via the table
    // kernel below (poly1 already handled small x). Get log(uhi) as (hi, lo_rest)
    // and inject the correction ulo/uhi BEFORE the single final rounding:
    //   log(1+x) = log(uhi + ulo) = log(uhi) + ulo/uhi  (d log/du = 1/u).
    let (hi, mut lo_rest) = log_general_hi_lo(asuint64(uhi));
    lo_rest += ulo / uhi;
    hi + lo_rest
}

/// `ln(1 + x)` as Excel's PMT computes it (dc-j9k.12): the atanh series
/// `ln(1 + x) = 2 atanh(s) = 2s (1 + s^2/3 + s^4/5 + ...)`, `s = x / (2 + x)` (Abramowitz and
/// Stegun 4.1.27 with `z = 1 + x`), all in binary64: `s` from a binary64 sum and quotient, the
/// tail `T = 1/3 + s^2/5 + ... + s^20/23` by Horner, and the leading `2s` added last:
/// `2s + 2s * (s^2 * T)`. Not correctly rounded: it differs from the correctly rounded log1p
/// by up to 2 ulp, and those differences are what Excel's PMT shows. Twelve terms converge for
/// `|s| <= 1/5` (`-1/3 <= x <= 1/2`), the range the caller uses it on.
pub(crate) fn excel_log1p_atanh(x: f64) -> f64 {
    let s = x / (2.0 + x);
    let s2 = s * s;
    let mut tail = 1.0 / 23.0;
    for k in (1..=10).rev() {
        tail = tail * s2 + 1.0 / f64::from(2 * k + 1);
    }
    let two_s = 2.0 * s;
    two_s + two_s * (s2 * tail)
}

/// Public worksheet SQRT candidate: extended square root followed by a
/// binary64 store. W111 discovery distinguishes this from a single binary64
/// rounding at 153 of 5,676 fresh normal/zero inputs. Other callers keep their
/// separately evidenced `excel_sqrt` graph until tested.
#[inline]
pub(crate) fn excel_sqrt_staged(x: f64) -> f64 {
    Native::store(&Native::sqrt(&load(x), CW), CW)
}