prospicio-math 0.0.4

//! Special functions.

use std::f64::consts::{FRAC_1_SQRT_2, LN_2, PI};

/// Standard normal density.
pub fn norm_pdf(x: f64) -> f64 {
    (-0.5 * x * x).exp() / (2.0 * PI).sqrt()
}

/// Standard normal distribution function, accurate in both tails.
pub fn norm_cdf(x: f64) -> f64 {
    0.5 * libm::erfc(-x * FRAC_1_SQRT_2)
}

/// Natural log of the gamma function, `ln Γ(x)` for `x > 0`.
///
/// ```
/// use prospicio_math::special::ln_gamma;
///
/// // Γ(5) = 4! = 24.
/// assert!((ln_gamma(5.0) - 24f64.ln()).abs() < 1e-14);
/// ```
pub fn ln_gamma(x: f64) -> f64 {
    libm::lgamma(x)
}

/// Digamma function `ψ(x) = d ln Γ(x) / dx` for `x > 0`: the recurrence
/// `ψ(x) = ψ(x + 1) - 1/x` up to `x ≥ 10`, then the asymptotic series.
///
/// ```
/// use prospicio_math::special::digamma;
///
/// // ψ(1) = -γ (Euler's constant).
/// assert!((digamma(1.0) + 0.5772156649015329).abs() < 1e-14);
/// ```
pub fn digamma(x: f64) -> f64 {
    if x.is_nan() || x <= 0.0 {
        return f64::NAN;
    }
    let mut x = x;
    let mut acc = 0.0;
    while x < 10.0 {
        acc -= 1.0 / x;
        x += 1.0;
    }
    let f = 1.0 / (x * x);
    let series = f
        * (-1.0 / 12.0
            + f * (1.0 / 120.0
                + f * (-1.0 / 252.0
                    + f * (1.0 / 240.0 + f * (-1.0 / 132.0 + f * 691.0 / 32760.0)))));
    acc + x.ln() - 0.5 / x + series
}

/// Trigamma function `ψ'(x)`, the derivative of [`digamma`], for `x > 0`:
/// the recurrence `ψ'(x) = ψ'(x + 1) + 1/x²` up to 20, then the
/// asymptotic series, whose first omitted term is below `1e-15` there.
///
/// ```
/// use prospicio_math::special::trigamma;
///
/// // ψ'(1) = π²/6.
/// assert!((trigamma(1.0) - std::f64::consts::PI.powi(2) / 6.0).abs() < 1e-14);
/// ```
pub fn trigamma(x: f64) -> f64 {
    if x.is_nan() || x <= 0.0 {
        return f64::NAN;
    }
    let mut x = x;
    let mut acc = 0.0;
    while x < 20.0 {
        acc += 1.0 / (x * x);
        x += 1.0;
    }
    let f = 1.0 / (x * x);
    let series = f / x * (1.0 / 6.0 + f * (-1.0 / 30.0 + f * (1.0 / 42.0 + f * (-1.0 / 30.0))));
    acc + 1.0 / x + 0.5 * f + series
}

/// Standard normal quantile (inverse of [`norm_cdf`]).
///
/// Returns `-inf` at 0, `+inf` at 1 and NaN outside `[0, 1]`. Starts from
/// Abramowitz & Stegun 26.2.23 (error below 4.5e-4) and polishes with Newton
/// steps on `ln Phi`, which keeps full relative accuracy deep in the tail.
///
/// ```
/// use prospicio_math::special::norm_quantile;
///
/// assert!((norm_quantile(0.975) - 1.959963984540054).abs() < 1e-14);
/// ```
pub fn norm_quantile(p: f64) -> f64 {
    if p.is_nan() || !(0.0..=1.0).contains(&p) {
        return f64::NAN;
    }
    if p == 0.0 {
        return f64::NEG_INFINITY;
    }
    if p == 1.0 {
        return f64::INFINITY;
    }
    if p > 0.5 {
        return -lower_quantile(1.0 - p);
    }
    lower_quantile(p)
}

/// Quantile for `0 < p <= 0.5`.
fn lower_quantile(p: f64) -> f64 {
    let t = (-2.0 * p.ln()).sqrt();
    let mut x = -(t
        - (2.515517 + 0.802853 * t + 0.010328 * t * t)
            / (1.0 + 1.432788 * t + 0.189269 * t * t + 0.001308 * t * t * t));
    let target = p.ln();
    for _ in 0..8 {
        let cdf = norm_cdf(x);
        let step = (cdf.ln() - target) * cdf / norm_pdf(x);
        x -= step;
        if step.abs() <= 1e-15 * x.abs().max(1.0) {
            break;
        }
    }
    x
}

