gam_math/

probability.rs

use libm::erfc;
use statrs::function::beta::inv_beta_reg;

const INV_SQRT_PI: f64 = 0.564_189_583_547_756_3;
const SQRT_2_OVER_PI: f64 = 0.797_884_560_802_865_4;

/// Quantile (inverse CDF) of a Beta distribution with shape parameters `a > 0`
/// and `b > 0` at probability `p`: the value `x in [0, 1]` with
/// `I_x(a, b) = p`, where `I` is the regularized incomplete beta.
///
/// `p <= 0` maps to the support floor and `p >= 1` to the support ceiling. A
/// non-finite or non-positive shape yields `NaN`.
pub fn beta_quantile(p: f64, a: f64, b: f64) -> f64 {
    if !(a.is_finite() && a > 0.0 && b.is_finite() && b > 0.0) {
        return f64::NAN;
    }
    if !p.is_finite() || p <= 0.0 {
        return 0.0;
    }
    if p >= 1.0 {
        return 1.0;
    }
    inv_beta_reg(a, b, p)
}

/// The part of `x·x` that `f64` cannot hold: `x² = x*x + square_residual(x)`,
/// exactly, for every `x` whose square neither overflows nor goes subnormal.
///
/// This exists because of what `exp` does to a squared argument. Rounding
/// `x*x` perturbs it by at most `ulp(x²)/2` — a RELATIVE perturbation of
/// `ε/2`, which is unremarkable on its own. But `exp` converts a relative
/// perturbation `δ` of its ARGUMENT into a relative perturbation `x²·δ` of
/// its RESULT, so `exp(x*x)` carries `x²·ε/2` relative error: `3.7e-14` at
/// `x = 26`, and `7.7e-14` at the `x ≈ 37` where `φ(x)` finally underflows.
/// That is two orders worse than the `exp` evaluation's own rounding, and it
/// is the error `erfcx` and `normal_pdf` were both actually delivering.
///
/// The residual is the whole of that discarded term and is itself exactly
/// representable (Dekker's two-product theorem, in its one-FMA form), so
/// `exp(x²) = exp(x*x)·exp(residual)` and `exp(residual) = 1 + residual` to
/// `O(residual²)` — below `1e-27` over the entire domain either caller uses.
/// One multiply by `1 + residual` therefore buys back every digit, and the
/// callers below apply it fused so the correction itself costs one more
/// rounding and nothing else.
///
/// `mul_add` is a single instruction wherever FMA is in the baseline ISA
/// (aarch64, and x86-64 built with `+fma`); on a baseline x86-64 build it is
/// a `glibc` call, measured at ~2.5 ns. Against `erfcx`'s 38 ns that is 9%;
/// against `normal_pdf`'s 6.2 ns it is 40% of a function that is nowhere the
/// bottleneck of a row loop that also assembles a design row and a Hessian
/// block. Both callers guard the pathological arguments BEFORE calling this,
/// so it never has to defend `±∞` (whose residual would be `NaN`).
#[inline]
fn square_residual(x: f64, rounded_square: f64) -> f64 {
    x.mul_add(x, -rounded_square)
}

/// Standard normal PDF phi(x).
///
/// The squared argument is carried exactly (see [`square_residual`]); without
/// that, `exp(-½·fl(x*x))` degrades like `x²·ε/2` and reaches `5.7e-14`
/// relative before `φ` underflows, against the `3.3e-16` it holds with.
#[inline]
pub fn normal_pdf(x: f64) -> f64 {
    const INV_SQRT_2PI: f64 = 0.398_942_280_401_432_7;
    let rounded_square = x * x;
    let head = INV_SQRT_2PI * (-0.5 * rounded_square).exp();
    if head == 0.0 || head.is_nan() {
        // The pdf underflowed or `x` was `±∞` (head `0`), or `x` was `NaN`.
        // Neither admits a relative correction, and `±∞` would feed the
        // residual an `∞ − ∞`; return the limit the plain form gives.
        return head;
    }
    let residual = square_residual(x, rounded_square);
    head.mul_add(-0.5 * residual, head)
}

/// Standard normal CDF Phi(x) evaluated via the exact special-function identity
///
///   Phi(x) = 0.5 * erfc(-x / sqrt(2)).
///
/// This is the exact Gaussian CDF semantics used throughout the codebase. The
/// numerical `erfc` implementation may use internal approximations, but the
/// returned function is the standard normal CDF itself rather than a separate
/// polynomial surrogate surface.
#[inline]
pub fn normal_cdf(x: f64) -> f64 {
    0.5 * erfc(-x / std::f64::consts::SQRT_2)
}

/// Scaled complementary error function `erfcx(x) = exp(x²) · erfc(x)`,
/// specialized to the closed domain `x ∈ [0, +∞]`.
///
/// `+∞` maps to the exact limiting value `0`; `NaN` and negative inputs map to
/// `NaN` because they violate this restricted kernel's domain. For
/// `0 ≤ x < 26` the direct `exp(x²)·erfc(x)` form is finite. Beyond that point
/// a six-correction asymptotic expansion avoids overflow while retaining the
/// representable subnormal tail. At the switch, the first omitted term is
/// below `2e-17` relative to the leading term.
///
/// The direct branch carries `x²` exactly (see [`square_residual`]). Without
/// that correction the branch degraded like `x²·ε/2` — `1.4e-14` at `x = 10`,
/// `5.7e-14` by the top of its range — while the asymptotic branch that takes
/// over at `26` was already delivering `3e-16`. The seam was therefore a
/// 190-fold step DOWN in error at the point where the code switches to what
/// reads like the fallback, and the whole `[0, 26)` interval, where every
/// probit / Mills / log-CDF consumer actually lives, was the inaccurate side.
/// Both branches now hold `< 5e-16`, so the crossover is invisible.
#[inline]
pub fn erfcx_nonnegative(x: f64) -> f64 {
    if x.is_nan() || x < 0.0 {
        return f64::NAN;
    }
    if x == f64::INFINITY {
        return 0.0;
    }
    if x < 26.0 {
        // `x` is finite and in `[0, 26)`, so the square is exact-splittable and
        // `head` is finite and strictly positive (`erfc(26⁻) ≈ 1e-295`).
        let rounded_square = x * x;
        let head = rounded_square.exp() * erfc(x);
        head.mul_add(square_residual(x, rounded_square), head)
    } else {
        let inv = 1.0 / x;
        let inv2 = inv * inv;
        // erfcx(x) ~ 1/(sqrt(pi)x) * sum_n (-1)^n (2n-1)!!/(2x^2)^n.
        // Horner form keeps the correction well scaled when `inv2` is tiny.
        let poly = 1.0
            + inv2
                * (-0.5
                    + inv2
                        * (0.75
                            + inv2
                                * (-1.875
                                    + inv2 * (6.5625 + inv2 * (-29.53125 + inv2 * 162.421875)))));
        inv * poly * INV_SQRT_PI
    }
}

/// Computes `log(1 - exp(-a))` for `a >= 0` without cancellation.
#[inline]
pub fn log1mexp_positive(a: f64) -> f64 {
    assert!(a >= 0.0, "log1mexp_positive requires a >= 0: a={a}");
    if a > core::f64::consts::LN_2 {
        (-(-a).exp()).ln_1p()
    } else if a > 0.0 {
        (-(-a).exp_m1()).ln()
    } else {
        f64::NEG_INFINITY
    }
}

/// Numerically stable signed log-sum-exp.  Given pairs
/// `(log|aⱼ|, sign(aⱼ))` (with `signs[j] ∈ {−1, 0, +1}`), returns
/// `(log|S|, sign(S))` for `S = Σⱼ signs[j]·exp(log_mags[j])`.  Positive
/// and negative magnitudes are reduced separately with the standard
/// log-sum-exp trick (subtract the max, sum, log, add back); the two
/// partial sums are then combined via `log(|p − n|) =
/// max(log p, log n) + log1mexp(|log p − log n|)`, preserving accuracy
/// even when `p ≈ n` (catastrophic cancellation regime).  When all
/// signs are zero or all magnitudes are `−∞`, returns
/// `(NEG_INFINITY, 0.0)`.
///
/// A `+∞` log-magnitude denotes an infinite-magnitude term (`exp(+∞) = +∞`)
/// and dominates the sum: if it appears only with positive sign the result
/// is `(+∞, +1)`; only with negative sign, `(+∞, −1)` (a log-magnitude of
/// `+∞` with sign `−1` encodes the value `−∞`); with both signs the sum is
/// the indeterminate `+∞ − ∞`, returned as `(NaN, 0.0)`.  A `−∞`
/// log-magnitude is `exp(−∞) = 0` and is correctly dropped.
pub fn signed_log_sum_exp(log_mags: &[f64], signs: &[f64]) -> (f64, f64) {
    // Infinite-magnitude terms dominate any finite contribution, so resolve
    // them before the finite log-sum-exp reduction below. `−∞` log-magnitudes
    // are `exp(−∞) = 0` and need no special handling.
    let mut has_pos_inf = false;
    let mut has_neg_inf = false;
    for (idx, &lm) in log_mags.iter().enumerate() {
        if lm == f64::INFINITY {
            if signs[idx] > 0.0 {
                has_pos_inf = true;
            } else if signs[idx] < 0.0 {
                has_neg_inf = true;
            }
        }
    }
    match (has_pos_inf, has_neg_inf) {
        // P = +∞, N = +∞ ⇒ indeterminate +∞ − ∞.
        (true, true) => return (f64::NAN, 0.0),
        // P = +∞, N < ∞ ⇒ S = +∞.
        (true, false) => return (f64::INFINITY, 1.0),
        // N = +∞, P < ∞ ⇒ S = −∞, encoded as log-magnitude +∞ with sign −1.
        (false, true) => return (f64::INFINITY, -1.0),
        (false, false) => {}
    }

    let mut pos_max = f64::NEG_INFINITY;
    let mut neg_max = f64::NEG_INFINITY;
    for (idx, &lm) in log_mags.iter().enumerate() {
        if signs[idx] > 0.0 {
            pos_max = pos_max.max(lm);
        } else if signs[idx] < 0.0 {
            neg_max = neg_max.max(lm);
        }
    }

    let mut pos_sum = 0.0_f64;
    let mut neg_sum = 0.0_f64;
    for (idx, &lm) in log_mags.iter().enumerate() {
        if !lm.is_finite() {
            continue;
        }
        if signs[idx] > 0.0 {
            pos_sum += (lm - pos_max).exp();
        } else if signs[idx] < 0.0 {
            neg_sum += (lm - neg_max).exp();
        }
    }

    let log_pos = if pos_sum > 0.0 {
        pos_max + pos_sum.ln()
    } else {
        f64::NEG_INFINITY
    };
    let log_neg = if neg_sum > 0.0 {
        neg_max + neg_sum.ln()
    } else {
        f64::NEG_INFINITY
    };

    if log_pos == f64::NEG_INFINITY && log_neg == f64::NEG_INFINITY {
        // Both partial sums are empty: no terms at all, all signs zero, or every
        // magnitude `−∞` (each `exp(−∞) = 0`). The signed sum is exactly `0`, so
        // the contract requires `(−∞, 0.0)` — NOT the positive-sum convention,
        // which would mislabel a zero as `+1` and corrupt any downstream cascade
        // that reads back the sign.
        return (f64::NEG_INFINITY, 0.0);
    }
    if log_neg == f64::NEG_INFINITY {
        return (log_pos, 1.0);
    }
    if log_pos == f64::NEG_INFINITY {
        return (log_neg, -1.0);
    }
    if log_pos > log_neg {
        let gap = log_pos - log_neg;
        (log_pos + log1mexp_positive(gap), 1.0)
    } else if log_neg > log_pos {
        let gap = log_neg - log_pos;
        (log_neg + log1mexp_positive(gap), -1.0)
    } else {
        (f64::NEG_INFINITY, 0.0)
    }
}

