wai-quantum 0.3.7

//! 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.
//!
//! # The codes verify their own claims
//!
//! A code is named `[[n, k, d]]`, and `d` is the half that matters: below it,
//! correction is guaranteed rather than likely. [`x_distance`] computes it exactly
//! by exhaustion, so `surface(3)` is *shown* to have distance 3 and `surface(5)`
//! distance 5 — every lighter logical operator ruled out — rather than being
//! assumed from the constructor's name. The decoder is held to the same standard:
//! every error of weight below `d` is decoded and checked, exhaustively.
//!
//! # Two decoders, and which code each is for
//!
//! **Relay-BP** is a belief-propagation decoder and 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 merely underperform — it inverts the
//! property that makes a code family worth having. Measured across `p = 0.001`
//! to `0.02`, a distance-5 surface code decodes *worse* than distance-3 under
//! Relay-BP, and the toric code shows the same inversion. That is the known
//! weakness of BP on degenerate codes with short cycles, not a defect in the code
//! constructions (both verify as valid CSS codes with the expected `k`).
//!
//! [`union_find_correct`] is the answer for those codes — a matching decoder that
//! pairs syndromes on the code's own graph. It restores the property BP loses
//! (surface code, logical failure rate, 2000 trials):
//!
//! | `p`   | d=3 BP | d=3 UF | d=5 BP | d=5 UF | d=7 BP | d=7 UF |
//! |-------|--------|--------|--------|--------|--------|--------|
//! | 0.005 | 0.0205 | 0.0010 | 0.0410 | 0.0000 | 0.0605 | 0.0000 |
//! | 0.01  | 0.0370 | 0.0020 | 0.0785 | 0.0000 | 0.1170 | 0.0000 |
//! | 0.02  | 0.0715 | 0.0090 | 0.1535 | 0.0020 | 0.2365 | 0.0005 |
//!
//! # Faulty measurement: decoding in space-time
//!
//! Everything above assumes the syndrome is read perfectly, which no device does.
//! [`decode_phenomenological`] runs `rounds` noisy measurement rounds plus a final
//! perfect readout and decodes the block as a **space-time** matching problem:
//! detectors are *differences* between consecutive rounds, so a data error lights
//! two checks at the same time while a lying measurement lights one check at two
//! adjacent times. Surface code, `q = p`, `d` rounds:
//!
//! | `p`   | d=3    | d=5    | d=7    |
//! |-------|--------|--------|--------|
//! | 0.003 | 0.0025 | 0.0000 | 0.0000 |
//! | 0.005 | 0.0050 | 0.0008 | 0.0008 |
//! | 0.01  | 0.0142 | 0.0042 | 0.0067 |
//! | 0.02  | 0.0475 | 0.0458 | 0.0417 |
//!
//! Distance pays below `p ≈ 1%` and the curves converge by `2%` — a threshold in
//! the region the literature puts the phenomenological surface code, which is the
//! outcome that would have exposed a broken implementation.
//!
//! Edges are weighted by `−ln(rate)`, so the decoder blames the likelier fault
//! rather than the nearer one. That matters more than it sounds: with data noise
//! off and measurement noise on, an unweighted decoder happily "explains" a lying
//! measurement by flipping data qubits, and enough of those in a line is a logical
//! error invented out of nothing.
//!
//! Distance now buys what it is supposed to buy, monotonically. Use union-find for
//! surface and toric codes, Relay-BP for the higher-degree qLDPC families where a
//! qubit touches more than two checks and the matching graph does not exist.
//!
//! 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.
/// A matching-style **union-find decoder** (Delfosse–Nickerson).
///
/// Belief propagation struggles on topological codes: they are degenerate and
/// their Tanner graphs are full of short cycles, and the measurements in this
/// module's tests show the consequence — under Relay-BP a distance-5 surface code
/// decodes *worse* than distance-3. Union-find decodes the syndrome the way the
/// code's structure actually wants: as a matching problem on a graph whose nodes
/// are checks and whose edges are qubits.
///
/// Grow every odd cluster until no odd cluster remains, then peel each cluster's
/// spanning forest from the leaves inward to read off a correction. Almost-linear
/// time, and unlike BP it makes the distance count.
///
/// # Applicability
///
/// This only makes sense when a qubit touches at most two checks — true for
/// surface and toric codes, false for the higher-degree qLDPC families (a
/// bivariate-bicycle qubit sits in three). Returns `None` rather than guessing
/// when the code is not a matching graph; use Relay-BP for those.
/// Advance `c` to the next `w`-combination of `0..n` in lexicographic order.
/// Returns `false` once the last combination has been passed.
fn next_combination(c: &mut [usize], n: usize) -> bool {
    let w = c.len();
    let mut i = w;
    while i > 0 {
        i -= 1;
        if c[i] < n - (w - i) {
            c[i] += 1;
            for j in i + 1..w {
                c[j] = c[j - 1] + 1;
            }
            return true;
        }
    }
    false
}