/// Regularized incomplete beta function `I_x(a, b)` for `a, b > 0` and
/// `x` in `[0, 1]`; NaN outside that domain.
///
/// Evaluated by the continued fraction for `I_x(a, b)` (modified Lentz),
/// switching to `1 - I_{1-x}(b, a)` where that converges faster; relative
/// accuracy is about `1e-14` away from underflow.
///
/// ```
/// use prospicio_math::special::beta_inc;
///
/// // I_x(1, 1) = x and I_x(a, 1) = x^a.
/// assert!((beta_inc(1.0, 1.0, 0.3) - 0.3).abs() < 1e-15);
/// assert!((beta_inc(2.5, 1.0, 0.3) - 0.3f64.powf(2.5)).abs() < 1e-15);
/// ```
pub fn beta_inc(a: f64, b: f64, x: f64) -> f64 {
    if !(a > 0.0 && b > 0.0) || x.is_nan() || !(0.0..=1.0).contains(&x) {
        return f64::NAN;
    }
    if x == 0.0 || x == 1.0 {
        return x;
    }
    // x^a (1 - x)^b / (a B(a, b)), in logs.
    let ln_front = a * x.ln() + b * (-x).ln_1p() + ln_gamma(a + b) - ln_gamma(a) - ln_gamma(b);
    if x < (a + 1.0) / (a + b + 2.0) {
        (ln_front.exp() / a) * beta_cf(a, b, x)
    } else {
        1.0 - (ln_front.exp() / b) * beta_cf(b, a, 1.0 - x)
    }
}

/// Continued fraction for the incomplete beta (Numerical Recipes `betacf`),
/// by the modified Lentz method.
fn beta_cf(a: f64, b: f64, x: f64) -> f64 {
    const TINY: f64 = 1e-300;
    let (qab, qap, qam) = (a + b, a + 1.0, a - 1.0);
    let mut c = 1.0;
    let mut d = 1.0 - qab * x / qap;
    if d.abs() < TINY {
        d = TINY;
    }
    d = 1.0 / d;
    let mut h = d;
    for m in 1..=10_000 {
        let m = f64::from(m);
        let m2 = 2.0 * m;
        let aa = m * (b - m) * x / ((qam + m2) * (a + m2));
        d = 1.0 + aa * d;
        if d.abs() < TINY {
            d = TINY;
        }
        c = 1.0 + aa / c;
        if c.abs() < TINY {
            c = TINY;
        }
        d = 1.0 / d;
        h *= d * c;
        let aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2));
        d = 1.0 + aa * d;
        if d.abs() < TINY {
            d = TINY;
        }
        c = 1.0 + aa / c;
        if c.abs() < TINY {
            c = TINY;
        }
        d = 1.0 / d;
        let delta = d * c;
        h *= delta;
        if (delta - 1.0).abs() < 1e-16 {
            break;
        }
    }
    h
}

/// Student's t distribution function with `nu > 0` degrees of freedom.
///
/// Uses `P(T > |t|) = I_{nu / (nu + t^2)}(nu / 2, 1 / 2) / 2`, so both
/// tails keep full relative accuracy.
///
/// ```
/// use prospicio_math::special::student_t_cdf;
///
/// // With one degree of freedom, t is Cauchy: F(1) = 3/4.
/// assert!((student_t_cdf(1.0, 1.0) - 0.75).abs() < 1e-15);
/// ```
pub fn student_t_cdf(t: f64, nu: f64) -> f64 {
    if t.is_nan() || nu.is_nan() || nu <= 0.0 {
        return f64::NAN;
    }
    if t.is_infinite() {
        return if t > 0.0 { 1.0 } else { 0.0 };
    }
    let x = nu / (nu + t * t);
    let tail = 0.5 * beta_inc(0.5 * nu, 0.5, x);
    if t > 0.0 { 1.0 - tail } else { tail }
}

/// Quantile of Student's t with `nu > 0` degrees of freedom: the inverse
/// of [`student_t_cdf`], by the Illinois method on a bracket doubled out
/// from the origin. NaN for `p` outside `[0, 1]`; `-inf` and `inf` at the
/// ends.
///
/// ```
/// use prospicio_math::special::student_t_quantile;
///
/// // With two degrees of freedom, (2p - 1) / sqrt(2 p (1 - p)).
/// let p: f64 = 0.9;
/// let want = (2.0 * p - 1.0) / (2.0 * p * (1.0 - p)).sqrt();
/// assert!((student_t_quantile(p, 2.0) - want).abs() < 1e-12);
/// ```
pub fn student_t_quantile(p: f64, nu: f64) -> f64 {
    if !(0.0..=1.0).contains(&p) || nu.is_nan() || nu <= 0.0 {
        return f64::NAN;
    }
    if p == 0.0 {
        return f64::NEG_INFINITY;
    }
    if p == 1.0 {
        return f64::INFINITY;
    }
    if p == 0.5 {
        return 0.0;
    }
    // Solve in the upper half and reflect: the quantile is odd about 1/2.
    let (q, sign) = if p > 0.5 { (p, 1.0) } else { (1.0 - p, -1.0) };
    let f = |x: f64| student_t_cdf(x, nu) - q;
    let mut hi = 1.0;
    while f(hi) < 0.0 && hi < 1e300 {
        hi *= 2.0;
    }
    let lo = hi / 2.0;
    let (lo, f_lo) = if hi == 1.0 {
        (0.0, f(0.0))
    } else {
        (lo, f(lo))
    };
    sign * crate::roots::illinois(lo, hi, f_lo, f(hi), f)
}

