wai-quantum 0.3.4

//! Deterministic quantum LDPC DECODING — `wai.quantum.decode`
//! (extensions/quantum-ops § Decoding).
//!
//! A fault-tolerant quantum computer runs a classical decoder in its innermost
//! loop: every syndrome-measurement round produces a syndrome, and a decoder must
//! turn it into a correction fast enough and cheaply enough to keep up with the
//! code cycle. For quantum LDPC (qLDPC) codes — the high-rate codes that make
//! fault tolerance affordable — the decoder of record is **belief propagation**,
//! and plain BP *fails* on quantum codes: the stabilizer degeneracy and the short
//! cycles in the Tanner graph trap it in symmetric non-solutions. The SOTA fix is
//! **Relay-BP**: an ensemble of BP passes with disordered per-node memory
//! strengths, relayed (warm-started) from one another, which breaks the symmetry.
//!
//! This module is Relay-BP done the WAI way: **normalized min-sum** message
//! passing in pure fixed-point integer arithmetic — no float, no `tanh`/`atanh` —
//! so a decode is **byte-identical on every machine** and the emitted correction
//! is a portable hash. Every decode seals a signed, joule-metered
//! [`crate::quantum_ops::DecodeReceipt`] carrying **joules-per-decode**, the metric
//! no decoder reports and none signs.
//!
//! Min-sum on a single uniform channel is scale-invariant (check node = min,
//! variable node = sum, both homogeneous), so the channel prior is a fixed unit
//! and the exact physical error rate never enters the decoder — only the
//! deterministic error *sampler* uses it. That keeps the decode fully integer.
//!
//! # Decoder scope, measured
//!
//! Relay-BP is a belief-propagation decoder, and it earns its name on the qLDPC
//! family it was designed for: on the bivariate-bicycle code the relay legs beat
//! plain BP (at `p = 0.06`, 0.30 → 0.27 logical failure), and more legs help.
//!
//! On **topological** codes it does not show a distance benefit — measured from
//! `p = 0.001` to `0.02`, distance-5 loses to distance-3 on the surface code *and*
//! on the toric code alike. This is the known weakness of BP-family decoding on
//! degenerate codes with short cycles, not a defect in the code constructions
//! (both are verified valid CSS codes with the expected `k`). Decoding a surface
//! code competitively wants a matching decoder — MWPM or union-find — which this
//! module does not yet provide. The constructions are here, correct and
//! exportable; the matching decoder is the honest next piece of work.
//!
//! Codes built in: the **rotated surface code** (`[[d², 1, d]]`, the planar code
//! hardware actually targets), the **toric code** (its periodic-boundary cousin)
//! and a
//! **bivariate-bicycle** qLDPC code (the high-rate family; the \[\[72,12\]\] preset is
//! the "gross"-style construction). Logical success — knowable only in simulation,
//! since it needs the true error — is reported by the benchmark harness, never put
//! in the receipt (a real decoder does not know it).

use crate::quantum_ops::{content_hash, DecodeReceipt, GrantRef};
use ed25519_dalek::SigningKey;

/// Fixed-point fractional bits for messages/damping. `1.0` = `1 << QEC_FRAC`.
pub const QEC_FRAC: u32 = 12;
const SCALE: i64 = 1 << QEC_FRAC;
/// Channel prior (uniform): min-sum is scale-invariant, so any positive unit works.
const L0: i64 = SCALE;
/// Message saturation, keeps fixed-point sums bounded without changing decisions.
const MSG_MAX: i64 = 64 * SCALE;
/// Normalized-min-sum scaling α = 7/8 (the standard hardware value).
const ALPHA_NUM: i64 = 7;
/// Per-node relay memory strength range. Relay-BP draws *disordered* memory,
/// including negative values (over-relaxation) that aggressively break the
/// symmetric trapping sets plain BP gets stuck in. γ ∈ [-0.6, 0.9].
const GAMMA_MIN: i64 = -(6 * SCALE) / 10;
const GAMMA_MAX: i64 = (9 * SCALE) / 10;

fn splitmix64(state: &mut u64) -> u64 {
    *state = state.wrapping_add(0x9E37_79B9_7F4A_7C15);
    let mut z = *state;
    z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
    z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
    z ^ (z >> 31)
}

// ===========================================================================
// GF(2) linear algebra — row space membership (logical-triviality test)
// ===========================================================================

/// A reduced GF(2) basis of a set of binary row vectors, supporting membership
/// testing. Rows are packed into `u64` words. Used to test whether a residual
/// error is a product of stabilizer generators (i.e. logically trivial).
#[derive(Clone, Debug)]
pub struct Gf2Basis {
    words: usize,
    /// Reduced rows, each with a distinct pivot column.
    rows: Vec<Vec<u64>>,
    /// Pivot column of each reduced row (parallel to `rows`).
    pivots: Vec<usize>,
}

