//! modp — modular fast path for the §0 covariant semi-symmetry discriminator.
//!
//! Same order-3 jet AD pipeline as the exact `main.rs`, but the scalar is a
//! finite field 𝔽_p (u64 mod a large prime) instead of BigRational. This kills
//! the denominator blowup that makes off-centre lattice points take 30+ minutes
//! in exact ℚ: 𝔽_p elements are single machine words, arithmetic is O(1), and
//! memory is flat. Typical speedup is 100–1000× per point.
//!
//! SOUNDNESS (the OTOM gate is preserved for the verdicts we actually get):
//!   * A curvature/∇R component ≢ 0 mod p is a *proof* it is ≠ 0 over ℚ
//!     (p ∤ numerator). So "∇R ≠ 0" (not locally symmetric) is a hard proof.
//!   * PROPER (not pseudosymmetric) is certified WITHOUT division: pick a
//!     reference pair (rr₀,qq₀) with qq₀ ≠ 0; if any pair has
//!     rr·qq₀ − rr₀·qq ≢ 0 mod p, then R·R ≠ L·Q(g,R) for the only possible
//!     L = rr₀/qq₀ — a hard proof over ℚ. If all Q-side entries vanish but some
//!     R·R entry is ≢ 0, that is also PROPER (no finite L works).
//!   * The ONLY verdicts needing an exact recheck are the all-zero ones
//!     (candidate locally-symmetric / semisymmetric / pseudosymmetric): there a
//!     modular zero could be a false zero (p | numerator). We have never hit one
//!     in the exact survey, and they escalate to the ℚ path automatically.
//!
//! Signature and g'(β,β) stay in exact ℚ (they need order/sign, meaningless in
//! 𝔽_p) — but they are an n×n value-only computation, not the bottleneck.
//!
//! A modular pivot may vanish (p | a nonzero-over-ℚ pivot) → probe returns a
//! Pole; the caller retries with another prime or falls back to exact.

use num_bigint::BigInt;
use num_rational::BigRational;
use num_traits::Zero;

use crate::{drift_beta, i2, i3_table, signature_opt, Q};

// ── 𝔽_p, p = 2³¹ − 1 (Mersenne prime; products fit in u64) ────────────

pub const DEFAULT_PRIME: u64 = 2_147_483_647; // 2^31 - 1

#[derive(Clone, Copy, PartialEq, Eq)]
pub struct Fp {
    v: u64,
    p: u64,
}

impl Fp {
    #[inline]
    fn new(v: u64, p: u64) -> Fp {
        Fp { v: v % p, p }
    }
    #[inline]
    fn zero(p: u64) -> Fp {
        Fp { v: 0, p }
    }
    #[inline]
    fn one(p: u64) -> Fp {
        Fp { v: 1, p }
    }
    #[inline]
    fn is_zero(self) -> bool {
        self.v == 0
    }
    #[inline]
    fn add(self, o: Fp) -> Fp {
        let s = self.v + o.v;
        Fp { v: if s >= self.p { s - self.p } else { s }, p: self.p }
    }
    #[inline]
    fn sub(self, o: Fp) -> Fp {
        Fp {
            v: if self.v >= o.v { self.v - o.v } else { self.v + self.p - o.v },
            p: self.p,
        }
    }
    #[inline]
    fn mul(self, o: Fp) -> Fp {
        // p < 2^31 ⇒ product < 2^62 < 2^64, no overflow
        Fp { v: (self.v * o.v) % self.p, p: self.p }
    }
    #[inline]
    fn neg(self) -> Fp {
        Fp { v: if self.v == 0 { 0 } else { self.p - self.v }, p: self.p }
    }
    fn pow(self, mut e: u64) -> Fp {
        let mut r = Fp::one(self.p);
        let mut b = self;
        while e > 0 {
            if e & 1 == 1 {
                r = r.mul(b);
            }
            b = b.mul(b);
            e >>= 1;
        }
        r
    }
    /// Fermat inverse; returns None at a modular pole (self ≡ 0).
    fn inv(self) -> Option<Fp> {
        if self.is_zero() {
            None
        } else {
            Some(self.pow(self.p - 2))
        }
    }
    fn from_i64(x: i64, p: u64) -> Fp {
        let m = x.rem_euclid(p as i64);
        Fp { v: m as u64, p }
    }
    /// Reduce a BigRational mod p; None if p | denominator.
    fn from_q(qr: &BigRational, p: u64) -> Option<Fp> {
        let n = bigint_mod(qr.numer(), p);
        let d = bigint_mod(qr.denom(), p);
        Some(Fp::new(n, p).mul(Fp::new(d, p).inv()?))
    }
}