/// Regularized incomplete gamma functions `(P(a, x), Q(a, x))`, the lower
/// and upper parts, for shape `a > 0` and `x >= 0`: `P(a, x)` is the
/// distribution function of a unit-scale gamma with shape `a` at `x`, and
/// `Q = 1 - P`.
///
/// The smaller of the two is computed directly, by the power series below
/// `x = a + 1` and by Lentz's continued fraction above, so both tails keep
/// their relative precision; the other is `1` minus it. The common factor
/// `x^a e^(-x) / Γ(a)` is formed in Stirling form for `a >= 10`, which
/// keeps full precision for shapes in the thousands. NaN for invalid
/// arguments.
///
/// ```
/// use prospicio_math::special::gamma_inc;
///
/// // Shape 1 is the exponential: P(1, x) = 1 - e^(-x).
/// let (p, q) = gamma_inc(1.0, 2.0);
/// assert!((p - (1.0 - (-2f64).exp())).abs() < 1e-15);
/// assert!((q - (-2f64).exp()).abs() < 1e-16);
/// ```
pub fn gamma_inc(a: f64, x: f64) -> (f64, f64) {
    if a.is_nan() || x.is_nan() || a <= 0.0 || x < 0.0 || a == f64::INFINITY {
        return (f64::NAN, f64::NAN);
    }
    if x == 0.0 {
        return (0.0, 1.0);
    }
    if x == f64::INFINITY {
        return (1.0, 0.0);
    }
    let factor = gamma_inc_factor(a, x).exp();
    if x < a + 1.0 {
        // P = factor × Σ_n x^n / (a (a+1) ⋯ (a+n)).
        let mut term = 1.0 / a;
        let mut sum = term;
        for n in 1..100_000 {
            term *= x / (a + f64::from(n));
            sum += term;
            if term <= sum * 1e-17 {
                break;
            }
        }
        let p = (factor * sum).min(1.0);
        (p, 1.0 - p)
    } else {
        // Q = factor / (x + 1 - a - 1(1 - a) / (x + 3 - a - 2(2 - a) / ⋯)).
        const TINY: f64 = 1e-300;
        let mut b = x + 1.0 - a;
        let mut c = 1.0 / TINY;
        let mut d = 1.0 / b;
        let mut h = d;
        for i in 1..100_000 {
            let an = -f64::from(i) * (f64::from(i) - a);
            b += 2.0;
            d = an * d + b;
            if d.abs() < TINY {
                d = TINY;
            }
            c = b + an / c;
            if c.abs() < TINY {
                c = TINY;
            }
            d = 1.0 / d;
            let delta = d * c;
            h *= delta;
            if (delta - 1.0).abs() <= 1e-16 {
                break;
            }
        }
        let q = (factor * h).min(1.0);
        (1.0 - q, q)
    }
}

/// Euler's constant `γ`.
const EULER_GAMMA: f64 = 0.577_215_664_901_532_9;

/// The unnormalized lower incomplete beta
/// `B(a, b; x) = ∫_0^x t^(a-1) (1-t)^(b-1) dt` for `a > 0`, any `b`, and
/// `0 <= x < 1` (finite for every `b`, since `x < 1`).
///
/// For `b > 0` it is `I_x(a, b) B(a, b)`. Otherwise it steps `b` up to a
/// positive value (or zero) and back down by
/// `B(a, b; x) = ((a + b) B(a, b + 1; x) - x^a (1 - x)^b) / b`. These are
/// the limited moments of Pareto-tailed severities beyond the moments that
/// exist (loglogistic, Burr).
///
/// ```
/// use prospicio_math::special::beta_lower;
///
/// // B(1, 0; x) = -ln(1 - x) and B(1, -1; x) = x / (1 - x).
/// assert!((beta_lower(1.0, 0.0, 0.5) / 2f64.ln() - 1.0).abs() < 1e-14);
/// assert!((beta_lower(1.0, -1.0, 0.75) - 3.0).abs() < 1e-14);
/// ```
pub fn beta_lower(a: f64, b: f64, x: f64) -> f64 {
    if x <= 0.0 {
        return 0.0;
    }
    if b > 0.0 {
        let ln_b = ln_gamma(a) + ln_gamma(b) - ln_gamma(a + b);
        return beta_inc(a, b, x) * ln_b.exp();
    }
    // The smallest c = b + n that is positive, or zero.
    let steps = (-b).floor();
    let c = b + steps;
    let (mut value, mut c) = if c == 0.0 {
        (beta_zero(a, x), 0.0)
    } else {
        let c = c + 1.0;
        let ln_b = ln_gamma(a) + ln_gamma(c) - ln_gamma(a + c);
        (beta_inc(a, c, x) * ln_b.exp(), c)
    };
    while c > b {
        c -= 1.0;
        value = ((a + c) * value - x.powf(a) * (c * (-x).ln_1p()).exp()) / c;
    }
    value
}

/// `B(a, 0; x) = ∫_0^x t^(a-1) / (1 - t) dt`: the series
/// `Σ x^(a+n) / (a+n)` below `x = 1/2`; above it, with `y = 1 - x`,
/// `-ln y - ψ(a) - γ - Σ_{k≥1} (-1)^k C(a-1, k) y^k / k`, from
/// `∫_0^1 (1 - t^(a-1)) / (1 - t) dt = ψ(a) + γ`. For integer `a` the sum
/// stops and this is `-ln(1 - x) - Σ_{k<a} x^k / k`.
fn beta_zero(a: f64, x: f64) -> f64 {
    if x < 0.5 {
        let mut sum = 0.0;
        let mut power = x.powf(a);
        for n in 0..4000 {
            let add = power / (a + f64::from(n));
            sum += add;
            if add < 1e-17 * sum {
                break;
            }
            power *= x;
        }
        sum
    } else {
        let y = 1.0 - x;
        let mut sum = 0.0;
        // c_k = (-1)^k C(a - 1, k) y^k.
        let mut c = 1.0;
        for k in 1..4000 {
            let kf = f64::from(k);
            c *= -(a - kf) / kf * y;
            let add = c / kf;
            sum += add;
            if c == 0.0 || add.abs() < 1e-17 * sum.abs() {
                break;
            }
        }
        -y.ln() - digamma(a) - EULER_GAMMA - sum
    }
}