fn bit_get(row: &[u64], i: usize) -> bool {
    (row[i >> 6] >> (i & 63)) & 1 == 1
}
fn bit_set(row: &mut [u64], i: usize) {
    row[i >> 6] |= 1u64 << (i & 63);
}
fn bit_xor(dst: &mut [u64], src: &[u64]) {
    for (d, s) in dst.iter_mut().zip(src) {
        *d ^= *s;
    }
}
fn is_zero(row: &[u64]) -> bool {
    row.iter().all(|&w| w == 0)
}

impl Gf2Basis {
    /// Build from a list of rows, each given as the sorted set of its 1-columns.
    pub fn from_rows(rows_cols: &[Vec<usize>], ncols: usize) -> Gf2Basis {
        let words = ncols.div_ceil(64).max(1);
        let mut basis: Vec<Vec<u64>> = Vec::new();
        let mut pivots: Vec<usize> = Vec::new();
        for cols in rows_cols {
            let mut v = vec![0u64; words];
            for &c in cols {
                bit_set(&mut v, c);
            }
            // reduce against current basis
            for (bi, br) in basis.iter().enumerate() {
                if bit_get(&v, pivots[bi]) {
                    bit_xor(&mut v, br);
                }
            }
            if let Some(p) = Self::lowest_set(&v) {
                // eliminate this pivot from existing basis rows for a full RREF
                for br in basis.iter_mut() {
                    if bit_get(br, p) {
                        bit_xor(br, &v);
                    }
                }
                basis.push(v);
                pivots.push(p);
            }
        }
        Gf2Basis { words, rows: basis, pivots }
    }

    fn lowest_set(v: &[u64]) -> Option<usize> {
        for (w, &word) in v.iter().enumerate() {
            if word != 0 {
                return Some(w * 64 + word.trailing_zeros() as usize);
            }
        }
        None
    }

    /// Rank of the row space.
    pub fn rank(&self) -> usize {
        self.rows.len()
    }

    /// Is the vector (given as its 1-columns) in the row space? Reduce it against
    /// the basis and check it vanishes.
    pub fn contains(&self, cols: &[usize]) -> bool {
        let mut v = vec![0u64; self.words];
        for &c in cols {
            bit_set(&mut v, c);
        }
        for (bi, br) in self.rows.iter().enumerate() {
            if bit_get(&v, self.pivots[bi]) {
                bit_xor(&mut v, br);
            }
        }
        is_zero(&v)
    }
}

// ===========================================================================
// CSS code — check matrix (for decoding) + stabilizer basis (for success test)
// ===========================================================================

/// A CSS code prepared for decoding one error sector. `checks` are the parity
/// checks whose syndrome we decode (rows as sorted qubit lists); `stab_basis` is
/// the GF(2) row space of the complementary stabilizers, against which a residual
/// error is tested for logical triviality.
#[derive(Clone, Debug)]
pub struct CssCode {
    pub id: String,
    pub n: usize,
    /// Parity checks (syndrome producers), each a sorted list of qubit indices.
    pub checks: Vec<Vec<usize>>,
    /// The complementary stabilizer family, kept so callers can inspect the code
    /// (a CSS code is defined by both families, not just the one that produces the
    /// syndrome being decoded).
    pub stabs: Vec<Vec<usize>>,
    stab_basis: Gf2Basis,
    /// Number of encoded logical qubits, `k = n − rank(checks) − rank(stab)`.
    pub k: usize,
}

impl CssCode {
    fn new(id: String, n: usize, checks: Vec<Vec<usize>>, stab: Vec<Vec<usize>>) -> CssCode {
        let stab_basis = Gf2Basis::from_rows(&stab, n);
        let check_basis = Gf2Basis::from_rows(&checks, n);
        let k = n.saturating_sub(check_basis.rank() + stab_basis.rank());
        CssCode { id, n, checks, stabs: stab, stab_basis, k }
    }

    /// Canonical bytes of the check matrix (content-addresses the code).
    pub fn check_matrix_bytes(&self) -> Vec<u8> {
        let mut b = Vec::new();
        b.extend_from_slice(b"wai:qec-checks\x01");
        b.extend_from_slice(&(self.n as u64).to_le_bytes());
        for row in &self.checks {
            b.extend_from_slice(&(row.len() as u32).to_le_bytes());
            for &c in row {
                b.extend_from_slice(&(c as u32).to_le_bytes());
            }
        }
        b
    }

    /// A residual error (as its 1-qubits) is logically trivial iff it lies in the
    /// stabilizer row space.
    pub fn is_trivial(&self, residual_cols: &[usize]) -> bool {
        self.stab_basis.contains(residual_cols)
    }

    // --- Toric code -------------------------------------------------------

