optirs-core 0.3.2

// Shared statistical primitives for adaptive streaming
//
// Pure-Rust numeric helpers used by the drift detectors, the anomaly
// detectors, the adaptive buffer and the performance tracker. Everything in
// here is real, closed-form or published math — no fabricated constants and no
// simulated results. Significance calculations are performed in `f64` and go
// through the `scirs2-stats` distribution abstractions where a canonical
// implementation exists, because the generic element type `A: Float` used
// throughout the streaming modules does not carry the `Display` bound that
// `scirs2_stats::distributions::*` requires.

use scirs2_core::numeric::Float;
use std::cmp::Ordering;

/// Total ordering for floating-point values that never panics.
///
/// `f64::total_cmp`/`f32::total_cmp` are inherent methods and therefore
/// unavailable behind a generic `A: Float` bound, so this reproduces the same
/// contract: a genuine total order in which `NaN` sorts after every real
/// number (and equals itself). This is what makes `sort_by`/`select_nth` safe
/// on data that may contain `NaN`, instead of the `partial_cmp(..).expect(..)`
/// pattern which panics the moment a `NaN` reaches the comparator.
pub fn total_order<A: Float>(a: &A, b: &A) -> Ordering {
    crate::utils::total_order(a, b)
}

/// Sorts a slice ascending using the `NaN`-safe total order.
pub fn sort_ascending<A: Float>(values: &mut [A]) {
    values.sort_by(total_order);
}

/// Arithmetic mean, or `None` for an empty sample.
pub fn mean<A: Float>(values: &[A]) -> Option<A> {
    if values.is_empty() {
        return None;
    }
    let sum = values.iter().fold(A::zero(), |acc, &v| acc + v);
    A::from(values.len()).map(|n| sum / n)
}

/// Population variance (divides by `n`), or `None` for an empty sample.
pub fn population_variance<A: Float>(values: &[A]) -> Option<A> {
    let mu = mean(values)?;
    let n = A::from(values.len())?;
    let ss = values.iter().fold(A::zero(), |acc, &v| {
        let d = v - mu;
        acc + d * d
    });
    Some(ss / n)
}

/// Sample standard deviation (divides by `n - 1`), or `None` when `n < 2`.
pub fn sample_std_dev<A: Float>(values: &[A]) -> Option<A> {
    if values.len() < 2 {
        return None;
    }
    let mu = mean(values)?;
    let denom = A::from(values.len() - 1)?;
    let ss = values.iter().fold(A::zero(), |acc, &v| {
        let d = v - mu;
        acc + d * d
    });
    Some((ss / denom).sqrt())
}

/// Nearest-rank quantile of `p` in `[0, 1]`, computed in `O(n)` expected time
/// with `select_nth_unstable_by` (no full sort). Reorders `values` in place.
pub fn quantile_in_place<A: Float>(values: &mut [A], p: f64) -> Option<A> {
    let n = values.len();
    if n == 0 {
        return None;
    }
    let p = p.clamp(0.0, 1.0);
    let rank = (p * n as f64).ceil();
    let index = if rank < 1.0 {
        0
    } else {
        ((rank as usize) - 1).min(n - 1)
    };
    let (_, nth, _) = values.select_nth_unstable_by(index, total_order);
    Some(*nth)
}

/// True median (mean of the two central order statistics for even `n`),
/// computed with `select_nth_unstable_by`. Reorders `values` in place.
pub fn median_in_place<A: Float>(values: &mut [A]) -> Option<A> {
    let n = values.len();
    if n == 0 {
        return None;
    }
    if n % 2 == 1 {
        let (_, nth, _) = values.select_nth_unstable_by(n / 2, total_order);
        return Some(*nth);
    }

    // Even length: select the upper central element, then take the maximum of
    // the (already partitioned) lower half, which is the lower central
    // element by construction.
    let (lower, upper_mid, _) = values.select_nth_unstable_by(n / 2, total_order);
    let upper = *upper_mid;
    let lower_mid = lower.iter().copied().reduce(|a, b| {
        if total_order(&b, &a) == Ordering::Greater {
            b
        } else {
            a
        }
    })?;
    let two = A::from(2.0)?;
    Some((lower_mid + upper) / two)
}