/// The unnormalized upper incomplete gamma `Γ(s, z) = ∫_z^∞ t^(s-1) e^(-t) dt`
/// for any real `s` and `z > 0` (for `s <= 0` it is finite only because
/// `z > 0`). NaN for `z <= 0` with `s <= 0`.
///
/// For `s > 0` it is `Γ(s) Q(s, z)`. Otherwise it starts from `s + n` in
/// `(0, 1)`, or from the exponential integral `E_1(z) = Γ(0, z)` when `s`
/// is an integer, and steps down by `Γ(c, z) = (Γ(c + 1, z) - z^c e^(-z)) / c`.
/// These are the limited moments of the inverse gamma beyond the moments
/// that exist.
///
/// ```
/// use prospicio_math::special::gamma_upper;
///
/// // Γ(1, z) = e^(-z); Γ(-1, z) = e^(-z)/z - E_1(z).
/// assert!((gamma_upper(1.0, 2.0) - (-2f64).exp()).abs() < 1e-16);
/// let e1 = gamma_upper(0.0, 2.0);
/// assert!((gamma_upper(-1.0, 2.0) - ((-2f64).exp() / 2.0 - e1)).abs() < 1e-15);
/// ```
pub fn gamma_upper(s: f64, z: f64) -> f64 {
    if s.is_nan() || z.is_nan() || z < 0.0 {
        return f64::NAN;
    }
    if s > 0.0 {
        return ln_gamma(s).exp() * gamma_inc(s, z).1;
    }
    if z == 0.0 {
        return f64::NAN;
    }
    let steps = (-s).floor();
    let c = s + steps;
    let (mut value, mut c) = if c == 0.0 {
        (expint_e1(z), 0.0)
    } else {
        let c = c + 1.0;
        (ln_gamma(c).exp() * gamma_inc(c, z).1, c)
    };
    while c > s {
        c -= 1.0;
        value = (value - (c * z.ln() - z).exp()) / c;
    }
    value
}

/// The exponential integral `E_1(z) = ∫_z^∞ e^(-t) / t dt` for `z > 0`:
/// the power series up to `z = 1`, Lentz's continued fraction above.
///
/// ```
/// use prospicio_math::special::expint_e1;
///
/// // SciPy: scipy.special.exp1(1.0).
/// assert!((expint_e1(1.0) / 0.21938393439552029 - 1.0).abs() < 1e-14);
/// ```
pub fn expint_e1(z: f64) -> f64 {
    if z.is_nan() || z < 0.0 {
        return f64::NAN;
    }
    if z == 0.0 {
        return f64::INFINITY;
    }
    if z == f64::INFINITY {
        return 0.0;
    }
    if z <= 1.0 {
        // -γ - ln z - Σ_{k≥1} (-z)^k / (k k!).
        let mut sum = 0.0;
        let mut term = 1.0;
        for k in 1..200 {
            let kf = f64::from(k);
            term *= -z / kf;
            let add = term / kf;
            sum += add;
            if add.abs() < 1e-17 * sum.abs() {
                break;
            }
        }
        -EULER_GAMMA - z.ln() - sum
    } else {
        const TINY: f64 = 1e-300;
        let mut b = z + 1.0;
        let mut c = 1.0 / TINY;
        let mut d = 1.0 / b;
        let mut h = d;
        for i in 1..10_000 {
            let an = -f64::from(i) * f64::from(i);
            b += 2.0;
            d = 1.0 / (an * d + b);
            c = b + an / c;
            let delta = c * d;
            h *= delta;
            if (delta - 1.0).abs() < 1e-16 {
                break;
            }
        }
        h * (-z).exp()
    }
}

/// The logarithm of the modified Bessel function of the second kind,
/// `ln K_ν(x)`, for real order `ν` and `x > 0`, from its integral
/// `K_ν(x) = ∫_0^∞ exp(-x cosh t) cosh(ν t) dt` by the trapezoidal rule.
/// The integrand is analytic in a strip about the real axis and decays
/// doubly exponentially, so the rule converges geometrically in the step,
/// which is sized to the integrand's peak; the sum is taken on a log
/// scale, so the value is finite where `K_ν(x)` itself overflows (small
/// `x`, large `ν`) or underflows (large `x`). `K_{-ν} = K_ν`.
///
/// ```
/// use prospicio_math::special::ln_bessel_k;
///
/// // K_{1/2}(x) = √(π / (2x)) e^(-x).
/// let x = 3.0_f64;
/// let exact = 0.5 * (std::f64::consts::PI / (2.0 * x)).ln() - x;
/// assert!((ln_bessel_k(0.5, x) - exact).abs() < 1e-14);
/// // mpmath: log(besselk(10, 0.2)).
/// assert!((ln_bessel_k(10.0, 0.2) / 35.13342019549504 - 1.0).abs() < 1e-14);
/// ```
pub fn ln_bessel_k(nu: f64, x: f64) -> f64 {
    if nu.is_nan() || x.is_nan() || x < 0.0 {
        return f64::NAN;
    }
    if x == 0.0 {
        return f64::INFINITY;
    }
    if x == f64::INFINITY {
        return f64::NEG_INFINITY;
    }
    let nu = nu.abs();
    // The exponent of the integrand relative to e^(-x):
    // -x (cosh t - 1) + ln cosh(ν t), which peaks at sinh t = ν / x with
    // curvature √(x² + ν²) there.
    let exponent = |t: f64| {
        let y = nu * t;
        -x * (2.0 * (0.5 * t).sinh().powi(2)) + y + (-2.0 * y).exp().ln_1p() - LN_2
    };
    let peak = (nu / x).asinh();
    let width = (x * x + nu * nu).powf(-0.25);
    let h = (0.5 * width).min(0.1);
    let top = exponent(peak);
    // Terms below e^-40 of the peak's change nothing at double precision.
    let mut sum = 0.5 * (exponent(0.0) - top).exp();
    let mut j = 1u32;
    loop {
        let t = f64::from(j) * h;
        let r = (exponent(t) - top).exp();
        sum += r;
        if t > peak && r < 1e-18 * sum {
            break;
        }
        j += 1;
        if j > 10_000_000 {
            break;
        }
    }
    -x + top + (h * sum).ln()
}

