feat(semisym): exact-rational Rust discriminator — jet AD, checkpointed survey, Δ₇ verdict

Order-3 BigRational jet AD (no symbolic algebra): Fisher–Rao + drift-flip
metric, Γ/R/∇R/R·R/Q(g,R) at rational points. Validated against the sympy
Δ₃ ground truth (exact match incl. (∇R)_(0,0,1,0,1)=1987584/1500625 and
the drift obstruction scalar). Signature via congruence elimination
(Sylvester's law of inertia) — fixes a row-aliasing bug that reported
(4,3) at the Δ₇ centroid; truth is (6,1), numpy-confirmed.

Δ₇ centroid results (31s single-core): baseline constant curvature 1/4
holds; g' signature (6,1) with drift time-like (g'(β,β)=−35); verdict
PROPER (11076/16807 nonzero ∇R, 591 inconsistent L ratios) — matches Δ₃.

survey subcommand = resumable DAG: one receipt JSON per lattice point
(skip-if-exists), rayon-parallel, --max-points deterministic stride cap;
receipts are ingestable scar-map events (obstruction ⇒ scar candidate).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

rust/semisym/Cargo.toml

@@ -0,0 +1,15 @@
+[package]
+name = "semisym"
+version = "0.1.0"
+edition = "2021"
+
+[dependencies]
+num-bigint = "0.4"
+num-rational = "0.4"
+num-traits = "0.2"
+rayon = "1.10"
+serde = { version = "1", features = ["derive"] }
+serde_json = "1"
+
+[profile.release]
+lto = true

rust/semisym/src/main.rs