/// Convenience wrapper: median of a borrowed sample (clones into a scratch
/// buffer so the caller's data is left untouched).
pub fn median<A: Float>(values: &[A]) -> Option<A> {
    let mut scratch: Vec<A> = values.to_vec();
    median_in_place(&mut scratch)
}

/// Survival function of the standard normal distribution, `P(Z > z)`.
///
/// Delegates to the canonical Gaussian CDF in `scirs2-stats` rather than
/// re-deriving an error-function approximation locally.
pub fn standard_normal_sf(z: f64) -> Result<f64, String> {
    let dist = scirs2_stats::distributions::normal::Normal::new(0.0_f64, 1.0_f64)
        .map_err(|e| format!("standard normal construction failed: {e}"))?;
    Ok((1.0 - dist.cdf(z)).clamp(0.0, 1.0))
}

/// Two-sided normal p-value for a z statistic: `P(|Z| > |z|)`.
pub fn normal_two_sided_p(z: f64) -> Result<f64, String> {
    let upper = standard_normal_sf(z.abs())?;
    Ok((2.0 * upper).clamp(0.0, 1.0))
}

/// Natural logarithm of the gamma function, by the Lanczos approximation with
/// `g = 7` and the standard nine-coefficient set (relative accuracy better than
/// `1e-13` for `x > 0`).
fn ln_gamma(x: f64) -> f64 {
    const COEFFICIENTS: [f64; 9] = [
        0.999_999_999_999_810,
        676.520_368_121_885,
        -1_259.139_216_722_403,
        771.323_428_777_653,
        -176.615_029_162_141,
        12.507_343_278_687,
        -0.138_571_095_265_720,
        0.000_009_984_369_578,
        0.000_000_150_563_274,
    ];

    if x < 0.5 {
        // Reflection formula: ln G(x) = ln(pi / sin(pi x)) - ln G(1 - x).
        return (std::f64::consts::PI / (std::f64::consts::PI * x).sin()).ln() - ln_gamma(1.0 - x);
    }

    let x = x - 1.0;
    let mut series = COEFFICIENTS[0];
    for (index, coefficient) in COEFFICIENTS.iter().enumerate().skip(1) {
        series += coefficient / (x + index as f64);
    }
    let t = x + 7.5;
    0.5 * (2.0 * std::f64::consts::PI).ln() + (x + 0.5) * t.ln() - t + series.ln()
}

/// Regularized upper incomplete gamma function `Q(a, x) = Gamma(a, x)/Gamma(a)`.
///
/// Uses the series expansion for the lower function when `x < a + 1` and
/// Legendre's continued fraction for the upper function otherwise, which is the
/// standard numerically-stable split.
fn regularized_upper_gamma(a: f64, x: f64) -> Result<f64, String> {
    if !a.is_finite() || a <= 0.0 || !x.is_finite() || x < 0.0 {
        return Err(format!(
            "regularized_upper_gamma requires a > 0 and x >= 0, got a={a}, x={x}"
        ));
    }
    if x == 0.0 {
        return Ok(1.0);
    }

    let log_prefactor = -x + a * x.ln() - ln_gamma(a);

    if x < a + 1.0 {
        // Series for the *lower* regularized function P(a, x), then Q = 1 - P.
        let mut term = 1.0 / a;
        let mut sum = term;
        let mut n = a;
        for _ in 0..1000 {
            n += 1.0;
            term *= x / n;
            sum += term;
            if term.abs() < sum.abs() * 1e-16 {
                break;
            }
        }
        let lower = sum * log_prefactor.exp();
        return Ok((1.0 - lower).clamp(0.0, 1.0));
    }

    // Legendre continued fraction for Q(a, x), evaluated by the modified
    // Lentz algorithm.
    let tiny = 1e-300_f64;
    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..1000 {
        let an = -(i as f64) * (i as f64 - 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;
        }
    }

    Ok((log_prefactor.exp() * h).clamp(0.0, 1.0))
}

