ballistics_engine/

special.rs

//! Special functions implemented in-crate rather than pulled from a dependency.
//!
//! MBA-1347 needs the mass of a bivariate normal over a target rectangle/circle, which
//! comes down to evaluating the standard normal CDF (and therefore the error function).
//! The crate ships to thirteen platforms including big-endian MIPS, RISC-V and wasm32
//! and is deliberately dependency-light -- it already hand-rolls its statistical
//! constants elsewhere -- so `erf`/`erfc`/`normal_cdf` are implemented and tested here
//! instead of pulling in a crate like `libm` or `statrs`.
//!
//! `ln_gamma` and `ln_beta` were added for Plan C (MBA-1352): the beta-binomial
//! mixture confidence sequence needs the log gamma and log beta functions, and the
//! same dependency-light posture applies -- implemented and tested here rather than
//! pulled in from `libm` or `statrs`.
//!
//! There are two independent approximations to validate: [`erfc`], of which [`erf`] and
//! [`normal_cdf`] are thin wrappers computed directly rather than as `1.0 -
//! erf(...)` so evaluating deep in a tail never cancels away the answer; and
//! [`ln_gamma`] (Lanczos, `g = 7`), of which [`ln_beta`] is in turn a thin wrapper.

/// Error function: `erf(x) = 2/sqrt(pi) * integral_0^x exp(-t^2) dt`.
///
/// Odd (`erf(-x) == -erf(x)`) and bounded in `[-1, 1]`. `erf(NAN)` is `NAN`;
/// `erf(+-INFINITY)` is `+-1.0` exactly.
pub fn erf(x: f64) -> f64 {
    if x.is_nan() {
        return f64::NAN;
    }
    if x < 0.0 {
        -(1.0 - erfc(-x))
    } else {
        1.0 - erfc(x)
    }
}

/// Complementary error function, `erfc(x) = 1 - erf(x)`.
///
/// Chebyshev-coefficient form (Numerical Recipes, 3rd ed., section 6.2.2, the `Erf`
/// struct's `erfccheb`), valid for any nonnegative `z` via the substitution `t = 2 / (2
/// + z)`, which maps `[0, infinity)` onto `(0, 1]` -- so the same 28-term expansion
/// covers the whole range without a separate small-`x`/large-`x` split. Negative
/// arguments use the exact identity `erfc(-x) = 2 - erfc(x)`.
///
/// Accurate to a couple of ULP of the true double-precision result across the whole
/// range (see this module's tests: the Abramowitz & Stegun Table 7.1 values, a denser
/// libm cross-check out to `x = 6`, and the far tail at `x = 6/sqrt(2)` used by
/// `normal_cdf(-6.0)`). `erfc(NAN)` is `NAN`; `erfc(INFINITY)` is `0.0`; `erfc(NEG_INFINITY)`
/// is `2.0`.
pub fn erfc(x: f64) -> f64 {
    if x < 0.0 {
        return 2.0 - erfc(-x);
    }
    let z = x.abs();
    let t = 2.0 / (2.0 + z);
    let ty = 4.0 * t - 2.0;
    const C: [f64; 28] = [
        -1.3026537197817094,
        6.419_697_923_564_902e-1,
        1.9476473204185836e-2,
        -9.561_514_786_808_63e-3,
        -9.46595344482036e-4,
        3.66839497852761e-4,
        4.2523324806907e-5,
        -2.0278578112534e-5,
        -1.624290004647e-6,
        1.303655835580e-6,
        1.5626441722e-8,
        -8.5238095915e-8,
        6.529054439e-9,
        5.059343495e-9,
        -9.91364156e-10,
        -2.27365122e-10,
        9.6467911e-11,
        2.394038e-12,
        -6.886027e-12,
        8.94487e-13,
        3.13092e-13,
        -1.12708e-13,
        3.81e-16,
        7.106e-15,
        -1.523e-15,
        -9.4e-17,
        1.21e-16,
        -2.8e-17,
    ];
    let mut d = 0.0f64;
    let mut dd = 0.0f64;
    for j in (1..C.len()).rev() {
        let tmp = d;
        d = ty * d - dd + C[j];
        dd = tmp;
    }
    t * (-z * z + 0.5 * (C[0] + ty * d) - dd).exp()
}