fn bigint_mod(b: &BigInt, p: u64) -> u64 {
    let r = (b % BigInt::from(p)).to_string();
    // r is a decimal in (-p, p); parse with sign
    let x: i64 = r.parse().unwrap_or(0);
    x.rem_euclid(p as i64) as u64
}

// ── order-3 jets over 𝔽_p ─────────────────────────────────────────────

#[derive(Clone)]
struct JetF {
    v: Fp,
    d1: Vec<Fp>,
    d2: Vec<Fp>,
    d3: Vec<Fp>,
}

struct CtxF {
    n: usize,
    n2: usize,
    n3: usize,
    p: u64,
    t3: Vec<[usize; 3]>,
    t3idx: Vec<Vec<Vec<usize>>>,
}

impl CtxF {
    fn new(n: usize, p: u64) -> Self {
        let (t3, t3idx) = i3_table(n);
        CtxF { n, n2: n * (n + 1) / 2, n3: t3.len(), p, t3, t3idx }
    }
    fn zero(&self) -> JetF {
        JetF {
            v: Fp::zero(self.p),
            d1: vec![Fp::zero(self.p); self.n],
            d2: vec![Fp::zero(self.p); self.n2],
            d3: vec![Fp::zero(self.p); self.n3],
        }
    }
    fn constant(&self, c: Fp) -> JetF {
        let mut j = self.zero();
        j.v = c;
        j
    }
    fn coord(&self, i: usize, val: Fp) -> JetF {
        let mut j = self.zero();
        j.v = val;
        j.d1[i] = Fp::one(self.p);
        j
    }
    fn add(&self, a: &JetF, b: &JetF) -> JetF {
        let mut c = a.clone();
        c.v = c.v.add(b.v);
        for i in 0..self.n {
            c.d1[i] = c.d1[i].add(b.d1[i]);
        }
        for i in 0..self.n2 {
            c.d2[i] = c.d2[i].add(b.d2[i]);
        }
        for i in 0..self.n3 {
            c.d3[i] = c.d3[i].add(b.d3[i]);
        }
        c
    }
    fn sub(&self, a: &JetF, b: &JetF) -> JetF {
        let mut c = a.clone();
        c.v = c.v.sub(b.v);
        for i in 0..self.n {
            c.d1[i] = c.d1[i].sub(b.d1[i]);
        }
        for i in 0..self.n2 {
            c.d2[i] = c.d2[i].sub(b.d2[i]);
        }
        for i in 0..self.n3 {
            c.d3[i] = c.d3[i].sub(b.d3[i]);
        }
        c
    }
    fn neg(&self, a: &JetF) -> JetF {
        self.sub(&self.zero(), a)
    }
    fn mul(&self, a: &JetF, b: &JetF) -> JetF {
        let n = self.n;
        let mut c = self.zero();
        c.v = a.v.mul(b.v);
        for i in 0..n {
            c.d1[i] = a.d1[i].mul(b.v).add(a.v.mul(b.d1[i]));
        }
        for i in 0..n {
            for j in i..n {
                let ij = i2(n, i, j);
                c.d2[ij] = a.d2[ij]
                    .mul(b.v)
                    .add(a.d1[i].mul(b.d1[j]))
                    .add(a.d1[j].mul(b.d1[i]))
                    .add(a.v.mul(b.d2[ij]));
            }
        }
        for (id, t) in self.t3.iter().enumerate() {
            let [i, j, k] = *t;
            c.d3[id] = a.d3[id]
                .mul(b.v)
                .add(a.d2[i2(n, i, j)].mul(b.d1[k]))
                .add(a.d2[i2(n, i, k)].mul(b.d1[j]))
                .add(a.d2[i2(n, j, k)].mul(b.d1[i]))
                .add(a.d1[i].mul(b.d2[i2(n, j, k)]))
                .add(a.d1[j].mul(b.d2[i2(n, i, k)]))
                .add(a.d1[k].mul(b.d2[i2(n, i, j)]))
                .add(a.v.mul(b.d3[id]));
        }
        c
    }
    /// reciprocal 1/b; None at a modular pole (b.v ≡ 0).
    fn inv(&self, b: &JetF) -> Option<JetF> {
        let n = self.n;
        let bv_inv = b.v.inv()?;
        let mut r = self.zero();
        r.v = bv_inv;
        for i in 0..n {
            r.d1[i] = b.d1[i].mul(r.v).neg().mul(bv_inv);
        }
        for i in 0..n {
            for j in i..n {
                let ij = i2(n, i, j);
                let s = b.d2[ij]
                    .mul(r.v)
                    .add(b.d1[i].mul(r.d1[j]))
                    .add(b.d1[j].mul(r.d1[i]));
                r.d2[ij] = s.neg().mul(bv_inv);
            }
        }
        for (id, t) in self.t3.iter().enumerate() {
            let [i, j, k] = *t;
            let s = b.d3[id]
                .mul(r.v)
                .add(b.d2[i2(n, i, j)].mul(r.d1[k]))
                .add(b.d2[i2(n, i, k)].mul(r.d1[j]))
                .add(b.d2[i2(n, j, k)].mul(r.d1[i]))
                .add(b.d1[i].mul(r.d2[i2(n, j, k)]))
                .add(b.d1[j].mul(r.d2[i2(n, i, k)]))
                .add(b.d1[k].mul(r.d2[i2(n, i, j)]));
            r.d3[id] = s.neg().mul(bv_inv);
        }
        Some(r)
    }
    fn div(&self, a: &JetF, b: &JetF) -> Option<JetF> {
        Some(self.mul(a, &self.inv(b)?))
    }
    fn shift(&self, a: &JetF, dir: usize) -> JetF {
        let n = self.n;
        let mut s = self.zero();
        s.v = a.d1[dir];
        for i in 0..n {
            s.d1[i] = a.d2[i2(n, dir, i)];
        }
        for i in 0..n {
            for j in i..n {
                s.d2[i2(n, i, j)] = a.d3[self.t3idx[dir][i][j]];
            }
        }
        s
    }
}