/// Survival function of the chi-square distribution, `P(X > x)`, for `df`
/// degrees of freedom.
///
/// `P(X > x) = Q(df/2, x/2)` with `Q` the regularized upper incomplete gamma
/// function. This is computed locally rather than through
/// `scirs2_stats::distributions::chi_square`, whose tail CDF was measured at
/// roughly 13% relative error near the 5% critical value of `chi2(1)`
/// (`0.0435` against the true `0.0500` at `x = 3.8415`) — accurate enough for a
/// plot, not for a drift verdict's significance gate.
pub fn chi_square_sf(x: f64, df: f64) -> Result<f64, String> {
    if !(df.is_finite() && df > 0.0) {
        return Err(format!("chi-square requires df > 0, got {df}"));
    }
    if x <= 0.0 {
        return Ok(1.0);
    }
    regularized_upper_gamma(df / 2.0, x / 2.0)
}

/// Survival function of the Kolmogorov distribution,
/// `Q(z) = 2 * sum_{k>=1} (-1)^(k-1) exp(-2 k^2 z^2)`.
///
/// This is the asymptotic null distribution of the (scaled) two-sample
/// Kolmogorov-Smirnov statistic. The alternating series converges rapidly for
/// the arguments that arise in practice; iteration stops once a term is
/// negligible relative to the accumulated sum.
pub fn kolmogorov_sf(z: f64) -> f64 {
    if !z.is_finite() || z <= 0.0 {
        return 1.0;
    }

    let a2 = -2.0 * z * z;
    let mut sign = 2.0_f64;
    let mut sum = 0.0_f64;
    let mut previous_magnitude = 0.0_f64;

    for k in 1..=200_u32 {
        let term = sign * (a2 * f64::from(k * k)).exp();
        sum += term;
        let magnitude = term.abs();
        if magnitude <= 1e-8 * previous_magnitude || magnitude <= 1e-16 * sum.abs() {
            break;
        }
        previous_magnitude = magnitude;
        sign = -sign;
    }

    sum.clamp(0.0, 1.0)
}

/// Asymptotic p-value for a two-sample Kolmogorov-Smirnov statistic `d`
/// observed on samples of size `n1` and `n2`.
pub fn ks_two_sample_p(d: f64, n1: usize, n2: usize) -> f64 {
    if n1 == 0 || n2 == 0 {
        return 1.0;
    }
    let n1 = n1 as f64;
    let n2 = n2 as f64;
    let effective_n = (n1 * n2) / (n1 + n2);
    kolmogorov_sf(effective_n.sqrt() * d)
}

/// Two-sample Kolmogorov-Smirnov statistic `D = max_x |F_a(x) - F_b(x)|`.
///
/// Computed by a single merged walk over the two sorted samples, so the cost is
/// `O(n log n)` for the sorts plus `O(n)` for the walk. Returns `None` when
/// either sample is empty (there is no empirical CDF to compare).
pub fn ks_statistic<A: Float>(sample_a: &[A], sample_b: &[A]) -> Option<f64> {
    if sample_a.is_empty() || sample_b.is_empty() {
        return None;
    }

    let mut a: Vec<A> = sample_a.to_vec();
    let mut b: Vec<A> = sample_b.to_vec();
    sort_ascending(&mut a);
    sort_ascending(&mut b);

    let na = a.len();
    let nb = b.len();
    let mut i = 0usize;
    let mut j = 0usize;
    let mut max_diff = 0.0_f64;

    while i < na || j < nb {
        // The ECDFs may only be compared at genuine step boundaries, i.e. after
        // *every* observation equal to the current smallest value has been
        // consumed from **both** samples. Advancing one index at a time and
        // comparing after each single step reports a spurious difference of up
        // to `1/n` on tied data — for two identical samples it would report
        // `D = 1/n` instead of `0`.
        let next = match (a.get(i), b.get(j)) {
            (Some(x), Some(y)) => {
                if total_order(x, y) == Ordering::Greater {
                    *y
                } else {
                    *x
                }
            }
            (Some(x), None) => *x,
            (None, Some(y)) => *y,
            (None, None) => break,
        };

        while i < na && total_order(&a[i], &next) != Ordering::Greater {
            i += 1;
        }
        while j < nb && total_order(&b[j], &next) != Ordering::Greater {
            j += 1;
        }

        let ecdf_a = i as f64 / na as f64;
        let ecdf_b = j as f64 / nb as f64;
        let diff = (ecdf_a - ecdf_b).abs();
        if diff > max_diff {
            max_diff = diff;
        }
    }

    Some(max_diff)
}