/// Standard normal (mean 0, variance 1) cumulative distribution function: `P(Z <= z)`.
///
/// Computed as `0.5 * erfc(-z / sqrt(2))` -- routing through [`erfc`] directly rather
/// than `0.5 * (1.0 + erf(z / sqrt(2)))` -- so the far lower tail (`z` very negative,
/// where `erf(z / sqrt(2))` is close to `-1.0`) is evaluated without cancellation.
/// `normal_cdf(NAN)` is `NAN`; `normal_cdf(INFINITY)` is `1.0`; `normal_cdf(NEG_INFINITY)`
/// is `0.0`.
pub fn normal_cdf(z: f64) -> f64 {
    0.5 * erfc(-z / std::f64::consts::SQRT_2)
}

/// Natural log of the gamma function on the positive reals.
///
/// Lanczos approximation (g = 7, 9 coefficients). The honest error claim is relative, not
/// absolute -- `ln_gamma`'s output ranges over many orders of magnitude, so a fixed absolute
/// bound is not meaningful at the extremes: measured relative error is <= 2.2e-15 (about 1
/// ULP) across the tested range, comfortably inside both the 1e-12 absolute bound the tests
/// pin there and the ~1e-10 the beta-binomial mixture confidence sequence actually needs.
/// Returns NaN for `x <= 0` and non-finite inputs -- this crate's special functions are total
/// functions with NaN as the out-of-domain signal rather than panicking, though unlike
/// [`erf`] (which returns a defined limit, `+-1.0`, at `+-infinity`), `ln_gamma` has no finite
/// limit there and returns NaN for every non-finite input, not just NaN itself.
pub fn ln_gamma(x: f64) -> f64 {
    if !x.is_finite() || x <= 0.0 {
        return f64::NAN;
    }
    const G: f64 = 7.0;
    const COEF: [f64; 9] = [
        0.999_999_999_999_809_9,
        676.520_368_121_885_1,
        -1_259.139_216_722_402_8,
        771.323_428_777_653_1,
        -176.615_029_162_140_6,
        12.507_343_278_686_905,
        -0.138_571_095_265_720_12,
        9.984_369_578_019_572e-6,
        1.505_632_735_149_311_6e-7,
    ];
    if x < 0.5 {
        // Reflection: ln Gamma(x) = ln(pi / sin(pi*x)) - ln Gamma(1 - x)
        let pi = std::f64::consts::PI;
        return (pi / (pi * x).sin()).ln() - ln_gamma(1.0 - x);
    }
    let x = x - 1.0;
    let mut acc = COEF[0];
    for (i, c) in COEF.iter().enumerate().skip(1) {
        acc += c / (x + i as f64);
    }
    let t = x + G + 0.5;
    0.5 * (2.0 * std::f64::consts::PI).ln() + (x + 0.5) * t.ln() - t + acc.ln()
}