/// The code's **X-distance**: the least weight of an undetectable error that is
/// not a stabilizer — a logical operator.
///
/// This is the `d` in `[[n, k, d]]`, and it is the claim a code family actually
/// rests on: any error of weight below `d` is correctable, because no error that
/// small can be logical. Asserting `k` while leaving `d` unchecked verifies the
/// less interesting half of the name.
///
/// Searched exactly by increasing weight, so a returned `d` is a proof: every
/// lighter operator has been ruled out by exhaustion. `None` means no logical of
/// weight `≤ max_weight` exists. The search is exponential in the cap — `d ≤ 5`
/// is instant, `d = 7` is not.
pub fn x_distance(code: &CssCode, max_weight: usize) -> Option<usize> {
    let n = code.n;
    for w in 1..=max_weight.min(n) {
        let mut combo: Vec<usize> = (0..w).collect();
        loop {
            let mut bits = vec![false; n];
            for &q in &combo {
                bits[q] = true;
            }
            // undetectable (zero syndrome) and not a stabilizer => logical
            if syndrome_of(&code.checks, &bits).iter().all(|b| !b) && !code.is_trivial(&combo) {
                return Some(w);
            }
            if !next_combination(&mut combo, n) {
                break;
            }
        }
    }
    None
}

/// The union-find core, on an abstract graph.
///
/// `edges` are `(u, v)` node pairs; `boundary` is the node that may absorb any
/// leftover parity (open boundaries, and the time boundary of a space-time graph).
/// Returns which edges the decoder selected.
///
/// Grow every odd cluster until none remain, then peel each cluster's spanning
/// forest from the leaves inward. Growth adds a **half-edge from each odd
/// endpoint**, so an edge joining two syndromes closes in one round while an edge
/// to a neutral node needs two — that asymmetry is what makes the decoder pair
/// syndromes with each other rather than draining them all to the boundary.
///
/// `weights` are per-edge costs: an edge closes at `2 × weight`, so cheap edges
/// close first. Uniform weights recover plain union-find; deriving them from the
/// error rates is what stops the decoder explaining a lying measurement with a
/// data correction when data errors are the less likely cause.
fn uf_decode_graph(
    nodes: usize,
    boundary: usize,
    edges: &[(usize, usize)],
    weights: &[u32],
    syndrome: &[bool],
) -> Vec<bool> {
    let mut parent: Vec<usize> = (0..nodes).collect();
    fn find(parent: &mut Vec<usize>, mut x: usize) -> usize {
        while parent[x] != x {
            parent[x] = parent[parent[x]];
            x = parent[x];
        }
        x
    }
    let mut syn: Vec<bool> = (0..nodes).map(|i| syndrome.get(i).copied().unwrap_or(false)).collect();
    syn[boundary] = false;
    let mut grown: Vec<u32> = vec![0; edges.len()];
    let maxw = weights.iter().copied().max().unwrap_or(1).max(1) as usize;

    for _ in 0..(4 * nodes * maxw + 8) {
        let mut parity = vec![false; nodes];
        let mut touches = vec![false; nodes];
        for v in 0..nodes {
            let r = find(&mut parent, v);
            if syn[v] {
                parity[r] ^= true;
            }
            if v == boundary {
                touches[r] = true;
            }
        }
        if !(0..nodes).any(|v| find(&mut parent, v) == v && parity[v] && !touches[v]) {
            break;
        }
        let mut to_union: Vec<usize> = Vec::new();
        for (ei, &(a, b)) in edges.iter().enumerate() {
            if grown[ei] >= 2 * weights.get(ei).copied().unwrap_or(1).max(1) {
                continue;
            }
            let (ra, rb) = (find(&mut parent, a), find(&mut parent, b));
            let odd_a = parity[ra] && !touches[ra];
            let odd_b = parity[rb] && !touches[rb];
            let inc = if ra == rb { u8::from(odd_a) } else { u8::from(odd_a) + u8::from(odd_b) };
            if inc > 0 {
                let cap = 2 * weights.get(ei).copied().unwrap_or(1).max(1);
                grown[ei] = grown[ei].saturating_add(inc as u32);
                if grown[ei] >= cap {
                    to_union.push(ei);
                }
            }
        }
        for ei in to_union {
            let (a, b) = edges[ei];
            let (ra, rb) = (find(&mut parent, a), find(&mut parent, b));
            if ra != rb {
                parent[ra] = rb;
            }
        }
    }

    let mut adj: Vec<Vec<(usize, usize)>> = vec![Vec::new(); nodes];
    for (ei, &(a, b)) in edges.iter().enumerate() {
        if grown[ei] >= 2 * weights.get(ei).copied().unwrap_or(1).max(1) {
            adj[a].push((b, ei));
            adj[b].push((a, ei));
        }
    }
    let mut picked = vec![false; edges.len()];
    let mut seen = vec![false; nodes];
    for root in std::iter::once(boundary).chain(0..nodes) {
        if seen[root] {
            continue;
        }
        let mut order = vec![root];
        let mut parent_edge: Vec<Option<(usize, usize)>> = vec![None; nodes];
        seen[root] = true;
        let mut i = 0;
        while i < order.len() {
            let v = order[i];
            i += 1;
            for &(w, ei) in &adj[v] {
                if !seen[w] {
                    seen[w] = true;
                    parent_edge[w] = Some((v, ei));
                    order.push(w);
                }
            }
        }
        for &v in order.iter().rev() {
            if v == root || !syn[v] {
                continue;
            }
            if let Some((u, ei)) = parent_edge[v] {
                picked[ei] ^= true;
                syn[v] = false;
                syn[u] ^= true;
            }
        }
    }
    picked
}