/// Numerically stable `ln Φ(x)` for the standard normal CDF. For `x ≥ 0`,
/// evaluates `ln(1 - 0.5 erfc(x/sqrt(2)))` with `ln_1p`, retaining the small
/// negative result after `Φ(x)` itself rounds to one. For `x < 0`, rewrites
/// `ln Φ(x) = −u² + ln(½·erfcx(u))`, `u = −x/√2`,
/// which preserves digits throughout the representable left tail without a
/// probability floor. Returns the corresponding IEEE limit at infinities and
/// propagates `NaN`.
#[inline]
pub fn normal_logcdf(x: f64) -> f64 {
    if x == f64::INFINITY {
        return 0.0;
    }
    if x == f64::NEG_INFINITY {
        return f64::NEG_INFINITY;
    }
    if x.is_nan() {
        return f64::NAN;
    }
    if x < 0.0 {
        let (u, scaled_tail) = negative_normal_tail_components(x);
        negative_normal_logcdf_from_scaled_tail(u, scaled_tail)
    } else {
        let upper_tail = 0.5 * erfc(x / std::f64::consts::SQRT_2);
        (-upper_tail).ln_1p()
    }
}

/// Numerically stable `ln(1 − Φ(x)) = ln Φ(−x)` for the standard normal
/// survival function.  Delegates to `normal_logcdf(-x)` so the deep-right
/// tail benefits from the same `erfcx`-based representation.
#[inline]
pub fn normal_logsf(x: f64) -> f64 {
    normal_logcdf(-x)
}

/// Joint evaluation of `ln Φ(x)` and the Mills-ratio analogue
/// `φ(x) / Φ(x)`, signed for the symmetric branch.  Used by the latent
/// probit families where the inverse-link gradient needs the ratio and
/// the likelihood needs the log-CDF on the same `x`; computing both in
/// one call shares the `erfcx` evaluation that dominates the cost in the
/// deep tail.
#[inline]
pub fn signed_probit_logcdf_and_mills_ratio(x: f64) -> (f64, f64) {
    if x == f64::INFINITY {
        return (0.0, 0.0);
    }
    if x == f64::NEG_INFINITY {
        return (f64::NEG_INFINITY, f64::INFINITY);
    }
    if x.is_nan() {
        return (f64::NAN, f64::NAN);
    }
    if x < 0.0 {
        let (u, scaled_tail) = negative_normal_tail_components(x);
        (
            negative_normal_logcdf_from_scaled_tail(u, scaled_tail),
            SQRT_2_OVER_PI / scaled_tail,
        )
    } else {
        let upper_tail = 0.5 * erfc(x / std::f64::consts::SQRT_2);
        let cdf = 1.0 - upper_tail;
        let lambda = normal_pdf(x) / cdf;
        ((-upper_tail).ln_1p(), lambda)
    }
}

#[inline]
fn negative_normal_tail_components(x: f64) -> (f64, f64) {
    assert!(x.is_finite() && x < 0.0);
    let u = -x / std::f64::consts::SQRT_2;
    (u, erfcx_nonnegative(u))
}

#[inline]
fn negative_normal_logcdf_from_scaled_tail(u: f64, scaled_tail: f64) -> f64 {
    -u * u + scaled_tail.ln() - std::f64::consts::LN_2
}

/// Stable value and first four derivatives of `ln Φ(x)`.
///
/// The moderate regime uses the exact Mills-ratio recurrence, with the brackets
/// collected in `q = λ + x` once `x < 0` so that they do not cancel as `λ`
/// closes on `−x`. In the deep left tail, differentiating the Laplace continued
/// fraction
///
/// `φ(t)/Φ(-t) = t + 1/(t + 2/(t + 3/(...)))`, `t = -x`,
///
/// carries the small correction to `t` independently, so `f'' -> -1` and the
/// higher derivatives approach zero without subtracting nearly equal `f64`s.
/// In the right tail, signed log-magnitude sums preserve polynomially weighted
/// derivatives even when `φ(x)/Φ(x)` itself has rounded to zero.
#[inline]
pub fn normal_logcdf_derivatives(x: f64) -> [f64; 5] {
    if x.is_nan() {
        return [f64::NAN; 5];
    }
    if x == f64::INFINITY {
        return [0.0; 5];
    }
    if x == f64::NEG_INFINITY {
        return [f64::NEG_INFINITY, f64::INFINITY, -1.0, 0.0, 0.0];
    }

    const RIGHT_LOG_MAGNITUDE_SWITCH: f64 = 8.0;
    if x <= LEFT_CONTINUED_FRACTION_SWITCH {
        return normal_logcdf_derivatives_left_tail(x);
    }
    if x >= RIGHT_LOG_MAGNITUDE_SWITCH {
        return normal_logcdf_derivatives_right_tail(x);
    }

    let (log_cdf, lambda) = signed_probit_logcdf_and_mills_ratio(x);
    let x2 = x * x;
    if x < 0.0 {
        // Left of the origin the brackets below are collected in the SAME Mills
        // correction `q = λ + x` the continued-fraction branch carries, because
        // written in `λ` they cancel catastrophically long before the branch
        // ends. `λ(x) → −x` as `x → −∞`, so every term of, say,
        // `(x³−3x) + (7x²−4)λ + 12xλ² + 6λ³` grows like `|x|³` while their sum
        // decays: at `x = −4` they are `−52`, `456`, `−857`, `453` and add to
        // `−0.0023`, a cancellation of 380000 that costs eleven digits. In `q`
        // the same bracket is `−6q³ + 6xq² + (4−x²)q − x`, whose terms are
        // `−0.069`, `−1.22`, `−2.71`, `4` — a cancellation of 1847, three
        // orders milder. The reformulation is exact (`λ = q − x` substituted and
        // re-collected), costs the same flops, and buys 16–34x across the whole
        // branch: worst over `x ∈ [−4, 0]` falls from `4.5e−11` to `2.8e−12`.
        //
        // `q` itself is safe to form here: `λ/2 ≤ |x| ≤ 2λ` holds over most of
        // the range, so `λ + x` is EXACT by Sterbenz, and where it is not (`x`
        // near 0) `q` is the same size as `λ` and nothing cancels. That is the
        // whole reason the rewrite works — it moves the cancellation out of the
        // brackets and into a subtraction that has none.
        //
        // Past the origin `q → x` is no longer small, the `λ` form has nothing
        // to cancel (`λ → 0` and `x² − 1` dominates), and it is the more
        // accurate of the two — hence the sign test rather than a blanket swap.
        let q = lambda + x;
        let q2 = q * q;
        return [
            log_cdf,
            lambda,
            -lambda * q,
            lambda * (2.0 * q2 - x * q - 1.0),
            lambda * (-6.0 * q2 * q + 6.0 * x * q2 + (4.0 - x2) * q - x),
        ];
    }
    let lambda2 = lambda * lambda;
    let lambda3 = lambda2 * lambda;
    [
        log_cdf,
        lambda,
        -lambda * (x + lambda),
        lambda * (x2 - 1.0 + 3.0 * x * lambda + 2.0 * lambda2),
        -lambda
            * ((x * x2 - 3.0 * x) + (7.0 * x2 - 4.0) * lambda + 12.0 * x * lambda2 + 6.0 * lambda3),
    ]
}

#[derive(Clone, Copy)]
struct MillsCorrectionDerivatives {
    value: f64,
    first: f64,
    second: f64,
    third: f64,
}

/// `x` at or below which the left-tail Mills ratio is taken from the Laplace
/// continued fraction rather than from `erfcx`. Equivalently `t = −x ≥ 4`.
const LEFT_CONTINUED_FRACTION_SWITCH: f64 = -4.0;

/// The Laplace continued-fraction **correction** to the left-tail Mills ratio,
///
/// `q(t) = λ(−t) − t = 1/(t + 2/(t + 3/(...)))`,   `λ(x) = φ(x)/Φ(x)`,
///
/// together with its first three derivatives in `t`. Requires `t ≥ 4`.
///
/// `q` is the whole content of the left tail that is NOT the leading `t`: it
/// decays like `1/t − 2/t³ + 10/t⁵ − ...`, and every operation building it is
/// a division or an addition of positive quantities, so it carries full
/// relative precision no matter how small it gets. That is the property its
/// two consumers need, and it is why the correction is returned separately
/// instead of pre-added to `t`:
///
/// * [`normal_logcdf_derivatives_left_tail`] needs `f'' = −(1 + q')` and the
///   higher derivatives, which tend to `−1` and `0` and would be destroyed by
///   differencing nearly equal `f64`s.
/// * [`cone_boundary_log_factor_and_derivatives`] needs `∂corr/∂a = b − q(t)`,
///   which is the same statement one substitution away (#2306 §4).
///
/// Recovering `q` from a separately computed `λ` — `q = λ − t` — is exactly the
/// cancellation this exists to avoid, and it is not a small effect: at `t = 1e8`
/// it costs every significant digit, and past `t ≈ 2e8` it returns the wrong
/// SIGN. The reference itself has to be carried at ~120 decimal digits before it
/// reproduces what this recursion gives in binary64.
#[inline]
fn mills_correction_continued_fraction(t: f64) -> MillsCorrectionDerivatives {
    assert!(t.is_finite() && t >= 4.0);
    let mut q = MillsCorrectionDerivatives {
        value: 0.0,
        first: 0.0,
        second: 0.0,
        third: 0.0,
    };
    // The truncation error is damped by a product of the continued-fraction
    // sensitivities `n/(t + q)^2`, so the depth must be sized at `t = 4` — the
    // LEAST converged point of the domain, and the one the log-CDF branch sits
    // exactly on. Each successive derivative converges roughly 15x slower than
    // the last, because differentiating the recursion multiplies each level's
    // contribution by another factor of that same sensitivity. Measured against
    // a 60-digit reference at `t = 4`:
    //
    // ```text
    //            q         q'        q''       q'''
    //   32   1.9e-15    7.0e-14    1.4e-12    2.1e-11
    //   64   2.3e-23    1.4e-21    4.4e-20    1.0e-18
    // ```
    //
    // 32 levels is enough for the VALUE and nothing else: it leaves `q'''` — the
    // fourth log-CDF derivative — wrong in its eleventh digit. The depths that
    // first reach `1e-17` at `t = 4` are 41, 47, 53 and 60 for the four
    // channels, so 64 covers the worst of them with ~200x of margin, and the
    // requirement falls off fast enough (33 levels at `t = 6`, 24 at `t = 8`,
    // 12 at `t = 20`) that one constant sized for the edge is safe everywhere
    // above it. The extra levels are pure convergence — every step divides
    // positive quantities — so they cannot destabilise a large `t`.
    for n in (1..=64).rev() {
        let denominator = t + q.value;
        let inv_denominator = denominator.recip();
        let value = f64::from(n) / denominator;
        let denominator_first = 1.0 + q.first;
        let a = denominator_first * inv_denominator;
        let b = q.second * inv_denominator;
        let c = q.third * inv_denominator;
        q = MillsCorrectionDerivatives {
            value,
            first: -value * denominator_first / denominator,
            second: value * (2.0 * a * a - b),
            third: value * (-6.0 * a * a * a + 6.0 * a * b - c),
        };
    }
    q
}