/// Mills ratio `R(t) = Φ(-t) / φ(t)` for `t >= 0`, the normal tail
/// relative to the density, finite where both underflow: directly below
/// `t = 10`, by the continued fraction `1 / (t + 1 / (t + 2 / (t + ⋯)))`
/// above.
///
/// ```
/// use prospicio_math::special::{mills_ratio, norm_cdf, norm_pdf};
///
/// assert!((mills_ratio(1.0) - norm_cdf(-1.0) / norm_pdf(1.0)).abs() < 1e-15);
/// // R(t) ~ 1/t far in the tail.
/// assert!((mills_ratio(1e6) * 1e6 - 1.0).abs() < 1e-11);
/// ```
pub fn mills_ratio(t: f64) -> f64 {
    if t.is_nan() {
        return f64::NAN;
    }
    if t < 10.0 {
        return norm_cdf(-t) / norm_pdf(t);
    }
    if t == f64::INFINITY {
        return 0.0;
    }
    // Backward evaluation of the continued fraction; 200 levels are far
    // more than t >= 10 needs.
    let mut v = 0.0;
    for k in (1..=200).rev() {
        v = f64::from(k) / (t + v);
    }
    1.0 / (t + v)
}

/// `ln(x^a e^(-x) / Γ(a))`. For `a >= 10`, as
/// `a (ln(1 + u) - u) + ln(a / 2π) / 2 - c(a)` with `u = (x - a) / a` and
/// `c` the Stirling series remainder of `ln Γ(a)`, so the large, nearly
/// cancelling terms `a ln x`, `x` and `ln Γ(a)` are never formed.
fn gamma_inc_factor(a: f64, x: f64) -> f64 {
    if a < 10.0 {
        return a * x.ln() - x - ln_gamma(a);
    }
    let u = (x - a) / a;
    let log1pmx = if u.abs() < 0.5 {
        // ln(1 + u) - u = Σ_{k≥2} (-1)^(k+1) u^k / k.
        let mut sum = 0.0;
        let mut power = u * u;
        for k in 2..200 {
            let term = power / f64::from(k);
            sum += if k % 2 == 0 { -term } else { term };
            if term.abs() <= 1e-17 * sum.abs() {
                break;
            }
            power *= u;
        }
        sum
    } else {
        // ln(x / a) directly: 1 + u would lose x's precision for x << a.
        (x / a).ln() - u
    };
    let inv = 1.0 / a;
    let inv2 = inv * inv;
    // ln Γ(a) - ((a - 1/2) ln a - a + ln(2π)/2), to about 1e-17 at a = 10.
    let stirling_remainder = inv
        * (1.0 / 12.0
            - inv2
                * (1.0 / 360.0
                    - inv2
                        * (1.0 / 1260.0
                            - inv2
                                * (1.0 / 1680.0
                                    - inv2
                                        * (1.0 / 1188.0
                                            - inv2 * (691.0 / 360_360.0 - inv2 / 156.0))))));
    a * log1pmx + 0.5 * (a / (2.0 * PI)).ln() - stirling_remainder
}

/// Gauss–Legendre nodes (negative half) and weights on `[-1, 1]` for 6, 12
/// and 20 points, as in Genz's BVND.
const GL6: ([f64; 3], [f64; 3]) = (
    [
        -0.932_469_514_203_152_2,
        -0.661_209_386_466_264_5,
        -0.238_619_186_083_197,
    ],
    [
        0.171_324_492_379_170_5,
        0.360_761_573_048_138_4,
        0.467_913_934_572_690_4,
    ],
);
const GL12: ([f64; 6], [f64; 6]) = (
    [
        -0.981_560_634_246_719_1,
        -0.904_117_256_370_475,
        -0.769_902_674_194_305,
        -0.587_317_954_286_617_1,
        -0.367_831_498_998_180_2,
        -0.125_233_408_511_469_2,
    ],
    [
        0.047_175_336_386_511_77,
        0.106_939_325_995_318_3,
        0.160_078_328_543_346_4,
        0.203_167_426_723_065_9,
        0.233_492_536_538_354_7,
        0.249_147_045_813_402_9,
    ],
);
const GL20: ([f64; 10], [f64; 10]) = (
    [
        -0.993_128_599_185_094_9,
        -0.963_971_927_277_913_8,
        -0.912_234_428_251_326,
        -0.839_116_971_822_218_8,
        -0.746_331_906_460_150_8,
        -0.636_053_680_726_515,
        -0.510_867_001_950_827_1,
        -0.373_706_088_715_419_6,
        -0.227_785_851_141_645_1,
        -0.076_526_521_133_497_33,
    ],
    [
        0.017_614_007_139_152_12,
        0.040_601_429_800_386_94,
        0.062_672_048_334_109_06,
        0.083_276_741_576_704_75,
        0.101_930_119_817_240_4,
        0.118_194_531_961_518_4,
        0.131_688_638_449_176_6,
        0.142_096_109_318_382_1,
        0.149_172_986_472_603_7,
        0.152_753_387_130_725_9,
    ],
);