/// A matching-style **union-find decoder** (Delfosse–Nickerson) for one shot of
/// syndrome, assuming perfect measurement.
///
/// Belief propagation struggles on topological codes: they are degenerate and
/// their Tanner graphs are full of short cycles. Union-find decodes the syndrome
/// the way the code's structure wants — as matching on a graph whose nodes are
/// checks and whose edges are qubits — and unlike BP it makes the distance count.
///
/// # Applicability
///
/// Only meaningful when a qubit touches at most two checks: true for surface and
/// toric codes, false for higher-degree qLDPC (a bivariate-bicycle qubit sits in
/// three). Returns `None` rather than guessing; use Relay-BP for those.
pub fn union_find_correct(code: &CssCode, syndrome: &[bool]) -> Option<Vec<bool>> {
    let m = code.checks.len();
    let boundary = m;
    let mut owners: Vec<Vec<usize>> = vec![Vec::new(); code.n];
    for (ci, chk) in code.checks.iter().enumerate() {
        for &q in chk {
            owners[q].push(ci);
        }
    }
    let mut edges: Vec<(usize, usize)> = Vec::new();
    let mut edge_qubit: Vec<usize> = Vec::new();
    for (q, own) in owners.iter().enumerate() {
        match own.len() {
            0 => {}
            1 => {
                edges.push((own[0], boundary));
                edge_qubit.push(q);
            }
            2 => {
                edges.push((own[0], own[1]));
                edge_qubit.push(q);
            }
            _ => return None,
        }
    }
    let weights = vec![1u32; edges.len()];
    let picked = uf_decode_graph(m + 1, boundary, &edges, &weights, syndrome);
    let mut correction = vec![false; code.n];
    for (ei, &on) in picked.iter().enumerate() {
        if on {
            correction[edge_qubit[ei]] ^= true;
        }
    }
    Some(correction)
}

/// Repeated-round noise with **faulty syndrome measurement** — the
/// phenomenological model.
///
/// Single-shot decoding assumes the syndrome itself is read perfectly, which no
/// device does. Here each of `rounds` measurement rounds can misreport a check,
/// so a lone flipped bit is ambiguous: a data error, or a lying measurement?
/// The resolution is to decode in **space-time**. Detectors are *differences*
/// between consecutive rounds, so a measurement error lights two detectors in the
/// same check at adjacent times, while a data error lights two checks at the same
/// time. A final perfect round closes the block.
#[derive(Clone, Copy, Debug)]
pub struct PhenomConfig {
    /// Noisy measurement rounds (a perfect readout round is appended).
    pub rounds: u32,
    /// Data error per qubit per round (fixed-point, scale 2^12).
    pub p_data_fx: i64,
    /// Measurement error per check per noisy round (fixed-point, scale 2^12).
    pub q_meas_fx: i64,
}