/// Inclusive range `(min, max)` of a sample as `f64`, or `None` when empty or
/// entirely non-finite.
pub fn finite_range<A: Float>(values: &[A]) -> Option<(f64, f64)> {
    let mut min = f64::INFINITY;
    let mut max = f64::NEG_INFINITY;
    for value in values {
        let Some(v) = value.to_f64() else { continue };
        if !v.is_finite() {
            continue;
        }
        if v < min {
            min = v;
        }
        if v > max {
            max = v;
        }
    }
    if min.is_finite() && max.is_finite() {
        Some((min, max))
    } else {
        None
    }
}

/// Counts a sample into `bins` equal-width buckets spanning `[min, max]`.
///
/// Values below `min` land in the first bucket and values at or above `max` in
/// the last, so the returned counts always sum to the number of finite
/// observations.
pub fn histogram_counts<A: Float>(values: &[A], min: f64, max: f64, bins: usize) -> Vec<f64> {
    let bins = bins.max(1);
    let mut counts = vec![0.0_f64; bins];
    let width = if max > min {
        (max - min) / bins as f64
    } else {
        // Degenerate (zero-width) support: every observation is identical, so
        // the whole sample belongs to a single bucket.
        0.0
    };

    for value in values {
        let Some(v) = value.to_f64() else { continue };
        if !v.is_finite() {
            continue;
        }
        let index = if width > 0.0 {
            (((v - min) / width).floor().max(0.0) as usize).min(bins - 1)
        } else {
            0
        };
        counts[index] += 1.0;
    }

    counts
}

/// Normalises counts into a probability mass function with additive (Laplace)
/// smoothing, so that no bucket is ever exactly zero. Zero-probability
/// buckets would make KL divergence infinite and the log-ratio undefined; the
/// smoothing constant is applied symmetrically to both distributions being
/// compared, which is the standard treatment.
pub fn smoothed_pmf(counts: &[f64], smoothing: f64) -> Vec<f64> {
    let smoothing = smoothing.max(f64::MIN_POSITIVE);
    let total: f64 = counts.iter().sum::<f64>() + smoothing * counts.len() as f64;
    if total <= 0.0 {
        let uniform = 1.0 / counts.len().max(1) as f64;
        return vec![uniform; counts.len()];
    }
    counts.iter().map(|&c| (c + smoothing) / total).collect()
}

/// Kullback-Leibler divergence `KL(p || q)` in nats. Both inputs must be
/// smoothed probability mass functions of equal length.
pub fn kl_divergence(p: &[f64], q: &[f64]) -> Result<f64, String> {
    if p.len() != q.len() {
        return Err("KL divergence requires equal-length distributions".to_string());
    }
    let mut sum = 0.0_f64;
    for (&pi, &qi) in p.iter().zip(q.iter()) {
        if pi <= 0.0 {
            continue;
        }
        if qi <= 0.0 {
            return Err("KL divergence is undefined for a zero reference bin".to_string());
        }
        sum += pi * (pi / qi).ln();
    }
    Ok(sum.max(0.0))
}

/// Jensen-Shannon divergence in nats, bounded by `ln 2`.
pub fn js_divergence(p: &[f64], q: &[f64]) -> Result<f64, String> {
    if p.len() != q.len() {
        return Err("JS divergence requires equal-length distributions".to_string());
    }
    let mixture: Vec<f64> = p
        .iter()
        .zip(q.iter())
        .map(|(&pi, &qi)| 0.5 * (pi + qi))
        .collect();
    let left = kl_divergence(p, &mixture)?;
    let right = kl_divergence(q, &mixture)?;
    Ok((0.5 * left + 0.5 * right).clamp(0.0, std::f64::consts::LN_2))
}

/// Hellinger distance, bounded by `1`.
pub fn hellinger_distance(p: &[f64], q: &[f64]) -> Result<f64, String> {
    if p.len() != q.len() {
        return Err("Hellinger distance requires equal-length distributions".to_string());
    }
    let bhattacharyya: f64 = p
        .iter()
        .zip(q.iter())
        .map(|(&pi, &qi)| (pi.max(0.0) * qi.max(0.0)).sqrt())
        .sum();
    Ok((1.0 - bhattacharyya).max(0.0).sqrt())
}