    /// The `L×L` toric code: `n = 2L²` qubits on edges, `L²` star (X) checks and
    /// `L²` plaquette (Z) checks, `k = 2`. We decode X-errors: `checks = Z`
    /// (plaquettes), triviality against the X-stabilizers (stars).
    pub fn toric(l: usize) -> CssCode {
        let n = 2 * l * l;
        let h = |i: usize, j: usize| (i % l) * l + (j % l); // horizontal edge id
        let v = |i: usize, j: usize| l * l + (i % l) * l + (j % l); // vertical edge id
        let mut plaq = Vec::with_capacity(l * l); // Z checks
        let mut star = Vec::with_capacity(l * l); // X stabilizers
        for i in 0..l {
            for j in 0..l {
                // plaquette (i,j): top/bottom horizontals, left/right verticals
                let mut p = vec![h(i, j), h(i + 1, j), v(i, j), v(i, j + 1)];
                p.sort_unstable();
                p.dedup();
                plaq.push(p);
                // star at vertex (i,j): incident horizontals + verticals
                let mut s = vec![h(i, j), h(i, j + l - 1), v(i, j), v(i + l - 1, j)];
                s.sort_unstable();
                s.dedup();
                star.push(s);
            }
        }
        CssCode::new(format!("toric:L{l}"), n, plaq, star)
    }

    // --- Rotated surface code ---------------------------------------------

    /// The distance-`d` **rotated surface code**, `[[d², 1, d]]` (`d` odd).
    ///
    /// The toric code above is this code's periodic-boundary cousin: elegant, but
    /// it needs a torus. The rotated surface code is the planar one — nearest-
    /// neighbour on a `d × d` patch with open boundaries — and it is what
    /// essentially every hardware roadmap actually targets.
    ///
    /// Data qubit `(r, c)` is index `r·d + c`. Bulk plaquettes cover each 2×2
    /// square and alternate in a checkerboard; the boundaries carry weight-2
    /// stabilizers, X along the top and bottom, Z along the left and right. The
    /// parities are the ones that make the two families commute — the neighbouring
    /// choice does not, which the `css_commutation` test would catch.
    pub fn surface(d: usize) -> CssCode {
        assert!(d >= 3 && d % 2 == 1, "rotated surface code needs an odd distance >= 3");
        let n = d * d;
        let q = |r: usize, c: usize| r * d + c;
        let mut z_checks: Vec<Vec<usize>> = Vec::new();
        let mut x_stabs: Vec<Vec<usize>> = Vec::new();

        // Bulk: every 2x2 square, checkerboarded.
        for r in 0..d - 1 {
            for c in 0..d - 1 {
                let face = vec![q(r, c), q(r, c + 1), q(r + 1, c), q(r + 1, c + 1)];
                if (r + c) % 2 == 0 {
                    z_checks.push(face);
                } else {
                    x_stabs.push(face);
                }
            }
        }
        // X boundary: top row on even columns, bottom row on odd columns.
        for c in (0..d - 1).step_by(2) {
            x_stabs.push(vec![q(0, c), q(0, c + 1)]);
        }
        for c in (1..d - 1).step_by(2) {
            x_stabs.push(vec![q(d - 1, c), q(d - 1, c + 1)]);
        }
        // Z boundary: left column on odd rows, right column on even rows.
        for r in (1..d - 1).step_by(2) {
            z_checks.push(vec![q(r, 0), q(r + 1, 0)]);
        }
        for r in (0..d - 1).step_by(2) {
            z_checks.push(vec![q(r, d - 1), q(r + 1, d - 1)]);
        }
        for v in z_checks.iter_mut().chain(x_stabs.iter_mut()) {
            v.sort_unstable();
            v.dedup();
        }
        CssCode::new(format!("surface:d{d}"), n, z_checks, x_stabs)
    }

    // --- Bivariate-bicycle qLDPC code -------------------------------------