/// The standard bivariate normal distribution function
/// `Φ₂(x, y; ρ) = P(X ≤ x, Y ≤ y)` for standard normals with correlation
/// `ρ` in `[-1, 1]`, by Genz's BVND (Drezner and Wesolowsky's method,
/// Gauss–Legendre on 6, 12 or 20 points as `|ρ|` grows): about `1e-15`
/// absolute.
///
/// ```
/// use prospicio_math::special::{bivariate_norm_cdf, norm_cdf};
///
/// // Independent at ρ = 0; P(X ≤ 0, Y ≤ 0) = 1/4 + asin(ρ) / (2π).
/// assert!((bivariate_norm_cdf(0.5, -1.0, 0.0) - norm_cdf(0.5) * norm_cdf(-1.0)).abs() < 1e-15);
/// let rho: f64 = 0.6;
/// let want = 0.25 + rho.asin() / (2.0 * std::f64::consts::PI);
/// assert!((bivariate_norm_cdf(0.0, 0.0, rho) - want).abs() < 1e-15);
/// ```
pub fn bivariate_norm_cdf(x: f64, y: f64, rho: f64) -> f64 {
    if x.is_nan() || y.is_nan() || rho.is_nan() {
        return f64::NAN;
    }
    if x == f64::NEG_INFINITY || y == f64::NEG_INFINITY {
        return 0.0;
    }
    if x == f64::INFINITY {
        return norm_cdf(y);
    }
    if y == f64::INFINITY {
        return norm_cdf(x);
    }
    bvn_upper(-x, -y, rho.clamp(-1.0, 1.0))
}