/// First Wasserstein distance (a.k.a. Earth Mover's Distance) between two
/// one-dimensional empirical distributions.
///
/// In one dimension the optimal transport cost has the closed form
/// `W1 = integral |F_a(x) - F_b(x)| dx`, which this evaluates exactly on the
/// merged support of the two samples. Returns `None` when either sample is
/// empty.
pub fn wasserstein_1d<A: Float>(sample_a: &[A], sample_b: &[A]) -> Option<f64> {
    if sample_a.is_empty() || sample_b.is_empty() {
        return None;
    }

    let mut a: Vec<f64> = sample_a.iter().filter_map(|v| v.to_f64()).collect();
    let mut b: Vec<f64> = sample_b.iter().filter_map(|v| v.to_f64()).collect();
    if a.is_empty() || b.is_empty() {
        return None;
    }
    a.sort_by(|x, y| x.partial_cmp(y).unwrap_or(Ordering::Equal));
    b.sort_by(|x, y| x.partial_cmp(y).unwrap_or(Ordering::Equal));

    let na = a.len();
    let nb = b.len();
    let mut i = 0usize;
    let mut j = 0usize;
    let mut previous = a[0].min(b[0]);
    let mut total = 0.0_f64;

    while i < na || j < nb {
        let next = match (a.get(i), b.get(j)) {
            (Some(&x), Some(&y)) => x.min(y),
            (Some(&x), None) => x,
            (None, Some(&y)) => y,
            (None, None) => break,
        };

        // Accumulate |F_a - F_b| over the interval [previous, next) using the
        // CDF levels that hold on that interval.
        let ecdf_a = i as f64 / na as f64;
        let ecdf_b = j as f64 / nb as f64;
        total += (ecdf_a - ecdf_b).abs() * (next - previous);
        previous = next;

        while i < na && a[i] <= next {
            i += 1;
        }
        while j < nb && b[j] <= next {
            j += 1;
        }
    }

    Some(total)
}