@@ -0,0 +1,757 @@
+//! semisym — Milestone §0 covariant semi-symmetry discriminator (exact ℚ).
+//!
+//! Rust port of python/semisymmetry_discriminator.py, validated against its
+//! Δ₃ ground truth (docs/reviews/SEMISYMMETRY_DISCRIMINATOR_DELTA3_LOG.txt).
+//! No floats anywhere: BigRational jet (Taylor) automatic differentiation to
+//! third order at a rational point — no symbolic algebra, no expression swell.
+//!
+//! Geometry: Fisher–Rao on Δ_{m−1}, reduced coords p₁..p_{m−1};
+//! drift-flip metric g' = g − 2 β♭⊗β♭ / g(β,β) with β from the
+//! rossbyDriftFromChirality weights (left +1, right −1, scarred +1/2,
+//! achiral 0), mean-centered (BraidStateN.lean).
+//!
+//! Derivative-order bookkeeping (why order-3 jets suffice):
+//!   g' jets carry exact ∂⁰..∂³  →  ∂g' order 2  →  Γ order 2
+//!   →  ∂Γ order 1  →  R order 1  →  ∂R order 0  →  ∇R exact values.
+//!
+//! Conventions (Kobayashi–Nomizu, identical to the Python reference):
+//!   R^ρ_{σμν} = ∂_μ Γ^ρ_{νσ} − ∂_ν Γ^ρ_{μσ} + Γ^ρ_{μλ}Γ^λ_{νσ} − Γ^ρ_{νλ}Γ^λ_{μσ}
+//!   (A·T)_{ijkl} = −A^m_i T_{mjkl} − A^m_j T_{imkl} − A^m_k T_{ijml} − A^m_l T_{ijkm}
+//!
+//! Usage:
+//!   semisym probe  --m 8 --point 1/8,1/8,1/8,1/8,1/8,1/8,1/8
+//!   semisym survey --m 8 --grid 3 --checkpoint ckpt/   (resumable DAG:
+//!       one receipt JSON per lattice point; existing receipts are skipped)
+//!
+//! Each receipt is an ingestable DAG-node event for the ENE scar map:
+//! obstruction present ⇒ scar candidate; see receipt fields.
+
+use num_bigint::BigInt;
+use num_rational::BigRational;
+use num_traits::{One, Zero};
+use rayon::prelude::*;
+use serde::Serialize;
+use std::fmt::Write as _;
+use std::fs;
+use std::path::PathBuf;
+use std::str::FromStr;
+
+type Q = BigRational;
+
+fn q(n: i64, d: i64) -> Q {
+    Q::new(BigInt::from(n), BigInt::from(d))
+}
+
+fn parse_q(s: &str) -> Q {
+    let s = s.trim();
+    match s.split_once('/') {
+        Some((a, b)) => Q::new(
+            BigInt::from_str(a.trim()).expect("numerator"),
+            BigInt::from_str(b.trim()).expect("denominator"),
+        ),
+        None => Q::new(BigInt::from_str(s).expect("integer"), BigInt::one()),
+    }
+}
+
+// ── symmetric index helpers (n vars) ─────────────────────────────────
+
+#[inline]
+fn i2(n: usize, i: usize, j: usize) -> usize {
+    let (i, j) = if i <= j { (i, j) } else { (j, i) };
+    i * n - i * (i + 1) / 2 + j
+}
+
+fn i3_table(n: usize) -> (Vec<[usize; 3]>, Vec<Vec<Vec<usize>>>) {
+    let mut list = Vec::new();
+    let mut idx = vec![vec![vec![0usize; n]; n]; n];
+    for i in 0..n {
+        for j in i..n {
+            for k in j..n {
+                let id = list.len();
+                list.push([i, j, k]);
+                // all permutations map to id
+                let perms = [
+                    [i, j, k], [i, k, j], [j, i, k], [j, k, i], [k, i, j], [k, j, i],
+                ];
+                for p in perms {
+                    idx[p[0]][p[1]][p[2]] = id;
+                }
+            }
+        }
+    }
+    (list, idx)
+}
+
+// ── order-3 jets over ℚ ──────────────────────────────────────────────
+
+#[derive(Clone)]
+struct Jet {
+    n: usize,
+    v: Q,
+    d1: Vec<Q>,      // n
+    d2: Vec<Q>,      // n(n+1)/2, raw ∂²
+    d3: Vec<Q>,      // C(n+2,3), raw ∂³
+}
+
+struct JetCtx {
+    n: usize,
+    n2: usize,
+    n3: usize,
+    t3: Vec<[usize; 3]>,
+    t3idx: Vec<Vec<Vec<usize>>>,
+}
+
+impl JetCtx {
+    fn new(n: usize) -> Self {
+        let (t3, t3idx) = i3_table(n);
+        JetCtx { n, n2: n * (n + 1) / 2, n3: t3.len(), t3, t3idx }
+    }
+    fn zero(&self) -> Jet {
+        Jet {
+            n: self.n,
+            v: Q::zero(),
+            d1: vec![Q::zero(); self.n],
+            d2: vec![Q::zero(); self.n2],
+            d3: vec![Q::zero(); self.n3],
+        }
+    }
+    fn constant(&self, c: Q) -> Jet {
+        let mut j = self.zero();
+        j.v = c;
+        j
+    }
+    fn coord(&self, i: usize, val: Q) -> Jet {
+        let mut j = self.zero();
+        j.v = val;
+        j.d1[i] = Q::one();
+        j
+    }
+    fn add(&self, a: &Jet, b: &Jet) -> Jet {
+        let mut c = a.clone();
+        c.v += &b.v;
+        for i in 0..self.n { c.d1[i] += &b.d1[i]; }
+        for i in 0..self.n2 { c.d2[i] += &b.d2[i]; }
+        for i in 0..self.n3 { c.d3[i] += &b.d3[i]; }
+        c
+    }
+    fn sub(&self, a: &Jet, b: &Jet) -> Jet {
+        let mut c = a.clone();
+        c.v -= &b.v;
+        for i in 0..self.n { c.d1[i] -= &b.d1[i]; }
+        for i in 0..self.n2 { c.d2[i] -= &b.d2[i]; }
+        for i in 0..self.n3 { c.d3[i] -= &b.d3[i]; }
+        c
+    }
+    fn neg(&self, a: &Jet) -> Jet {
+        self.sub(&self.zero(), a)
+    }
+    fn mul(&self, a: &Jet, b: &Jet) -> Jet {
+        let n = self.n;
+        let mut c = self.zero();
+        c.v = &a.v * &b.v;
+        for i in 0..n {
+            c.d1[i] = &a.d1[i] * &b.v + &a.v * &b.d1[i];
+        }
+        for i in 0..n {
+            for j in i..n {
+                c.d2[i2(n, i, j)] = &a.d2[i2(n, i, j)] * &b.v
+                    + &a.d1[i] * &b.d1[j]
+                    + &a.d1[j] * &b.d1[i]
+                    + &a.v * &b.d2[i2(n, i, j)];
+            }
+        }
+        for (id, t) in self.t3.iter().enumerate() {
+            let [i, j, k] = *t;
+            c.d3[id] = &a.d3[id] * &b.v
+                + &a.d2[i2(n, i, j)] * &b.d1[k]
+                + &a.d2[i2(n, i, k)] * &b.d1[j]
+                + &a.d2[i2(n, j, k)] * &b.d1[i]
+                + &a.d1[i] * &b.d2[i2(n, j, k)]
+                + &a.d1[j] * &b.d2[i2(n, i, k)]
+                + &a.d1[k] * &b.d2[i2(n, i, j)]
+                + &a.v * &b.d3[id];
+        }
+        c
+    }
+    /// reciprocal 1/b, solving b·r = 1 order by order
+    fn inv(&self, b: &Jet) -> Jet {
+        let n = self.n;
+        assert!(!b.v.is_zero(), "jet reciprocal at pole");
+        let mut r = self.zero();
+        r.v = Q::one() / &b.v;
+        for i in 0..n {
+            r.d1[i] = -(&b.d1[i] * &r.v) / &b.v;
+        }
+        for i in 0..n {
+            for j in i..n {
+                let s = &b.d2[i2(n, i, j)] * &r.v
+                    + &b.d1[i] * &r.d1[j]
+                    + &b.d1[j] * &r.d1[i];
+                r.d2[i2(n, i, j)] = -s / &b.v;
+            }
+        }
+        for (id, t) in self.t3.iter().enumerate() {
+            let [i, j, k] = *t;
+            let s = &b.d3[id] * &r.v
+                + &b.d2[i2(n, i, j)] * &r.d1[k]
+                + &b.d2[i2(n, i, k)] * &r.d1[j]
+                + &b.d2[i2(n, j, k)] * &r.d1[i]
+                + &b.d1[i] * &r.d2[i2(n, j, k)]
+                + &b.d1[j] * &r.d2[i2(n, i, k)]
+                + &b.d1[k] * &r.d2[i2(n, i, j)];
+            r.d3[id] = -s / &b.v;
+        }
+        r
+    }
+    fn div(&self, a: &Jet, b: &Jet) -> Jet {
+        self.mul(a, &self.inv(b))
+    }
+    /// ∂_dir shift: returns jet valid one order lower (d3 of result zeroed)
+    fn shift(&self, a: &Jet, dir: usize) -> Jet {
+        let n = self.n;
+        let mut s = self.zero();
+        s.v = a.d1[dir].clone();
+        for i in 0..n {
+            s.d1[i] = a.d2[i2(n, dir, i)].clone();
+        }
+        for i in 0..n {
+            for j in i..n {
+                s.d2[i2(n, i, j)] = a.d3[self.t3idx[dir][i][j]].clone();
+            }
+        }
+        s
+    }
+}
+
+// ── receipts ─────────────────────────────────────────────────────────
+
+#[derive(Serialize)]
+struct Receipt {
+    schema: &'static str,
+    m: usize,
+    geometrization: String,
+    point: Vec<String>,
+    beta_reduced: Vec<String>,
+    baseline_constant_curvature_quarter: bool,
+    gprime_signature: (usize, usize),
+    gprime_beta_norm: String,
+    nabla_r_nonzero: usize,
+    nabla_r_total: usize,
+    nabla_r_sample: Option<(Vec<usize>, String)>,
+    rr_nonzero: usize,
+    pseudo_l_candidates: Vec<String>,
+    verdict: String,
+    obstruction_along_drift: String,
+}
+
+fn qs(x: &Q) -> String {
+    x.to_string()
+}
+
+// ── geometry pipeline at a rational point ────────────────────────────
+
+fn drift_beta(m: usize) -> Vec<Q> {
+    // BraidStateN.lean rossbyDriftFromChirality weights, cycled
+    let base = [q(1, 1), q(-1, 1), q(1, 2), q(0, 1)];
+    let w: Vec<Q> = (0..m).map(|i| base[i % 4].clone()).collect();
+    let mean: Q = w.iter().fold(Q::zero(), |a, b| a + b) / Q::from(BigInt::from(m as i64));
+    w.iter().map(|wi| wi - &mean).collect()
+}
+
+fn probe(m: usize, point: &[Q]) -> Receipt {
+    let n = m - 1;
+    let cx = JetCtx::new(n);
+    let ps: Vec<Jet> = (0..n).map(|i| cx.coord(i, point[i].clone())).collect();
+    let mut pm = cx.constant(Q::one());
+    for p in &ps {
+        pm = cx.sub(&pm, p);
+    }
+    assert!(pm.v > Q::zero(), "point outside simplex");
+
+    // Fisher–Rao g
+    let inv_pm = cx.inv(&pm);
+    let mut g: Vec<Vec<Jet>> = vec![vec![cx.zero(); n]; n];
+    for i in 0..n {
+        for j in 0..n {
+            g[i][j] = if i == j {
+                cx.add(&cx.inv(&ps[i]), &inv_pm)
+            } else {
+                inv_pm.clone()
+            };
+        }
+    }
+
+    // baseline check at this point: R = 1/4 (g∧g) — via full pipeline on g
+    let base_ok = check_constant_curvature(&cx, &g, n);
+
+    // drift-flip g'
+    let beta_full = drift_beta(m);
+    let beta: Vec<Q> = beta_full[..n].to_vec();
+    let bjets: Vec<Jet> = beta.iter().map(|b| cx.constant(b.clone())).collect();
+    let mut beta_flat: Vec<Jet> = vec![cx.zero(); n];
+    for i in 0..n {
+        for j in 0..n {
+            beta_flat[i] = cx.add(&beta_flat[i], &cx.mul(&g[i][j], &bjets[j]));
+        }
+    }
+    let mut bnorm = cx.zero();
+    for i in 0..n {
+        bnorm = cx.add(&bnorm, &cx.mul(&bjets[i], &beta_flat[i]));
+    }
+    let two = cx.constant(q(2, 1));
+    let mut gp: Vec<Vec<Jet>> = vec![vec![cx.zero(); n]; n];
+    for i in 0..n {
+        for j in 0..n {
+            let corr = cx.div(
+                &cx.mul(&two, &cx.mul(&beta_flat[i], &beta_flat[j])),
+                &bnorm,
+            );
+            gp[i][j] = cx.sub(&g[i][j], &corr);
+        }
+    }
+
+    // signature via leading principal minors (values only)
+    let gm: Vec<Vec<Q>> = (0..n)
+        .map(|i| (0..n).map(|j| gp[i][j].v.clone()).collect())
+        .collect();
+    let sig = signature(&gm);
+
+    // g'(β,β)
+    let mut bb = Q::zero();
+    for i in 0..n {
+        for j in 0..n {
+            bb += &gm[i][j] * &beta[i] * &beta[j];
+        }
+    }
+
+    // curvature pipeline on g'
+    let (rup, rdn, gamma) = curvature(&cx, &gp, n);
+
+    // ∇R (values, exact)
+    let mut nabla_nz = 0usize;
+    let mut nabla_sample: Option<(Vec<usize>, String)> = None;
+    let mut drift_contract = Q::zero();
+    for a in 0..n {
+        for i in 0..n {
+            for jj in 0..n {
+                for k in 0..n {
+                    for l in 0..n {
+                        // ∂_a R_ijkl from R jets (order ≥1) minus Γ corrections
+                        let mut v = rdn[i][jj][k][l].d1[a].clone();
+                        for mm in 0..n {
+                            v -= &gamma[mm][a][i].v * &rdn[mm][jj][k][l].v;
+                            v -= &gamma[mm][a][jj].v * &rdn[i][mm][k][l].v;
+                            v -= &gamma[mm][a][k].v * &rdn[i][jj][mm][l].v;
+                            v -= &gamma[mm][a][l].v * &rdn[i][jj][k][mm].v;
+                        }
+                        if !v.is_zero() {
+                            nabla_nz += 1;
+                            if nabla_sample.is_none() {
+                                nabla_sample =
+                                    Some((vec![a, i, jj, k, l], qs(&v)));
+                            }
+                        }
+                        let term = &beta[a] * &v;
+                        drift_contract += &term * &term;
+                    }
+                }
+            }
+        }
+    }
+
+    // R·R and Q(g,R), pseudosymmetry ratios
+    let mut rr_nz = 0usize;
+    let mut ratios: Vec<Q> = Vec::new();
+    let mut proper = false;
+    for a in 0..n {
+        for b in (a + 1)..n {
+            for i in 0..n {
+                for jj in 0..n {
+                    for k in 0..n {
+                        for l in 0..n {
+                            let mut rr = Q::zero();
+                            let mut qq = Q::zero();
+                            for mm in 0..n {
+                                let a_r = &rup[mm][i][a][b].v;
+                                rr -= a_r * &rdn[mm][jj][k][l].v;
+                                let a_r = &rup[mm][jj][a][b].v;
+                                rr -= a_r * &rdn[i][mm][k][l].v;
+                                let a_r = &rup[mm][k][a][b].v;
+                                rr -= a_r * &rdn[i][jj][mm][l].v;
+                                let a_r = &rup[mm][l][a][b].v;
+                                rr -= a_r * &rdn[i][jj][k][mm].v;
+                            }
+                            let wedge = |mm: usize, idx: usize| -> Q {
+                                let mut w = Q::zero();
+                                if mm == a { w += &gm[b][idx]; }
+                                if mm == b { w -= &gm[a][idx]; }
+                                w
+                            };
+                            for mm in 0..n {
+                                qq -= wedge(mm, i) * &rdn[mm][jj][k][l].v;
+                                qq -= wedge(mm, jj) * &rdn[i][mm][k][l].v;
+                                qq -= wedge(mm, k) * &rdn[i][jj][mm][l].v;
+                                qq -= wedge(mm, l) * &rdn[i][jj][k][mm].v;
+                            }
+                            if !rr.is_zero() {
+                                rr_nz += 1;
+                                if qq.is_zero() {
+                                    proper = true;
+                                } else {
+                                    let r = rr / qq;
+                                    if !ratios.contains(&r) {
+                                        ratios.push(r);
+                                    }
+                                }
+                            } else if !qq.is_zero() {
+                                let z = Q::zero();
+                                if !ratios.contains(&z) {
+                                    ratios.push(z);
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    let verdict = if nabla_nz == 0 {
+        "LOCALLY_SYMMETRIC".to_string()
+    } else if rr_nz == 0 {
+        "SEMISYMMETRIC_AT_POINT".to_string()
+    } else if !proper && ratios.len() == 1 {
+        format!("PSEUDOSYMMETRIC_AT_POINT L={}", qs(&ratios[0]))
+    } else {
+        "PROPER".to_string()
+    };
+
+    Receipt {
+        schema: "semisym_receipt_v1",
+        m,
+        geometrization: "drift_flip".into(),
+        point: point.iter().map(qs).collect(),
+        beta_reduced: beta.iter().map(qs).collect(),
+        baseline_constant_curvature_quarter: base_ok,
+        gprime_signature: sig,
+        gprime_beta_norm: qs(&bb),
+        nabla_r_nonzero: nabla_nz,
+        nabla_r_total: n.pow(5),
+        nabla_r_sample: nabla_sample,
+        rr_nonzero: rr_nz,
+        pseudo_l_candidates: ratios.iter().map(qs).collect(),
+        verdict,
+        obstruction_along_drift: qs(&drift_contract),
+    }
+}
+
+/// Γ (order-2 jets), R^ρ_σμν and R_ρσμν (order-1 jets)
+#[allow(clippy::type_complexity)]
+fn curvature(
+    cx: &JetCtx,
+    g: &[Vec<Jet>],
+    n: usize,
+) -> (Vec<Vec<Vec<Vec<Jet>>>>, Vec<Vec<Vec<Vec<Jet>>>>, Vec<Vec<Vec<Jet>>>) {
+    let ginv = mat_inv(cx, g, n);
+    let half = cx.constant(q(1, 2));
+    // Γ^l_ij
+    let mut gamma = vec![vec![vec![cx.zero(); n]; n]; n];
+    for l in 0..n {
+        for i in 0..n {
+            for j in i..n {
+                let mut s = cx.zero();
+                for k in 0..n {
+                    let t = cx.sub(
+                        &cx.add(&cx.shift(&g[k][j], i), &cx.shift(&g[k][i], j)),
+                        &cx.shift(&g[i][j], k),
+                    );
+                    s = cx.add(&s, &cx.mul(&ginv[l][k], &t));
+                }
+                let v = cx.mul(&half, &s);
+                gamma[l][i][j] = v.clone();
+                gamma[l][j][i] = v;
+            }
+        }
+    }
+    // R^ρ_σμν
+    let mut rup = vec![vec![vec![vec![cx.zero(); n]; n]; n]; n];
+    for rho in 0..n {
+        for sg in 0..n {
+            for mu in 0..n {
+                for nu in (mu + 1)..n {
+                    let mut t = cx.sub(
+                        &cx.shift(&gamma[rho][nu][sg], mu),
+                        &cx.shift(&gamma[rho][mu][sg], nu),
+                    );
+                    for lam in 0..n {
+                        t = cx.add(&t, &cx.mul(&gamma[rho][mu][lam], &gamma[lam][nu][sg]));
+                        t = cx.sub(&t, &cx.mul(&gamma[rho][nu][lam], &gamma[lam][mu][sg]));
+                    }
+                    rup[rho][sg][mu][nu] = t.clone();
+                    rup[rho][sg][nu][mu] = cx.neg(&t);
+                }
+            }
+        }
+    }
+    // lower first index
+    let mut rdn = vec![vec![vec![vec![cx.zero(); n]; n]; n]; n];
+    for rho in 0..n {
+        for sg in 0..n {
+            for mu in 0..n {
+                for nu in (mu + 1)..n {
+                    let mut s = cx.zero();
+                    for lam in 0..n {
+                        s = cx.add(&s, &cx.mul(&g[rho][lam], &rup[lam][sg][mu][nu]));
+                    }
+                    rdn[rho][sg][mu][nu] = s.clone();
+                    rdn[rho][sg][nu][mu] = cx.neg(&s);
+                }
+            }
+        }
+    }
+    (rup, rdn, gamma)
+}
+
+fn check_constant_curvature(cx: &JetCtx, g: &[Vec<Jet>], n: usize) -> bool {
+    let (_, rdn, _) = curvature(cx, g, n);
+    let quarter = q(1, 4);
+    for i in 0..n {
+        for j in 0..n {
+            for k in 0..n {
+                for l in 0..n {
+                    let expect =
+                        &quarter * (&g[i][k].v * &g[j][l].v - &g[i][l].v * &g[j][k].v);
+                    if rdn[i][j][k][l].v != expect {
+                        return false;
+                    }
+                }
+            }
+        }
+    }
+    true
+}
+
+fn mat_inv(cx: &JetCtx, a: &[Vec<Jet>], n: usize) -> Vec<Vec<Jet>> {
+    // Gauss-Jordan with jet arithmetic; exact, pivots nonzero for metrics
+    let mut m: Vec<Vec<Jet>> = a.to_vec();
+    let mut inv: Vec<Vec<Jet>> = (0..n)
+        .map(|i| {
+            (0..n)
+                .map(|j| if i == j { cx.constant(Q::one()) } else { cx.zero() })
+                .collect()
+        })
+        .collect();
+    for col in 0..n {
+        // pivot: find row with nonzero value
+        let piv = (col..n)
+            .find(|&r| !m[r][col].v.is_zero())
+            .expect("singular jet matrix");
+        m.swap(col, piv);
+        inv.swap(col, piv);
+        let pinv = cx.inv(&m[col][col]);
+        for j in 0..n {
+            m[col][j] = cx.mul(&m[col][j], &pinv);
+            inv[col][j] = cx.mul(&inv[col][j], &pinv);
+        }
+        for r in 0..n {
+            if r != col && !(m[r][col].v.is_zero()
+                && m[r][col].d1.iter().all(|x| x.is_zero())
+                && m[r][col].d2.iter().all(|x| x.is_zero())
+                && m[r][col].d3.iter().all(|x| x.is_zero()))
+            {
+                let f = m[r][col].clone();
+                for j in 0..n {
+                    let t = cx.mul(&f, &m[col][j]);
+                    m[r][j] = cx.sub(&m[r][j], &t);
+                    let t = cx.mul(&f, &inv[col][j]);
+                    inv[r][j] = cx.sub(&inv[r][j], &t);
+                }
+            }
+        }
+    }
+    inv
+}
+
+fn signature(gm: &[Vec<Q>]) -> (usize, usize) {
+    // Exact signature by symmetric (congruence) elimination — robust to
+    // vanishing leading minors, which Sylvester's criterion is not.
+    // Congruences A ↦ EᵀAE preserve signature (Sylvester's law of inertia).
+    let n = gm.len();
+    let mut a: Vec<Vec<Q>> = gm.to_vec();
+    let (mut pos, mut neg) = (0usize, 0usize);
+    for k in 0..n {
+        if a[k][k].is_zero() {
+            // symmetric pivot repair: prefer swapping in a nonzero diagonal
+            if let Some(j) = ((k + 1)..n).find(|&j| !a[j][j].is_zero()) {
+                a.swap(k, j);
+                for row in a.iter_mut() {
+                    row.swap(k, j);
+                }
+            } else if let Some(j) = ((k + 1)..n).find(|&j| !a[k][j].is_zero()) {
+                // row/col addition: (k) += (j) symmetrically ⇒ a[k][k] = 2a[k][j] ≠ 0
+                for t in k..n {
+                    let v = a[j][t].clone();
+                    a[k][t] += v;
+                }
+                for t in k..n {
+                    let v = a[t][j].clone();
+                    a[t][k] += v;
+                }
+            } else {
+                panic!("degenerate metric block (rank deficiency)");
+            }
+        }
+        if a[k][k] > Q::zero() {
+            pos += 1;
+        } else {
+            neg += 1;
+        }
+        let pv = a[k][k].clone();
+        // snapshot row k: zeroing a[k][i] mid-loop must not corrupt the
+        // updates for later rows (this exact aliasing bug produced a wrong
+        // (4,3) signature at the Δ₇ centroid; truth is (6,1))
+        let row_k: Vec<Q> = a[k].clone();
+        for i in (k + 1)..n {
+            let f = &a[i][k] / &pv;
+            for j in (k + 1)..n {
+                let t = &f * &row_k[j];
+                a[i][j] -= t;
+            }
+            a[i][k] = Q::zero();
+            a[k][i] = Q::zero();
+        }
+    }
+    (pos, neg)
+}
+
+fn det(a: &[Vec<Q>]) -> Q {
+    let n = a.len();
+    let mut m = a.to_vec();
+    let mut d = Q::one();
+    for col in 0..n {
+        let piv = match (col..n).find(|&r| !m[r][col].is_zero()) {
+            Some(p) => p,
+            None => return Q::zero(),
+        };
+        if piv != col {
+            m.swap(col, piv);
+            d = -d;
+        }
+        d *= m[col][col].clone();
+        let pv = m[col][col].clone();
+        for r in (col + 1)..n {
+            let f = &m[r][col] / &pv;
+            for j in col..n {
+                let t = &f * &m[col][j];
+                m[r][j] -= t;
+            }
+        }
+    }
+    d
+}
+
+// ── lattice survey with checkpoint/resume ────────────────────────────
+
+fn lattice_points(m: usize, grid: u64) -> Vec<Vec<Q>> {
+    // interior compositions: p_i = c_i / (grid*m) with c_i ≥ 1, Σ c = grid*m
+    // sampled coarsely: c_i ∈ {1..} via stars-and-bars over a small grid
+    let n = m - 1;
+    let total = grid * m as u64;
+    let mut pts = Vec::new();
+    fn rec(n: usize, left: u64, cur: &mut Vec<u64>, out: &mut Vec<Vec<u64>>) {
+        if n == 0 {
+            if left >= 1 {
+                let mut c = cur.clone();
+                c.push(left);
+                out.push(c);
+            }
+            return;
+        }
+        for c in 1..=(left.saturating_sub(n as u64)) {
+            cur.push(c);
+            rec(n - 1, left - c, cur, out);
+            cur.pop();
+        }
+    }
+    let mut combos = Vec::new();
+    rec(n, total, &mut Vec::new(), &mut combos);
+    for c in combos {
+        pts.push(
+            c[..n]
+                .iter()
+                .map(|&ci| Q::new(BigInt::from(ci), BigInt::from(total)))
+                .collect(),
+        );
+    }
+    pts
+}
+
+fn main() {
+    let args: Vec<String> = std::env::args().collect();
+    let cmd = args.get(1).map(String::as_str).unwrap_or("probe");
+    let mut m = 4usize;
+    let mut point: Option<Vec<Q>> = None;
+    let mut grid = 2u64;
+    let mut ckpt = PathBuf::from("semisym_ckpt");
+    let mut max_points = 200usize;
+    let mut i = 2;
+    while i < args.len() {
+        match args[i].as_str() {
+            "--m" => { m = args[i + 1].parse().unwrap(); i += 2; }
+            "--point" => {
+                point = Some(args[i + 1].split(',').map(parse_q).collect());
+                i += 2;
+            }
+            "--grid" => { grid = args[i + 1].parse().unwrap(); i += 2; }
+            "--checkpoint" => { ckpt = PathBuf::from(&args[i + 1]); i += 2; }
+            "--max-points" => { max_points = args[i + 1].parse().unwrap(); i += 2; }
+            other => panic!("unknown arg {other}"),
+        }
+    }
+
+    match cmd {
+        "probe" => {
+            let n = m - 1;
+            let pt = point.unwrap_or_else(|| {
+                (0..n).map(|_| Q::new(BigInt::one(), BigInt::from(m as i64))).collect()
+            });
+            let r = probe(m, &pt);
+            println!("{}", serde_json::to_string_pretty(&r).unwrap());
+        }
+        "survey" => {
+            fs::create_dir_all(&ckpt).unwrap();
+            let mut pts = lattice_points(m, grid);
+            if pts.len() > max_points {
+                // deterministic stride subsample — same DAG nodes every run,
+                // so resume-by-receipt stays stable
+                let stride = pts.len() / max_points;
+                pts = pts.into_iter().step_by(stride.max(1)).take(max_points).collect();
+            }
+            eprintln!("survey: {} lattice points, checkpoint {}", pts.len(), ckpt.display());
+            let done: Vec<(usize, String)> = pts
+                .par_iter()
+                .enumerate()
+                .map(|(idx, pt)| {
+                    let mut name = format!("m{m}_g{grid}_");
+                    for x in pt {
+                        let _ = write!(name, "{}_", x.numer());
+                    }
+                    let file = ckpt.join(format!("{name}.json"));
+                    if file.exists() {
+                        return (idx, "skip".to_string());
+                    }
+                    let r = probe(m, pt);
+                    let v = r.verdict.clone();
+                    fs::write(&file, serde_json::to_string_pretty(&r).unwrap()).unwrap();
+                    (idx, v)
+                })
+                .collect();
+            let mut counts = std::collections::BTreeMap::new();
+            for (_, v) in &done {
+                *counts.entry(v.clone()).or_insert(0usize) += 1;
+            }
+            println!("{counts:?}");
+        }
+        other => panic!("unknown command {other}"),
+    }
+}