#[inline]
fn normal_logcdf_derivatives_left_tail(x: f64) -> [f64; 5] {
    assert!(x.is_finite() && x <= LEFT_CONTINUED_FRACTION_SWITCH);
    let t = -x;
    let q = mills_correction_continued_fraction(t);
    [
        normal_logcdf(x),
        t + q.value,
        -(1.0 + q.first),
        q.second,
        -q.third,
    ]
}

#[inline]
fn normal_logcdf_derivatives_right_tail(x: f64) -> [f64; 5] {
    assert!(x.is_finite() && x >= 8.0);
    const LOG_SQRT_2PI: f64 = 0.918_938_533_204_672_7;
    let log_cdf = normal_logcdf(x);
    let u = x / std::f64::consts::SQRT_2;
    let log_lambda = -u * u - LOG_SQRT_2PI - log_cdf;
    let log_x = x.ln();
    let inv_x2 = x.recip() * x.recip();

    let first = log_lambda.exp();
    let second = signed_exp_sum(&[log_x + log_lambda, 2.0 * log_lambda], &[-1.0, -1.0]);
    let third = signed_exp_sum(
        &[
            2.0 * log_x + (-inv_x2).ln_1p() + log_lambda,
            3.0_f64.ln() + log_x + 2.0 * log_lambda,
            2.0_f64.ln() + 3.0 * log_lambda,
        ],
        &[1.0, 1.0, 1.0],
    );
    let fourth = signed_exp_sum(
        &[
            3.0 * log_x + (-3.0 * inv_x2).ln_1p() + log_lambda,
            7.0_f64.ln() + 2.0 * log_x + (-(4.0 / 7.0) * inv_x2).ln_1p() + 2.0 * log_lambda,
            12.0_f64.ln() + log_x + 3.0 * log_lambda,
            6.0_f64.ln() + 4.0 * log_lambda,
        ],
        &[-1.0, -1.0, -1.0, -1.0],
    );
    [log_cdf, first, second, third, fourth]
}

#[inline]
fn signed_exp_sum(log_magnitudes: &[f64], signs: &[f64]) -> f64 {
    let (log_magnitude, sign) = signed_log_sum_exp(log_magnitudes, signs);
    if sign == 0.0 {
        0.0
    } else {
        sign * log_magnitude.exp()
    }
}

#[inline]
fn acklam_lower_tail_quantile_from_log_probability(log_p: f64) -> f64 {
    const C: [f64; 6] = [
        -7.784_894_002_430_293e-3,
        -3.223_964_580_411_365e-1,
        -2.400_758_277_161_838,
        -2.549_732_539_343_734,
        4.374_664_141_464_968,
        2.938_163_982_698_783,
    ];
    const D: [f64; 4] = [
        7.784_695_709_041_462e-3,
        3.224_671_290_700_398e-1,
        2.445_134_137_142_996,
        3.754_408_661_907_416,
    ];
    let q = (-2.0 * log_p).sqrt();
    (((((C[0] * q + C[1]) * q + C[2]) * q + C[3]) * q + C[4]) * q + C[5])
        / ((((D[0] * q + D[1]) * q + D[2]) * q + D[3]) * q + 1.0)
}

/// Standard normal quantile Φ⁻¹(p) using Acklam's rational approximation.
#[inline]
pub fn standard_normal_quantile(p: f64) -> Result<f64, String> {
    if !(p.is_finite() && p > 0.0 && p < 1.0) {
        return Err(format!("normal quantile requires p in (0,1), got {p}"));
    }

    const A: [f64; 6] = [
        -3.969_683_028_665_376e1,
        2.209_460_984_245_205e2,
        -2.759_285_104_469_687e2,
        1.383_577_518_672_69e2,
        -3.066_479_806_614_716e1,
        2.506_628_277_459_239,
    ];
    const B: [f64; 5] = [
        -5.447_609_879_822_406e1,
        1.615_858_368_580_409e2,
        -1.556_989_798_598_866e2,
        6.680_131_188_771_972e1,
        -1.328_068_155_288_572e1,
    ];
    const P_LOW: f64 = 0.02425;
    const P_HIGH: f64 = 1.0 - P_LOW;

    let mut x = if p < P_LOW {
        acklam_lower_tail_quantile_from_log_probability(p.ln())
    } else if p <= P_HIGH {
        let q = p - 0.5;
        let r = q * q;
        (((((A[0] * r + A[1]) * r + A[2]) * r + A[3]) * r + A[4]) * r + A[5]) * q
            / (((((B[0] * r + B[1]) * r + B[2]) * r + B[3]) * r + B[4]) * r + 1.0)
    } else {
        -acklam_lower_tail_quantile_from_log_probability((1.0 - p).ln())
    };
    for _ in 0..2 {
        let density = normal_pdf(x);
        if !(density.is_finite() && density > 0.0) {
            break;
        }
        // Residual F(x) − p, formed without catastrophic cancellation in
        // either tail. For an upper-tail iterate `x > 0`, `normal_cdf(x)`
        // saturates to ~1, so the direct `normal_cdf(x) − p` annihilates the
        // tiny residual the polish must act on; instead use the upper-tail
        // complement `F(x) − p = (1 − p) − 0.5·erfc(x/√2)`, where both terms
        // are the small upper-tail quantities (`1 − p` is exact by Sterbenz
        // for `p ∈ [½,1)`). For `x ≤ 0`, `normal_cdf(x) = 0.5·erfc(|x|/√2)` is
        // itself the faithfully carried small lower-tail value, so the direct
        // form is already cancellation-free.
        let residual = if x > 0.0 {
            (1.0 - p) - 0.5 * erfc(x / std::f64::consts::SQRT_2)
        } else {
            normal_cdf(x) - p
        };
        let correction = residual / density;
        let denominator = 1.0 + 0.5 * x * correction;
        if !(correction.is_finite() && denominator.is_finite() && denominator != 0.0) {
            break;
        }
        let step = correction / denominator;
        if !step.is_finite() {
            break;
        }
        x -= step;
        if step.abs() <= 2.0 * f64::EPSILON * x.abs().max(1.0) {
            break;
        }
    }
    Ok(x)
}

/// Standard normal quantile from `log_p = ln Φ(x)`.
///
/// Unlike [`standard_normal_quantile`], this remains defined when `Φ(x)` is
/// smaller than the least positive `f64`, and when `Φ(x)` is so close to one
/// that exponentiating `log_p` rounds to exactly one. Acklam's lower-tail
/// approximation supplies the initial point; Newton polishing solves
/// `ln Φ(x) = log_p` with the stable log-CDF and Mills ratio, so neither tail
/// forms a probability-space subtraction.
#[inline]
pub fn standard_normal_quantile_from_log_cdf(log_p: f64) -> Result<f64, String> {
    if !(log_p.is_finite() && log_p < 0.0) {
        return Err(format!(
            "normal log-quantile requires finite log_p < 0, got {log_p}"
        ));
    }

    if log_p > -std::f64::consts::LN_2 {
        // Reflect through the upper tail without forming `1 - exp(log_p)`.
        let log_q = (-log_p.exp_m1()).ln();
        return standard_normal_quantile_from_log_cdf(log_q).map(|x| -x);
    }

    let p = log_p.exp();
    let mut x = if p > 0.0 {
        standard_normal_quantile(p)?
    } else {
        acklam_lower_tail_quantile_from_log_probability(log_p)
    };
    for _ in 0..4 {
        let (current_log_p, mills_ratio) = signed_probit_logcdf_and_mills_ratio(x);
        if !(current_log_p.is_finite() && mills_ratio.is_finite() && mills_ratio > 0.0) {
            break;
        }
        let step = (current_log_p - log_p) / mills_ratio;
        if !step.is_finite() {
            break;
        }
        x -= step;
        if step.abs() <= 2.0 * f64::EPSILON * x.abs().max(1.0) {
            break;
        }
    }
    Ok(x)
}

/// Log of the standardized one-sided truncated-Gaussian boundary factor for
/// the constrained-LAML cone correction (gam#2306 §4).
///
/// For a constraint coordinate with Lagrange multiplier `μ ≥ 0`, normal
/// curvature `h > 0`, and signed interior slack `s ≥ 0`, the exact 1-D
/// boundary integral is
///
/// ```text
///   ∫_{−s}^{∞} exp(−μ·u − ½·h·u²) du
///     = √(2π/h) · exp(μ²/(2h)) · Φ(s·√h − μ/√h),
/// ```
///
/// and the correction of the Laplace criterion RELATIVE to the unrestricted
/// Gaussian factor `√(2π/h)` is, in the standardized arguments
/// `a = μ/√h ≥ 0`, `b = s·√h ≥ 0`:
///
/// ```text
///   corr(a, b) = a²/2 + ln Φ(b − a).
/// ```
///
/// Key limits (the #2306 derivation's continuity contract): an activating
/// row (`a = 0`, `b = 0`) contributes exactly `ln ½` (the half-Gaussian); a
/// deep interior row (`b − a → ∞`) contributes `→ 0`, reducing byte-exactly
/// to the unrestricted LAML; a hard-pushed active row (`a → ∞`, `b = 0`)
/// follows the exact linear-decay limit `corr → −ln(a·√(2π))`.
///
/// Evaluated FUSED: computing `a²/2` and `ln Φ(b−a)` as two separate f64
/// terms cancels catastrophically once `a ≳ 10⁴` (both grow like `±a²/2`).
/// On the `b < a` branch the sum collapses analytically to
/// `a·b − b²/2 + ln(erfcx((a−b)/√2)/2)`, which is cancellation-free (for an
/// active row, `b = 0`, it is a single `erfcx` evaluation).
#[must_use]
pub fn cone_boundary_log_factor(mu_over_sqrt_h: f64, slack_times_sqrt_h: f64) -> f64 {
    let a = mu_over_sqrt_h;
    let b = slack_times_sqrt_h;
    if !(a.is_finite() && b.is_finite()) || a < 0.0 || b < 0.0 {
        return f64::NAN;
    }
    let xi = b - a;
    if xi >= 0.0 {
        // Interior-dominant: ln Φ(ξ) is a small negative number and a²/2 is
        // exact; no cancellation between them (a ≤ b here, so a²/2 ≤ ab −
        // b²/2 + O(1) stays modest whenever the factor itself is modest).
        0.5 * a * a + normal_logcdf(xi)
    } else {
        // Active-dominant: fused analytic collapse of a²/2 + ln Φ(−(a−b)).
        let u = (a - b) / std::f64::consts::SQRT_2;
        a * b - 0.5 * b * b + (0.5 * erfcx_nonnegative(u)).ln()
    }
}