    /// A bivariate-bicycle code on `Z_l × Z_m`. `a` and `b` are lists of monomials
    /// `(px, py)`; `A = Σ x^px y^py`, `B` likewise (mod-2 sums of shift matrices).
    /// `H_X = [A | B]`, `H_Z = [Bᵀ | Aᵀ]`, `n = 2·l·m`. We decode X-errors:
    /// `checks = H_Z`, triviality against `H_X`.
    pub fn bivariate_bicycle(l: usize, m: usize, a: &[(usize, usize)], b: &[(usize, usize)]) -> CssCode {
        let lm = l * m;
        let n = 2 * lm;
        let cell = |i: usize, j: usize| (i % l) * m + (j % m);
        // A[row][col]: monomial (px,py) sends col-cell (i,j) → row-cell (i+px, j+py).
        // Build row-wise supports for A and B (over the lm cells).
        let build = |monos: &[(usize, usize)]| -> Vec<Vec<usize>> {
            let mut rows = vec![Vec::new(); lm];
            for i in 0..l {
                for j in 0..m {
                    let col = cell(i, j);
                    for &(px, py) in monos {
                        let row = cell(i + px, j + py);
                        rows[row].push(col);
                    }
                }
            }
            for r in rows.iter_mut() {
                r.sort_unstable();
                r.dedup();
            }
            rows
        };
        let a_rows = build(a);
        let b_rows = build(b);
        // transpose helper (for Bᵀ, Aᵀ)
        let transpose = |rows: &[Vec<usize>]| -> Vec<Vec<usize>> {
            let mut t = vec![Vec::new(); lm];
            for (r, cols) in rows.iter().enumerate() {
                for &c in cols {
                    t[c].push(r);
                }
            }
            for x in t.iter_mut() {
                x.sort_unstable();
            }
            t
        };
        let at = transpose(&a_rows);
        let bt = transpose(&b_rows);

        // H_X = [A | B]  (columns 0..lm = block L, lm..2lm = block R)
        let hx: Vec<Vec<usize>> = (0..lm)
            .map(|r| {
                let mut row: Vec<usize> = a_rows[r].clone();
                row.extend(b_rows[r].iter().map(|&c| c + lm));
                row.sort_unstable();
                row
            })
            .collect();
        // H_Z = [Bᵀ | Aᵀ]
        let hz: Vec<Vec<usize>> = (0..lm)
            .map(|r| {
                let mut row: Vec<usize> = bt[r].clone();
                row.extend(at[r].iter().map(|&c| c + lm));
                row.sort_unstable();
                row
            })
            .collect();

        // decode X-errors: checks = H_Z, triviality against H_X
        let code = CssCode::new(String::new(), n, hz, hx);
        CssCode { id: format!("bb:[[{},{}]]", code.n, code.k), ..code }
    }

    /// The `[[72,12]]` bivariate-bicycle preset (l=m=6, A=x³+y+y², B=y³+x+x²).
    pub fn gross() -> CssCode {
        CssCode::bivariate_bicycle(6, 6, &[(3, 0), (0, 1), (0, 2)], &[(0, 3), (1, 0), (2, 0)])
    }
}

// ===========================================================================
// Tanner graph + normalized min-sum belief propagation
// ===========================================================================

struct Tanner {
    n: usize,
    n_edges: usize,
    e_var: Vec<u32>,          // variable (qubit) of each edge
    chk_range: Vec<(u32, u32)>, // contiguous edge range per check
    var_edges: Vec<Vec<u32>>, // edges incident to each variable
}

impl Tanner {
    fn build(checks: &[Vec<usize>], n: usize) -> Tanner {
        let mut e_var = Vec::new();
        let mut chk_range = Vec::with_capacity(checks.len());
        let mut var_edges = vec![Vec::new(); n];
        for cols in checks {
            let start = e_var.len() as u32;
            for &v in cols {
                var_edges[v].push(e_var.len() as u32);
                e_var.push(v as u32);
            }
            chk_range.push((start, e_var.len() as u32));
        }
        let n_edges = e_var.len();
        Tanner { n, n_edges, e_var, chk_range, var_edges }
    }
}

/// One min-sum BP leg over the current messages. `gamma` is the per-variable relay
/// memory strength (fixed-point, 0 = no memory). Returns `(hard_decision,
/// converged, iters, messages)`. Messages are warm-started from `mu_v2c`.
#[allow(clippy::too_many_arguments)]
fn bp_leg(
    t: &Tanner,
    syndrome: &[bool],
    gamma: &[i64],
    max_iter: u32,
    mu_v2c: &mut [i64],
    mu_c2v: &mut [i64],
) -> (Vec<bool>, bool, u32, u64) {
    let mut hard = vec![false; t.n];
    let mut messages = 0u64;
    let mut iters = 0;
    for _ in 0..max_iter {
        iters += 1;
        // --- check update (min-sum) ---
        for (&(a, b), &s) in t.chk_range.iter().zip(syndrome) {
            let (mut min1, mut min2) = (i64::MAX, i64::MAX);
            let mut arg = a;
            let mut neg_parity = false;
            for e in a..b {
                let x = mu_v2c[e as usize];
                if x < 0 {
                    neg_parity = !neg_parity;
                }
                let mag = x.abs();
                if mag < min1 {
                    min2 = min1;
                    min1 = mag;
                    arg = e;
                } else if mag < min2 {
                    min2 = mag;
                }
            }
            for e in a..b {
                let x = mu_v2c[e as usize];
                let base = if e == arg { min2 } else { min1 };
                let mag = (base * ALPHA_NUM) >> 3;
                let is_neg = neg_parity ^ (x < 0) ^ s;
                mu_c2v[e as usize] = if is_neg { -mag } else { mag };
            }
        }
        // --- variable update + hard decision ---
        for v in 0..t.n {
            let mut total = L0;
            for &e in &t.var_edges[v] {
                total += mu_c2v[e as usize];
            }
            hard[v] = total < 0;
            let g = gamma[v];
            for &e in &t.var_edges[v] {
                let target = total - mu_c2v[e as usize];
                let nv = if g == 0 {
                    target
                } else {
                    (g * mu_v2c[e as usize] + (SCALE - g) * target) >> QEC_FRAC
                };
                mu_v2c[e as usize] = nv.clamp(-MSG_MAX, MSG_MAX);
            }
        }
        messages += 2 * t.n_edges as u64;
        // --- convergence: H·hard == syndrome ---
        let mut ok = true;
        for (&(a, b), &s) in t.chk_range.iter().zip(syndrome) {
            let mut par = false;
            for e in a..b {
                par ^= hard[t.e_var[e as usize] as usize];
            }
            if par != s {
                ok = false;
                break;
            }
        }
        if ok {
            return (hard, true, iters, messages);
        }
    }
    (hard, false, iters, messages)
}

