diff --git a/arborist/attest/__init__.py b/arborist/attest/__init__.py
index 0d8d82c..d2d4d0a 100644
--- a/arborist/attest/__init__.py
+++ b/arborist/attest/__init__.py
@@ -23,6 +23,14 @@ from arborist.attest.chain import (
     entity_chain_iter,
     verify_entity_chain,
 )
+from arborist.attest.consistency import (
+    KSResult,
+    estimate_cov,
+    estimate_mean,
+    invert_matrix,
+    ks_two_sample,
+    mahalanobis,
+)
 from arborist.attest.fingerprint import (
     Fingerprint,
     canonical_bytes,
@@ -32,11 +40,17 @@ from arborist.attest.fingerprint import (
 
 __all__ = [
     "Fingerprint",
+    "KSResult",
     "canonical_bytes",
     "commit_fingerprint",
     "entity_chain_hash",
     "entity_chain_iter",
+    "estimate_cov",
+    "estimate_mean",
     "from_samples",
+    "invert_matrix",
+    "ks_two_sample",
     "leaf_hash",
+    "mahalanobis",
     "verify_entity_chain",
 ]
diff --git a/arborist/attest/consistency.py b/arborist/attest/consistency.py
new file mode 100644
index 0000000..60adcc3
--- /dev/null
+++ b/arborist/attest/consistency.py
@@ -0,0 +1,196 @@
+"""Pure-stats consistency checks for Fingerprint sequences.
+
+Two deterministic, stdlib-only statistics:
+
+- **Two-sample Kolmogorov-Smirnov** on raw latency samples. Returns
+  the K-S statistic ``D`` and an asymptotic ``p`` from the
+  Kolmogorov distribution. ``D=0, p=1`` for identical CDFs; ``D=1,
+  p≈0`` for disjoint supports.
+- **Mahalanobis squared distance** of a moment vector from a
+  historical baseline ``(mean, inv_cov)``. Returns the raw
+  squared distance.
+
+Both functions return raw numbers. The caller picks the
+accept/reject threshold (e.g. ``reject if p < 0.05`` for K-S, or
+``reject if mahalanobis > chi2_95(k)`` for k-degree moment
+vectors). Tuning a threshold by training a discriminator against
+the corpus = sneaking a soft classifier in the side door —
+forbidden by the project's verifier-stays-binary rule.
+"""
+
+from __future__ import annotations
+
+import math
+from dataclasses import dataclass
+from typing import Sequence
+
+
+@dataclass(frozen=True)
+class KSResult:
+    """Two-sample Kolmogorov-Smirnov result."""
+
+    statistic: float
+    p_value: float
+    n1: int
+    n2: int
+
+
+def _empirical_cdf_step(
+    sorted_samples: Sequence[int | float], x: float
+) -> float:
+    """Fraction of sorted_samples <= x."""
+    # Binary search for rightmost index with value <= x.
+    lo, hi = 0, len(sorted_samples)
+    while lo < hi:
+        mid = (lo + hi) // 2
+        if sorted_samples[mid] <= x:
+            lo = mid + 1
+        else:
+            hi = mid
+    return lo / len(sorted_samples)
+
+
+def _ks_p_value(d: float, n1: int, n2: int, terms: int = 100) -> float:
+    """Asymptotic p-value for two-sample K-S via the Kolmogorov series.
+
+    Q(x) = 2 * sum_{k=1..inf} (-1)^(k-1) * exp(-2 * k^2 * x^2)
+    with x = D * sqrt(n_eff) and n_eff = n1*n2 / (n1+n2). Series
+    converges fast; 100 terms is overkill for any practical D > 0.
+    """
+    if n1 == 0 or n2 == 0:
+        return 1.0
+    if d == 0.0:
+        return 1.0
+    n_eff = (n1 * n2) / (n1 + n2)
+    x = d * math.sqrt(n_eff)
+    s = 0.0
+    sign = 1.0
+    for k in range(1, terms + 1):
+        term = sign * math.exp(-2.0 * k * k * x * x)
+        s += term
+        sign = -sign
+        if abs(term) < 1e-20:
+            break
+    p = 2.0 * s
+    return max(0.0, min(1.0, p))
+
+
+def ks_two_sample(
+    samples_a: Sequence[int | float], samples_b: Sequence[int | float]
+) -> KSResult:
+    """Two-sample Kolmogorov-Smirnov.
+
+    Raises ``ValueError`` if either sample is empty.
+    """
+    n1, n2 = len(samples_a), len(samples_b)
+    if n1 == 0 or n2 == 0:
+        raise ValueError("samples must be non-empty")
+    a = sorted(samples_a)
+    b = sorted(samples_b)
+    # Evaluate D at every distinct value present in either sample.
+    # The maximum |F1 - F2| is always reached at one of these points.
+    grid = sorted(set(a) | set(b))
+    d = 0.0
+    for x in grid:
+        f1 = _empirical_cdf_step(a, x)
+        f2 = _empirical_cdf_step(b, x)
+        gap = abs(f1 - f2)
+        if gap > d:
+            d = gap
+    p = _ks_p_value(d, n1, n2)
+    return KSResult(statistic=d, p_value=p, n1=n1, n2=n2)
+
+
+def estimate_mean(vectors: Sequence[Sequence[float]]) -> list[float]:
+    """Sample mean of N row-vectors of length k."""
+    if not vectors:
+        raise ValueError("vectors is empty")
+    k = len(vectors[0])
+    if any(len(v) != k for v in vectors):
+        raise ValueError("vectors must be same length")
+    n = len(vectors)
+    return [sum(v[i] for v in vectors) / n for i in range(k)]
+
+
+def estimate_cov(
+    vectors: Sequence[Sequence[float]], mean: Sequence[float] | None = None
+) -> list[list[float]]:
+    """Sample covariance (n-1 normalization) of N row-vectors of length k.
+
+    Returns a kxk matrix as a list of lists. Raises ``ValueError``
+    if ``n < 2`` (covariance undefined).
+    """
+    n = len(vectors)
+    if n < 2:
+        raise ValueError("need n>=2 vectors for sample covariance")
+    k = len(vectors[0])
+    mu = list(mean) if mean is not None else estimate_mean(vectors)
+    cov = [[0.0] * k for _ in range(k)]
+    for v in vectors:
+        for i in range(k):
+            di = v[i] - mu[i]
+            for j in range(k):
+                cov[i][j] += di * (v[j] - mu[j])
+    denom = n - 1
+    for i in range(k):
+        for j in range(k):
+            cov[i][j] /= denom
+    return cov
+
+
+def invert_matrix(m: Sequence[Sequence[float]]) -> list[list[float]]:
+    """Gauss-Jordan inverse of a square matrix.
+
+    Raises ``ValueError`` on singular input (pivot below 1e-12).
+    """
+    n = len(m)
+    if any(len(row) != n for row in m):
+        raise ValueError("matrix must be square")
+    # Build augmented [m | I] and reduce.
+    aug = [[float(m[i][j]) for j in range(n)] + [1.0 if i == j else 0.0 for j in range(n)] for i in range(n)]
+    for col in range(n):
+        # Partial pivot.
+        pivot = col
+        for r in range(col + 1, n):
+            if abs(aug[r][col]) > abs(aug[pivot][col]):
+                pivot = r
+        if abs(aug[pivot][col]) < 1e-12:
+            raise ValueError("matrix is singular")
+        if pivot != col:
+            aug[col], aug[pivot] = aug[pivot], aug[col]
+        # Scale pivot row.
+        p = aug[col][col]
+        for j in range(2 * n):
+            aug[col][j] /= p
+        # Eliminate other rows.
+        for r in range(n):
+            if r == col:
+                continue
+            factor = aug[r][col]
+            if factor == 0.0:
+                continue
+            for j in range(2 * n):
+                aug[r][j] -= factor * aug[col][j]
+    return [row[n:] for row in aug]
+
+
+def mahalanobis(
+    point: Sequence[float],
+    mean: Sequence[float],
+    inv_cov: Sequence[Sequence[float]],
+) -> float:
+    """Mahalanobis SQUARED distance from ``point`` to ``mean`` under ``inv_cov``.
+
+    ``return = (x - mu)^T * inv_cov * (x - mu)``. Always >= 0 for a
+    positive-definite ``inv_cov``. Caller compares against
+    chi-squared critical values for the dimensionality.
+    """
+    k = len(point)
+    if len(mean) != k:
+        raise ValueError("point and mean dimensionalities differ")
+    if len(inv_cov) != k or any(len(row) != k for row in inv_cov):
+        raise ValueError("inv_cov shape mismatch")
+    diff = [point[i] - mean[i] for i in range(k)]
+    # (M v)_i = sum_j M[i][j] * v[j]
+    mv = [sum(inv_cov[i][j] * diff[j] for j in range(k)) for i in range(k)]
+    return sum(diff[i] * mv[i] for i in range(k))
diff --git a/tests/test_attest_consistency.py b/tests/test_attest_consistency.py
new file mode 100644
index 0000000..1672063
--- /dev/null
+++ b/tests/test_attest_consistency.py
@@ -0,0 +1,224 @@
+"""Tests for `arborist.attest.consistency`."""
+
+from __future__ import annotations
+
+import math
+
+import pytest
+
+from arborist.attest import (
+    KSResult,
+    estimate_cov,
+    estimate_mean,
+    invert_matrix,
+    ks_two_sample,
+    mahalanobis,
+)
+from arborist.attest.consistency import _empirical_cdf_step, _ks_p_value
+
+
+# ---------------------------------------------------------------- K-S
+
+
+def test_ks_identical_samples():
+    r = ks_two_sample([1, 2, 3, 4, 5], [1, 2, 3, 4, 5])
+    assert isinstance(r, KSResult)
+    assert r.statistic == 0.0
+    assert r.p_value == 1.0
+    assert r.n1 == 5 and r.n2 == 5
+
+
+def test_ks_disjoint_supports():
+    """D=1 needs n>=~6 per side for the asymptotic p to clear 0.05.
+
+    With n1=n2=10 and D=1, n_eff=5, x=sqrt(5)≈2.236,
+    Q ≈ 2*exp(-2*5) ≈ 9e-5 — well below 0.01.
+    """
+    r = ks_two_sample(list(range(1, 11)), list(range(100, 110)))
+    assert r.statistic == 1.0
+    assert r.p_value < 0.01
+
+
+def test_ks_monotone_as_samples_diverge():
+    """As B drifts further from A, K-S statistic should not decrease."""
+    a = [1000, 1100, 1200, 1300, 1400]
+    d_near = ks_two_sample(a, [1050, 1150, 1250, 1350, 1450]).statistic
+    d_far = ks_two_sample(a, [2000, 2100, 2200, 2300, 2400]).statistic
+    assert d_far >= d_near
+
+
+def test_ks_rejects_empty():
+    with pytest.raises(ValueError):
+        ks_two_sample([], [1, 2, 3])
+    with pytest.raises(ValueError):
+        ks_two_sample([1, 2, 3], [])
+
+
+def test_ks_p_value_bounded():
+    """p in [0, 1] for any D in [0, 1]."""
+    for d in (0.0, 0.1, 0.5, 0.9, 1.0):
+        p = _ks_p_value(d, 10, 10)
+        assert 0.0 <= p <= 1.0
+
+
+def test_ks_p_value_zero_D_is_one():
+    assert _ks_p_value(0.0, 100, 100) == 1.0
+
+
+def test_ks_handles_overlap_partial():
+    """Partial overlap → 0 < D < 1, 0 < p < 1."""
+    a = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
+    b = [6, 7, 8, 9, 10, 11, 12, 13, 14, 15]
+    r = ks_two_sample(a, b)
+    assert 0.0 < r.statistic < 1.0
+    # Half-overlap with n=10 each → D=0.5; p moderate, not 0 or 1.
+    assert 0.0 < r.p_value < 1.0
+
+
+def test_empirical_cdf_step_basic():
+    s = sorted([10, 20, 30, 40, 50])
+    assert _empirical_cdf_step(s, 5) == 0.0
+    assert _empirical_cdf_step(s, 10) == 0.2
+    assert _empirical_cdf_step(s, 30) == 0.6
+    assert _empirical_cdf_step(s, 50) == 1.0
+    assert _empirical_cdf_step(s, 100) == 1.0
+
+
+# ---------------------------------------------------------------- mean / cov
+
+
+def test_estimate_mean_1d():
+    assert estimate_mean([[1.0], [2.0], [3.0]]) == [2.0]
+
+
+def test_estimate_mean_3d():
+    vecs = [[1.0, 10.0, 100.0], [2.0, 20.0, 200.0], [3.0, 30.0, 300.0]]
+    assert estimate_mean(vecs) == [2.0, 20.0, 200.0]
+
+
+def test_estimate_mean_rejects_ragged():
+    with pytest.raises(ValueError):
+        estimate_mean([[1.0, 2.0], [1.0]])
+
+
+def test_estimate_mean_rejects_empty():
+    with pytest.raises(ValueError):
+        estimate_mean([])
+
+
+def test_estimate_cov_independent_dimensions():
+    """Independent dimensions → diagonal covariance (within float epsilon)."""
+    # x ~ {1, 2, 3, 4, 5}; y constant; covariance off-diagonals = 0.
+    vecs = [[i, 0.0] for i in (1.0, 2.0, 3.0, 4.0, 5.0)]
+    cov = estimate_cov(vecs)
+    assert cov[0][1] == 0.0
+    assert cov[1][0] == 0.0
+
+
+def test_estimate_cov_rejects_n_lt_2():
+    with pytest.raises(ValueError):
+        estimate_cov([[1.0, 2.0]])
+
+
+def test_estimate_cov_symmetric():
+    vecs = [[1.0, 2.0], [3.0, 1.0], [5.0, 4.0], [2.0, 7.0]]
+    cov = estimate_cov(vecs)
+    assert math.isclose(cov[0][1], cov[1][0])
+
+
+# ---------------------------------------------------------------- invert
+
+
+def test_invert_identity_is_identity():
+    inv = invert_matrix([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]])
+    for i in range(3):
+        for j in range(3):
+            expected = 1.0 if i == j else 0.0
+            assert math.isclose(inv[i][j], expected, abs_tol=1e-12)
+
+
+def test_invert_diagonal():
+    inv = invert_matrix([[2.0, 0.0], [0.0, 4.0]])
+    assert math.isclose(inv[0][0], 0.5)
+    assert math.isclose(inv[1][1], 0.25)
+    assert math.isclose(inv[0][1], 0.0)
+    assert math.isclose(inv[1][0], 0.0)
+
+
+def test_invert_round_trip_3x3():
+    m = [[4.0, 7.0, 2.0], [3.0, 5.0, 1.0], [2.0, 1.0, 3.0]]
+    inv = invert_matrix(m)
+    # m * inv ≈ I
+    for i in range(3):
+        for j in range(3):
+            s = sum(m[i][k] * inv[k][j] for k in range(3))
+            expected = 1.0 if i == j else 0.0
+            assert math.isclose(s, expected, abs_tol=1e-10)
+
+
+def test_invert_rejects_singular():
+    with pytest.raises(ValueError):
+        invert_matrix([[1.0, 2.0], [2.0, 4.0]])
+
+
+def test_invert_rejects_non_square():
+    with pytest.raises(ValueError):
+        invert_matrix([[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]])
+
+
+# ---------------------------------------------------------------- mahalanobis
+
+
+def test_mahalanobis_at_mean_is_zero():
+    mean = [10.0, 20.0]
+    inv_cov = [[1.0, 0.0], [0.0, 1.0]]
+    assert mahalanobis(mean, mean, inv_cov) == 0.0
+
+
+def test_mahalanobis_identity_cov_is_squared_euclidean():
+    inv_cov = [[1.0, 0.0], [0.0, 1.0]]
+    d2 = mahalanobis([3.0, 4.0], [0.0, 0.0], inv_cov)
+    assert math.isclose(d2, 25.0)  # 3^2 + 4^2
+
+
+def test_mahalanobis_uses_inverse_covariance():
+    """Stretching inv_cov along an axis should scale that axis's contribution."""
+    inv_cov_iso = [[1.0, 0.0], [0.0, 1.0]]
+    inv_cov_x_heavy = [[4.0, 0.0], [0.0, 1.0]]  # x distance counts 4x more
+    d_iso = mahalanobis([1.0, 0.0], [0.0, 0.0], inv_cov_iso)
+    d_xh = mahalanobis([1.0, 0.0], [0.0, 0.0], inv_cov_x_heavy)
+    assert math.isclose(d_iso, 1.0)
+    assert math.isclose(d_xh, 4.0)
+
+
+def test_mahalanobis_dimension_mismatch():
+    with pytest.raises(ValueError):
+        mahalanobis([1.0, 2.0], [0.0, 0.0, 0.0], [[1.0, 0.0], [0.0, 1.0]])
+
+
+def test_mahalanobis_inv_cov_shape_mismatch():
+    with pytest.raises(ValueError):
+        mahalanobis([1.0, 2.0], [0.0, 0.0], [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0]])
+
+
+def test_mahalanobis_end_to_end_with_estimate():
+    """Estimate mean+cov from history → invert → score a new point.
+
+    Dimensions are chosen to vary independently so the sample
+    covariance is well away from singular. Callers must supply
+    non-singular history; Mahalanobis is undefined otherwise.
+    """
+    history = [
+        [1000.0, 100.0, 0.0],
+        [1100.0, 95.0, 10.0],
+        [950.0, 120.0, -5.0],
+        [1050.0, 85.0, 5.0],
+        [1020.0, 110.0, -2.0],
+    ]
+    mu = estimate_mean(history)
+    cov = estimate_cov(history, mean=mu)
+    inv_cov = invert_matrix(cov)
+    d_near = mahalanobis([1020.0, 100.0, 2.0], mu, inv_cov)
+    d_far = mahalanobis([5000.0, 500.0, 100.0], mu, inv_cov)
+    assert d_far > d_near
+    assert d_near >= 0.0