impl PhenomConfig {
    pub fn from_p(rounds: u32, p_data: f64, q_meas: f64) -> PhenomConfig {
        PhenomConfig {
            rounds,
            p_data_fx: (p_data * 4096.0) as i64,
            q_meas_fx: (q_meas * 4096.0) as i64,
        }
    }
}

/// One space-time decode.
#[derive(Clone, Debug)]
pub struct PhenomResult {
    pub rounds: u32,
    /// Detectors that fired, as `round * n_checks + check`.
    pub detectors: Vec<u32>,
    /// Accumulated data error over the whole block.
    pub data_error: Vec<u32>,
    pub correction: Vec<u32>,
    pub logical_success: bool,
}

/// Decode a phenomenological-noise block with union-find over the space-time
/// graph. Returns `None` when the code is not a matching graph.
pub fn decode_phenomenological(
    code: &CssCode,
    cfg: &PhenomConfig,
    seed: u64,
) -> Option<PhenomResult> {
    let m = code.checks.len();
    let t_rounds = cfg.rounds.max(1) as usize;
    let layers = t_rounds + 1; // + the final perfect readout

    // qubit -> its checks (matching graph precondition)
    let mut owners: Vec<Vec<usize>> = vec![Vec::new(); code.n];
    for (ci, chk) in code.checks.iter().enumerate() {
        for &q in chk {
            owners[q].push(ci);
        }
    }
    if owners.iter().any(|o| o.len() > 2) {
        return None;
    }

    // --- simulate the block ---
    let mut st = seed ^ 0x51ED_2701_A55A_1234;
    let draw = |rate: i64, st: &mut u64| -> bool {
        *st = st
            .wrapping_mul(6364136223846793005)
            .wrapping_add(1442695040888963407);
        (((*st >> 40) & 0xFFF) as i64) < rate
    };
    let mut cum_err = vec![false; code.n];
    let mut s_prev = vec![false; m];
    let mut detector_bits = vec![false; layers * m];
    for t in 0..layers {
        if t < t_rounds {
            for q in 0..code.n {
                if draw(cfg.p_data_fx, &mut st) {
                    cum_err[q] ^= true;
                }
            }
        }
        let s_true = syndrome_of(&code.checks, &cum_err);
        let mut s_obs = s_true.clone();
        if t < t_rounds {
            for c in 0..m {
                if draw(cfg.q_meas_fx, &mut st) {
                    s_obs[c] ^= true;
                }
            }
        }
        for c in 0..m {
            detector_bits[t * m + c] = s_obs[c] ^ s_prev[c];
        }
        s_prev = s_obs;
    }

    // --- space-time matching graph ---
    let boundary = layers * m;
    let node = |c: usize, t: usize| t * m + c;
    // Cost of an edge is -ln(rate): the less likely a fault, the more the decoder
    // must pay to blame it. A rate of zero means the fault cannot happen, so the
    // edge is omitted rather than merely made expensive.
    let cost = |rate_fx: i64| -> u32 {
        let r = rate_fx as f64 / 4096.0;
        if r <= 0.0 { 0 } else { (-r.ln()).round().clamp(1.0, 12.0) as u32 }
    };
    let w_space = cost(cfg.p_data_fx);
    let w_time = cost(cfg.q_meas_fx);
    let mut edges: Vec<(usize, usize)> = Vec::new();
    let mut weights: Vec<u32> = Vec::new();
    // space edges: a data error at layer t lights its two checks at layer t
    let mut space_qubit: Vec<Option<usize>> = Vec::new();
    if w_space > 0 {
        for t in 0..layers {
            for (q, own) in owners.iter().enumerate() {
                match own.len() {
                    1 => {
                        edges.push((node(own[0], t), boundary));
                        weights.push(w_space);
                        space_qubit.push(Some(q));
                    }
                    2 => {
                        edges.push((node(own[0], t), node(own[1], t)));
                        weights.push(w_space);
                        space_qubit.push(Some(q));
                    }
                    _ => {}
                }
            }
        }
    }
    // time edges: a measurement error at round t lights the same check at t and t+1
    if w_time > 0 {
        for t in 0..layers - 1 {
            for c in 0..m {
                edges.push((node(c, t), node(c, t + 1)));
                weights.push(w_time);
                space_qubit.push(None);
            }
        }
    }

    let picked = uf_decode_graph(boundary + 1, boundary, &edges, &weights, &detector_bits);
    let mut correction = vec![false; code.n];
    for (ei, &on) in picked.iter().enumerate() {
        if on && let Some(q) = space_qubit[ei] {
            correction[q] ^= true;
        }
    }

    let residual: Vec<usize> = (0..code.n).filter(|&i| cum_err[i] ^ correction[i]).collect();
    Some(PhenomResult {
        rounds: cfg.rounds,
        detectors: bits_to_ids(&detector_bits),
        data_error: bits_to_ids(&cum_err),
        correction: bits_to_ids(&correction),
        logical_success: code.is_trivial(&residual),
    })
}