/// Relay-BP: leg 0 is plain min-sum; each subsequent leg draws disordered
/// per-variable memory strengths (deterministic from `seed`) and warm-starts
/// (relays) from the previous leg. Stops at the first leg that satisfies the
/// syndrome. `max_legs = 1` recovers plain BP.
struct DecodeCore {
    correction: Vec<bool>,
    converged: bool,
    legs: u32,
    iters: u32,
    messages: u64,
}

fn relay_bp(
    t: &Tanner,
    syndrome: &[bool],
    seed: u64,
    max_iter: u32,
    max_legs: u32,
) -> DecodeCore {
    let mut mu_v2c = vec![L0; t.n_edges];
    let mut mu_c2v = vec![0i64; t.n_edges];
    let gamma0 = vec![0i64; t.n];
    let (mut hard, mut conv, mut it, mut msg) =
        bp_leg(t, syndrome, &gamma0, max_iter, &mut mu_v2c, &mut mu_c2v);
    let mut legs = 1u32;
    let (mut tot_it, mut tot_msg) = (it, msg);
    while !conv && legs < max_legs {
        let mut gamma = vec![0i64; t.n];
        let mut st = seed
            .wrapping_mul(0x100_0001)
            .wrapping_add(legs as u64)
            .wrapping_add(0xD15EA5E);
        let span = (GAMMA_MAX - GAMMA_MIN) as u64 + 1;
        for gv in gamma.iter_mut() {
            let r = splitmix64(&mut st);
            *gv = GAMMA_MIN + (r % span) as i64;
        }
        let out = bp_leg(t, syndrome, &gamma, max_iter, &mut mu_v2c, &mut mu_c2v);
        hard = out.0;
        conv = out.1;
        it = out.2;
        msg = out.3;
        tot_it += it;
        tot_msg += msg;
        legs += 1;
    }
    let _ = it;
    let _ = msg;
    DecodeCore { correction: hard, converged: conv, legs, iters: tot_it, messages: tot_msg }
}

// ===========================================================================
// Simulation front-end + decode result
// ===========================================================================

/// Deterministic decoder configuration.
#[derive(Clone, Copy, Debug)]
pub struct DecodeConfig {
    /// Physical error rate per qubit (fixed-point, scale 2^12).
    pub p_fx: i64,
    /// Max BP iterations per leg.
    pub max_iter: u32,
    /// Max relay legs (1 = plain BP).
    pub max_legs: u32,
}

impl DecodeConfig {
    pub fn from_p(p: f64) -> DecodeConfig {
        DecodeConfig { p_fx: (p * SCALE as f64) as i64, max_iter: 32, max_legs: 12 }
    }
}

/// One decode: sample an error at rate `p`, decode its syndrome, and check
/// logical success against the true error. Deterministic given `(code, cfg, seed)`.
#[derive(Clone, Debug)]
pub struct DecodeResult {
    pub code_id: String,
    pub n: usize,
    pub error: Vec<u32>,      // true error qubits (simulation only)
    pub correction: Vec<u32>, // decoder output qubits
    pub syndrome: Vec<u32>,   // unsatisfied check indices
    pub converged: bool,
    pub logical_success: bool,
    pub legs: u32,
    pub iters: u32,
    pub messages: u64,
    syndrome_bits: Vec<bool>,
    correction_bits: Vec<bool>,
    check_matrix_hash: [u8; 32],
}

/// Sample a bit-flip error deterministically at rate `p_fx`.
fn sample_error(n: usize, p_fx: i64, seed: u64) -> Vec<bool> {
    let mut e = vec![false; n];
    let mut st = seed.wrapping_mul(0x2545_F491_4F6C_DD1D).wrapping_add(1);
    for ei in e.iter_mut() {
        let r = (splitmix64(&mut st) >> (64 - QEC_FRAC)) as i64; // [0, SCALE)
        if r < p_fx {
            *ei = true;
        }
    }
    e
}

fn syndrome_of(checks: &[Vec<usize>], e: &[bool]) -> Vec<bool> {
    checks
        .iter()
        .map(|row| row.iter().fold(false, |acc, &q| acc ^ e[q]))
        .collect()
}