/// `P(X > h, Y > k)` for standard normals with correlation `r` (Genz's
/// BVNU).
fn bvn_upper(h: f64, k: f64, r: f64) -> f64 {
    let (nodes, weights): (&[f64], &[f64]) = if r.abs() < 0.3 {
        (&GL6.0, &GL6.1)
    } else if r.abs() < 0.75 {
        (&GL12.0, &GL12.1)
    } else {
        (&GL20.0, &GL20.1)
    };
    let two_pi = 2.0 * PI;
    let mut hk = h * k;
    let mut bvn = 0.0;
    if r.abs() < 0.925 {
        let hs = (h * h + k * k) / 2.0;
        let asr = r.asin();
        for (&x, &w) in nodes.iter().zip(weights) {
            for side in [-1.0, 1.0] {
                let sn = (asr * (side * x + 1.0) / 2.0).sin();
                bvn += w * ((sn * hk - hs) / (1.0 - sn * sn)).exp();
            }
        }
        return bvn * asr / (2.0 * two_pi) + norm_cdf(-h) * norm_cdf(-k);
    }
    let mut k = k;
    if r < 0.0 {
        k = -k;
        hk = -hk;
    }
    if r.abs() < 1.0 {
        let a2 = (1.0 - r) * (1.0 + r);
        let mut a = a2.sqrt();
        let bs = (h - k) * (h - k);
        let c = (4.0 - hk) / 8.0;
        let d = (12.0 - hk) / 16.0;
        let asr = -(bs / a2 + hk) / 2.0;
        if asr > -100.0 {
            bvn = a
                * asr.exp()
                * (1.0 - c * (bs - a2) * (1.0 - d * bs / 5.0) / 3.0 + c * d * a2 * a2 / 5.0);
        }
        if hk > -100.0 {
            let b = bs.sqrt();
            bvn -= (-hk / 2.0).exp()
                * two_pi.sqrt()
                * norm_cdf(-b / a)
                * b
                * (1.0 - c * bs * (1.0 - d * bs / 5.0) / 3.0);
        }
        a /= 2.0;
        for (&x, &w) in nodes.iter().zip(weights) {
            for side in [-1.0, 1.0] {
                let xs = (a * (side * x + 1.0)).powi(2);
                let rs = (1.0 - xs).sqrt();
                let asr = -(bs / xs + hk) / 2.0;
                if asr > -100.0 {
                    bvn += a
                        * w
                        * asr.exp()
                        * ((-hk * (1.0 - rs) / (2.0 * (1.0 + rs))).exp() / rs
                            - (1.0 + c * xs * (1.0 + d * xs)));
                }
            }
        }
        bvn = -bvn / two_pi;
    }
    if r > 0.0 {
        bvn + norm_cdf(-h.max(k))
    } else {
        let mut v = -bvn;
        if k > h {
            v += norm_cdf(k) - norm_cdf(h);
        }
        v.max(0.0)
    }
}

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

    #[test]
    fn ln_bessel_k_matches_mpmath() {
        // mpmath 1.3, 30 digits: log(besselk(nu, x)).
        let cases = [
            (0.0, 1.0, -0.8650643989067881),
            (0.3, 0.01, 1.930085981618933),
            (2.5, 3.0, -2.476216931302124),
            (-1.7, 0.5, 1.4915900467181742),
            (10.0, 0.2, 35.13342019549504),
            (0.5, 1000.0, -1003.2280862868463),
            (4.2, 2500.0, -2503.6827543489107),
            (30.0, 5.0, 42.860473318170534),
            (0.1, 1e-6, 2.9467464507698966),
        ];
        for (nu, x, want) in cases {
            let got = ln_bessel_k(nu, x);
            assert!(
                (got - want).abs() < 2e-15 * want.abs().max(1.0),
                "{nu} {x}: {got} {want}"
            );
        }
        // The recurrence K_{ν+1} = K_{ν-1} + (2ν / x) K_ν.
        let (nu, x) = (1.3, 0.7);
        let k = |n: f64| ln_bessel_k(n, x).exp();
        assert!((k(nu + 1.0) / (k(nu - 1.0) + 2.0 * nu / x * k(nu)) - 1.0).abs() < 1e-14);
    }

    #[test]
    fn digamma_matches_the_derivative_of_ln_gamma() {
        // ψ(1/2) = -γ - 2 ln 2; ψ(n) = H_{n-1} - γ.
        let gamma = 0.5772156649015329;
        assert!((digamma(0.5) + gamma + 2.0 * 2f64.ln()).abs() < 1e-14);
        assert!((digamma(10.0) - (7129.0 / 2520.0 - gamma)).abs() < 1e-14);
        for x in [0.01, 0.3, 2.7, 45.0, 1e4] {
            let h = 1e-5 * x;
            let numeric = (ln_gamma(x + h) - ln_gamma(x - h)) / (2.0 * h);
            assert!(
                (digamma(x) - numeric).abs() < 1e-7 * numeric.abs().max(1.0),
                "{x}"
            );
        }
        assert!(digamma(0.0).is_nan());
    }

    #[test]
    fn quantile_inverts_cdf() {
        for &p in &[
            1e-300,
            1e-100,
            1e-20,
            1e-8,
            0.001,
            0.1,
            0.3,
            0.5,
            0.7,
            0.99,
            1.0 - 1e-12,
        ] {
            let x = norm_quantile(p);
            let back = norm_cdf(x);
            assert!(
                ((back - p) / p).abs() < 1e-13,
                "p = {p}, x = {x}, back = {back}"
            );
        }
    }

    #[test]
    fn quantile_known_values() {
        // SciPy 1.x: scipy.stats.norm.ppf.
        let cases = [
            (0.5, 0.0),
            (0.995, 2.5758293035489004),
            (0.001, -3.090232306167813),
            (1e-10, -6.361340902404056),
        ];
        for (p, want) in cases {
            let got = norm_quantile(p);
            assert!(
                (got - want).abs() < 1e-13,
                "p = {p}: got {got}, want {want}"
            );
        }
    }

    #[test]
    fn quantile_edges() {
        assert_eq!(norm_quantile(0.0), f64::NEG_INFINITY);
        assert_eq!(norm_quantile(1.0), f64::INFINITY);
        assert!(norm_quantile(1.5).is_nan());
        assert!(norm_quantile(f64::NAN).is_nan());
    }

    #[test]
    fn beta_inc_closed_forms_and_symmetry() {
        for &x in &[1e-10, 0.01, 0.2, 0.5, 0.8, 0.999] {
            // I_x(a, b) = 1 - I_{1-x}(b, a); skipped where rounding 1 - x
            // alone moves the result by more than the tolerance.
            for &(a, b) in &[(0.5, 0.5), (2.0, 3.0), (10.0, 0.7), (50.0, 40.0)] {
                if x < 1e-3 {
                    continue;
                }
                let s = beta_inc(a, b, x) + beta_inc(b, a, 1.0 - x);
                assert!((s - 1.0).abs() < 1e-13, "a = {a}, b = {b}, x = {x}");
            }
            // I_x(1/2, 1/2) = (2 / pi) asin(sqrt(x)).
            let arcsine = 2.0 / PI * x.sqrt().asin();
            assert!((beta_inc(0.5, 0.5, x) - arcsine).abs() < 1e-14 * arcsine.max(1e-300));
        }
        assert!(beta_inc(0.0, 1.0, 0.5).is_nan());
        assert!(beta_inc(1.0, 1.0, 1.5).is_nan());
        assert_eq!(beta_inc(2.0, 3.0, 0.0), 0.0);
        assert_eq!(beta_inc(2.0, 3.0, 1.0), 1.0);
    }

    #[test]
    fn student_t_quantile_inverts_the_cdf() {
        for &nu in &[1.0, 2.0, 3.5, 30.0] {
            for &p in &[1e-6, 0.01, 0.3, 0.5, 0.77, 0.999_999] {
                let x = student_t_quantile(p, nu);
                assert!(
                    (student_t_cdf(x, nu) - p).abs() < 1e-12 * p.max(1e-3),
                    "{nu} {p}"
                );
            }
        }
        // Cauchy: tan(pi (p - 1/2)).
        let p: f64 = 0.01;
        assert!((student_t_quantile(p, 1.0) / (PI * (p - 0.5)).tan() - 1.0).abs() < 1e-10);
        assert!(student_t_quantile(1.5, 2.0).is_nan());
        assert_eq!(student_t_quantile(0.0, 2.0), f64::NEG_INFINITY);
    }

    #[test]
    fn trigamma_matches_known_values() {
        // ψ'(1/2) = π²/2, ψ'(2) = π²/6 − 1, and ψ' is the slope of ψ.
        assert!((trigamma(0.5) - PI * PI / 2.0).abs() < 1e-13);
        assert!((trigamma(2.0) - (PI * PI / 6.0 - 1.0)).abs() < 1e-14);
        let (x, h) = (3.7, 1e-5);
        let slope = (digamma(x + h) - digamma(x - h)) / (2.0 * h);
        assert!((trigamma(x) - slope).abs() < 1e-8);
        assert!(trigamma(0.0).is_nan());
    }

    #[test]
    fn student_t_limits() {
        // nu = 2 has a closed form: F(t) = 1/2 + t / (2 sqrt(2 + t^2)).
        for &t in &[-30.0, -2.0, -0.1, 0.0, 0.7, 5.0] {
            let want = 0.5 + t / (2.0 * (2.0f64 + t * t).sqrt());
            assert!((student_t_cdf(t, 2.0) - want).abs() < 1e-15, "t = {t}");
        }
        // Large nu approaches the normal.
        assert!((student_t_cdf(1.3, 1e7) - norm_cdf(1.3)).abs() < 1e-7);
        assert_eq!(student_t_cdf(f64::INFINITY, 3.0), 1.0);
        assert_eq!(student_t_cdf(f64::NEG_INFINITY, 3.0), 0.0);
        assert!(student_t_cdf(1.0, 0.0).is_nan());
    }

    /// mpmath 1.4.1 at 30 digits: `quad` of the integrand for `B(a, b; x)`,
    /// `gammainc(s, z)` for `Γ(s, z)`, `e1(z)` for `E_1(z)`.
    #[test]
    fn beta_lower_matches_mpmath() {
        for (a, b, x, want) in [
            (0.7, -2.0, 0.6, 3.598744771047808),
            (1.5, 0.5, 0.7, 0.5328990169356083),
            (2.0, -1.0, 0.9, 6.697414907005956),
            (2.25, -0.25, 0.4, 0.08655486451119004),
            (3.0, -1.0, 0.2, 0.003712897371580489),
            (2.5, 0.0, 0.3, 0.02525438280078752),
            (2.5, 0.0, 0.8, 0.6213887331578339),
            (3.0, 0.0, 0.9, 0.9975850929940459),
        ] {
            let got = beta_lower(a, b, x);
            assert!((got / want - 1.0).abs() < 1e-13, "{a} {b} {x}: {got}");
        }
    }

    #[test]
    fn gamma_upper_matches_mpmath() {
        assert!((expint_e1(0.1) / 1.8229239584193908 - 1.0).abs() < 1e-14);
        assert!((expint_e1(5.0) / 0.0011482955912753257 - 1.0).abs() < 1e-14);
        assert!((expint_e1(50.0) / 3.783264029550459e-24 - 1.0).abs() < 1e-13);
        for (s, z, want) in [
            (-0.5, 0.3, 1.1503670473551644),
            (-1.5, 2.0, 0.011832994103345996),
            (-2.0, 0.7, 0.3389003309406555),
            (0.4, 1.5, 0.13628632343383457),
            (-0.25, 4.0, 0.0025577114691076545),
        ] {
            let got = gamma_upper(s, z);
            assert!((got / want - 1.0).abs() < 1e-13, "{s} {z}: {got}");
        }
    }

    #[test]
    fn mills_ratio_is_continuous_at_the_switch() {
        let below = norm_cdf(-10.0) / norm_pdf(10.0);
        assert!((mills_ratio(10.0) / below - 1.0).abs() < 1e-13);
        assert!((mills_ratio(0.0) - (PI / 2.0).sqrt()).abs() < 1e-15);
    }
    #[test]
    fn bivariate_normal_edges_and_symmetry() {
        let (x, y) = (0.7, -0.4);
        // ρ = 1: P(X ≤ min(x, y)); ρ = −1: P(−y ≤ X ≤ x) or 0.
        assert!((bivariate_norm_cdf(x, y, 1.0) - norm_cdf(y)).abs() < 1e-15);
        assert!((bivariate_norm_cdf(x, y, -1.0) - (norm_cdf(x) - norm_cdf(-y))).abs() < 1e-15);
        assert_eq!(bivariate_norm_cdf(-1.0, 0.5, -1.0), 0.0);
        assert_eq!(bivariate_norm_cdf(f64::NEG_INFINITY, y, 0.3), 0.0);
        assert_eq!(bivariate_norm_cdf(f64::INFINITY, y, 0.3), norm_cdf(y));
        assert!(bivariate_norm_cdf(f64::NAN, y, 0.3).is_nan());
        for r in [-0.99, -0.6, -0.2, 0.0, 0.4, 0.8, 0.95, 0.999] {
            let c = bivariate_norm_cdf(x, y, r);
            assert!((c - bivariate_norm_cdf(y, x, r)).abs() < 1e-15, "{r}");
            // P(X ≤ x, Y ≤ y) + P(X ≤ x, Y > y) = Φ(x), with −Y.
            let rest = bivariate_norm_cdf(x, -y, -r);
            assert!((c + rest - norm_cdf(x)).abs() < 1e-14, "{r}");
            assert!((0.0..=1.0).contains(&c));
        }
    }
}