/// Logical failure rate under phenomenological noise.
pub fn benchmark_phenomenological(
    code: &CssCode,
    cfg: &PhenomConfig,
    trials: u32,
    seed: u64,
) -> (u32, u32) {
    let mut fails = 0;
    for tix in 0..trials {
        let s = seed.wrapping_add(tix as u64).wrapping_mul(0x9E37_79B9_7F4A_7C15);
        match decode_phenomenological(code, cfg, s) {
            Some(r) if r.logical_success => {}
            _ => fails += 1,
        }
    }
    (trials, fails)
}

/// One union-find decode of a sampled error — the matching counterpart to
/// [`decode_once`].
pub fn decode_once_uf(code: &CssCode, cfg: &DecodeConfig, seed: u64) -> DecodeResult {
    let error = sample_error(code.n, cfg.p_fx, seed);
    let syndrome_bits = syndrome_of(&code.checks, &error);
    let (correction, converged) = match union_find_correct(code, &syndrome_bits) {
        Some(c) => (c, true),
        None => (vec![false; code.n], false),
    };
    let residual: Vec<usize> = (0..code.n).filter(|&i| error[i] ^ correction[i]).collect();
    let logical_success = converged && code.is_trivial(&residual);
    DecodeResult {
        code_id: code.id.clone(),
        n: code.n,
        error: bits_to_ids(&error),
        correction: bits_to_ids(&correction),
        syndrome: bits_to_ids(&syndrome_bits),
        converged,
        logical_success,
        legs: 0,
        iters: 0,
        messages: 0,
        syndrome_bits,
        correction_bits: correction,
        check_matrix_hash: content_hash(&code.check_matrix_bytes()),
    }
}