/// [`cone_boundary_log_factor`] together with its exact partial derivatives
/// in the standardized arguments — the pieces the outer ρ-gradient chains
/// through `(μ̃, h̃, s)(ρ)` (gam#2306 §4 "the g-factors differentiate in
/// closed form"). With `ξ = b − a` and the Mills ratio `λ(ξ) = φ(ξ)/Φ(ξ)`:
///
/// ```text
///   ∂corr/∂a = a − λ(ξ),      ∂corr/∂b = λ(ξ).
/// ```
///
/// `∂corr/∂b = λ(ξ)` is a single `erfcx` evaluation and needs nothing further.
///
/// `∂corr/∂a` does. Written literally as `a − λ(ξ)` it is a subtraction of two
/// quantities that both grow like `a`, because `λ(−t) = t + q(t)` with
/// `q(t) ~ 1/t`: the answer is the SMALL correction `q`, and forming it by
/// subtraction destroys `log₁₀(a²·ε)` digits of it. The value
/// [`cone_boundary_log_factor`] is fused precisely to dodge the twin of this
/// cancellation, and the gradient has to be fused the same way rather than
/// re-derived from a `λ` that has already lost the digits.
///
/// So on the active branch the correction is taken directly from the Laplace
/// continued fraction ([`mills_correction_continued_fraction`]), the same one
/// the left-tail log-CDF derivatives use, under the substitution
///
/// ```text
///   ξ = b − a,  t = −ξ = a − b  ⇒  ∂corr/∂a = a − λ(ξ) = a − (t + q(t)) = b − q(t),
/// ```
///
/// which is cancellation-free for every `a`: `b ≥ 0` and `q(t) ∈ (0, ¼]`. The
/// deep-active limit `∂corr/∂a → −1/a` then holds to full relative precision
/// instead of to none, and the sign is right (the factor is strictly decreasing
/// in `a`, so `∂corr/∂a < 0` whenever `b = 0`).
///
/// Measured on `a ∈ [10⁻², 10¹⁴] × b ∈ {0, …, 10³}` against a 250-digit
/// reference: the subtractive form reaches `6.1e6` relative error and turns
/// positive past `a ≈ 2e8`; this form is within `7.5e-16` — about 3 ulp — of
/// the truth, measured against the magnitudes entering the subtraction rather
/// than against the result. That is the right denominator because `∂corr/∂a`
/// genuinely passes through ZERO along the curve `b = q(a − b)` (the factor is
/// increasing in `a` for slack rows and decreasing for active ones), and no
/// representation carries relative precision across its own root; near it the
/// error is bounded in absolute terms by `ε·b`, which is what a gradient
/// consumer needs.
#[must_use]
pub fn cone_boundary_log_factor_and_derivatives(
    mu_over_sqrt_h: f64,
    slack_times_sqrt_h: f64,
) -> (f64, f64, f64) {
    let a = mu_over_sqrt_h;
    let b = slack_times_sqrt_h;
    let value = cone_boundary_log_factor(a, b);
    if value.is_nan() {
        // The value's domain guard (finite, non-negative `a` and `b`) is the
        // function's domain; a gradient off it is not defined either, and
        // returning a finite one next to a NaN value would read as usable.
        return (value, f64::NAN, f64::NAN);
    }
    let xi = b - a;
    let (_, mills) = signed_probit_logcdf_and_mills_ratio(xi);
    let d_a = if xi <= LEFT_CONTINUED_FRACTION_SWITCH {
        b - mills_correction_continued_fraction(-xi).value
    } else {
        // `|ξ| < 4`, so `λ(ξ) < λ(−4) ≈ 4.26` and `a = b − ξ` is bounded by it:
        // the subtraction is between two `O(1)` quantities and loses nothing
        // that matters.
        a - mills
    };
    (value, d_a, mills)
}

#[cfg(test)]
mod cone_boundary_factor_tests {
    use super::*;

    /// Adaptive-free Simpson quadrature of the exact 1-D boundary integral
    /// `∫_{−s}^{U} exp(−μu − ½hu²) du` on a truncation `U` chosen so the
    /// discarded tail is below 1e-18 of the mass.
    fn quadrature_log_relative_factor(mu: f64, h: f64, s: f64) -> f64 {
        let upper = ((-mu / h) + 12.0 / h.sqrt()).max(-s + 12.0 / h.sqrt());
        let lower = -s;
        let n = 40_000usize;
        let step = (upper - lower) / n as f64;
        let f = |u: f64| (-mu * u - 0.5 * h * u * u).exp();
        let mut acc = f(lower) + f(upper);
        for i in 1..n {
            let u = lower + step * i as f64;
            acc += if i % 2 == 1 { 4.0 } else { 2.0 } * f(u);
        }
        let integral = acc * step / 3.0;
        (integral / (2.0 * std::f64::consts::PI / h).sqrt()).ln()
    }

    /// The closed form must match direct quadrature of the defining integral
    /// across active (s=0), interior (μ=0), and mixed regimes (gam#2306 §4).
    #[test]
    fn boundary_factor_matches_quadrature_across_regimes() {
        let cases: [(f64, f64, f64); 8] = [
            (0.0, 1.0, 0.0),  // activating row: exactly ln ½
            (0.0, 4.0, 0.0),  // curvature does not move the standardized value
            (2.5, 1.0, 0.0),  // active with a real multiplier
            (30.0, 9.0, 0.0), // deep linear-decay limit
            (0.0, 1.0, 0.7),  // interior near-boundary
            (0.0, 2.0, 4.0),  // interior far: → 0
            (1.5, 0.5, 2.0),  // mixed multiplier + slack
            (4.0, 2.0, 1.0),  // active-dominant mixed
        ];
        for &(mu, h, s) in &cases {
            let a = mu / h.sqrt();
            let b = s * h.sqrt();
            let closed = cone_boundary_log_factor(a, b);
            let quad = quadrature_log_relative_factor(mu, h, s);
            assert!(
                (closed - quad).abs() <= 1e-9 * (1.0 + quad.abs()),
                "(μ={mu}, h={h}, s={s}): closed {closed} vs quadrature {quad}"
            );
        }
        assert!(
            (cone_boundary_log_factor(0.0, 0.0) - 0.5_f64.ln()).abs() < 1e-15,
            "an activating row must contribute exactly the half-Gaussian ln ½"
        );
    }

    /// The deep-active limit is the exact linear decay `corr → −ln(a·√(2π))`,
    /// and the deep-interior limit vanishes — the two continuity anchors that
    /// make the constrained criterion reduce to the unrestricted LAML away
    /// from the boundary.
    #[test]
    fn boundary_factor_limits_are_exact() {
        let a = 1.0e6;
        let expected = -(a * (2.0 * std::f64::consts::PI).sqrt()).ln();
        let got = cone_boundary_log_factor(a, 0.0);
        assert!(
            (got - expected).abs() <= 1e-9 * expected.abs(),
            "deep-active: got {got}, expected {expected}"
        );
        let interior = cone_boundary_log_factor(0.0, 40.0);
        assert!(
            interior.abs() < 1e-300 || interior > -1e-12,
            "deep-interior must vanish; got {interior}"
        );
    }

    /// The deep-active GRADIENT has to survive as far as the deep-active VALUE
    /// does. `∂corr/∂a = a − λ(−a)` is the small residual left by two terms
    /// that both grow like `a`, so writing it as that subtraction loses
    /// `log₁₀(a²·ε)` digits: at `a = 1e6` it was already 4 digits down, at
    /// `a = 2e8` it came back POSITIVE, and past `a = 5e8` it was flat zero
    /// while the true value is `−2e-9`. The value alongside it was correct to
    /// 15 digits the whole way, which is what made the defect quiet.
    ///
    /// The reference here is the asymptotic series of the Mills correction,
    /// `λ(−a) = a + 1/a − 2/a³ + 10/a⁵ − 74/a⁷ + …` (so `∂corr/∂a = −1/a +
    /// 2/a³ − …`), which is the cheapest exact statement of the limit and is
    /// good to well past f64 from `a = 100` up. Finite differences cannot gate
    /// this: the quantity under test is smaller than any usable FD step's own
    /// truncation error.
    #[test]
    fn boundary_factor_active_gradient_holds_to_the_representable_limit() {
        let mut a = 100.0_f64;
        while a <= 1.0e14 {
            let (_, d_a, _) = cone_boundary_log_factor_and_derivatives(a, 0.0);
            let inv = 1.0 / a;
            let expected = -inv + 2.0 * inv.powi(3) - 10.0 * inv.powi(5) + 74.0 * inv.powi(7);
            assert!(
                d_a < 0.0,
                "corr is strictly decreasing in a at b=0, so ∂a must stay negative; \
                 got {d_a} at a={a}"
            );
            assert!(
                (d_a - expected).abs() <= 1.0e-13 * expected.abs(),
                "deep-active ∂a at a={a}: got {d_a}, expected {expected} \
                 (rel {:.3e})",
                (d_a - expected).abs() / expected.abs()
            );
            a *= 10.0;
        }
    }

    /// `∂corr/∂a + ∂corr/∂b = a` identically, since the two partials are
    /// `a − λ(ξ)` and `λ(ξ)` for the same `ξ`. The two are now computed by
    /// different routes in the active branch — a continued fraction and an
    /// `erfcx` — so this is the gate that they still describe one function.
    #[test]
    fn boundary_factor_partials_sum_to_a() {
        for &a in &[0.0_f64, 0.5, 3.0, 4.0, 12.0, 1.0e3, 1.0e7, 1.0e12] {
            for &b in &[0.0_f64, 1.0e-3, 0.9, 5.0, 1.0e3] {
                let (_, d_a, d_b) = cone_boundary_log_factor_and_derivatives(a, b);
                assert!(
                    (d_a + d_b - a).abs() <= 1.0e-14 * a.max(d_b).max(1.0),
                    "(a={a}, b={b}): ∂a {d_a} + ∂b {d_b} = {} ≠ a",
                    d_a + d_b
                );
            }
        }
    }

