//! 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;

mod modp;

pub(crate) 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]
pub(crate) 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
}

pub(crate) 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)]
pub(crate) 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 ────────────────────────────

pub(crate) 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()
}

pub(crate) 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
}

pub(crate) fn signature(gm: &[Vec<Q>]) -> (usize, usize) {
    // Panicking wrapper for the exact/oracle path, where a degenerate g'
    // block is not expected. The robust modular path uses signature_opt and
    // records an edge receipt instead of crashing.
    signature_opt(gm).expect("degenerate metric block (rank deficiency)")
}

/// 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). Returns None (rather than
/// panicking) when the block is rank-deficient / degenerate — a genuine
/// geometric edge (e.g. a null drift direction) that callers record explicitly.
pub(crate) fn signature_opt(gm: &[Vec<Q>]) -> Option<(usize, usize)> {
    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 {
                return None; // rank-deficient block (e.g. null drift) — edge case
            }
        }
        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();
        }
    }
    Some((pos, neg))
}

#[allow(dead_code)]
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 primes: Option<Vec<u64>> = None; // Some ⇒ modular fast path
    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; }
            "--mod" => {
                // bare `--mod` uses the default multi-prime set; `--mod p1,p2,..`
                // overrides (each must be prime — guarded below).
                let list = args.get(i + 1).and_then(|s| {
                    s.split(',').map(|x| x.parse::<u64>()).collect::<Result<Vec<u64>, _>>().ok()
                });
                match list {
                    Some(l) if !l.is_empty() => { primes = Some(l); i += 2; }
                    _ => { primes = Some(modp::DEFAULT_PRIMES.to_vec()); i += 1; }
                }
            }
            other => panic!("unknown arg {other}"),
        }
    }

    // primality guard: a composite modulus makes the Fermat inverse — and every
    // modular "proof" — silently unsound, so reject it up front (edge accounted
    // for, not assumed away).
    if let Some(ref ps) = primes {
        for &p in ps {
            assert!(
                modp::is_prime(p),
                "--mod modulus {p} is not prime; modular soundness requires a prime field"
            );
        }
    }

    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()
            });
            match &primes {
                Some(ps) => {
                    let (_verdict, json) = modp::probe_survey(m, &pt, ps);
                    println!("{json}");
                }
                None => {
                    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();
            }
            let tag = match &primes {
                Some(ps) => format!("modular (primes={ps:?})"),
                None => "exact ℚ".into(),
            };
            eprintln!(
                "survey [{tag}]: {} 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 {
                        // encode numer AND denom: reduced numerators alone
                        // collide (1/16,1/8,1/4,1/2 all → "1"), silently folding
                        // distinct lattice points onto one checkpoint file
                        let _ = write!(name, "{}-{}_", x.numer(), x.denom());
                    }
                    let file = ckpt.join(format!("{name}.json"));
                    if file.exists() {
                        return (idx, "skip".to_string());
                    }
                    // catch_unwind backstop: an UNFORESEEN panic on one point
                    // must never abort the whole survey (which would discard all
                    // in-flight rayon work). Assume it happens; contain it.
                    let compute = std::panic::AssertUnwindSafe(|| match &primes {
                        Some(ps) => modp::probe_survey(m, pt, ps),
                        None => {
                            let r = probe(m, pt);
                            (r.verdict.clone(), serde_json::to_string_pretty(&r).unwrap())
                        }
                    });
                    let (v, json) = match std::panic::catch_unwind(compute) {
                        Ok(x) => x,
                        Err(_) => {
                            let pt_s: Vec<String> = pt.iter().map(|x| x.to_string()).collect();
                            let e = serde_json::json!({
                                "schema": "semisym_receipt_edge_v1",
                                "m": m,
                                "point": pt_s,
                                "verdict": "PANIC_CAUGHT",
                                "reason": "unforeseen_panic",
                                "note": "contained by survey backstop; needs manual review",
                            });
                            ("PANIC_CAUGHT".to_string(), serde_json::to_string_pretty(&e).unwrap())
                        }
                    };
                    fs::write(&file, json).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:?}");
        }
        // Cross-check the modular path against the exact-ℚ receipt oracle:
        // for every exact receipt in --checkpoint, rerun modp at that point and
        // assert method agreement (different arithmetic, so this is a genuine
        // independent check, not the same-algorithm trap that hid the old
        // signature bug). Prints any mismatch; exits nonzero if any found.
        "validate-modp" => {
            let p = primes
                .as_ref()
                .and_then(|v| v.first().copied())
                .unwrap_or(modp::DEFAULT_PRIMES[0]);
            let entries: Vec<PathBuf> = fs::read_dir(&ckpt)
                .unwrap_or_else(|_| panic!("no checkpoint dir {}", ckpt.display()))
                .filter_map(|e| e.ok().map(|e| e.path()))
                .filter(|p| p.extension().is_some_and(|x| x == "json"))
                .collect();
            let mut checked = 0usize;
            let mut mismatches = 0usize;
            let mut poles = 0usize;
            for path in &entries {
                let txt = fs::read_to_string(path).unwrap();
                let v: serde_json::Value = serde_json::from_str(&txt).unwrap();
                // only compare against EXACT receipts
                if v.get("schema").and_then(|s| s.as_str()) != Some("semisym_receipt_v1") {
                    continue;
                }
                let mm = v["m"].as_u64().unwrap() as usize;
                let pt: Vec<Q> = v["point"].as_array().unwrap().iter()
                    .map(|s| parse_q(s.as_str().unwrap())).collect();
                let exact_verdict = v["verdict"].as_str().unwrap();
                let exact_nabla = v["nabla_r_nonzero"].as_u64().unwrap() as usize;
                let exact_sig = (
                    v["gprime_signature"][0].as_u64().unwrap() as usize,
                    v["gprime_signature"][1].as_u64().unwrap() as usize,
                );
                let r = match modp::probe_modp(mm, &pt, p) {
                    Ok(r) => r,
                    Err(reason) => {
                        // exact receipt exists ⇒ point is interior/non-degenerate,
                        // so a geometric edge here is a genuine disagreement
                        mismatches += 1;
                        eprintln!("MISMATCH {}: modp edge '{reason}' but exact receipt present", path.display());
                        continue;
                    }
                };
                checked += 1;
                if r.modular_pole {
                    // expected edge for a single prime, not a correctness failure
                    poles += 1;
                    continue;
                }
                let mut problems = Vec::new();
                if r.gprime_signature != exact_sig {
                    problems.push(format!(
                        "signature exact {exact_sig:?} vs modp {:?}",
                        r.gprime_signature
                    ));
                }
                // modp certified-nonzero must be a subset of exact nonzero
                if r.nabla_r_certified_nonzero > exact_nabla {
                    problems.push(format!(
                        "∇R certified {} > exact nonzero {}",
                        r.nabla_r_certified_nonzero, exact_nabla
                    ));
                }
                // if exact says PROPER, modp (sound) must not disagree
                if exact_verdict == "PROPER" && r.sound && r.verdict != "PROPER" {
                    problems.push(format!("exact PROPER vs modp {}", r.verdict));
                }
                if !problems.is_empty() {
                    mismatches += 1;
                    eprintln!("MISMATCH {}: {}", path.display(), problems.join("; "));
                }
            }
            println!("validate-modp: checked {checked} exact receipts, {mismatches} mismatches, {poles} single-prime poles (p={p})");
            if mismatches > 0 {
                std::process::exit(1);
            }
        }
        other => panic!("unknown command {other}"),
    }
}