fn bits_to_ids(bits: &[bool]) -> Vec<u32> {
    bits.iter()
        .enumerate()
        .filter(|&(_, &b)| b)
        .map(|(i, _)| i as u32)
        .collect()
}

fn canonical_bits(tag: &[u8], bits: &[bool]) -> Vec<u8> {
    let mut out = tag.to_vec();
    out.extend_from_slice(&(bits.len() as u64).to_le_bytes());
    for chunk in bits.chunks(8) {
        let mut byte = 0u8;
        for (i, &b) in chunk.iter().enumerate() {
            if b {
                byte |= 1 << i;
            }
        }
        out.push(byte);
    }
    out
}

/// Run one full decode on `code` with `cfg`, sampling error stream `seed`.
pub fn decode_once(code: &CssCode, cfg: &DecodeConfig, seed: u64) -> DecodeResult {
    let t = Tanner::build(&code.checks, code.n);
    let error = sample_error(code.n, cfg.p_fx, seed);
    let syndrome_bits = syndrome_of(&code.checks, &error);
    let core = relay_bp(&t, &syndrome_bits, seed, cfg.max_iter, cfg.max_legs);
    // residual = error XOR correction
    let residual: Vec<usize> = (0..code.n)
        .filter(|&i| error[i] ^ core.correction[i])
        .collect();
    let logical_success = core.converged && code.is_trivial(&residual);
    DecodeResult {
        code_id: code.id.clone(),
        n: code.n,
        error: bits_to_ids(&error),
        correction: bits_to_ids(&core.correction),
        syndrome: bits_to_ids(&syndrome_bits),
        converged: core.converged,
        logical_success,
        legs: core.legs,
        iters: core.iters,
        messages: core.messages,
        syndrome_bits,
        correction_bits: core.correction.clone(),
        check_matrix_hash: content_hash(&code.check_matrix_bytes()),
    }
}

impl DecodeResult {
    /// Seal this decode into a signed, joule-metered receipt. `joules_micro` is the
    /// decode's energy: on real hardware a RAPL/NVML meter; here a caller-supplied
    /// modeled figure (the `messages` work is the exact, deterministic basis).
    pub fn seal(
        &self,
        signer: &SigningKey,
        signer_id: impl Into<String>,
        joules_micro: u64,
        grant: GrantRef,
    ) -> DecodeReceipt {
        DecodeReceipt::seal(
            signer,
            signer_id,
            self.code_id.clone(),
            "relay-bp.min-sum",
            self.check_matrix_hash,
            content_hash(&canonical_bits(b"wai:qec-syndrome\x01", &self.syndrome_bits)),
            content_hash(&canonical_bits(b"wai:qec-correction\x01", &self.correction_bits)),
            self.converged,
            self.legs,
            self.messages,
            joules_micro,
            grant,
            None,
        )
    }

    /// Canonical correction bytes (what the receipt's `correction_hash` binds).
    pub fn correction_bytes(&self) -> Vec<u8> {
        canonical_bits(b"wai:qec-correction\x01", &self.correction_bits)
    }
}

/// A Monte-Carlo comparison of plain BP vs Relay-BP on the same error samples —
/// the demonstration that relay is what makes BP work on quantum codes.
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct BenchResult {
    pub trials: u32,
    pub plain_logical_failures: u32,
    pub relay_logical_failures: u32,
    pub relay_total_messages: u64,
}

/// Run `trials` decodes at rate `p`, each with plain BP (`max_legs=1`) and
/// Relay-BP, over identical error samples, counting logical failures of each.
pub fn benchmark(code: &CssCode, p: f64, trials: u32, seed: u64, max_legs: u32) -> BenchResult {
    let base = DecodeConfig::from_p(p);
    let plain = DecodeConfig { max_legs: 1, ..base };
    let relay = DecodeConfig { max_legs, ..base };
    let mut pf = 0;
    let mut rf = 0;
    let mut msg = 0u64;
    for tix in 0..trials {
        let s = seed.wrapping_add(tix as u64).wrapping_mul(0x9E37_79B9_7F4A_7C15);
        if !decode_once(code, &plain, s).logical_success {
            pf += 1;
        }
        let r = decode_once(code, &relay, s);
        if !r.logical_success {
            rf += 1;
        }
        msg += r.messages;
    }
    BenchResult { trials, plain_logical_failures: pf, relay_logical_failures: rf, relay_total_messages: msg }
}

#[cfg(test)]
mod tests {
    use super::*;

    fn key(s: u8) -> SigningKey {
        SigningKey::from_bytes(&[s; 32])
    }