/// ln B(a, b) = ln Gamma(a) + ln Gamma(b) - ln Gamma(a + b), same domain rule per argument.
pub fn ln_beta(a: f64, b: f64) -> f64 {
    ln_gamma(a) + ln_gamma(b) - ln_gamma(a + b)
}

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

    /// Reference values from Abramowitz & Stegun Table 7.1.
    #[test]
    fn erf_matches_published_reference_values() {
        for (x, want) in [
            (0.0, 0.0),
            (0.5, 0.5204998778),
            (1.0, 0.8427007929),
            (2.0, 0.9953222650),
            (3.0, 0.9999779095),
        ] {
            assert!(
                (erf(x) - want).abs() < 1e-9,
                "erf({x}) = {} want {want}",
                erf(x)
            );
        }
    }

    #[test]
    fn erf_is_odd_and_bounded() {
        for x in [0.3, 1.7, 4.2, 9.0] {
            assert!((erf(-x) + erf(x)).abs() < 1e-12);
            assert!(erf(x) <= 1.0 && erf(x) >= -1.0);
        }
    }

    #[test]
    fn normal_cdf_hits_known_points() {
        assert!((normal_cdf(0.0) - 0.5).abs() < 1e-12);
        assert!((normal_cdf(1.959_963_984_540_054) - 0.975).abs() < 1e-9);
        assert!((normal_cdf(-6.0)).abs() < 1e-8);
    }

    /// Denser cross-check against libm `erf` (Python `math.erf`, which on this platform
    /// wraps the system libm), x = 0.00..=6.00 step 0.25. The A&S table above only pins
    /// five points; this catches a bug confined to a narrow band the table would miss.
    /// Reference values are `repr()`-precision `f64` literals, i.e. exact round-trips.
    #[test]
    fn erf_matches_libm_cross_check_over_a_range() {
        let cases = [
            (0.0, 0.0),
            (0.25, 0.2763263901682369),
            (0.5, 0.5204998778130465),
            (0.75, 0.7111556336535152),
            (1.0, 0.8427007929497148),
            (1.25, 0.9229001282564582),
            (1.5, 0.9661051464753108),
            (1.75, 0.9866716712191824),
            (2.0, 0.9953222650189527),
            (2.25, 0.9985372834133188),
            (2.5, 0.999593047982555),
            (2.75, 0.9998993780778804),
            (3.0, 0.9999779095030015),
            (3.25, 0.9999956972205364),
            (3.5, 0.9999992569016276),
            (3.75, 0.9999998862727435),
            (4.0, 0.9999999845827421),
            (4.25, 0.9999999981494259),
            (4.5, 0.9999999998033839),
            (4.75, 0.999999999981515),
            (5.0, 0.9999999999984626),
            (5.25, 0.999999999999887),
            (5.5, 0.9999999999999927),
            (5.75, 0.9999999999999996),
            (6.0, 1.0),
        ];
        let mut worst: f64 = 0.0;
        for (x, want) in cases {
            let got = erf(x);
            let diff = (got - want).abs();
            worst = worst.max(diff);
            assert!(
                diff < 1e-12,
                "erf({x}) = {got:.17e} want {want:.17e} diff {diff:.3e}"
            );
        }
        assert!(worst < 1e-12, "worst diff over the sweep was {worst:.3e}");
    }

    /// A hit-probability integral truncated at roughly +-6 sigma depends on the tail
    /// being a real (if tiny) number rather than a flushed-to-zero sentinel. Reference
    /// value cross-checked against libm `erfc` (Python `math.erfc(6/sqrt(2))/2`):
    /// `normal_cdf(-6.0) == 9.865876450377014e-10`.
    #[test]
    fn normal_cdf_tail_is_real_not_flushed_to_zero() {
        let far_tail = normal_cdf(-6.0);
        assert!(
            far_tail > 0.0,
            "normal_cdf(-6.0) must be a real positive number, got {far_tail}"
        );
        assert_ne!(far_tail, 0.0, "tail was flushed to exactly 0.0");
        assert!(
            (far_tail - 9.865876450377014e-10).abs() < 1e-16,
            "normal_cdf(-6.0) = {far_tail:.17e}, want ~9.865876450377014e-10"
        );

        // erfc itself must stay strictly positive and strictly decreasing well past the
        // x = 6 the brief asks for (checked here to x = 8, still nowhere near the ~x >
        // 27 point where exp() genuinely underflows to zero).
        //
        // erfc(0.0) is mathematically exactly 1.0, but this is a general Chebyshev
        // evaluation with no z=0 special case, so it lands 1 ULP away
        // (0.9999999999999999) rather than bit-exact -- that is excellent accuracy, not
        // a defect, so this checks it approximately rather than with assert_eq!.
        let mut prev = erfc(0.0);
        assert!(
            (prev - 1.0).abs() < 1e-14,
            "erfc(0.0) = {prev:.17}, want ~1.0"
        );
        for i in 1..=32 {
            let x = i as f64 * 0.25;
            let v = erfc(x);
            assert!(v > 0.0, "erfc({x}) = {v} is not strictly positive");
            assert!(
                v < prev,
                "erfc not strictly decreasing at x={x}: prev={prev:e} v={v:e}"
            );
            prev = v;
        }
    }

    /// Across a sweep, `erf` must be strictly increasing and stay within `[-1, 1]`;
    /// `normal_cdf` strictly increasing within `[0, 1]`.
    ///
    /// The sweep stops at |x| = 5.0, short of where `erf(x) = 1.0 - erfc(x)`
    /// mathematically saturates to the nearest representable `f64` (empirically just
    /// past x = 5.9, once `erfc(x)` drops below half a ULP of 1.0): once two distinct
    /// real inputs round to the identical double there is no such thing as "strictly
    /// increasing" between them, in any correct double-precision implementation (libm's
    /// `erf` saturates at exactly the same place). That is a property of IEEE-754
    /// binary64, not a defect, so it is checked separately below instead of folded into
    /// the strict-monotonicity sweep.
    #[test]
    fn erf_and_normal_cdf_are_monotonic_and_bounded_over_a_sweep() {
        let xs = [
            -5.0, -4.0, -3.0, -2.5, -2.0, -1.5, -1.0, -0.75, -0.5, -0.25, -0.1, 0.0, 0.1, 0.25,
            0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0,
        ];

        let mut prev_erf = f64::NEG_INFINITY;
        let mut prev_cdf = f64::NEG_INFINITY;
        for x in xs {
            let e = erf(x);
            let c = normal_cdf(x);

            assert!((-1.0..=1.0).contains(&e), "erf({x}) = {e} out of [-1, 1]");
            assert!(
                (0.0..=1.0).contains(&c),
                "normal_cdf({x}) = {c} out of [0, 1]"
            );
            assert!(
                e > prev_erf,
                "erf not strictly increasing at x={x}: prev={prev_erf} now={e}"
            );
            assert!(
                c > prev_cdf,
                "normal_cdf not strictly increasing at x={x}: prev={prev_cdf} now={c}"
            );

            prev_erf = e;
            prev_cdf = c;
        }

        // Deep in the tail, multiple distinct real inputs necessarily collapse onto the
        // same double; what must still hold is that the bound is never exceeded. These
        // x values (10, 20, 50, 1e10) sit many orders of magnitude past the saturation
        // boundary on either side, so this is robust regardless of exactly which ULP
        // `erfc`'s Chebyshev approximation saturates on.
        for x in [10.0, 20.0, 50.0, 1.0e10] {
            let e = erf(x);
            let ne = erf(-x);
            let c = normal_cdf(x);
            let nc = normal_cdf(-x);

            assert_eq!(e, 1.0, "erf({x}) should have saturated to exactly 1.0");
            assert_eq!(ne, -1.0, "erf(-{x}) should have saturated to exactly -1.0");
            assert_eq!(
                c, 1.0,
                "normal_cdf({x}) should have saturated to exactly 1.0"
            );
            assert!(
                (0.0..=1.0).contains(&nc),
                "normal_cdf(-{x}) = {nc} out of [0, 1]"
            );
            assert!(!nc.is_nan());
        }
    }

    /// `normal_cdf(z) + normal_cdf(-z) == 1` for any `z`: the identity is exact in the
    /// implementation (`erfc(-w) = 2.0 - erfc(w)` is computed directly, not derived by
    /// cancellation), so this should hold to within a handful of ULP, not just the
    /// requested 1e-12.
    #[test]
    fn normal_cdf_symmetry_around_zero() {
        for z in [
            0.0,
            0.05,
            0.1,
            0.5,
            1.0,
            1.959_963_984_540_054,
            2.0,
            3.0,
            4.5,
            6.0,
            100.0,
        ] {
            let sum = normal_cdf(z) + normal_cdf(-z);
            assert!(
                (sum - 1.0).abs() < 1e-12,
                "normal_cdf({z}) + normal_cdf(-{z}) = {sum:.17}, want 1.0"
            );
        }
    }

    /// Non-finite inputs must behave predictably rather than accidentally (e.g. via a
    /// stray NaN comparison or an unguarded subtraction of infinities).
    #[test]
    fn non_finite_inputs_are_well_defined() {
        assert_eq!(erf(f64::INFINITY), 1.0);
        assert_eq!(erf(f64::NEG_INFINITY), -1.0);
        assert!(erf(f64::NAN).is_nan());

        assert_eq!(erfc(f64::INFINITY), 0.0);
        assert_eq!(erfc(f64::NEG_INFINITY), 2.0);
        assert!(erfc(f64::NAN).is_nan());

        assert_eq!(normal_cdf(f64::INFINITY), 1.0);
        assert_eq!(normal_cdf(f64::NEG_INFINITY), 0.0);
        assert!(normal_cdf(f64::NAN).is_nan());
    }

    // Reference values: DLMF 5.4 / Abramowitz & Stegun 6.1; ln Γ(0.5) = ln √π.
    #[test]
    fn ln_gamma_matches_published_reference_values() {
        let cases = [
            (0.5_f64, 0.572_364_942_924_700_1_f64), // ln sqrt(pi)
            (1.0, 0.0),
            (1.5, -0.120_782_237_635_245_22),
            (2.0, 0.0),
            (3.0, std::f64::consts::LN_2),  // ln 2! = ln 2
            (10.0, 12.801_827_480_081_469), // ln 9!
            (100.0, 359.134_205_369_575_4), // ln 99!
            (0.1, 2.252_712_651_734_206),   // reflection-region case
        ];
        for (x, expected) in cases {
            let got = ln_gamma(x);
            assert!(
                (got - expected).abs() < 1e-12,
                "ln_gamma({x}) = {got}, expected {expected}"
            );
        }
    }

    #[test]
    fn ln_gamma_rejects_nonpositive_and_nonfinite() {
        assert!(ln_gamma(0.0).is_nan());
        assert!(ln_gamma(-1.0).is_nan());
        assert!(ln_gamma(f64::NAN).is_nan());
        assert!(ln_gamma(f64::INFINITY).is_nan());
    }

    #[test]
    fn ln_gamma_satisfies_the_recurrence_identity() {
        // ln Γ(x+1) = ln Γ(x) + ln x — a property pin independent of the table above.
        for x in [0.3_f64, 0.7, 1.2, 4.5, 25.0, 170.0] {
            let lhs = ln_gamma(x + 1.0);
            let rhs = ln_gamma(x) + x.ln();
            assert!(
                (lhs - rhs).abs() < 1e-10,
                "recurrence failed at x={x}: {lhs} vs {rhs}"
            );
        }
    }

    #[test]
    fn ln_beta_is_the_gamma_combination_and_is_symmetric() {
        // B(1,1) = 1 → ln = 0; B(2,3) = 1/12 → ln = -ln 12.
        assert!((ln_beta(1.0, 1.0) - 0.0).abs() < 1e-14);
        assert!((ln_beta(2.0, 3.0) - (-(12.0_f64).ln())).abs() < 1e-12);
        assert!((ln_beta(7.3, 0.4) - ln_beta(0.4, 7.3)).abs() < 1e-13);
        assert!(ln_beta(0.0, 1.0).is_nan());
    }
}