From 6c58feaf53643ef136a70083348609ae90b9fcaa Mon Sep 17 00:00:00 2001 From: allaun Date: Thu, 2 Jul 2026 02:43:54 -0500 Subject: [PATCH] =?UTF-8?q?feat(semisym):=20exact-rational=20Rust=20discri?= =?UTF-8?q?minator=20=E2=80=94=20jet=20AD,=20checkpointed=20survey,=20?= =?UTF-8?q?=CE=94=E2=82=87=20verdict?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 --- rust/semisym/Cargo.toml | 15 + rust/semisym/src/main.rs | 757 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 772 insertions(+) create mode 100644 rust/semisym/Cargo.toml create mode 100644 rust/semisym/src/main.rs diff --git a/rust/semisym/Cargo.toml b/rust/semisym/Cargo.toml new file mode 100644 index 00000000..33386939 --- /dev/null +++ b/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 diff --git a/rust/semisym/src/main.rs b/rust/semisym/src/main.rs new file mode 100644 index 00000000..8d470ce0 --- /dev/null +++ b/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>>) { + 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, // n + d2: Vec, // n(n+1)/2, raw ∂² + d3: Vec, // C(n+2,3), raw ∂³ +} + +struct JetCtx { + n: usize, + n2: usize, + n3: usize, + t3: Vec<[usize; 3]>, + t3idx: Vec>>, +} + +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, + beta_reduced: Vec, + 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, String)>, + rr_nonzero: usize, + pseudo_l_candidates: Vec, + 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 { + // BraidStateN.lean rossbyDriftFromChirality weights, cycled + let base = [q(1, 1), q(-1, 1), q(1, 2), q(0, 1)]; + let w: Vec = (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 = (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![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 = beta_full[..n].to_vec(); + let bjets: Vec = beta.iter().map(|b| cx.constant(b.clone())).collect(); + let mut beta_flat: Vec = 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![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> = (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, 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 = 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], + n: usize, +) -> (Vec>>>, Vec>>>, Vec>>) { + 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], 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], n: usize) -> Vec> { + // Gauss-Jordan with jet arithmetic; exact, pivots nonzero for metrics + let mut m: Vec> = a.to_vec(); + let mut inv: Vec> = (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]) -> (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> = 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 = 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 { + 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> { + // 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, out: &mut Vec>) { + 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 = std::env::args().collect(); + let cmd = args.get(1).map(String::as_str).unwrap_or("probe"); + let mut m = 4usize; + let mut point: Option> = 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}"), + } +}