// ── modular receipt ───────────────────────────────────────────────────

#[derive(serde::Serialize)]
pub struct ReceiptModp {
    pub schema: &'static str,
    pub m: usize,
    pub prime: u64,
    pub geometrization: String,
    pub point: Vec<String>,
    pub baseline_constant_curvature_quarter: bool,
    pub gprime_signature: (usize, usize),
    pub gprime_beta_norm: String,
    /// components of ∇R that are ≢ 0 mod p — each a PROOF of ≠ 0 over ℚ
    pub nabla_r_certified_nonzero: usize,
    pub nabla_r_total: usize,
    pub rr_nonzero_modp: usize,
    /// index tuple + nonzero residue certifying PROPER, if found
    pub proper_certificate: Option<(Vec<usize>, String)>,
    pub verdict: String,
    /// true ⇒ the verdict is a hard proof; false ⇒ needs exact ℚ recheck
    pub sound: bool,
    pub modular_pole: bool,
}

/// Exact ℚ value-only metric block at the point → g', signature input & β-norm.
/// Returns a named edge instead of a bare None so the caller can record WHY a
/// point has no valid drift-flip geometry (these edges are assumed to occur on
/// every large run, not treated as unreachable):
///   "outside_simplex" — Σp_i ≥ 1 (p_m ≤ 0), point not interior
///   "null_drift"      — g(β,β) = 0, the reflection g − 2β♭⊗β♭/g(β,β) is
///                       undefined (drift direction is g-null); itself a
///                       meaningful scar-map feature, not a failure.
fn exact_gprime(m: usize, point: &[Q]) -> Result<(Vec<Vec<Q>>, Q), &'static str> {
    use num_traits::One;
    let n = m - 1;
    if point.len() < n {
        return Err("outside_simplex");
    }
    let mut pm = Q::one();
    for pv in point.iter().take(n) {
        if pv <= &Q::zero() {
            return Err("outside_simplex");
        }
        pm -= pv;
    }
    if pm <= Q::zero() {
        return Err("outside_simplex");
    }
    let inv_pm = Q::one() / &pm;
    let mut g = vec![vec![Q::zero(); n]; n];
    for i in 0..n {
        for j in 0..n {
            g[i][j] = if i == j {
                Q::one() / &point[i] + &inv_pm
            } else {
                inv_pm.clone()
            };
        }
    }
    let beta_full = drift_beta(m);
    let beta: Vec<Q> = beta_full[..n].to_vec();
    let mut beta_flat = vec![Q::zero(); n];
    for i in 0..n {
        for j in 0..n {
            beta_flat[i] += &g[i][j] * &beta[j];
        }
    }
    let mut bnorm = Q::zero();
    for i in 0..n {
        bnorm += &beta[i] * &beta_flat[i];
    }
    if bnorm.is_zero() {
        return Err("null_drift"); // g-null drift: reflection undefined (edge)
    }
    let mut gp = vec![vec![Q::zero(); n]; n];
    for i in 0..n {
        for j in 0..n {
            gp[i][j] = &g[i][j] - &(Q::new(BigInt::from(2), BigInt::from(1)) * &beta_flat[i] * &beta_flat[j] / &bnorm);
        }
    }
    let mut bb = Q::zero();
    for i in 0..n {
        for j in 0..n {
            bb += &gp[i][j] * &beta[i] * &beta[j];
        }
    }
    Ok((gp, bb))
}