/// Union-find logical failure rate over `trials` sampled errors.
pub fn benchmark_uf(code: &CssCode, p: f64, trials: u32, seed: u64) -> (u32, u32) {
    let cfg = DecodeConfig::from_p(p);
    let mut fails = 0;
    for tix in 0..trials {
        let s = seed.wrapping_add(tix as u64).wrapping_mul(0x9E37_79B9_7F4A_7C15);
        if !decode_once_uf(code, &cfg, s).logical_success {
            fails += 1;
        }
    }
    (trials, fails)
}

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));
        }
    }



    /// The property BP could not deliver: with a matching decoder, distance counts.
    /// Under Relay-BP a distance-5 surface code decoded WORSE than distance-3;
    /// union-find inverts that decisively.
    #[test]
    fn union_find_makes_distance_count() {
        let p = 0.01;
        let rate = |d: usize| {
            let (t, f) = benchmark_uf(&CssCode::surface(d), p, 1500, 11);
            f as f64 / t as f64
        };
        let (r3, r5) = (rate(3), rate(5));
        assert!(r5 < r3, "d=5 ({r5}) must beat d=3 ({r3}) below threshold");
        assert!(r5 <= 0.005, "d=5 at p={p} should almost never fail: {r5}");
    }

    /// And it must beat belief propagation on the codes BP is bad at.
    #[test]
    fn union_find_beats_bp_on_topological_codes() {
        let p = 0.02;
        for d in [5usize, 7] {
            let code = CssCode::surface(d);
            let bp = {
                let r = benchmark(&code, p, 1000, 11, 16);
                r.relay_logical_failures as f64 / r.trials as f64
            };
            let (t, f) = benchmark_uf(&code, p, 1000, 11);
            let uf = f as f64 / t as f64;
            assert!(uf < bp, "d={d}: union-find {uf} should beat Relay-BP {bp}");
        }
    }

    /// A qubit in three checks is not a matching graph. Refuse, do not guess —
    /// Relay-BP is the decoder for those families.
    #[test]
    fn union_find_refuses_non_matching_codes() {
        assert!(
            union_find_correct(&CssCode::gross(), &vec![false; CssCode::gross().checks.len()]).is_none(),
            "a bivariate-bicycle qubit sits in three checks; UF must decline"
        );
        assert!(union_find_correct(&CssCode::surface(5), &vec![false; 12]).is_some());
    }

    /// Ground truth: every single-qubit error must be corrected exactly.
    #[test]
    fn union_find_corrects_every_single_qubit_error() {
        for d in [3usize, 5] {
            let code = CssCode::surface(d);
            for q in 0..code.n {
                let mut err = vec![false; code.n];
                err[q] = true;
                let syn = syndrome_of(&code.checks, &err);
                let corr = union_find_correct(&code, &syn).expect("matching graph");
                let residual: Vec<usize> = (0..code.n).filter(|&i| err[i] ^ corr[i]).collect();
                assert!(
                    code.is_trivial(&residual),
                    "d={d}: single error on qubit {q} left a logical residual {residual:?}"
                );
            }
        }
    }



    /// Measurement errors alone must never become a logical error. With no data
    /// noise the true error is empty, so any space edge the decoder picks IS the
    /// residual — this asserts it explains lying measurements with time edges
    /// instead of inventing data corrections.
    #[test]
    fn measurement_noise_alone_never_causes_a_logical_error() {
        for d in [3usize, 5] {
            let code = CssCode::surface(d);
            let cfg = PhenomConfig::from_p(d as u32, 0.0, 0.05);
            for t in 0..300u64 {
                let r = decode_phenomenological(&code, &cfg, t.wrapping_mul(0x9E37_79B9))
                    .expect("matching graph");
                assert!(
                    r.logical_success,
                    "d={d} seed={t}: pure measurement noise produced a logical failure \
                     (correction {:?} on an empty data error)",
                    r.correction
                );
            }
        }
    }

    /// Below threshold, distance must still buy protection when the syndrome
    /// itself is unreliable — the property that matters for a real device.
    #[test]
    fn phenomenological_distance_helps() {
        let p = 0.005;
        let rate = |d: usize| {
            let cfg = PhenomConfig::from_p(d as u32, p, p);
            let (t, f) = benchmark_phenomenological(&CssCode::surface(d), &cfg, 1200, 11);
            f as f64 / t as f64
        };
        let (r3, r5) = (rate(3), rate(5));
        assert!(r5 < r3, "d=5 ({r5}) must beat d=3 ({r3}) under noisy measurement");
    }

    /// With perfect measurement the space-time decoder must reduce to the
    /// single-shot one.
    #[test]
    fn perfect_measurement_reduces_to_single_shot() {
        for d in [3usize, 5] {
            let code = CssCode::surface(d);
            let cfg = PhenomConfig::from_p(1, 0.02, 0.0);
            let (t1, f1) = benchmark_phenomenological(&code, &cfg, 1200, 11);
            let (t2, f2) = benchmark_uf(&code, 0.02, 1200, 11);
            let (a, b) = (f1 as f64 / t1 as f64, f2 as f64 / t2 as f64);
            assert!((a - b).abs() < 0.02, "d={d}: spacetime {a} vs single-shot {b}");
        }
    }

    #[test]
    fn phenomenological_refuses_non_matching_codes() {
        let cfg = PhenomConfig::from_p(3, 0.01, 0.01);
        assert!(decode_phenomenological(&CssCode::gross(), &cfg, 1).is_none());
    }


    /// The `d` in `[[d², 1, d]]` — verified, not asserted by naming. A code whose
    /// distance is smaller than advertised silently corrects less than it claims.
    #[test]
    fn surface_code_distance_is_d() {
        for d in [3usize, 5] {
            let code = CssCode::surface(d);
            assert_eq!(
                x_distance(&code, d + 1),
                Some(d),
                "surface d={d} should have X-distance exactly {d}"
            );
        }
    }

    /// The toric code's distance is L, and its k is 2 — a different point in the
    /// same family, so this also guards the distance search itself.
    #[test]
    fn toric_code_distance_is_l() {
        let code = CssCode::toric(3);
        assert_eq!(code.k, 2, "toric encodes two logical qubits");
        assert_eq!(x_distance(&code, 4), Some(3));
    }

    /// Below the distance, correction must be guaranteed — not merely likely.
    /// Every error of weight < d, exhaustively, must decode to a trivial residual.
    #[test]
    fn union_find_corrects_every_error_below_the_distance() {
        for d in [3usize, 5] {
            let code = CssCode::surface(d);
            let w = (d - 1) / 2; // guaranteed-correctable weight
            let mut combo: Vec<usize> = (0..w).collect();
            let mut checked = 0u32;
            loop {
                let mut err = vec![false; code.n];
                for &q in &combo {
                    err[q] = true;
                }
                let syn = syndrome_of(&code.checks, &err);
                let corr = union_find_correct(&code, &syn).expect("matching graph");
                let residual: Vec<usize> =
                    (0..code.n).filter(|&i| err[i] ^ corr[i]).collect();
                assert!(
                    code.is_trivial(&residual),
                    "d={d}: weight-{w} error {combo:?} decoded to a logical residual {residual:?}"
                );
                checked += 1;
                if !next_combination(&mut combo, code.n) {
                    break;
                }
            }
            assert!(checked > 0, "d={d}: nothing checked");
        }
    }

    #[test]
    #[ignore]
    fn probe_phenomenological() {
        println!("\n  surface code, phenomenological noise (data p, measurement q=p), d rounds");
        println!("  p        d=3       d=5       d=7");
        for &pp in &[0.001f64, 0.003, 0.005, 0.01, 0.02] {
            let r = |d: usize| {
                let cfg = PhenomConfig::from_p(d as u32, pp, pp);
                let (t, f) = benchmark_phenomenological(&CssCode::surface(d), &cfg, 1200, 11);
                f as f64 / t as f64
            };
            println!("  {pp:<8} {:<9.4} {:<9.4} {:<9.4}", r(3), r(5), r(7));
        }
        println!("\n  sanity: q=0 (perfect measurement) should track the single-shot decoder");
        for &pp in &[0.005f64, 0.02] {
            let st = |d: usize| {
                let cfg = PhenomConfig::from_p(1, pp, 0.0);
                let (t, f) = benchmark_phenomenological(&CssCode::surface(d), &cfg, 1200, 11);
                f as f64 / t as f64
            };
            let ss = |d: usize| { let (t,f)=benchmark_uf(&CssCode::surface(d), pp, 1200, 11); f as f64/t as f64 };
            println!("  p={pp:<7} d=3 spacetime {:<8.4} single-shot {:<8.4} | d=5 spacetime {:<8.4} single-shot {:<8.4}",
                st(3), ss(3), st(5), ss(5));
        }
    }

    #[test]
    #[ignore]
    fn probe_union_find_vs_bp() {
        println!("\n  surface code — logical failure rate");
        println!("  p        d=3 BP    d=3 UF    d=5 BP    d=5 UF    d=7 BP    d=7 UF");
        for &pp in &[0.001f64, 0.005, 0.01, 0.02, 0.05] {
            let bp = |d: usize| {
                let r = benchmark(&CssCode::surface(d), pp, 2000, 11, 16);
                r.relay_logical_failures as f64 / r.trials as f64
            };
            let uf = |d: usize| {
                let (t, f) = benchmark_uf(&CssCode::surface(d), pp, 2000, 11);
                f as f64 / t as f64
            };
            println!("  {pp:<8} {:<9.4} {:<9.4} {:<9.4} {:<9.4} {:<9.4} {:<9.4}",
                bp(3), uf(3), bp(5), uf(5), bp(7), uf(7));
        }
        println!("\n  toric code — UF");
        for &pp in &[0.005f64, 0.02] {
            let uf = |l: usize| { let (t,f)=benchmark_uf(&CssCode::toric(l), pp, 2000, 11); f as f64/t as f64 };
            println!("  p={pp:<7} L=3 {:<8.4} L=5 {:<8.4}", uf(3), uf(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);
    }
}