    #[test]
    fn gf2_membership() {
        // rows {0,1}, {1,2} span {0,1},{1,2},{0,2}; not {0},{2}.
        let basis = Gf2Basis::from_rows(&[vec![0, 1], vec![1, 2]], 3);
        assert_eq!(basis.rank(), 2);
        assert!(basis.contains(&[0, 1]));
        assert!(basis.contains(&[1, 2]));
        assert!(basis.contains(&[0, 2]));
        assert!(basis.contains(&[])); // zero vector always in the space
        assert!(!basis.contains(&[0]));
        assert!(!basis.contains(&[2]));
    }


    /// A CSS code is only a code if the two families commute: every X stabilizer
    /// must meet every Z check in an even number of qubits. This is the test that
    /// decides whether a boundary convention is right — the first parities I tried
    /// for the surface code failed it.
    #[test]
    fn surface_code_is_a_valid_css_code() {
        for d in [3usize, 5, 7] {
            let code = CssCode::surface(d);
            assert_eq!(code.n, d * d, "d={d} qubit count");
            assert_eq!(code.k, 1, "d={d} must encode exactly one logical qubit");
            let x: Vec<&Vec<usize>> = code.stabs.iter().collect();
            for zc in &code.checks {
                for xs in &x {
                    let overlap = zc.iter().filter(|q| xs.contains(q)).count();
                    assert_eq!(overlap % 2, 0, "d={d}: X{xs:?} and Z{zc:?} anticommute");
                }
            }
            // stabilizer counts: (d^2-1)/2 of each family
            assert_eq!(code.checks.len(), (d * d - 1) / 2, "d={d} Z-check count");
            assert_eq!(x.len(), (d * d - 1) / 2, "d={d} X-stabilizer count");
            // planarity: every stabilizer is weight 4 (bulk) or weight 2 (boundary)
            for v in code.checks.iter().chain(x.iter().copied()) {
                assert!(v.len() == 4 || v.len() == 2, "d={d}: weight {} stabilizer", v.len());
            }
        }
    }

    /// The decoder must actually protect the code, and Relay-BP must earn its
    /// name: plain belief propagation degrades on quantum codes (degeneracy plus
    /// short cycles), which is the whole reason the relay legs exist.
    #[test]
    fn surface_code_decodes() {
        let d3 = CssCode::surface(3);
        let r = benchmark(&d3, 0.01, 500, 7, 4);
        let relay = r.relay_logical_failures as f64 / r.trials as f64;
        assert!(relay <= 0.05, "d=3 at p=0.01 should rarely fail: {relay}");
        assert!(
            r.relay_logical_failures <= r.plain_logical_failures,
            "relay {} must not do worse than plain BP {}",
            r.relay_logical_failures,
            r.plain_logical_failures
        );
    }

    /// Relay-BP's actual claim, on the code family it was designed for: the relay
    /// legs must beat plain belief propagation on a qLDPC (bivariate-bicycle) code.
    ///
    /// Note what is deliberately NOT asserted here: that a larger *topological*
    /// code decodes better. Measured across p = 0.001…0.02, distance-5 loses to
    /// distance-3 on both the surface code and the toric code — see
    /// `probe_surface_threshold`. That is a property of BP-family decoding on
    /// degenerate topological codes, not of the code constructions, and inventing
    /// a passing threshold for it would be asserting something untrue.
    #[test]
    fn relay_beats_plain_bp_on_qldpc() {
        let bb = CssCode::gross();
        let r = benchmark(&bb, 0.06, 400, 11, 16);
        assert!(
            r.relay_logical_failures < r.plain_logical_failures,
            "relay {} should beat plain BP {} on a qLDPC code",
            r.relay_logical_failures,
            r.plain_logical_failures
        );
    }

    #[test]
    #[ignore]
    fn probe_surface_threshold() {
        println!("\n-- relay legs on the code Relay-BP was designed for: BB [[72,12]] --");
        let bb = CssCode::gross();
        println!(" p       legs=1(plain)  legs=4    legs=16");
        for &pp in &[0.02f64, 0.04, 0.06] {
            let f = |legs: u32| {
                let r = benchmark(&bb, pp, 400, 11, legs);
                r.relay_logical_failures as f64 / r.trials as f64
            };
            let plain = {
                let r = benchmark(&bb, pp, 400, 11, 4);
                r.plain_logical_failures as f64 / r.trials as f64
            };
            println!(" {pp:<8} {plain:<13.4} {:<9.4} {:<9.4}", f(4), f(16));
        }
        println!("\n-- the same sweep on topological codes --");
        println!(" p       surf d=3  surf d=5  toric L=3 toric L=5");
        for &pp in &[0.001f64, 0.005, 0.01] {
            let sf = |d: usize| {
                let r = benchmark(&CssCode::surface(d), pp, 1500, 11, 16);
                r.relay_logical_failures as f64 / r.trials as f64
            };
            let tf = |l: usize| {
                let r = benchmark(&CssCode::toric(l), pp, 1500, 11, 16);
                r.relay_logical_failures as f64 / r.trials as f64
            };
            println!(" {pp:<8} {:<9.4} {:<9.4} {:<9.4} {:<9.4}", sf(3), sf(5), tf(3), tf(5));
        }
    }