fn curvature_f(
    cx: &CtxF,
    g: &[Vec<JetF>],
    n: usize,
) -> Option<(Vec<Vec<Vec<Vec<JetF>>>>, Vec<Vec<Vec<Vec<JetF>>>>, Vec<Vec<Vec<JetF>>>)> {
    let ginv = mat_inv_f(cx, g, n)?;
    let half = cx.constant(Fp::new(1, cx.p).mul(Fp::new(2, cx.p).inv().unwrap()));
    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;
            }
        }
    }
    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);
                }
            }
        }
    }
    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);
                }
            }
        }
    }
    Some((rup, rdn, gamma))
}

fn mat_inv_f(cx: &CtxF, a: &[Vec<JetF>], n: usize) -> Option<Vec<Vec<JetF>>> {
    let mut m: Vec<Vec<JetF>> = a.to_vec();
    let mut inv: Vec<Vec<JetF>> = (0..n)
        .map(|i| {
            (0..n)
                .map(|j| if i == j { cx.constant(Fp::one(cx.p)) } else { cx.zero() })
                .collect()
        })
        .collect();
    for col in 0..n {
        let piv = (col..n).find(|&r| !m[r][col].v.is_zero())?; // None ⇒ modular pole
        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 {
            let entry_zero = 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());
            if r != col && !entry_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);
                }
            }
        }
    }
    Some(inv)
}

fn check_constant_curvature_f(cx: &CtxF, g: &[Vec<JetF>], n: usize) -> bool {
    let (_, rdn, _) = match curvature_f(cx, g, n) {
        Some(x) => x,
        None => return false,
    };
    let quarter = Fp::new(1, cx.p).mul(Fp::new(4, cx.p).inv().unwrap());
    for i in 0..n {
        for j in 0..n {
            for k in 0..n {
                for l in 0..n {
                    let expect = quarter.mul(g[i][k].v.mul(g[j][l].v).sub(g[i][l].v.mul(g[j][k].v)));
                    if rdn[i][j][k][l].v != expect {
                        return false;
                    }
                }
            }
        }
    }
    true
}