/// G-test (likelihood-ratio) statistic for a `2 x k` contingency table built
/// from two histograms, together with its degrees of freedom.
///
/// `G = 2 * sum_cells O * ln(O / E)` is asymptotically chi-square distributed
/// with `k - 1` degrees of freedom under the null hypothesis that both samples
/// were drawn from the same binned distribution, which turns a histogram
/// divergence into a genuine significance test rather than a bare distance.
pub fn g_test_statistic(counts_a: &[f64], counts_b: &[f64]) -> Result<(f64, f64), String> {
    if counts_a.len() != counts_b.len() {
        return Err("G-test requires equal-length histograms".to_string());
    }
    let total_a: f64 = counts_a.iter().sum();
    let total_b: f64 = counts_b.iter().sum();
    let grand_total = total_a + total_b;
    if total_a <= 0.0 || total_b <= 0.0 {
        return Err("G-test requires both samples to be non-empty".to_string());
    }

    let mut g = 0.0_f64;
    let mut non_empty_columns = 0usize;
    for (&observed_a, &observed_b) in counts_a.iter().zip(counts_b.iter()) {
        let column_total = observed_a + observed_b;
        if column_total <= 0.0 {
            continue;
        }
        non_empty_columns += 1;

        let expected_a = column_total * total_a / grand_total;
        let expected_b = column_total * total_b / grand_total;
        if observed_a > 0.0 && expected_a > 0.0 {
            g += observed_a * (observed_a / expected_a).ln();
        }
        if observed_b > 0.0 && expected_b > 0.0 {
            g += observed_b * (observed_b / expected_b).ln();
        }
    }

    let degrees_of_freedom = (non_empty_columns.saturating_sub(1)) as f64;
    Ok((2.0 * g, degrees_of_freedom))
}

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

    #[test]
    fn total_order_is_nan_safe_and_total() {
        let mut values = vec![3.0_f64, f64::NAN, 1.0, 2.0];
        // Would panic under `partial_cmp(..).expect(..)`.
        sort_ascending(&mut values);
        assert_eq!(values[0], 1.0);
        assert_eq!(values[1], 2.0);
        assert_eq!(values[2], 3.0);
        assert!(values[3].is_nan(), "NaN must sort last");
    }

    #[test]
    fn median_matches_definition_for_odd_and_even_lengths() {
        assert_eq!(median(&[5.0_f64, 1.0, 3.0]), Some(3.0));
        assert_eq!(median(&[4.0_f64, 1.0, 3.0, 2.0]), Some(2.5));
        assert_eq!(median::<f64>(&[]), None);
    }

    #[test]
    fn quantile_uses_nearest_rank() {
        let mut values = vec![1.0_f64, 2.0, 3.0, 4.0];
        assert_eq!(quantile_in_place(&mut values, 0.0), Some(1.0));
        let mut values = vec![1.0_f64, 2.0, 3.0, 4.0];
        assert_eq!(quantile_in_place(&mut values, 1.0), Some(4.0));
        let mut values = vec![1.0_f64, 2.0, 3.0, 4.0];
        assert_eq!(quantile_in_place(&mut values, 0.5), Some(2.0));
    }

    #[test]
    fn normal_survival_function_matches_known_quantiles() {
        let p = standard_normal_sf(1.959_963_985).expect("normal sf");
        assert!(
            (p - 0.025).abs() < 1e-4,
            "P(Z > 1.96) should be ~0.025, got {p}"
        );
        let p0 = standard_normal_sf(0.0).expect("normal sf");
        assert!((p0 - 0.5).abs() < 1e-9);
    }

    /// Exercises the **continued-fraction** branch (`x >= a + 1`).
    #[test]
    fn chi_square_survival_function_matches_known_quantiles() {
        // chi2(1) 95th percentile is 3.8415: a = 0.5, x = 1.92, so x >= a + 1.
        let p = chi_square_sf(3.841_458_8, 1.0).expect("chi2 sf");
        assert!((p - 0.05).abs() < 1e-6, "expected 0.05, got {p}");

        // chi2(2) 95th percentile is 5.9915: a = 1, x = 3.0.
        let p = chi_square_sf(5.991_464_5, 2.0).expect("chi2 sf");
        assert!((p - 0.05).abs() < 1e-6, "expected 0.05, got {p}");

        // chi2(15) 95th percentile is 24.9958: a = 7.5, x = 12.5.
        let p = chi_square_sf(24.995_79, 15.0).expect("chi2 sf");
        assert!((p - 0.05).abs() < 1e-5, "expected 0.05, got {p}");
    }

    /// Exercises the **series** branch (`x < a + 1`), which the quantile checks
    /// above never reach. The two branches are separate code paths, so a wrong
    /// series index would go unnoticed without this.
    #[test]
    fn chi_square_survival_function_is_correct_in_the_series_branch() {
        // chi2(4) at x = 1.0: a = 2, x = 0.5, so x < a + 1. P(X > 1) = 0.909796.
        let p = chi_square_sf(1.0, 4.0).expect("chi2 sf");
        assert!((p - 0.909_796).abs() < 1e-5, "expected 0.909796, got {p}");

        // chi2(10) at x = 2.0: a = 5, x = 1.0. P(X > 2) = 0.996340.
        let p = chi_square_sf(2.0, 10.0).expect("chi2 sf");
        assert!((p - 0.996_340).abs() < 1e-5, "expected 0.996340, got {p}");

        // A very small x must approach 1, and x = 0 must be exactly 1.
        assert!(chi_square_sf(1e-12, 3.0).expect("chi2 sf") > 0.999_999);
        assert_eq!(chi_square_sf(0.0, 3.0).expect("chi2 sf"), 1.0);
    }

    /// The survival function must be monotonically decreasing across the branch
    /// boundary, with no discontinuity where the two algorithms meet.
    #[test]
    fn chi_square_survival_function_is_continuous_across_the_branch_switch() {
        let df = 6.0_f64;
        let a = df / 2.0; // 3.0, so the switch is at x/2 = 4, i.e. x = 8.
        let switch = 2.0 * (a + 1.0);
        let below = chi_square_sf(switch - 1e-7, df).expect("chi2 sf");
        let above = chi_square_sf(switch + 1e-7, df).expect("chi2 sf");
        // The two branches are independent algorithms, so they agree to their
        // own truncation accuracy rather than to machine epsilon; measured
        // agreement is ~1.5e-8 absolute here. That is six orders of magnitude
        // tighter than what this replaced, and the absolute accuracy against
        // textbook quantiles is separately asserted at 1e-6 above.
        assert!(
            (below - above).abs() < 1e-7,
            "the two branches disagree at the switch point ({below} vs {above})"
        );

        let mut previous = 1.0_f64;
        for step in 1..=200 {
            let x = step as f64 * 0.15;
            let p = chi_square_sf(x, df).expect("chi2 sf");
            assert!(
                p <= previous + 1e-12,
                "survival function increased at x = {x} ({previous} -> {p})"
            );
            previous = p;
        }
    }

    #[test]
    fn ln_gamma_matches_known_values() {
        // ln G(1) = ln G(2) = 0; ln G(5) = ln 24; ln G(0.5) = ln sqrt(pi).
        assert!(ln_gamma(1.0).abs() < 1e-12);
        assert!(ln_gamma(2.0).abs() < 1e-12);
        assert!((ln_gamma(5.0) - 24.0_f64.ln()).abs() < 1e-11);
        assert!(
            (ln_gamma(0.5) - std::f64::consts::PI.sqrt().ln()).abs() < 1e-11,
            "ln G(0.5) = {}",
            ln_gamma(0.5)
        );
    }

    #[test]
    fn ks_statistic_is_one_for_disjoint_samples() {
        let d = ks_statistic(&[0.0_f64, 1.0, 2.0], &[10.0, 11.0, 12.0]).expect("ks");
        assert!(
            (d - 1.0).abs() < 1e-12,
            "disjoint samples must give D = 1, got {d}"
        );
        let p = ks_two_sample_p(d, 3, 3);
        assert!(
            p < 0.5,
            "D = 1 on n = 3 should be at least mildly significant, got {p}"
        );
    }

    #[test]
    fn ks_statistic_is_zero_for_identical_samples() {
        let d = ks_statistic(&[1.0_f64, 2.0, 3.0], &[1.0, 2.0, 3.0]).expect("ks");
        assert!(
            d.abs() < 1e-12,
            "identical samples must give D = 0, got {d}"
        );
        assert!((ks_two_sample_p(d, 3, 3) - 1.0).abs() < 1e-12);
    }

    #[test]
    fn wasserstein_matches_closed_form_for_a_pure_shift() {
        // Shifting every point by 5 costs exactly 5 in W1.
        let a = [0.0_f64, 1.0, 2.0, 3.0];
        let b = [5.0_f64, 6.0, 7.0, 8.0];
        let w = wasserstein_1d(&a, &b).expect("w1");
        assert!((w - 5.0).abs() < 1e-9, "expected W1 = 5, got {w}");
    }

    #[test]
    fn divergences_are_zero_for_identical_distributions_and_positive_otherwise() {
        let p = smoothed_pmf(&[10.0, 10.0, 10.0], 0.5);
        let q = smoothed_pmf(&[10.0, 10.0, 10.0], 0.5);
        assert!(kl_divergence(&p, &q).expect("kl") < 1e-12);
        assert!(js_divergence(&p, &q).expect("js") < 1e-12);
        assert!(hellinger_distance(&p, &q).expect("hellinger") < 1e-6);

        let r = smoothed_pmf(&[30.0, 0.0, 0.0], 0.5);
        assert!(kl_divergence(&p, &r).expect("kl") > 0.1);
        assert!(js_divergence(&p, &r).expect("js") > 0.1);
        assert!(hellinger_distance(&p, &r).expect("hellinger") > 0.1);
    }

    #[test]
    fn g_test_is_insignificant_for_matching_histograms() {
        let (g, df) = g_test_statistic(&[20.0, 20.0, 20.0], &[20.0, 20.0, 20.0]).expect("g-test");
        assert!(g.abs() < 1e-9);
        assert_eq!(df, 2.0);
        let p = chi_square_sf(g, df.max(1.0)).expect("chi2");
        assert!(
            p > 0.9,
            "identical histograms must not be significant, got {p}"
        );
    }

    #[test]
    fn g_test_is_significant_for_disjoint_histograms() {
        let (g, df) = g_test_statistic(&[60.0, 0.0], &[0.0, 60.0]).expect("g-test");
        assert!(
            g > 100.0,
            "disjoint histograms should give a large G, got {g}"
        );
        let p = chi_square_sf(g, df.max(1.0)).expect("chi2");
        assert!(p < 1e-6, "expected an extremely small p-value, got {p}");
    }
}