    /// The continued-fraction branch and the direct `a − λ` form must agree
    /// just inside the `ξ ≤ −4` switch, where the subtraction still has most of
    /// its digits. Without this, the branch could be precise and WRONG — the
    /// accuracy gate above pins a limit the continued fraction could hit while
    /// disagreeing with the function it is supposed to be differentiating.
    ///
    /// The band is set by the instrument being compared against, not by taste.
    /// `direct` is `a − λ` with `λ` from the `erfcx` route, whose measured
    /// relative accuracy is `~5e-14` (libm `erfc` plus the `exp(x²)` multiply);
    /// its absolute error is therefore `~5e-14·λ`, and the subtraction cannot
    /// remove it. Note how little room that leaves already: the amplification
    /// `λ/|a−λ|` is 18x at `a = 4` and 403x at `a = 20`, so at the top of this
    /// range the direct form is down to ~11 correct digits — five short — while
    /// the continued fraction still matches a 250-digit reference to 16. This
    /// test is deliberately capped at `a = 20` for that reason; it is the last
    /// place the two CAN be compared.
    #[test]
    fn boundary_factor_active_branch_agrees_with_the_direct_form_where_both_are_valid() {
        const LAMBDA_REL_ACCURACY: f64 = 5.0e-14;
        for &a in &[4.0_f64, 4.5, 6.0, 9.0, 20.0] {
            for &b in &[0.0_f64, 0.25, 1.5] {
                if b - a > LEFT_CONTINUED_FRACTION_SWITCH {
                    continue; // not on the continued-fraction branch
                }
                let (_, d_a, _) = cone_boundary_log_factor_and_derivatives(a, b);
                let (_, mills) = signed_probit_logcdf_and_mills_ratio(b - a);
                let direct = a - mills;
                assert!(
                    (d_a - direct).abs() <= LAMBDA_REL_ACCURACY * mills,
                    "(a={a}, b={b}): continued fraction {d_a} vs direct {direct} \
                     (gap {:.3e}, budget {:.3e})",
                    (d_a - direct).abs(),
                    LAMBDA_REL_ACCURACY * mills
                );
            }
        }
    }

    /// A gradient off the domain must not read as usable next to a NaN value.
    #[test]
    fn boundary_factor_derivatives_are_nan_off_the_domain() {
        for &(a, b) in &[
            (-1.0_f64, 0.0_f64),
            (1.0, -1.0),
            (f64::NAN, 1.0),
            (f64::INFINITY, 0.0),
        ] {
            let (v, d_a, d_b) = cone_boundary_log_factor_and_derivatives(a, b);
            assert!(
                v.is_nan() && d_a.is_nan() && d_b.is_nan(),
                "(a={a}, b={b}) is off-domain: got value {v}, ∂a {d_a}, ∂b {d_b}"
            );
        }
    }