/// Single-prime modular probe.
///   Err(reason)          — a prime-independent geometric edge (outside_simplex
///                          / null_drift / degenerate_signature_block); the
///                          caller records it once, does not retry other primes.
///   Ok(r), r.modular_pole — this prime divides a nonzero-over-ℚ quantity; the
///                          caller retries another prime.
///   Ok(r), r.sound        — PROPER is a hard proof (see module header).
///   Ok(r), !r.sound       — a candidate all-zero verdict; caller escalates ℚ.
pub fn probe_modp(m: usize, point: &[Q], p: u64) -> Result<ReceiptModp, &'static str> {
    let n = m - 1;

    // exact ℚ side: signature + β-norm (needs sign/order). Named edges bubble up.
    let (gm_q, bb) = exact_gprime(m, point)?;
    let sig = signature_opt(&gm_q).ok_or("degenerate_signature_block")?;

    let cx = CtxF::new(n, p);
    let pole_receipt = |base_ok: bool| ReceiptModp {
        schema: "semisym_receipt_modp_v1",
        m,
        prime: p,
        geometrization: "drift_flip".into(),
        point: point.iter().map(|x| x.to_string()).collect(),
        baseline_constant_curvature_quarter: base_ok,
        gprime_signature: sig,
        gprime_beta_norm: bb.to_string(),
        nabla_r_certified_nonzero: 0,
        nabla_r_total: n.pow(5),
        rr_nonzero_modp: 0,
        proper_certificate: None,
        verdict: "MODULAR_POLE".into(),
        sound: false,
        modular_pole: true,
    };

    // build 𝔽_p coord jets (p | a coord denominator ⇒ modular pole, retry prime)
    let ps: Vec<JetF> = match (0..n)
        .map(|i| Fp::from_q(&point[i], p).map(|c| cx.coord(i, c)))
        .collect::<Option<Vec<_>>>()
    {
        Some(v) => v,
        None => return Ok(pole_receipt(false)),
    };
    let mut pm = cx.constant(Fp::one(p));
    for pj in &ps {
        pm = cx.sub(&pm, pj);
    }
    let inv_pm = match cx.inv(&pm) {
        Some(x) => x,
        None => return Ok(pole_receipt(false)),
    };
    let mut g = vec![vec![cx.zero(); n]; n];
    for i in 0..n {
        let inv_pi = match cx.inv(&ps[i]) {
            Some(x) => x,
            None => return Ok(pole_receipt(false)),
        };
        for j in 0..n {
            g[i][j] = if i == j { cx.add(&inv_pi, &inv_pm) } else { inv_pm.clone() };
        }
    }
    let base_ok = check_constant_curvature_f(&cx, &g, n);

    // drift-flip g'
    let beta_full = drift_beta(m);
    let bjets: Vec<JetF> = match (0..n)
        .map(|i| Fp::from_q(&beta_full[i], p).map(|b| cx.constant(b)))
        .collect::<Option<Vec<_>>>()
    {
        Some(v) => v,
        None => return Ok(pole_receipt(base_ok)),
    };
    let mut beta_flat = 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(Fp::new(2, p));
    let mut gp = vec![vec![cx.zero(); n]; n];
    for i in 0..n {
        for j in 0..n {
            let corr = match cx.div(&cx.mul(&two, &cx.mul(&beta_flat[i], &beta_flat[j])), &bnorm) {
                Some(x) => x,
                None => return Ok(pole_receipt(base_ok)),
            };
            gp[i][j] = cx.sub(&g[i][j], &corr);
        }
    }

    let (rup, rdn, gamma) = match curvature_f(&cx, &gp, n) {
        Some(x) => x,
        None => return Ok(pole_receipt(base_ok)),
    };

    // ∇R — count certified-nonzero components
    let mut nabla_nz = 0usize;
    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 {
                        let mut v = rdn[i][jj][k][l].d1[a];
                        for mm in 0..n {
                            v = v.sub(gamma[mm][a][i].v.mul(rdn[mm][jj][k][l].v));
                            v = v.sub(gamma[mm][a][jj].v.mul(rdn[i][mm][k][l].v));
                            v = v.sub(gamma[mm][a][k].v.mul(rdn[i][jj][mm][l].v));
                            v = v.sub(gamma[mm][a][l].v.mul(rdn[i][jj][k][mm].v));
                        }
                        if !v.is_zero() {
                            nabla_nz += 1;
                        }
                    }
                }
            }
        }
    }

    // R·R and Q(g,R): collect pairs, certify PROPER by cross-multiplication
    let mut rr_nz = 0usize;
    let mut pairs: Vec<(Fp, Fp, Vec<usize>)> = Vec::new();
    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 = Fp::zero(p);
                            let mut qq = Fp::zero(p);
                            for mm in 0..n {
                                rr = rr.sub(rup[mm][i][a][b].v.mul(rdn[mm][jj][k][l].v));
                                rr = rr.sub(rup[mm][jj][a][b].v.mul(rdn[i][mm][k][l].v));
                                rr = rr.sub(rup[mm][k][a][b].v.mul(rdn[i][jj][mm][l].v));
                                rr = rr.sub(rup[mm][l][a][b].v.mul(rdn[i][jj][k][mm].v));
                            }
                            let gv = |r: usize, c: usize| -> Fp { gp[r][c].v };
                            let wedge = |mm: usize, idx: usize| -> Fp {
                                let mut w = Fp::zero(p);
                                if mm == a {
                                    w = w.add(gv(b, idx));
                                }
                                if mm == b {
                                    w = w.sub(gv(a, idx));
                                }
                                w
                            };
                            for mm in 0..n {
                                qq = qq.sub(wedge(mm, i).mul(rdn[mm][jj][k][l].v));
                                qq = qq.sub(wedge(mm, jj).mul(rdn[i][mm][k][l].v));
                                qq = qq.sub(wedge(mm, k).mul(rdn[i][jj][mm][l].v));
                                qq = qq.sub(wedge(mm, l).mul(rdn[i][jj][k][mm].v));
                            }
                            if !rr.is_zero() {
                                rr_nz += 1;
                            }
                            if !(rr.is_zero() && qq.is_zero()) {
                                pairs.push((rr, qq, vec![a, b, i, jj, k, l]));
                            }
                        }
                    }
                }
            }
        }
    }

    // PROPER certificate: reference pair with qq ≠ 0, then any pair breaking
    // rr·qq₀ − rr₀·qq ≡ 0.  If all qq ≡ 0 but some rr ≠ 0 ⇒ also PROPER.
    let mut proper_cert: Option<(Vec<usize>, String)> = None;
    if let Some((rr0, qq0, _)) = pairs.iter().find(|(_, q, _)| !q.is_zero()).cloned() {
        for (rr, qq, idx) in &pairs {
            let cross = rr.mul(qq0).sub(rr0.mul(*qq));
            if !cross.is_zero() {
                proper_cert = Some((idx.clone(), cross.v.to_string()));
                break;
            }
        }
    } else if let Some((_, _, idx)) = pairs.iter().find(|(r, _, _)| !r.is_zero()) {
        // no Q-side entry nonzero but R·R nonzero ⇒ no finite L works
        proper_cert = Some((idx.clone(), "rr_nonzero_q_all_zero".into()));
    }

    let (verdict, sound) = if nabla_nz == 0 {
        ("CANDIDATE_LOCALLY_SYMMETRIC(needs exact)".to_string(), false)
    } else if rr_nz == 0 {
        ("CANDIDATE_SEMISYMMETRIC(needs exact)".to_string(), false)
    } else if proper_cert.is_some() {
        ("PROPER".to_string(), true) // hard proof: ∇R≠0 (certified) & not pseudosym
    } else {
        ("CANDIDATE_PSEUDOSYMMETRIC(needs exact)".to_string(), false)
    };

    Ok(ReceiptModp {
        schema: "semisym_receipt_modp_v1",
        m,
        prime: p,
        geometrization: "drift_flip".into(),
        point: point.iter().map(|x| x.to_string()).collect(),
        baseline_constant_curvature_quarter: base_ok,
        gprime_signature: sig,
        gprime_beta_norm: bb.to_string(),
        nabla_r_certified_nonzero: nabla_nz,
        nabla_r_total: n.pow(5),
        rr_nonzero_modp: rr_nz,
        proper_certificate: proper_cert,
        verdict,
        sound,
        modular_pole: false,
    })
}