    #[test]
    fn toric_code_shape() {
        let c = CssCode::toric(5);
        assert_eq!(c.n, 50);
        assert_eq!(c.checks.len(), 25);
        assert_eq!(c.k, 2, "toric code encodes 2 logical qubits");
        // every check has weight 4
        assert!(c.checks.iter().all(|r| r.len() == 4));
    }

    #[test]
    fn gross_code_is_72_12() {
        let c = CssCode::gross();
        assert_eq!(c.n, 72);
        assert_eq!(c.k, 12, "the gross-style BB code encodes 12 logical qubits");
        // qLDPC: bounded check weight (6 here), independent of n
        assert!(c.checks.iter().all(|r| r.len() == 6));
    }

    #[test]
    fn css_orthogonality_gross() {
        // H_X · H_Zᵀ = 0: every Z-check commutes with every X-stabilizer, i.e.
        // every decode check overlaps every stabilizer in an even number of qubits.
        let c = CssCode::gross();
        // stab_basis holds H_X's row space; reconstruct raw H_X rows to test.
        let hx = CssCode::bivariate_bicycle(6, 6, &[(3, 0), (0, 1), (0, 2)], &[(0, 3), (1, 0), (2, 0)]);
        for zc in &c.checks {
            let zset: std::collections::HashSet<usize> = zc.iter().copied().collect();
            for xs in raw_hx(&hx) {
                let overlap = xs.iter().filter(|q| zset.contains(q)).count();
                assert_eq!(overlap % 2, 0, "CSS orthogonality violated");
            }
        }
    }

    // helper: recover raw H_X rows (stab side) for the orthogonality test
    fn raw_hx(_c: &CssCode) -> Vec<Vec<usize>> {
        // rebuild H_X = [A|B] directly
        let (l, m) = (6usize, 6usize);
        let a = [(3usize, 0usize), (0, 1), (0, 2)];
        let b = [(0usize, 3usize), (1, 0), (2, 0)];
        let lm = l * m;
        let cell = |i: usize, j: usize| (i % l) * m + (j % m);
        let build = |monos: &[(usize, usize)]| -> Vec<Vec<usize>> {
            let mut rows = vec![Vec::new(); lm];
            for i in 0..l {
                for j in 0..m {
                    let col = cell(i, j);
                    for &(px, py) in monos {
                        rows[cell(i + px, j + py)].push(col);
                    }
                }
            }
            rows
        };
        let ar = build(&a);
        let br = build(&b);
        (0..lm)
            .map(|r| {
                let mut row = ar[r].clone();
                row.extend(br[r].iter().map(|&c| c + lm));
                row.sort_unstable();
                row
            })
            .collect()
    }

    #[test]
    fn decode_recovers_low_weight_errors() {
        // At a low error rate, decoding should almost always succeed logically.
        let c = CssCode::toric(5);
        let cfg = DecodeConfig::from_p(0.02);
        let mut fails = 0;
        for s in 0..200u64 {
            if !decode_once(&c, &cfg, s.wrapping_mul(0x9E37_79B9)).logical_success {
                fails += 1;
            }
        }
        assert!(fails < 20, "too many logical failures at p=0.02: {fails}/200");
    }

    #[test]
    fn decode_is_deterministic() {
        let c = CssCode::gross();
        let cfg = DecodeConfig::from_p(0.03);
        let a = decode_once(&c, &cfg, 12345);
        let b = decode_once(&c, &cfg, 12345);
        assert_eq!(a.correction, b.correction);
        assert_eq!(a.messages, b.messages);
        assert_eq!(a.legs, b.legs);
    }

    #[test]
    fn relay_beats_plain_bp() {
        // The whole point of Relay-BP: on a quantum code, relay legs recover
        // decodes that plain BP leaves failed. Relay failures ≤ plain failures.
        let c = CssCode::gross();
        let bench = benchmark(&c, 0.04, 150, 0xBEEF, 12);
        assert!(
            bench.relay_logical_failures <= bench.plain_logical_failures,
            "relay {} must not exceed plain {}",
            bench.relay_logical_failures, bench.plain_logical_failures
        );
    }

    #[test]
    fn decode_seals_verifying_receipt() {
        let c = CssCode::gross();
        let cfg = DecodeConfig::from_p(0.03);
        let r = decode_once(&c, &cfg, 7);
        let rec = r.seal(&key(1), "did:key:lab", 500_000, GrantRef::unbounded("quantum.decode"));
        assert!(rec.verify());
        assert!(rec.correction_matches(&r.correction_bytes()));
        assert_eq!(rec.code_id, "bb:[[72,12]]");
        assert_eq!(rec.work_messages, r.messages);
    }
}