    /// Closed-form partials against finite differences of the value
    /// (test-only FD; the production gradient consumes the analytic form).
    /// The domain is `a, b ≥ 0`, so a coordinate sitting exactly on the
    /// boundary uses a one-sided forward difference instead of stepping
    /// outside the domain (where the factor is deliberately NaN).
    #[test]
    fn boundary_factor_derivatives_match_finite_differences() {
        let cases: [(f64, f64); 5] = [(0.3, 0.0), (2.0, 0.5), (0.0, 1.2), (5.0, 0.2), (0.7, 3.0)];
        let step = 1e-6;
        let fd = |lo: f64, mid: f64, hi: f64, coord: f64| -> f64 {
            if coord >= step {
                (hi - lo) / (2.0 * step)
            } else {
                (hi - mid) / step
            }
        };
        for &(a, b) in &cases {
            let (_, d_a, d_b) = cone_boundary_log_factor_and_derivatives(a, b);
            let fd_a = fd(
                cone_boundary_log_factor((a - step).max(0.0), b),
                cone_boundary_log_factor(a, b),
                cone_boundary_log_factor(a + step, b),
                a,
            );
            let fd_b = fd(
                cone_boundary_log_factor(a, (b - step).max(0.0)),
                cone_boundary_log_factor(a, b),
                cone_boundary_log_factor(a, b + step),
                b,
            );
            // One-sided differences on boundary coordinates carry O(step)
            // truncation error, so the band is a few multiples of step.
            assert!(
                (d_a - fd_a).abs() <= 5e-6 * (1.0 + fd_a.abs()),
                "(a={a}, b={b}): ∂a analytic {d_a} vs FD {fd_a}"
            );
            assert!(
                (d_b - fd_b).abs() <= 5e-6 * (1.0 + fd_b.abs()),
                "(a={a}, b={b}): ∂b analytic {d_b} vs FD {fd_b}"
            );
        }
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    const TOL: f64 = 1e-12;

    fn rel_err(got: f64, expected: f64) -> f64 {
        (got - expected).abs() / expected.abs().max(1e-300)
    }

    #[test]
    fn beta_quantile_matches_known_reference_values() {
        let cases: [(f64, f64, f64, f64); 8] = [
            (0.025, 2.0, 2.0, 0.094_299_3),
            (0.975, 2.0, 2.0, 0.905_700_7),
            (0.5, 2.0, 2.0, 0.5),
            (0.025, 0.8, 4.0, 0.002_339_1),
            (0.975, 0.8, 4.0, 0.564_717_3),
            (0.025, 5.0, 1.5, 0.408_549_1),
            (0.5, 20.0, 80.0, 0.197_994_8),
            (0.975, 20.0, 80.0, 0.283_367_6),
        ];
        for (p, a, b, expected) in cases {
            let got = beta_quantile(p, a, b);
            let abs = (got - expected).abs();
            assert!(
                abs < 1e-5,
                "beta_quantile(p={p}, a={a}, b={b}) = {got}, expected ≈ {expected} (abs err {abs})"
            );
        }
    }

    #[test]
    fn beta_quantile_boundaries_and_degeneracy() {
        assert_eq!(beta_quantile(0.0, 2.0, 3.0), 0.0);
        assert_eq!(beta_quantile(-0.5, 2.0, 3.0), 0.0);
        assert_eq!(beta_quantile(1.0, 2.0, 3.0), 1.0);
        assert_eq!(beta_quantile(1.5, 2.0, 3.0), 1.0);
        assert!(beta_quantile(0.5, -1.0, 3.0).is_nan());
        assert!(beta_quantile(0.5, 2.0, 0.0).is_nan());
        assert!(beta_quantile(0.5, f64::NAN, 3.0).is_nan());
        let mut prev = 0.0;
        for i in 1..100 {
            let p = i as f64 / 100.0;
            let q = beta_quantile(p, 3.0, 5.0);
            assert!(q > prev, "beta quantile not increasing at p={p}");
            prev = q;
        }
    }

    // ── normal_pdf ────────────────────────────────────────────────────────────

    #[test]
    fn normal_pdf_at_zero() {
        let expected = 1.0 / (2.0 * std::f64::consts::PI).sqrt();
        assert!((normal_pdf(0.0) - expected).abs() < TOL);
    }

    #[test]
    fn normal_pdf_symmetry() {
        for &x in &[0.5, 1.0, 2.0, 3.0, 5.0] {
            assert_eq!(normal_pdf(x), normal_pdf(-x), "symmetry failed at x={x}");
        }
    }

    /// `x*x` is exact-splittable and the split is what `exp` needs.
    ///
    /// Two independent statements, because the correction is only worth what
    /// its residual is worth. First, `x*x + residual` is `x²` EXACTLY: checked
    /// against a Veltkamp/Dekker split, which reaches the same residual through
    /// pure multiplies and adds and shares no code path with the `mul_add`
    /// route. Second, the residual is not decorative — for these arguments it
    /// is a relative perturbation of `x²` big enough that `exp` amplifies it
    /// past a single ulp of the result.
    #[test]
    fn square_residual_completes_the_rounded_square_exactly() {
        // 2^27 + 1: Veltkamp's splitting factor, exact for any `x` whose
        // scaled form does not overflow.
        const SPLIT: f64 = 134_217_729.0;
        let mut saw_amplified = false;
        for &x in &[
            0.1, 0.7, 1.3, 2.9, 6.1, 10.5, 14.3, 19.7, 23.9, 25.9999, 34.7,
        ] {
            let rounded = x * x;
            let residual = square_residual(x, rounded);

            let c = x * SPLIT;
            let head = c - (c - x);
            let tail = x - head;
            let dekker = ((head * head - rounded) + 2.0 * head * tail) + tail * tail;
            assert_eq!(
                residual, dekker,
                "x={x}: mul_add residual {residual:e} != Dekker residual {dekker:e}"
            );

            // `exp` multiplies a relative argument perturbation by the argument.
            let amplified = (residual / rounded).abs() * rounded;
            if amplified > f64::EPSILON {
                saw_amplified = true;
            }
        }
        assert!(
            saw_amplified,
            "no test argument had a residual `exp` could amplify past one ulp; \
             the correction under test would be untested"
        );
    }

    /// `φ(x)` against an EXTERNAL high-precision reference (mpmath, dps=60).
    ///
    /// Every argument here has an INEXACT square, which is the whole point.
    /// `exp(−½·fl(x*x))` misplaces the argument by `x²·ε/2` RELATIVE, and `exp`
    /// hands that straight back as relative error in the result: `1.4e-14` at
    /// `x ≈ 17`, `5.7e-14` by `x ≈ 35`, where `φ` is still a normal `f64`. Only
    /// the top of the range makes that visible, so the table has to reach it —
    /// a `φ` table that stops at `x = 5` cannot tell the two forms apart.
    ///
    /// `1.5e-15` (≈7 ulp) is the portability allowance: `f64::exp` is the
    /// platform libm and the only part of this that is not fixed by the crate
    /// graph, and it is worth ~1 ulp on the implementations in use. That still
    /// leaves 38x of margin against the defect at the top of the table.
    #[test]
    fn normal_pdf_matches_high_precision_reference() {
        const TOLERANCE: f64 = 1.5e-15;
        let refs: &[(f64, f64)] = &[
            (0.5, 0.35206532676429947),
            (1.0, 0.24197072451914334),
            (2.5, 0.017528300493568537),
            (4.0, 0.00013383022576488534),
            (7.3, 1.0693837871541648e-12),
            (11.9, 7.090702668428078e-32),
            (17.4, 7.201308152719057e-67),
            (23.6, 4.555989824112156e-122),
            (29.1, 5.229437243665329e-185),
            (34.7, 1.368008224488383e-262),
        ];
        for &(x, reference) in refs {
            // The small arguments anchor the ordinary range; the large ones are
            // where the defect lives, and every one of THOSE has to have a
            // square `f64` cannot hold or it exercises nothing.
            assert!(
                x <= 5.0 || square_residual(x, x * x) != 0.0,
                "x={x} squares exactly, so it cannot exercise the correction"
            );
            let rel = rel_err(normal_pdf(x), reference);
            assert!(
                rel < TOLERANCE,
                "normal_pdf({x}) = {:.17e}, reference {reference:.17e}, rel {rel:.3e}",
                normal_pdf(x)
            );
        }
    }

    /// `φ` off the ordinary domain, where the square has no usable residual:
    /// `±∞` squares to `∞` and would hand the correction an `∞ − ∞`.
    #[test]
    fn normal_pdf_nonfinite_and_underflowed_arguments() {
        assert_eq!(normal_pdf(f64::INFINITY), 0.0);
        assert_eq!(normal_pdf(f64::NEG_INFINITY), 0.0);
        assert!(normal_pdf(f64::NAN).is_nan());
        // Past ~38.6 the pdf underflows; it must reach zero, not NaN.
        assert_eq!(normal_pdf(40.0), 0.0);
        assert_eq!(normal_pdf(-40.0), 0.0);
        assert_eq!(normal_pdf(f64::MAX), 0.0);
        // Just inside the underflow edge the result is subnormal but positive.
        let edge = normal_pdf(38.0);
        assert!(edge > 0.0 && edge.is_subnormal(), "phi(38) = {edge:e}");
    }

    #[test]
    fn normal_pdf_positive() {
        for &x in &[-5.0, -1.0, 0.0, 1.0, 5.0] {
            assert!(normal_pdf(x) > 0.0, "pdf should be positive at x={x}");
        }
    }

    // ── normal_cdf ────────────────────────────────────────────────────────────

    #[test]
    fn normal_cdf_at_zero_is_half() {
        assert!((normal_cdf(0.0) - 0.5).abs() < TOL);
    }

    #[test]
    fn normal_cdf_symmetry() {
        for &x in &[0.5, 1.0, 2.0, 3.0] {
            let sum = normal_cdf(x) + normal_cdf(-x);
            assert!(
                (sum - 1.0).abs() < TOL,
                "cdf symmetry failed at x={x}: sum={sum}"
            );
        }
    }

    #[test]
    fn normal_cdf_bounds() {
        assert!(normal_cdf(10.0) > 0.9999);
        assert!(normal_cdf(-10.0) < 1e-22);
        assert!(normal_cdf(0.0) > 0.0);
        assert!(normal_cdf(0.0) < 1.0);
    }

    #[test]
    fn normal_cdf_at_1_96_near_0975() {
        // Phi(1.96) ≈ 0.975 — canonical two-sided 5% critical value.
        let p = normal_cdf(1.959_963_985);
        assert!((p - 0.975).abs() < 1e-8, "p={p}");
    }

    // ── erfcx_nonnegative ─────────────────────────────────────────────────────

    #[test]
    fn erfcx_zero_is_one_and_negative_domain_is_rejected() {
        assert_eq!(erfcx_nonnegative(0.0), 1.0);
        assert!(erfcx_nonnegative(-f64::MIN_POSITIVE).is_nan());
        assert!(erfcx_nonnegative(-1.0).is_nan());
        assert!(erfcx_nonnegative(f64::NEG_INFINITY).is_nan());
    }

    #[test]
    fn erfcx_positive_inf_returns_zero() {
        assert_eq!(erfcx_nonnegative(f64::INFINITY), 0.0);
    }

    #[test]
    fn erfcx_nan_propagates() {
        assert!(erfcx_nonnegative(f64::NAN).is_nan());
    }

    #[test]
    fn erfcx_small_positive_matches_direct() {
        use libm::erfc;
        for &x in &[0.1_f64, 0.5, 1.0, 5.0, 10.0, 25.0] {
            let got = erfcx_nonnegative(x);
            let expected = (x * x).exp() * erfc(x);
            let err = rel_err(got, expected);
            assert!(
                err < 1e-10,
                "x={x}: got={got} expected={expected} rel={err}"
            );
        }
    }

    #[test]
    fn erfcx_large_x_positive_and_finite() {
        // For x >= 26 the asymptotic branch must remain positive and finite.
        let got = erfcx_nonnegative(50.0);
        assert!(got.is_finite() && got > 0.0, "erfcx(50)={got}");
        // Leading asymptotic term: 1/(x*sqrt(pi)).
        let asymptotic = 1.0 / (50.0 * std::f64::consts::PI.sqrt());
        assert!(
            rel_err(got, asymptotic) < 1e-3,
            "got={got} asymptotic={asymptotic}"
        );
    }

    /// The two branches must describe one function across `x = 26`.
    ///
    /// Note WHY the plain `exp(x*x)·erfc(x)` below is a legitimate oracle at
    /// this particular argument and nowhere else: `26² = 676` is exactly
    /// representable, so the rounded square carries no residual and the direct
    /// form is momentarily as good as the corrected one. That is also exactly
    /// why this check was blind to the `x²·ε/2` defect it looks like it should
    /// have caught — at `25.9` the same comparison would have failed by
    /// `5.7e-14`, but the seam was only ever probed at the one point in the
    /// neighbourhood where the defect vanishes. The bit-adjacent step below
    /// cannot substitute for it either: `d(ln erfcx)/dx ≈ −2x` at the switch,
    /// so one ulp of `x` moves the true value by `1.8e-13`, three times the
    /// defect. It takes a reference at a DISTANCE from the seam — the table in
    /// `erfcx_matches_high_precision_reference` — to see the defect at all.
    #[test]
    fn erfcx_asymptotic_switch_matches_finite_direct_identity() {
        let switch = 26.0_f64;
        assert_eq!(
            square_residual(switch, switch * switch),
            0.0,
            "676 must be exact for the direct form below to be an oracle"
        );
        let direct = (switch * switch).exp() * erfc(switch);
        let asymptotic = erfcx_nonnegative(switch);
        assert!(
            rel_err(asymptotic, direct) < 1.0e-15,
            "switch mismatch: asymptotic={asymptotic:.17e}, direct={direct:.17e}"
        );

        // Continuity across the branch cut, up to how fast the function itself
        // moves over one ulp of `x` (`|d ln erfcx/dx| ≈ 2x` ⇒ ~1.9e-13 here).
        let immediately_below = f64::from_bits(switch.to_bits() - 1);
        let below = erfcx_nonnegative(immediately_below);
        let step = 2.0 * switch * (switch - immediately_below);
        assert!(
            rel_err(asymptotic, below) < 2.0 * step,
            "discontinuous switch: below={below:.17e}, at={asymptotic:.17e}, \
             one-ulp travel {step:.3e}"
        );
    }

    #[test]
    fn erfcx_preserves_representable_subnormal_tail() {
        let tail = erfcx_nonnegative(f64::MAX);
        assert!(tail > 0.0 && tail.is_subnormal(), "erfcx(MAX)={tail:e}");
    }

    /// Absolute-accuracy pin against an EXTERNAL high-precision reference
    /// (mpmath, dps=60) spanning the direct branch `[0.1, 26)`. This is the
    /// root-cause guard: the previous `exp(x²)·erfc(x)` direct form was built on
    /// `statrs::erfc`, whose ~1e-10 relative accuracy silently poisoned every
    /// downstream probit / Mills / log-CDF derivative.
    ///
    /// The table had a SECOND job it was not doing. Of its twelve arguments,
    /// eleven — `0.5`, `2`, `3.5`, `6`, `9`, `13`, `18`, `22`, `25.5`, and the
    /// two whose squares are far too small to matter — square EXACTLY in `f64`,
    /// so `fl(x*x) = x²` and the `x²·ε/2` error the rounded square feeds `exp`
    /// was identically zero at every one of them. The twelfth, `25.9999`, does
    /// not square exactly; it was the one point in the table where the defect
    /// was live, and its literal had been recorded WITH the defect in it —
    /// `0.021683668126370212` against a true `0.021683668126369115`, off by
    /// `5.1e-14`. Three independent high-precision routes (`exp(x²)·erfc(x)`,
    /// the 12-term asymptotic series, and a 400-level Laplace continued
    /// fraction) and `scipy.special.erfcx` all agree on the corrected value.
    /// A `1e-13` tolerance then accepted a reference that was itself wrong by
    /// half the tolerance, which is how a 190x accuracy defect sat under a
    /// test named for high precision.
    ///
    /// So the table now RUNS ON arguments with inexact squares (`10.5`,
    /// `14.3`, `19.7`, `23.9` alongside the original grid) and the tolerance is
    /// `1.5e-15` — 38x below the defect at the top of the range, and still ~7
    /// ulp of headroom for the platform `f64::exp` (the only part of this path
    /// not pinned by the crate graph; `erfc` comes from the `libm` crate and is
    /// identical everywhere).
    #[test]
    fn erfcx_matches_high_precision_reference() {
        const TOLERANCE: f64 = 1.5e-15;
        // (x, mpmath exp(x²)·erfc(x) at dps=60, rounded to f64).
        let refs: &[(f64, f64)] = &[
            (0.1, 0.8964569799691267),
            (0.5, 0.6156903441929259),
            (1.0, 0.427583576155807),
            (2.0, 0.25539567631050575),
            (3.5, 0.1552936556088943),
            (6.0, 0.09277656780053835),
            (9.0, 0.06230772403777468),
            (10.5, 0.05349189974656412),
            (13.0, 0.043271921864609694),
            (14.3, 0.0393580473372741),
            (18.0, 0.03129571781590521),
            (19.7, 0.028602309402825203),
            (22.0, 0.025618570005879453),
            (23.9, 0.023585649371803793),
            (25.5, 0.022108108052519827),
            (25.9999, 0.021683668126369115),
        ];
        for &(x, reference) in refs {
            let got = erfcx_nonnegative(x);
            let rel = rel_err(got, reference);
            assert!(
                rel < TOLERANCE,
                "erfcx({x}) = {got:.17e}, reference {reference:.17e}, rel {rel:.3e}"
            );
        }
        // The point of the added arguments: at least four of them must have a
        // square `f64` cannot hold, or the table is back to testing nothing.
        let inexact = refs
            .iter()
            .filter(|&&(x, _)| square_residual(x, x * x) != 0.0)
            .count();
        assert!(
            inexact >= 4,
            "only {inexact} of {} reference arguments have an inexact square",
            refs.len()
        );
    }

    // ── log1mexp_positive ─────────────────────────────────────────────────────

    #[test]
    fn log1mexp_at_zero_is_neg_inf() {
        assert_eq!(log1mexp_positive(0.0), f64::NEG_INFINITY);
    }

    #[test]
    fn log1mexp_recovers_log_one_minus_exp() {
        // Verify exp(log1mexp(a)) + exp(-a) ≈ 1 for several a > 0. This
        // roundtrip avoids computing `(1 - exp(-a)).ln()` directly, which
        // suffers catastrophic cancellation for large a (e.g. a=20 where
        // `1.0 - exp(-20)` loses 9 decimal digits from the subtraction).
        for &a in &[0.001_f64, 0.5, std::f64::consts::LN_2, 1.0, 5.0, 20.0] {
            let lm = log1mexp_positive(a);
            let roundtrip = lm.exp() + (-a).exp();
            assert!(
                (roundtrip - 1.0).abs() < 1e-14,
                "a={a}: exp(log1mexp(a)) + exp(-a) = {roundtrip}, expected 1.0"
            );
        }
    }

    #[test]
    fn log1mexp_at_ln2_is_neg_ln2() {
        let ln2 = std::f64::consts::LN_2;
        let got = log1mexp_positive(ln2);
        assert!((got - (-ln2)).abs() < TOL, "got={got}");
    }

    // ── signed_log_sum_exp ────────────────────────────────────────────────────

    #[test]
    fn slse_all_positive_single() {
        let (lm, sg) = signed_log_sum_exp(&[2.0], &[1.0]);
        assert!((lm - 2.0).abs() < TOL);
        assert!((sg - 1.0).abs() < TOL);
    }

    #[test]
    fn slse_difference_recovers_log2() {
        // 3 - 1 = 2 → log|2| = ln(2), sign = +1.
        let log3 = 3.0_f64.ln();
        let log1 = 0.0_f64; // ln(1)
        let (lm, sg) = signed_log_sum_exp(&[log3, log1], &[1.0, -1.0]);
        assert!((lm - 2.0_f64.ln()).abs() < TOL, "lm={lm}");
        assert!((sg - 1.0).abs() < TOL, "sg={sg}");
    }

    #[test]
    fn slse_cancellation_gives_neg_inf() {
        // a - a = 0 → log|0| = -∞.
        let ln2 = 2.0_f64.ln();
        let (lm, sg) = signed_log_sum_exp(&[ln2, ln2], &[1.0, -1.0]);
        assert_eq!(lm, f64::NEG_INFINITY);
        assert_eq!(sg, 0.0);
    }

    #[test]
    fn slse_empty_returns_neg_inf_with_zero_sign() {
        // With no terms the sum is exactly 0, so the docstring contract is
        // `(−∞, 0.0)`. (This test previously encoded the buggy `+1.0` positive-sum
        // convention, which contradicted both the docstring and the cancellation
        // test below; rewritten to the correct zero sign.)
        let (lm, sg) = signed_log_sum_exp(&[], &[]);
        assert_eq!(lm, f64::NEG_INFINITY);
        assert_eq!(sg, 0.0);
    }

    #[test]
    fn slse_all_zero_signs_return_zero_sign() {
        // A single term whose sign is 0 contributes nothing; S = 0 ⇒ (−∞, 0.0).
        let (lm, sg) = signed_log_sum_exp(&[0.0], &[0.0]);
        assert_eq!(lm, f64::NEG_INFINITY);
        assert_eq!(sg, 0.0);
    }

    #[test]
    fn slse_all_neg_inf_magnitudes_return_zero_sign() {
        // Every magnitude is exp(−∞) = 0 regardless of sign, so the sum is 0 and
        // the reported sign must be 0.0, not +1.0.
        let (lm, sg) = signed_log_sum_exp(&[f64::NEG_INFINITY, f64::NEG_INFINITY], &[1.0, -1.0]);
        assert_eq!(lm, f64::NEG_INFINITY);
        assert_eq!(sg, 0.0);
    }

    #[test]
    fn slse_pos_inf_dominates() {
        let (lm, sg) = signed_log_sum_exp(&[f64::INFINITY, 1.0], &[1.0, -1.0]);
        assert_eq!(lm, f64::INFINITY);
        assert_eq!(sg, 1.0);
    }

    #[test]
    fn slse_neg_inf_dominates() {
        let (lm, sg) = signed_log_sum_exp(&[f64::INFINITY, 1.0], &[-1.0, 1.0]);
        assert_eq!(lm, f64::INFINITY);
        assert_eq!(sg, -1.0);
    }

    #[test]
    fn slse_both_inf_signs_gives_nan() {
        let (lm, sg) = signed_log_sum_exp(&[f64::INFINITY, f64::INFINITY], &[1.0, -1.0]);
        assert!(lm.is_nan());
        assert_eq!(sg, 0.0);
    }

    // ── normal_logcdf ─────────────────────────────────────────────────────────

    #[test]
    fn logcdf_at_zero_is_log_half() {
        let got = normal_logcdf(0.0);
        let expected = 0.5_f64.ln();
        assert!((got - expected).abs() < TOL, "got={got}");
    }

    #[test]
    fn logcdf_pos_inf_is_zero() {
        assert_eq!(normal_logcdf(f64::INFINITY), 0.0);
    }

    #[test]
    fn logcdf_neg_inf_is_neg_inf() {
        assert_eq!(normal_logcdf(f64::NEG_INFINITY), f64::NEG_INFINITY);
    }

    #[test]
    fn logcdf_nan_is_nan() {
        assert!(normal_logcdf(f64::NAN).is_nan());
    }

    #[test]
    fn logcdf_matches_log_cdf_for_moderate_x() {
        for &x in &[-2.0_f64, -1.0, 0.0, 1.0, 2.0, 3.0] {
            let got = normal_logcdf(x);
            let expected = normal_cdf(x).ln();
            assert!(
                (got - expected).abs() < 1e-10,
                "x={x}: got={got} expected={expected}"
            );
        }
    }

    #[test]
    fn logcdf_deep_left_tail_stays_finite() {
        // For very negative x, normal_cdf(x) underflows to 0, but logcdf should
        // remain finite and large-negative.
        let got = normal_logcdf(-20.0);
        assert!(got.is_finite() && got < -100.0, "logcdf(-20)={got}");
    }

    #[test]
    fn logcdf_positive_tail_does_not_round_through_unit_cdf() {
        let x = 10.0_f64;
        let got = normal_logcdf(x);
        let expected = (-0.5 * erfc(x / std::f64::consts::SQRT_2)).ln_1p();
        assert!(
            got < 0.0,
            "logcdf(10) must retain its negative tail: {got:e}"
        );
        assert_eq!(got.to_bits(), expected.to_bits());
    }

    #[test]
    fn log_cdf_quantile_round_trips_both_unrepresentable_tails() {
        for x in [-1.0e6, -40.0, -10.0, -2.0, 0.0, 2.0, 10.0] {
            let log_p = normal_logcdf(x);
            let recovered = standard_normal_quantile_from_log_cdf(log_p)
                .expect("finite strict log-CDF has a quantile");
            assert!(
                (recovered - x).abs() <= 2.0e-12 * x.abs().max(1.0),
                "log-quantile round trip at x={x}: log_p={log_p}, recovered={recovered}"
            );
        }
    }

    // ── normal_logsf ─────────────────────────────────────────────────────────

    #[test]
    fn logsf_at_zero_is_log_half() {
        let got = normal_logsf(0.0);
        let expected = 0.5_f64.ln();
        assert!((got - expected).abs() < TOL, "got={got}");
    }

    #[test]
    fn logsf_mirrors_logcdf() {
        // logsf(x) = logcdf(-x) by definition.
        for &x in &[-3.0_f64, -1.0, 0.0, 1.0, 3.0] {
            assert_eq!(normal_logsf(x), normal_logcdf(-x));
        }
    }

    // ── signed_probit_logcdf_and_mills_ratio ──────────────────────────────────

    #[test]
    fn probit_at_pos_inf() {
        let (lc, mr) = signed_probit_logcdf_and_mills_ratio(f64::INFINITY);
        assert_eq!(lc, 0.0);
        assert_eq!(mr, 0.0);
    }

    #[test]
    fn probit_at_neg_inf() {
        let (lc, mr) = signed_probit_logcdf_and_mills_ratio(f64::NEG_INFINITY);
        assert_eq!(lc, f64::NEG_INFINITY);
        assert_eq!(mr, f64::INFINITY);
    }

    #[test]
    fn probit_nan_propagates() {
        let (lc, mr) = signed_probit_logcdf_and_mills_ratio(f64::NAN);
        assert!(lc.is_nan() && mr.is_nan());
    }

    #[test]
    fn probit_at_zero_logcdf_and_mills() {
        let (lc, mr) = signed_probit_logcdf_and_mills_ratio(0.0);
        assert!((lc - 0.5_f64.ln()).abs() < TOL, "lc={lc}");
        // phi(0)/Phi(0) = 0.3989.../0.5 ≈ 0.7979.
        assert!((mr - 0.797_884_560_802_865).abs() < 1e-10, "mr={mr}");
    }

    #[test]
    fn probit_positive_branch_matches_logcdf() {
        for &x in &[0.5_f64, 1.0, 2.0, 3.0] {
            let (lc, mr) = signed_probit_logcdf_and_mills_ratio(x);
            let lc_ref = normal_logcdf(x);
            let mr_ref = normal_pdf(x) / normal_cdf(x);
            assert!(
                (lc - lc_ref).abs() < 1e-10,
                "x={x}: lc={lc} lc_ref={lc_ref}"
            );
            assert!(
                (mr - mr_ref).abs() < 1e-10,
                "x={x}: mr={mr} mr_ref={mr_ref}"
            );
        }
    }

    #[test]
    fn probit_negative_branch_matches_logcdf() {
        for &x in &[-0.5_f64, -1.0, -2.0, -5.0] {
            let (lc, mr) = signed_probit_logcdf_and_mills_ratio(x);
            let lc_ref = normal_logcdf(x);
            assert!(
                (lc - lc_ref).abs() < 1e-10,
                "x={x}: lc={lc} lc_ref={lc_ref}"
            );
            assert!(mr.is_finite() && mr > 0.0, "x={x}: mr={mr}");
        }
    }

    #[test]
    fn probit_mills_ratio_has_no_deep_tail_floor() {
        let x = -1.0e305_f64;
        let (log_cdf, mills_ratio) = signed_probit_logcdf_and_mills_ratio(x);
        assert_eq!(log_cdf, f64::NEG_INFINITY);
        assert!(mills_ratio.is_finite());
        assert!(
            ((mills_ratio / -x) - 1.0).abs() < 5.0e-15,
            "mills({x:e})={mills_ratio:e}"
        );
    }

    #[test]
    fn normal_logcdf_derivative_stack_has_honest_infinite_limits() {
        assert_eq!(normal_logcdf_derivatives(f64::INFINITY), [0.0; 5]);
        assert_eq!(
            normal_logcdf_derivatives(f64::NEG_INFINITY),
            [f64::NEG_INFINITY, f64::INFINITY, -1.0, 0.0, 0.0]
        );
        assert!(
            normal_logcdf_derivatives(f64::NAN)
                .into_iter()
                .all(f64::is_nan)
        );

        for x in [-1.0e200_f64, 1.0e200_f64] {
            let derivatives = normal_logcdf_derivatives(x);
            assert!(
                derivatives.into_iter().all(|value| !value.is_nan()),
                "NaN derivative at x={x:e}: {derivatives:?}"
            );
        }
    }

    #[test]
    fn normal_logcdf_left_tail_derivatives_do_not_cancel() {
        let x = -1.0e100_f64;
        let derivatives = normal_logcdf_derivatives(x);
        assert_eq!(derivatives[2], -1.0);
        assert!(derivatives[3] > 0.0 && derivatives[3].is_finite());
        assert!(
            (derivatives[3] / 2.0e-300 - 1.0).abs() < 2.0e-14,
            "third derivative={:e}",
            derivatives[3]
        );
        assert_eq!(derivatives[4], 0.0);
    }

    #[test]
    fn normal_logcdf_right_tail_preserves_weighted_subnormal_derivatives() {
        let derivatives = normal_logcdf_derivatives(38.6);
        assert_eq!(derivatives[1], 0.0);
        assert!(derivatives[2] < 0.0 && derivatives[2].is_subnormal());
        assert!(derivatives[3] > 0.0 && derivatives[3].is_subnormal());
        assert!(derivatives[4] < 0.0 && derivatives[4].is_subnormal());
    }

    #[test]
    fn normal_logcdf_tail_stack_is_finite_difference_consistent() {
        let h = 1.0e-4_f64;
        for x in [-8.0_f64, -4.0, 8.0, 20.0] {
            let center = normal_logcdf_derivatives(x);
            let left = normal_logcdf_derivatives(x - h);
            let right = normal_logcdf_derivatives(x + h);
            for order in 1..=3 {
                let finite_difference = (right[order] - left[order]) / (2.0 * h);
                let expected = center[order + 1];
                let relative = (finite_difference - expected).abs() / expected.abs().max(1.0e-300);
                assert!(
                    relative < 2.0e-5,
                    "x={x}, order={order}: fd={finite_difference:e}, expected={expected:e}, rel={relative:e}"
                );
            }
        }
    }

    /// Absolute-accuracy pin of the full `ln Φ(x)` derivative tower against an
    /// EXTERNAL high-precision reference (mpmath, dps=60), covering all three
    /// branches (continued-fraction left tail at x=−4, the moderate Mills
    /// recurrence for x∈(−4, 8), and both signs). Before the `erfc` root-cause
    /// fix the moderate branch's `λ = φ/Φ` inherited `statrs::erfc`'s ~1e-10
    /// error, so `f''` was wrong by ~1e-9 near the −4 seam; this pins every
    /// entry to `2e-11` relative, catching that regression head-on rather than
    /// through a seam-straddling finite difference.
    #[test]
    fn normal_logcdf_derivative_tower_matches_high_precision_reference() {
        // (x, [value, f', f'', f''', f''''] from mpmath at dps=60).
        let refs: &[(f64, [f64; 5])] = &[
            (
                -4.0,
                [
                    -10.360101486527291,
                    4.2256071444894711,
                    -0.95332716160257737,
                    0.017856339307658426,
                    0.0095065764315958691,
                ],
            ),
            // Two points well inside the continued-fraction branch, where the
            // truncation the depth controls is the ONLY error source: at -4 the
            // branch is at its least converged, and these confirm it stays put.
            (
                -10.0,
                [
                    -53.231285150512471,
                    10.098093233962512,
                    -0.99055462217434374,
                    0.0017864003921165069,
                    0.00049785382237944016,
                ],
            ),
            (
                -6.0,
                [
                    -20.736768949974706,
                    6.1584826045445989,
                    -0.97601236321083323,
                    0.0069535374991643118,
                    0.0028992056785575027,
                ],
            ),
            (
                -2.0,
                [
                    -3.7831843336820319,
                    2.3732155328228409,
                    -0.88572089958591874,
                    0.059355861291565813,
                    0.039421993865946813,
                ],
            ),
            (
                -1.0,
                [
                    -1.8410216450092635,
                    1.5251352761609812,
                    -0.80090233442965121,
                    0.11693119540604883,
                    0.07917498368074563,
                ],
            ),
            (
                -0.3,
                [
                    -0.96210281816885066,
                    0.99816596885848332,
                    -0.69688551072964971,
                    0.18398317992442132,
                    0.11037564722092704,
                ],
            ),
            (
                0.5,
                [
                    -0.36894641528865639,
                    0.50916043383703349,
                    -0.5138245643036329,
                    0.27099012446870783,
                    0.088167801929197554,
                ],
            ),
            (
                2.0,
                [
                    -0.023012909328963488,
                    0.055247862678989959,
                    -0.11354805168857645,
                    0.18439481503247759,
                    -0.18785468561160969,
                ],
            ),
        ];
        // The moderate-branch statrs regression produced ~1e-9 errors in f''.
        // The bound used to sit at 1e-10 to respect what was called the
        // continued-fraction branch's "inherent" ~2e-11 in f''''; that was not
        // inherent but a depth, and at 64 levels the branch reproduces this
        // 60-digit reference EXACTLY at x = -4, -6 and -10. What remains is the
        // moderate branch, where the brackets are already collected in `q` and
        // the floor is `λ`'s own relative error amplified by `λ/q` (18.7 at the
        // switch): 1.8e-13 at x = -2, the worst point here. 1e-11 keeps 55x of
        // headroom over that while still failing the 32-level truncation head-on.
        for &(x, reference) in refs {
            let got = normal_logcdf_derivatives(x);
            for (order, (&g, &r)) in got.iter().zip(reference.iter()).enumerate() {
                let rel = (g - r).abs() / r.abs().max(1.0e-3);
                assert!(
                    rel < 1.0e-11,
                    "normal_logcdf_derivatives({x})[{order}] = {g:.17e}, reference {r:.17e}, \
                     rel {rel:.3e} >= 1e-11"
                );
            }
        }
    }

    // ── standard_normal_quantile ──────────────────────────────────────────────

    #[test]
    fn quantile_rejects_out_of_range() {
        assert!(standard_normal_quantile(0.0).is_err());
        assert!(standard_normal_quantile(1.0).is_err());
        assert!(standard_normal_quantile(-0.1).is_err());
        assert!(standard_normal_quantile(1.1).is_err());
        assert!(standard_normal_quantile(f64::NAN).is_err());
    }

    #[test]
    fn quantile_at_half_is_near_zero() {
        let q = standard_normal_quantile(0.5).unwrap();
        assert!(q.abs() < 1e-10, "quantile(0.5)={q}");
    }

    #[test]
    fn quantile_at_0975_is_near_196() {
        let q = standard_normal_quantile(0.975).unwrap();
        assert!((q - 1.959_963_984_540_054).abs() < 1e-14, "q={q}");
    }

    /// `standard_normal_quantile` and its log-CDF sibling, against a 120-digit
    /// root of `Φ(x) = p` (respectively `ln Φ(x) = log_p`).
    ///
    /// The seed is Acklam's rational approximation, whose accuracy is `1.15e-9`
    /// relative; the two Halley steps after it are what make the result
    /// ulp-accurate. Deleting the polish loop entirely leaves EVERY other
    /// quantile test in this module green except `quantile_roundtrip_cdf`, and
    /// that one only by a factor of 1.9 — so the polish had no real gate. This
    /// table is that gate: it fails by six orders if the seed ships unpolished.
    ///
    /// The grid straddles Acklam's own `P_LOW = 0.02425` branch on both sides,
    /// runs out to `p = 1e-300` where the seed is far from the root, and covers
    /// the reflected upper tail where the residual must be formed from
    /// `(1 − p) − ½erfc(x/√2)` rather than `Φ(x) − p`.
    #[test]
    fn normal_quantiles_match_independent_high_precision_reference() {
        const QUANTILE_REFERENCE: [[f64; 2]; 22] = [
            [1e-300, -37.0470962993612],
            [1e-100, -21.273453560965326],
            [1e-20, -9.262340089798407],
            [1e-08, -5.612001244174789],
            [0.001, -3.0902323061678136],
            [0.02424, -1.9731366119445441],
            [0.02425, -1.972961051311885],
            [0.02426, -1.9727855514678605],
            [0.05, -1.6448536269514726],
            [0.1, -1.2815515655446004],
            [0.25, -0.6744897501960817],
            [0.4, -0.2533471031357997],
            [0.5, 0.0],
            [0.6, 0.2533471031357997],
            [0.75, 0.6744897501960817],
            [0.9, 1.2815515655446006],
            [0.95, 1.6448536269514722],
            [0.975, 1.9599639845400538],
            [0.99, 2.3263478740408408],
            [0.999, 3.090232306167813],
            [0.99999999, 5.612001243305505],
            [0.9999999999999999, 8.209536151601387],
        ];
        for [p, want] in QUANTILE_REFERENCE {
            let got = standard_normal_quantile(p).expect("p in (0,1) has a quantile");
            let error = (got - want).abs();
            // `Φ⁻¹(½) = 0` exactly, so it is the one absolute comparison.
            let budget = if want == 0.0 {
                1e-16
            } else {
                4e-15 * want.abs()
            };
            assert!(
                error <= budget,
                "Φ⁻¹({p}): got {got:.17e}, want {want:.17e} (error {error:.3e} > {budget:.3e})"
            );
        }

        const LOG_CDF_QUANTILE_REFERENCE: [[f64; 2]; 9] = [
            [-0.7, -0.008559478582480282],
            [-2.0, -1.1015196284987503],
            [-10.0, -3.913946240531893],
            [-50.0, -9.674825283612357],
            [-200.0, -19.803669380301212],
            [-1000.0, -44.6157477319694],
            [-10000.0, -141.37983987312717],
            [-100000.0, -447.1978936785251],
            [-1000000.0, -1414.2077829910174],
        ];
        for [log_p, want] in LOG_CDF_QUANTILE_REFERENCE {
            let got =
                standard_normal_quantile_from_log_cdf(log_p).expect("finite log_p < 0 has a root");
            let error = (got - want).abs();
            // Rounding `log_p` itself to `f64` already moves the root by
            // `ulp(log_p)·dx/d(log_p)`, and `dx/d(log_p) = Φ/φ = 1/λ` — about
            // `1.25` near `p = ½` and `≈ 1/|x|` in the deep tail. That input
            // conditioning, not the solver, is what limits `log_p = −0.7`,
            // where the root sits at `−0.00856` and one ulp of `0.7` is already
            // `1.4e-16` of it.
            let conditioning = 8.0 * f64::EPSILON * log_p.abs() / want.abs().max(0.8);
            let budget = 4e-15 * want.abs() + conditioning;
            assert!(
                error <= budget,
                "Φ⁻¹(exp({log_p})): got {got:.17e}, want {want:.17e} \
                 (error {error:.3e} > {budget:.3e})"
            );
        }
    }

    #[test]
    fn quantile_antisymmetry() {
        let q_lo = standard_normal_quantile(0.1).unwrap();
        let q_hi = standard_normal_quantile(0.9).unwrap();
        assert!((q_lo + q_hi).abs() < 1e-10, "q_lo={q_lo} q_hi={q_hi}");
    }

    #[test]
    fn quantile_roundtrip_cdf() {
        for &p in &[
            0.001, 0.01, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.99, 0.999,
        ] {
            let q = standard_normal_quantile(p).unwrap();
            let p_back = normal_cdf(q);
            // RELATIVE, and sized by what the round trip can cost: a few ulp of
            // `q` propagated through `φ(q)`, plus a couple of ulp from `erfc`
            // itself. The former absolute `1e-10` bar was two orders looser than
            // an unpolished Acklam seed at its worst point.
            assert!(
                (p_back - p).abs() <= 1e-14 * p,
                "roundtrip failed at p={p}: q={q} p_back={p_back}"
            );
        }
    }
}