// ── robust orchestration: assume edges happen on every run ────────────
//
// GUARANTEE: probe_survey resolves EVERY point to exactly one of
//   (a) a proven-PROPER modular receipt   (fast, hard proof over ℚ),
//   (b) a sound exact-ℚ receipt           (auto-escalation), or
//   (c) a recorded edge receipt           (geometric edge, no crash).
// It never panics on a known edge and never emits an unsound verdict.

/// A geometric / arithmetic edge recorded in-band (never a panic).
#[derive(serde::Serialize)]
pub struct EdgeReceipt {
    pub schema: &'static str,
    pub m: usize,
    pub point: Vec<String>,
    pub verdict: String,
    pub reason: String,
    pub note: String,
}

/// Default multi-prime set: the three largest primes below 2³¹. Multi-prime
/// makes a modular pole (one prime dividing a nonzero-over-ℚ quantity) a
/// non-event — another prime almost surely resolves it — and guards a
/// single-prime false-zero from hiding a PROPER certificate.
pub const DEFAULT_PRIMES: [u64; 3] = [2_147_483_647, 2_147_483_629, 2_147_483_587];

fn mod_mul(a: u64, b: u64, m: u64) -> u64 {
    ((a as u128 * b as u128) % m as u128) as u64
}
fn mod_pow(mut a: u64, mut e: u64, m: u64) -> u64 {
    let mut r = 1u64 % m;
    while e > 0 {
        if e & 1 == 1 {
            r = mod_mul(r, a, m);
        }
        a = mod_mul(a, a, m);
        e >>= 1;
    }
    r
}

/// Deterministic Miller–Rabin; bases {2,3,5,7} are a proven witness set for all
/// n < 3.2·10⁹ ⊃ our p < 2³¹. Guards against a composite modulus, which would
/// make the Fermat inverse (and every "proof") silently unsound.
pub fn is_prime(n: u64) -> bool {
    if n < 2 {
        return false;
    }
    for &sp in &[2u64, 3, 5, 7, 11, 13, 17, 19, 23] {
        if n % sp == 0 {
            return n == sp;
        }
    }
    let mut d = n - 1;
    let mut r = 0u32;
    while d % 2 == 0 {
        d /= 2;
        r += 1;
    }
    'witness: for &a in &[2u64, 3, 5, 7] {
        let mut x = mod_pow(a % n, d, n);
        if x == 1 || x == n - 1 {
            continue;
        }
        for _ in 0..r.saturating_sub(1) {
            x = mod_mul(x, x, n);
            if x == n - 1 {
                continue 'witness;
            }
        }
        return false;
    }
    true
}

fn edge(m: usize, point: &[Q], reason: &str, note: &str) -> (String, String) {
    let e = EdgeReceipt {
        schema: "semisym_receipt_edge_v1",
        m,
        point: point.iter().map(|x| x.to_string()).collect(),
        verdict: format!("EDGE_{}", reason.to_uppercase()),
        reason: reason.into(),
        note: note.into(),
    };
    let j = serde_json::to_string_pretty(&e).unwrap();
    (e.verdict, j)
}

/// Robust per-point resolution used by surveys. Returns (verdict_for_tally,
/// pretty_json). See the module GUARANTEE above.
pub fn probe_survey(m: usize, point: &[Q], primes: &[u64]) -> (String, String) {
    let mut saw_candidate = false;
    let mut n_poles = 0usize;
    for &p in primes {
        match probe_modp(m, point, p) {
            // prime-independent geometric edge — record once, no retry
            Err(reason) => {
                return edge(m, point, reason, "no valid drift-flip geometry at this point");
            }
            Ok(r) if r.modular_pole => {
                n_poles += 1;
                continue; // this prime divides something; try the next
            }
            Ok(r) if r.sound => {
                // PROPER proven over ℚ by this prime — fast path, done
                let j = serde_json::to_string_pretty(&r).unwrap();
                return (r.verdict, j);
            }
            Ok(_) => {
                // candidate all-zero on this prime; another prime may still
                // expose a PROPER certificate hidden by a false-zero
                saw_candidate = true;
            }
        }
    }
    // No prime proved PROPER ⇒ escalate to exact ℚ for a sound verdict. Covers
    // genuine candidate verdicts, single-prime false-zeros, and all-poled.
    // (Reached only for interior, non-null-drift points, so crate::probe's
    //  interior/non-degenerate preconditions hold; the survey also wraps this
    //  in catch_unwind as a final backstop.)
    let ex = crate::probe(m, point);
    let why = if saw_candidate {
        "no_modular_PROPER_proof(candidate all-zero)"
    } else {
        "all_primes_poled"
    };
    let mut v: serde_json::Value = serde_json::to_value(&ex).unwrap();
    v["_escalated_from_modular"] = serde_json::json!({
        "reason": why,
        "primes_tried": primes,
        "poles": n_poles,
    });
    (ex.verdict.clone(), serde_json::to_string_pretty(&v).unwrap())
}
