wai/wai-quantum-serve/Cargo.lock

@@ -589,7 +589,7 @@ dependencies = [
 
 [[package]]
 name = "wai-quantum"
-version = "0.3.29"
+version = "0.3.30"
 dependencies = [
  "blake3",
  "ed25519-dalek",

wai/wai-quantum/Cargo.lock

@@ -415,7 +415,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a"
 
 [[package]]
 name = "wai-quantum"
-version = "0.3.29"
+version = "0.3.30"
 dependencies = [
  "blake3",
  "ed25519-dalek",

wai/wai-quantum/Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "wai-quantum"
-version = "0.3.29"
+version = "0.3.30"
 edition = "2024"
 rust-version = "1.97"
 license = "Apache-2.0"
@@ -68,6 +68,9 @@ quantum_sv = []
 # Architecture-level estimation: yoked storage, magic-state cultivation, CCZ
 # factories and reaction-limited time, assembled into whole-machine estimates.
 quantum_arch = []
+# qLDPC memory at circuit level: bivariate-bicycle codes, their depth-8 syndrome
+# cycle, and BP-OSD decoding of the circuit's error model.
+quantum_ldpc = ["quantum_frame"]
 # Mixed states and the observables that distinguish quantum theory from a classical
 # account: density matrices, the partial trace, exact Pauli expectations, CHSH.
 quantum_info = ["quantum"]
@@ -107,13 +110,17 @@ quantum_phasor_meter = ["_signing"]
 
 # Everything.
 full = [
-    "quantum", "quantum_toolchain", "quantum_mps", "quantum_stabilizer", "quantum_pauli", "quantum_spd", "quantum_spd_receipt", "quantum_bptn", "quantum_resource", "quantum_resource_receipt", "quantum_frame", "quantum_tn", "quantum_gbs", "quantum_sv", "quantum_arch", "quantum_qir", "quantum_info", "quantum_qudit", "quantum_comm", "quantum_source", "quantum_channel",
+    "quantum", "quantum_toolchain", "quantum_mps", "quantum_stabilizer", "quantum_pauli", "quantum_spd", "quantum_spd_receipt", "quantum_bptn", "quantum_resource", "quantum_resource_receipt", "quantum_frame", "quantum_tn", "quantum_gbs", "quantum_sv", "quantum_arch", "quantum_ldpc", "quantum_qir", "quantum_info", "quantum_qudit", "quantum_comm", "quantum_source", "quantum_channel",
     "quantum_receipt", "quantum_ops", "quantum_cal", "quantum_mitigate", "quantum_qec",
     "quantum_control", "quantum_noise", "quantum_compile", "quantum_atom", "quantum_qbom",
     "quantum_vml", "quantum_phasor", "quantum_qfhrr", "quantum_kernel", "quantum_qdata",
     "quantum_phasor_meter",
 ]
 
+[[example]]
+name = "ldpc_memory"
+required-features = ["quantum_ldpc"]
+
 [[example]]
 name = "arch_factoring"
 required-features = ["quantum_arch"]

wai/wai-quantum/README.md

@@ -64,6 +64,7 @@ never creates information.
 | `quantum_bptn` | **tensor-network states on the hardware graph**: simple-update evolution with belief-propagation environments and re-gauging, the device's couplings applied a layer at a time in parallel (bit-identical to one at a time). Exact on trees; on the 127-qubit lattice it gives the exactly computed 5-step magnetization to 4e-7 at bond dimension 32 — an independent check on sparse Pauli dynamics, whose errors it does not share |
 | `quantum_resource` | **fault-tolerant resource estimation**: from a program's logical counts (qubits, T, Toffoli and rotation gates, measurements) to code distance, T-state factories found by exhaustive search, physical qubits and runtime, on six reference qubits and three codes, with the qubits-versus-time frontier. Reproduces the published worked examples to their printed precision |
 | `quantum_resource_receipt` | a signed, energy-accounted receipt over an estimate, re-checkable bit for bit on any machine |
+| `quantum_ldpc` | **qLDPC memory at circuit level**: bivariate-bicycle codes from their two polynomials (`k` by rank, logical bases), the depth-8 syndrome cycle with circuit-level noise, its detector error model, and BP-OSD decoding (min-sum with adaptive scaling, ordered statistics with a combination sweep). The error model has the same columns as the source's own decoding matrices, and on identical syndromes the decoder makes the source's decisions shot for shot |
 | `quantum_arch` | **architecture-level estimation**: yoked surface-code storage (row and grid parity checks over patches, hot and cold), magic-state cultivation read off the published simulation curves, CCZ factories, and a clock set by whichever is slowest of lattice surgery, the magic-state supply and the control loop. Assembles whole-machine estimates; its storage search matches the source's footprint tool layout for layout |
 | `quantum_receipt` | signed, energy-accounted receipt over a reconstruction |
 | `quantum_ops` | the operations/attestation layer |
@@ -153,6 +154,30 @@ search limits, to the bit). Re-running it anywhere must reproduce every figure,
 down to each factory round, and the bytes are the same natively and as a WASI
 component.
 
+## qLDPC memory
+
+```sh
+wai-quantum ldpc --code 144 --p 0.004 --cycles 12 --shots 1000   # failures per cycle, both halves
+wai-quantum ldpc --circuit --code 72 --cycles 6 --p 0.003        # the circuit itself
+```
+
+High-rate codes are where most of the projected qubit savings in fault
+tolerance come from, and their headline numbers are logical failure rates per
+syndrome cycle under circuit noise. `quantum_ldpc` re-derives them for the
+bivariate-bicycle family (arXiv:2308.07915):
+- the codes are built from their polynomials;
+- the depth-8 syndrome cycle is written out gate by gate;
+- the circuit is sampled with the Pauli-frame simulator;
+- its error model is decoded with belief propagation plus ordered statistics.
+
+The source published its simulation software, and running it gives the
+oracle. On the 72-qubit code the derived error model has the source's columns,
+detector for detector, and on identical syndromes the decoder makes the
+source's decision shot for shot. Logical failure rates track the source's;
+the one open residual, 2.3σ on one half at the highest noise, is documented
+in the module. On the 144-qubit code the rates land on the source's published
+curve.
+
 ## Architecture
 
 ```sh

wai/wai-quantum/examples/ldpc_memory.rs

@@ -0,0 +1,45 @@
+//! Run a bivariate-bicycle memory experiment at circuit level and decode it
+//! with BP-OSD: logical failures per shot and per syndrome cycle.
+//!
+//! ```text
+//! cargo run --release --example ldpc_memory --features quantum_ldpc -- CODE P CYCLES SHOTS [THREADS] [MAX_ITER]
+//! ```
+//! CODE is 72, 90, 108, 144 or 288.
+use wai_quantum::quantum_frame::Basis;
+use wai_quantum::quantum_ldpc::{memory_experiment, BbCode, BpOsdConfig};
+
+fn main() {
+    let a: Vec<String> = std::env::args().collect();
+    let code = match a[1].as_str() {
+        "72" => BbCode::bb72(),
+        "90" => BbCode::bb90(),
+        "108" => BbCode::bb108(),
+        "144" => BbCode::bb144(),
+        "288" => BbCode::bb288(),
+        other => panic!("unknown code {other}"),
+    };
+    let p: f64 = a[2].parse().unwrap();
+    let cycles: u32 = a[3].parse().unwrap();
+    let shots: usize = a[4].parse().unwrap();
+    let threads: usize = a.get(5).map_or(8, |x| x.parse().unwrap());
+    let max_iter: u32 = a.get(6).map_or(10_000, |x| x.parse().unwrap());
+    let cfg = BpOsdConfig { max_iter, ..BpOsdConfig::default() };
+    let mut both_ok = 1.0;
+    for (name, basis, seed) in [("Z", Basis::Z, 1u64 << 32), ("X", Basis::X, 2u64 << 32)] {
+        let t0 = std::time::Instant::now();
+        let r = memory_experiment(&code, cycles, p, basis, shots, seed, cfg, threads).unwrap();
+        let dt = t0.elapsed().as_secs_f64();
+        both_ok *= 1.0 - r.per_shot();
+        println!(
+            "{name}-basis  shots {}  failures {}  per shot {:.5}  per cycle {:.3e}  bp-converged {:.3}  {:.1} ms/shot",
+            r.shots,
+            r.failures,
+            r.per_shot(),
+            r.per_cycle(),
+            r.converged as f64 / r.shots as f64,
+            1e3 * dt * threads as f64 / shots as f64
+        );
+    }
+    let p_shot = 1.0 - both_ok;
+    println!("either     per shot {:.5}  per cycle {:.3e}", p_shot, 1.0 - (1.0 - p_shot).powf(1.0 / f64::from(cycles)));
+}

wai/wai-quantum/src/lib.rs

@@ -106,6 +106,9 @@ pub mod quantum_arch;
 #[cfg(feature = "quantum_frame")]
 pub mod quantum_frame;
 
+#[cfg(feature = "quantum_ldpc")]
+pub mod quantum_ldpc;
+
 #[cfg(feature = "quantum_frame")]
 pub mod quantum_match;
 

wai/wai-quantum/src/quantum_ldpc.rs

@@ -0,0 +1,978 @@
+//! qLDPC MEMORY at circuit level — `wai.quantum.ldpc`.
+//!
+//! High-rate quantum LDPC codes promise many logical qubits per physical
+//! qubit. The bivariate-bicycle family (arXiv:2308.07915) is the
+//! best-documented case: 12 logical qubits in 144 data and 144 check qubits at
+//! distance 12. The claim that matters is how often such a memory fails per
+//! syndrome cycle when every gate, idle, preparation and measurement is noisy.
+//! That number comes from a circuit simulation and a decoder.
+//!
+//! This module builds both, so the claim can be re-derived rather than quoted:
+//!
+//! - **The codes.** [`BbCode`] builds the check matrices from the two
+//!   polynomials and counts the logical qubits by rank. It also finds a basis
+//!   of logical operators. Five codes from the source are presets.
+//! - **The circuit.** [`memory_circuit`] writes the source's depth-8 syndrome
+//!   cycle in the ecosystem's circuit text:
+//!   - each check qubit is coupled to its six neighbours in the published
+//!     order;
+//!   - X checks are measured in the X basis and Z checks in the Z basis, all
+//!     interleaved;
+//!   - circuit-level depolarizing noise of strength `p` follows every CNOT,
+//!     idle, preparation and measurement;
+//!   - two noiseless cycles at the end read the final syndrome.
+//!
+//!   `quantum_frame` samples it and derives its detector error model.
+//! - **The decoder.** [`BpOsd`] is belief propagation with ordered-statistics
+//!   post-processing:
+//!   - min-sum belief propagation in the log domain, with the adaptive
+//!     scaling `α = 1 − 2^(−t)`;
+//!   - when it fails to converge, ordered statistics: columns are ranked by
+//!     the posterior, an information set is taken in that order, and the
+//!     combination sweep is scored by `Σ ln(1/p)`.
+//!
+//!   It runs in `f64` with in-crate logarithms, so a decode is the same bits on
+//!   every machine.
+//! - **The experiment.** [`memory_experiment`] runs Z-basis or X-basis memory,
+//!   sampled in fixed chunks across threads, and counts logical failures.
+//!   [`LdpcPoint::per_cycle`] turns them into a rate per syndrome cycle.
+//!
+//! # Checked against the source's own software
+//!
+//! The source's simulation software is Apache-2.0. Running it gives the
+//! oracle:
+//! - **The error model.** Its decoding matrices for `[[72,12,6]]` over six
+//!   cycles at `p = 0.003` have the same columns as the ones derived here
+//!   from the circuit, detector for detector: 2,268 for X faults and 2,232
+//!   for Z faults. Half the probabilities agree to `10⁻¹³`. The largest gap
+//!   is 0.94%, on columns where several faults merge: the source adds their
+//!   probabilities, while this crate combines them as independent events.
+//! - **Logical failure rates.** These agree with the source's decoder on the
+//!   same codes and noise. On `[[72,12,6]]` over six cycles, the failure
+//!   rates per shot are:
+//!
+//!   | `p`   | here (8,000 shots) | source's software (2,000 trials) |
+//!   |-------|--------------------|----------------------------------|
+//!   | 0.003 | 0.028              | 0.034                            |
+//!   | 0.004 | 0.090              | 0.086                            |
+//!   | 0.005 | 0.230              | 0.209                            |
+//!
+//!   "Here" combines the two halves as independent. The source decodes both
+//!   halves of each run, and a Y fault can fail both together: in its software
+//!   at `p = 0.005` both fail in 2.9% of runs, against 1.3% if they were
+//!   independent. That accounts for most of the gap in the last row.
+//!
+//!   Taken half by half at `p = 0.005`, the source's software (9,000 runs for
+//!   the X-type half, 3,000 for the Z-type) gives 0.111 and 0.118. This
+//!   crate gives 0.121 (20,000 shots) and 0.122 (4,000). The Z-type half
+//!   agrees. The X-type half sits 2.3σ above.
+//!
+//!   Each piece of the chain was checked separately:
+//!   - the error models match column for column;
+//!   - the source's decoder, given this crate's syndromes and error model,
+//!     makes the same decision as [`BpOsd`] on every one of 1,500 shots;
+//!   - decoding with the source's priors instead of these moves 476
+//!     failures out of 4,000 to 478;
+//!   - both samplers' detector rates match the exact expectation of the
+//!     error model.
+//!
+//!   The residual lies in sampling, at a level not yet resolved.
+//! - **The 144-qubit code.** On `[[144,12,12]]` over twelve cycles, the
+//!   failure rate per run is 0.034 at `p = 0.004` and 0.180 at `p = 0.005`
+//!   (8,000 and 2,000 shots). The source's software gives 0.039 (1,522 runs)
+//!   and 0.193 (519 runs). Per cycle, the rates land on the source's
+//!   published curve: 2.90·10⁻³ at `p = 0.004` (fit 3.15·10⁻³) and
+//!   1.64·10⁻² at `p = 0.005` (fit 1.54·10⁻²).
+//!
+//! # Honest boundaries
+//!
+//! - **Memory only.** Logical gates on these codes (automorphisms, the
+//!   ZX-duality, measurements through an ancilla system) are not simulated.
+//! - **Two experiments, not one.** A Z-basis memory sees X-type failures and an
+//!   X-basis memory sees Z-type failures. Their rates are combined as
+//!   independent. The source decodes both halves of one noisy run instead,
+//!   which differs at second order.
+//! - **The decoder is the source's.** It is not the best decoder known for
+//!   these codes. A better one would make the code look better. The point here
+//!   is to reproduce the claim on its own terms.
+
+use crate::quantum_frame::{error_model, parse, sample, Basis, ErrorModel, FrameError};
+use crate::repro::ln;
+
+// ---------------------------------------------------------------------------
+// GF(2)
+// ---------------------------------------------------------------------------
+
+/// A set of GF(2) rows over `cols` columns, kept in reduced echelon form as
+/// they arrive, so membership and rank are cheap.
+#[derive(Clone, Debug)]
+struct Span {
+    /// Reduced rows with the column of their leading one.
+    rows: Vec<(usize, Vec<u64>)>,
+}
+
+impl Span {
+    fn new(cols: usize) -> Span {
+        let _ = cols;
+        Span { rows: Vec::new() }
+    }
+    fn reduce(&self, v: &mut [u64]) {
+        for (lead, r) in &self.rows {
+            if (v[lead / 64] >> (lead % 64)) & 1 == 1 {
+                for (a, b) in v.iter_mut().zip(r) {
+                    *a ^= b;
+                }
+            }
+        }
+    }
+    /// Add `v` unless it is already in the span; whether it was new.
+    fn insert(&mut self, v: &[u64]) -> bool {
+        let mut v = v.to_vec();
+        self.reduce(&mut v);
+        let Some(w) = v.iter().position(|&x| x != 0) else { return false };
+        let lead = w * 64 + v[w].trailing_zeros() as usize;
+        // Keep every row reduced at the new leading column.
+        for (_, r) in self.rows.iter_mut() {
+            if (r[lead / 64] >> (lead % 64)) & 1 == 1 {
+                for (a, b) in r.iter_mut().zip(&v) {
+                    *a ^= b;
+                }
+            }
+        }
+        self.rows.push((lead, v));
+        true
+    }
+    fn rank(&self) -> usize {
+        self.rows.len()
+    }
+}
+
+fn bits(cols: usize, ones: impl IntoIterator<Item = usize>) -> Vec<u64> {
+    let mut v = vec![0u64; cols.div_ceil(64)];
+    for c in ones {
+        v[c / 64] ^= 1 << (c % 64);
+    }
+    v
+}
+
+fn ones(v: &[u64]) -> Vec<usize> {
+    let mut out = Vec::new();
+    for (w, &x) in v.iter().enumerate() {
+        let mut x = x;
+        while x != 0 {
+            out.push(w * 64 + x.trailing_zeros() as usize);
+            x &= x - 1;
+        }
+    }
+    out
+}
+
+/// A basis of the null space of `rows` (each a list of columns) over `cols`
+/// columns.
+fn kernel(rows: &[Vec<usize>], cols: usize) -> Vec<Vec<u64>> {
+    // Row-reduce, then read the kernel off the free columns.
+    let mut span = Span::new(cols);
+    for r in rows {
+        span.insert(&bits(cols, r.iter().copied()));
+    }
+    let leads: Vec<usize> = span.rows.iter().map(|(l, _)| *l).collect();
+    let mut is_lead = vec![false; cols];
+    for &l in &leads {
+        is_lead[l] = true;
+    }
+    let mut out = Vec::new();
+    for free in (0..cols).filter(|&c| !is_lead[c]) {
+        let mut v = bits(cols, [free]);
+        for (lead, r) in &span.rows {
+            if (r[free / 64] >> (free % 64)) & 1 == 1 {
+                v[lead / 64] ^= 1 << (lead % 64);
+            }
+        }
+        out.push(v);
+    }
+    out
+}
+
+// ---------------------------------------------------------------------------
+// Bivariate-bicycle codes
+// ---------------------------------------------------------------------------
+
+/// A bivariate-bicycle code: `A = x^a₁ + y^a₂ + y^a₃` and
+/// `B = y^b₁ + x^b₂ + x^b₃` over `x^ℓ = y^m = 1`, with `H_X = [A | B]` and
+/// `H_Z = [Bᵀ | Aᵀ]`.
+#[derive(Clone, Copy, Debug, PartialEq, Eq)]
+pub struct BbCode {
+    pub l: u32,
+    pub m: u32,
+    pub a: [u32; 3],
+    pub b: [u32; 3],
+}
+
+/// Which monomial of `A` or `B` a check's neighbour comes from.
+#[derive(Clone, Copy)]
+enum Mono {
+    X(u32),
+    Y(u32),
+}
+
+impl BbCode {
+    /// `[[72,12,6]]`.
+    pub fn bb72() -> BbCode {
+        BbCode { l: 6, m: 6, a: [3, 1, 2], b: [3, 1, 2] }
+    }
+    /// `[[90,8,10]]`.
+    pub fn bb90() -> BbCode {
+        BbCode { l: 15, m: 3, a: [9, 1, 2], b: [0, 2, 7] }
+    }
+    /// `[[108,8,10]]`.
+    pub fn bb108() -> BbCode {
+        BbCode { l: 9, m: 6, a: [3, 1, 2], b: [3, 1, 2] }
+    }
+    /// `[[144,12,12]]`.
+    pub fn bb144() -> BbCode {
+        BbCode { l: 12, m: 6, a: [3, 1, 2], b: [3, 1, 2] }
+    }
+    /// `[[288,12,18]]`.
+    pub fn bb288() -> BbCode {
+        BbCode { l: 12, m: 12, a: [3, 2, 7], b: [3, 1, 2] }
+    }
+
+    /// Checks of each type, `ℓm`.
+    pub fn half(&self) -> usize {
+        (self.l * self.m) as usize
+    }
+    /// Data qubits, `2ℓm`.
+    pub fn n(&self) -> usize {
+        2 * self.half()
+    }
+
+    fn apply(&self, mono: Mono, i: usize) -> usize {
+        let (l, m) = (self.l as usize, self.m as usize);
+        let (r, c) = (i / m, i % m);
+        match mono {
+            Mono::X(s) => ((r + s as usize) % l) * m + c,
+            Mono::Y(s) => r * m + (c + s as usize) % m,
+        }
+    }
+    fn apply_t(&self, mono: Mono, i: usize) -> usize {
+        let (l, m) = (self.l as usize, self.m as usize);
+        let (r, c) = (i / m, i % m);
+        match mono {
+            Mono::X(s) => ((r + l - s as usize % l) % l) * m + c,
+            Mono::Y(s) => r * m + (c + m - s as usize % m) % m,
+        }
+    }
+    fn a_mono(&self) -> [Mono; 3] {
+        [Mono::X(self.a[0]), Mono::Y(self.a[1]), Mono::Y(self.a[2])]
+    }
+    fn b_mono(&self) -> [Mono; 3] {
+        [Mono::Y(self.b[0]), Mono::X(self.b[1]), Mono::X(self.b[2])]
+    }
+
+    /// The six data qubits of X check `i`, by direction 0 to 5: `A₁(i), A₂(i),
+    /// A₃(i)` on the left half, then `B₁(i), B₂(i), B₃(i)` on the right.
+    pub fn x_neighbours(&self, i: usize) -> [usize; 6] {
+        let (a, b, h) = (self.a_mono(), self.b_mono(), self.half());
+        [self.apply(a[0], i), self.apply(a[1], i), self.apply(a[2], i), h + self.apply(b[0], i), h + self.apply(b[1], i), h + self.apply(b[2], i)]
+    }
+    /// The six data qubits of Z check `i`: `B₁ᵀ(i), B₂ᵀ(i), B₃ᵀ(i)` on the
+    /// left, then `A₁ᵀ(i), A₂ᵀ(i), A₃ᵀ(i)` on the right.
+    pub fn z_neighbours(&self, i: usize) -> [usize; 6] {
+        let (a, b, h) = (self.a_mono(), self.b_mono(), self.half());
+        [self.apply_t(b[0], i), self.apply_t(b[1], i), self.apply_t(b[2], i), h + self.apply_t(a[0], i), h + self.apply_t(a[1], i), h + self.apply_t(a[2], i)]
+    }
+
+    /// X checks as lists of data qubits.
+    pub fn hx(&self) -> Vec<Vec<usize>> {
+        (0..self.half()).map(|i| self.x_neighbours(i).to_vec()).collect()
+    }
+    /// Z checks as lists of data qubits.
+    pub fn hz(&self) -> Vec<Vec<usize>> {
+        (0..self.half()).map(|i| self.z_neighbours(i).to_vec()).collect()
+    }
+
+    /// Logical qubits, `n − rank H_X − rank H_Z`.
+    pub fn k(&self) -> usize {
+        let rank = |rows: Vec<Vec<usize>>| {
+            let mut s = Span::new(self.n());
+            for r in rows {
+                s.insert(&bits(self.n(), r));
+            }
+            s.rank()
+        };
+        self.n() - rank(self.hx()) - rank(self.hz())
+    }
+
+    /// Whether every X check commutes with every Z check.
+    pub fn is_valid(&self) -> bool {
+        let hz = self.hz();
+        self.hx().iter().all(|x| hz.iter().all(|z| x.iter().filter(|q| z.contains(q)).count() % 2 == 0))
+    }
+
+    /// A basis of `k` logical operators of one type. `Basis::Z` gives Z-type
+    /// operators, which commute with the X checks and are not products of Z
+    /// checks; `Basis::X` gives the reverse. Each is a list of data qubits.
+    pub fn logicals(&self, basis: Basis) -> Vec<Vec<usize>> {
+        let (commute, stabilizers) = match basis {
+            Basis::Z => (self.hx(), self.hz()),
+            Basis::X => (self.hz(), self.hx()),
+        };
+        let n = self.n();
+        let mut span = Span::new(n);
+        for s in stabilizers {
+            span.insert(&bits(n, s));
+        }
+        let mut out = Vec::new();
+        for v in kernel(&commute, n) {
+            if span.insert(&v) {
+                out.push(ones(&v));
+            }
+        }
+        out
+    }
+}
+
+// ---------------------------------------------------------------------------
+// The syndrome-measurement circuit
+// ---------------------------------------------------------------------------
+
+/// The order in which X checks (`SX`) and Z checks (`SZ`) visit their six
+/// neighbours over the seven CNOT rounds of a cycle (arXiv:2308.07915, its
+/// table of the depth-8 cycle). `None` is a round without a CNOT.
+const SX: [Option<usize>; 7] = [None, Some(1), Some(4), Some(3), Some(5), Some(0), Some(2)];
+const SZ: [Option<usize>; 7] = [Some(3), Some(5), Some(0), Some(1), Some(2), Some(4), None];
+
+fn line(out: &mut String, name: &str, arg: Option<f64>, targets: &[usize]) {
+    if targets.is_empty() {
+        return;
+    }
+    out.push_str(name);
+    if let Some(a) = arg {
+        out.push_str(&format!("({a})"));
+    }
+    for t in targets {
+        out.push_str(&format!(" {t}"));
+    }
+    out.push('\n');
+}
+
+/// Text of a bivariate-bicycle memory experiment in `basis`.
+///
+/// - **Layout.** Qubits are X checks `0..ℓm`, data `ℓm..3ℓm` (left half, then
+///   right) and Z checks `3ℓm..4ℓm`.
+/// - **Start.** Data and Z checks start in `basis`, noiselessly.
+/// - **Cycles.** `cycles` noisy syndrome cycles at strength `p`, then two
+///   noiseless ones.
+/// - **Detectors.** Each check of `basis`, cycle by cycle: its first outcome,
+///   then each change.
+/// - **Observables.** The logical operators of [`BbCode::logicals`], read from
+///   a final noiseless measurement of the data.
+///
+/// Noise:
+/// - two-qubit depolarizing noise after each CNOT;
+/// - one-qubit depolarizing noise on idle data;
+/// - a flip after each preparation;
+/// - a flipped result on each measurement.
+pub fn memory_circuit(code: &BbCode, cycles: u32, p: f64, basis: Basis) -> String {
+    let h = code.half();
+    let xc: Vec<usize> = (0..h).collect();
+    let data: Vec<usize> = (h..3 * h).collect();
+    let zc: Vec<usize> = (3 * h..4 * h).collect();
+    let xn: Vec<[usize; 6]> = (0..h).map(|i| code.x_neighbours(i).map(|q| q + h)).collect();
+    let zn: Vec<[usize; 6]> = (0..h).map(|i| code.z_neighbours(i).map(|q| q + h)).collect();
+    let mut out = String::new();
+    match basis {
+        Basis::Z => line(&mut out, "R", None, &data),
+        Basis::X => line(&mut out, "RX", None, &data),
+    }
+    line(&mut out, "R", None, &zc);
+    let mut measured = 0usize;
+    // Absolute record of each check's outcome in the previous cycle.
+    let mut previous: Option<Vec<usize>> = None;
+    for cycle in 0..cycles + 2 {
+        let noise = if cycle < cycles { Some(p) } else { None };
+        let idle_cnots = |out: &mut String, t: usize, x_rounds: bool, z_rounds: bool| {
+            let mut pairs: Vec<usize> = Vec::new();
+            let mut busy = vec![false; 4 * h];
+            if let (true, Some(dir)) = (z_rounds, SZ[t]) {
+                for i in 0..h {
+                    pairs.extend([zn[i][dir], zc[i]]);
+                    busy[zn[i][dir]] = true;
+                }
+            }
+            if let (true, Some(dir)) = (x_rounds, SX[t]) {
+                for i in 0..h {
+                    pairs.extend([xc[i], xn[i][dir]]);
+                    busy[xn[i][dir]] = true;
+                }
+            }
+            line(out, "CX", None, &pairs);
+            if let Some(p) = noise {
+                line(out, "DEPOLARIZE2", Some(p), &pairs);
+                let idle: Vec<usize> = data.iter().copied().filter(|&q| !busy[q]).collect();
+                line(out, "DEPOLARIZE1", Some(p), &idle);
+            }
+        };
+        // Round 1: prepare the X checks; the Z checks take their first CNOT.
+        line(&mut out, "RX", None, &xc);
+        if let Some(p) = noise {
+            line(&mut out, "Z_ERROR", Some(p), &xc);
+        }
+        idle_cnots(&mut out, 0, false, true);
+        out.push_str("TICK\n");
+        // Rounds 2 to 6: both kinds of check.
+        for t in 1..6 {
+            idle_cnots(&mut out, t, true, true);
+            out.push_str("TICK\n");
+        }
+        // Round 7: measure the Z checks; the X checks take their last CNOT.
+        line(&mut out, "M", noise, &zc);
+        let z_records: Vec<usize> = (measured..measured + h).collect();
+        measured += h;
+        idle_cnots(&mut out, 6, true, false);
+        out.push_str("TICK\n");
+        // Round 8: data idle; measure the X checks; prepare the Z checks.
+        if let Some(p) = noise {
+            line(&mut out, "DEPOLARIZE1", Some(p), &data);
+        }
+        line(&mut out, "MX", noise, &xc);
+        let x_records: Vec<usize> = (measured..measured + h).collect();
+        measured += h;
+        line(&mut out, "R", None, &zc);
+        if let Some(p) = noise {
+            line(&mut out, "X_ERROR", Some(p), &zc);
+        }
+        out.push_str("TICK\n");
+        let current = match basis {
+            Basis::Z => z_records,
+            Basis::X => x_records,
+        };
+        for (i, &r) in current.iter().enumerate() {
+            match &previous {
+                None => out.push_str(&format!("DETECTOR rec[-{}]\n", measured - r)),
+                Some(prev) => out.push_str(&format!("DETECTOR rec[-{}] rec[-{}]\n", measured - r, measured - prev[i])),
+            }
+        }
+        previous = Some(current);
+    }
+    match basis {
+        Basis::Z => line(&mut out, "M", None, &data),
+        Basis::X => line(&mut out, "MX", None, &data),
+    }
+    let first = measured;
+    measured += data.len();
+    for (k, op) in code.logicals(basis).iter().enumerate() {
+        out.push_str(&format!("OBSERVABLE_INCLUDE({k})"));
+        for &q in op {
+            out.push_str(&format!(" rec[-{}]", measured - (first + q)));
+        }
+        out.push('\n');
+    }
+    out
+}
+
+// ---------------------------------------------------------------------------
+// BP-OSD
+// ---------------------------------------------------------------------------
+
+/// How belief propagation scales its min-sum messages.
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub enum Scaling {
+    /// `α = 1 − 2^(−t)` at iteration `t`.
+    Adaptive,
+    Fixed(f64),
+}
+
+/// Decoder settings. The default is the source's: min-sum, adaptive scaling,
+/// up to 10,000 iterations, combination sweep of order 7.
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub struct BpOsdConfig {
+    pub max_iter: u32,
+    pub scaling: Scaling,
+    /// Order of the combination sweep; 0 is OSD-0 alone.
+    pub osd_order: u32,
+}
+
+impl Default for BpOsdConfig {
+    fn default() -> BpOsdConfig {
+        BpOsdConfig { max_iter: 10_000, scaling: Scaling::Adaptive, osd_order: 7 }
+    }
+}
+
+/// A decode: the observables the correction flips, and how it was found.
+#[derive(Clone, Copy, Debug, PartialEq, Eq)]
+pub struct Decoded {
+    pub observables: u64,
+    /// Whether belief propagation alone explained the syndrome.
+    pub converged: bool,
+    pub iterations: u32,
+}
+
+/// Belief propagation with ordered-statistics post-processing over a
+/// detector error model.
+#[derive(Clone, Debug)]
+pub struct BpOsd {
+    cfg: BpOsdConfig,
+    checks: usize,
+    vars: usize,
+    /// Edges by check: `(var)`, CSR.
+    row_start: Vec<usize>,
+    row_var: Vec<usize>,
+    /// Edges by variable: indices into the check-ordered edge list, CSR.
+    col_start: Vec<usize>,
+    col_edge: Vec<usize>,
+    prior: Vec<f64>,
+    /// `ln(1/p)`, the OSD cost of each fault.
+    cost: Vec<f64>,
+    observables: Vec<u64>,
+    rank: usize,
+}
+
+impl BpOsd {
+    /// Build the decoder for `model`. Faults with probability 0 are dropped.
+    pub fn new(model: &ErrorModel, cfg: BpOsdConfig) -> BpOsd {
+        let mechs: Vec<_> = model.mechanisms.iter().filter(|m| m.probability > 0.0 && !m.detectors.is_empty()).collect();
+        let checks = model.detectors as usize;
+        let vars = mechs.len();
+        let mut by_row: Vec<Vec<usize>> = vec![Vec::new(); checks];
+        for (j, m) in mechs.iter().enumerate() {
+            for &d in &m.detectors {
+                by_row[d as usize].push(j);
+            }
+        }
+        let mut row_start = vec![0];
+        let mut row_var: Vec<usize> = Vec::new();
+        for r in &by_row {
+            row_var.extend(r);
+            row_start.push(row_var.len());
+        }
+        let mut by_col: Vec<Vec<usize>> = vec![Vec::new(); vars];
+        for (e, &j) in row_var.iter().enumerate() {
+            by_col[j].push(e);
+        }
+        let mut col_start = vec![0];
+        let mut col_edge: Vec<usize> = Vec::new();
+        for c in &by_col {
+            col_edge.extend(c);
+            col_start.push(col_edge.len());
+        }
+        let prior: Vec<f64> = mechs.iter().map(|m| ln((1.0 - m.probability) / m.probability)).collect();
+        let cost: Vec<f64> = mechs.iter().map(|m| ln(1.0 / m.probability)).collect();
+        let observables = mechs.iter().map(|m| m.observables).collect();
+        let mut span = Span::new(vars);
+        // Rank of H: row-reduce the checks over the faults.
+        for r in &by_row {
+            span.insert(&bits(vars, r.iter().copied()));
+        }
+        let rank = span.rank();
+        BpOsd { cfg, checks, vars, row_start, row_var, col_start, col_edge, prior, cost, observables, rank }
+    }
+
+    /// Faults (columns) in the model the decoder uses.
+    pub fn faults(&self) -> usize {
+        self.vars
+    }
+
+    /// Decode the detectors that fired.
+    pub fn decode(&self, fired: &[u32]) -> Decoded {
+        let mut synd = vec![false; self.checks];
+        for &d in fired {
+            synd[d as usize] = true;
+        }
+        let edges = self.row_var.len();
+        let mut v2c = vec![0.0f64; edges];
+        let mut c2v = vec![0.0f64; edges];
+        let mut sgn = vec![0u32; edges];
+        let mut post = vec![0.0f64; self.vars];
+        let mut hard = vec![false; self.vars];
+        for j in 0..self.vars {
+            for &e in &self.col_edge[self.col_start[j]..self.col_start[j + 1]] {
+                v2c[e] = self.prior[j];
+            }
+        }
+        let mut iterations = 0;
+        for t in 1..=self.cfg.max_iter {
+            iterations = t;
+            let alpha = match self.cfg.scaling {
+                Scaling::Adaptive => 1.0 - pow2_neg(t),
+                Scaling::Fixed(a) => a,
+            };
+            // Check to variable: the smallest other magnitude, the parity of
+            // the other signs and the syndrome; zero counts as negative.
+            for (i, &fired) in synd.iter().enumerate() {
+                let (s, e_end) = (self.row_start[i], self.row_start[i + 1]);
+                let mut min = 1e308f64;
+                let mut negatives = u32::from(fired);
+                for e in s..e_end {
+                    c2v[e] = min;
+                    sgn[e] = negatives;
+                    if v2c[e].abs() < min {
+                        min = v2c[e].abs();
+                    }
+                    if v2c[e] <= 0.0 {
+                        negatives += 1;
+                    }
+                }
+                let (mut min, mut negatives) = (1e308f64, 0u32);
+                for e in (s..e_end).rev() {
+                    if min < c2v[e] {
+                        c2v[e] = min;
+                    }
+                    sgn[e] += negatives;
+                    c2v[e] *= if sgn[e] % 2 == 1 { -alpha } else { alpha };
+                    if v2c[e].abs() < min {
+                        min = v2c[e].abs();
+                    }
+                    if v2c[e] <= 0.0 {
+                        negatives += 1;
+                    }
+                }
+            }
+            // Variable to check: the prior plus every other incoming message,
+            // summed forward then backward.
+            for j in 0..self.vars {
+                let col = &self.col_edge[self.col_start[j]..self.col_start[j + 1]];
+                let mut acc = self.prior[j];
+                for &e in col {
+                    v2c[e] = acc;
+                    acc += c2v[e];
+                }
+                post[j] = acc;
+                hard[j] = acc <= 0.0;
+                let mut acc = 0.0;
+                for &e in col.iter().rev() {
+                    v2c[e] += acc;
+                    acc += c2v[e];
+                }
+            }
+            if self.explains(&hard, &synd) {
+                return Decoded { observables: self.flips(&hard), converged: true, iterations };
+            }
+        }
+        let x = self.osd(&post, &synd);
+        Decoded { observables: self.flips(&x), converged: false, iterations }
+    }
+
+    fn explains(&self, x: &[bool], synd: &[bool]) -> bool {
+        (0..self.checks).all(|i| {
+            let parity = self.row_var[self.row_start[i]..self.row_start[i + 1]].iter().filter(|&&j| x[j]).count() % 2 == 1;
+            parity == synd[i]
+        })
+    }
+
+    fn flips(&self, x: &[bool]) -> u64 {
+        x.iter().zip(&self.observables).filter(|(b, _)| **b).fold(0, |m, (_, o)| m ^ o)
+    }
+
+    /// Ordered-statistics decoding: rank faults by posterior (likeliest
+    /// first), take the first independent columns as the information set, and
+    /// solve; then sweep single flips of every other column and pairs among
+    /// the first `osd_order`, keeping the cheapest solution.
+    fn osd(&self, post: &[f64], synd: &[bool]) -> Vec<bool> {
+        let mut order: Vec<usize> = (0..self.vars).collect();
+        order.sort_by(|&a, &b| post[a].total_cmp(&post[b]).then(a.cmp(&b)));
+        // Row-major matrix over the sorted columns, the syndrome appended.
+        let cols = self.vars + 1;
+        let words = cols.div_ceil(64);
+        let mut m = vec![0u64; self.checks * words];
+        let mut position = vec![0usize; self.vars];
+        for (p, &j) in order.iter().enumerate() {
+            position[j] = p;
+        }
+        for i in 0..self.checks {
+            for &j in &self.row_var[self.row_start[i]..self.row_start[i + 1]] {
+                let p = position[j];
+                m[i * words + p / 64] ^= 1 << (p % 64);
+            }
+            if synd[i] {
+                let p = self.vars;
+                m[i * words + p / 64] ^= 1 << (p % 64);
+            }
+        }
+        // Reduced echelon form, pivots in sorted order.
+        let mut pivots: Vec<usize> = Vec::with_capacity(self.rank);
+        let mut row = 0;
+        for p in 0..self.vars {
+            if row == self.checks || pivots.len() == self.rank {
+                break;
+            }
+            let (w, b) = (p / 64, p % 64);
+            let Some(r) = (row..self.checks).find(|&r| (m[r * words + w] >> b) & 1 == 1) else { continue };
+            if r != row {
+                for k in 0..words {
+                    m.swap(r * words + k, row * words + k);
+                }
+            }
+            let (head, tail) = m.split_at_mut(row * words);
+            let (pivot, rest) = tail.split_at_mut(words);
+            for other in head.chunks_exact_mut(words).chain(rest.chunks_exact_mut(words)) {
+                // Whole rows: the free columns before `p` are read later.
+                if (other[w] >> b) & 1 == 1 {
+                    for (a, c) in other.iter_mut().zip(pivot.iter()) {
+                        *a ^= c;
+                    }
+                }
+            }
+            pivots.push(p);
+            row += 1;
+        }
+        let bit = |r: usize, p: usize| (m[r * words + p / 64] >> (p % 64)) & 1 == 1;
+        let s = self.vars;
+        // OSD-0: the pivots take the reduced syndrome; everything else is 0.
+        let mut base = vec![false; self.vars];
+        for (r, &p) in pivots.iter().enumerate() {
+            base[p] = bit(r, s);
+        }
+        let cost_at = |p: usize| self.cost[order[p]];
+        let base_cost: f64 = (0..self.vars).filter(|&p| base[p]).map(cost_at).sum();
+        let mut best = (base_cost, Vec::new());
+        if self.cfg.osd_order > 0 {
+            let mut is_pivot = vec![false; self.vars];
+            for &p in &pivots {
+                is_pivot[p] = true;
+            }
+            let free: Vec<usize> = (0..self.vars).filter(|&p| !is_pivot[p]).collect();
+            // The change in cost from flipping a set of free columns.
+            let delta = |set: &[usize]| {
+                let mut d: f64 = set.iter().map(|&p| cost_at(p)).sum();
+                for (r, &pv) in pivots.iter().enumerate() {
+                    let flip = set.iter().fold(false, |a, &p| a ^ bit(r, p));
+                    if flip {
+                        d += if base[pv] { -cost_at(pv) } else { cost_at(pv) };
+                    }
+                }
+                d
+            };
+            let order_k = (self.cfg.osd_order as usize).min(free.len());
+            let mut candidates: Vec<Vec<usize>> = free.iter().map(|&p| vec![p]).collect();
+            for i in 0..order_k {
+                for j in i + 1..order_k {
+                    candidates.push(vec![free[i], free[j]]);
+                }
+            }
+            for set in candidates {
+                let c = base_cost + delta(&set);
+                if c < best.0 {
+                    best = (c, set);
+                }
+            }
+        }
+        let mut x = base;
+        for &p in &best.1 {
+            x[p] = true;
+            for (r, &pv) in pivots.iter().enumerate() {
+                if bit(r, p) {
+                    x[pv] = !x[pv];
+                }
+            }
+        }
+        let mut out = vec![false; self.vars];
+        for (p, &j) in order.iter().enumerate() {
+            out[j] = x[p];
+        }
+        out
+    }
+}
+
+fn pow2_neg(t: u32) -> f64 {
+    if t >= 1075 {
+        return 0.0;
+    }
+    let mut x = 1.0;
+    for _ in 0..t {
+        x /= 2.0;
+    }
+    x
+}
+
+// ---------------------------------------------------------------------------
+// Memory experiments
+// ---------------------------------------------------------------------------
+
+/// A memory experiment's tally.
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub struct LdpcPoint {
+    pub p: f64,
+    pub cycles: u32,
+    pub shots: u64,
+    /// Shots whose decoded observables were wrong (any logical flipped).
+    pub failures: u64,
+    /// Shots belief propagation alone decoded.
+    pub converged: u64,
+}
+
+impl LdpcPoint {
+    /// Failure rate per shot.
+    pub fn per_shot(&self) -> f64 {
+        self.failures as f64 / self.shots as f64
+    }
+    /// Failure rate per syndrome cycle, `1 − (1 − P)^(1/N)`.
+    pub fn per_cycle(&self) -> f64 {
+        1.0 - crate::repro::exp(ln(1.0 - self.per_shot()) / f64::from(self.cycles))
+    }
+}
+
+/// Shots per independently seeded chunk: results do not depend on how many
+/// threads ran them.
+const CHUNK: usize = 64;
+
+/// Run `shots` shots of the memory experiment in `basis` and decode each.
+/// Chunk `c` is sampled with seed `seed + c`. The tally is the same at any
+/// `threads` (one on wasm32, which has none).
+#[allow(clippy::too_many_arguments)]
+pub fn memory_experiment(code: &BbCode, cycles: u32, p: f64, basis: Basis, shots: usize, seed: u64, cfg: BpOsdConfig, threads: usize) -> Result<LdpcPoint, FrameError> {
+    let circuit = parse(&memory_circuit(code, cycles, p, basis))?;
+    let model = error_model(&circuit)?;
+    let decoder = BpOsd::new(&model, cfg);
+    let chunks = shots.div_ceil(CHUNK);
+    let run = |c: usize| {
+        let n = CHUNK.min(shots - c * CHUNK);
+        let det = sample(&circuit, n, seed.wrapping_add(c as u64));
+        let (mut fail, mut conv) = (0u64, 0u64);
+        for s in 0..n {
+            let d = decoder.decode(&det.fired(s));
+            fail += u64::from(d.observables != det.flips(s));
+            conv += u64::from(d.converged);
+        }
+        (fail, conv)
+    };
+    // wasm32 has no threads; the tally is the same either way.
+    let threads = if cfg!(target_arch = "wasm32") { 1 } else { threads.max(1).min(chunks.max(1)) };
+    let mut tallies = vec![(0u64, 0u64); chunks];
+    if threads == 1 {
+        for (c, t) in tallies.iter_mut().enumerate() {
+            *t = run(c);
+        }
+    } else {
+        let next = std::sync::atomic::AtomicUsize::new(0);
+        let results = std::sync::Mutex::new(&mut tallies);
+        std::thread::scope(|scope| {
+            for _ in 0..threads {
+                scope.spawn(|| loop {
+                    let c = next.fetch_add(1, std::sync::atomic::Ordering::Relaxed);
+                    if c >= chunks {
+                        break;
+                    }
+                    let r = run(c);
+                    results.lock().unwrap()[c] = r;
+                });
+            }
+        });
+    }
+    let (failures, converged) = tallies.iter().fold((0, 0), |(f, c), t| (f + t.0, c + t.1));
+    Ok(LdpcPoint { p, cycles, shots: shots as u64, failures, converged })
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn the_presets_have_their_published_parameters() {
+        for (code, n, k) in [(BbCode::bb72(), 72, 12), (BbCode::bb90(), 90, 8), (BbCode::bb108(), 108, 8), (BbCode::bb144(), 144, 12), (BbCode::bb288(), 288, 12)] {
+            assert!(code.is_valid());
+            assert_eq!((code.n(), code.k()), (n, k));
+            for basis in [Basis::X, Basis::Z] {
+                assert_eq!(code.logicals(basis).len(), k);
+            }
+        }
+    }
+
+    #[test]
+    fn logical_operators_commute_with_checks_and_pair_up() {
+        let code = BbCode::bb72();
+        let (lx, lz) = (code.logicals(Basis::X), code.logicals(Basis::Z));
+        let overlap = |a: &[usize], b: &[usize]| a.iter().filter(|q| b.contains(q)).count() % 2;
+        for z in &lz {
+            assert!(code.hx().iter().all(|c| overlap(c, z) == 0));
+        }
+        for x in &lx {
+            assert!(code.hz().iter().all(|c| overlap(c, x) == 0));
+        }
+        // The pairing matrix between the two bases is invertible.
+        let mut span = Span::new(lz.len());
+        for x in &lx {
+            span.insert(&bits(lz.len(), lz.iter().enumerate().filter(|(_, z)| overlap(x, z) == 1).map(|(i, _)| i)));
+        }
+        assert_eq!(span.rank(), lz.len());
+    }
+
+    #[test]
+    fn the_noiseless_cycle_is_deterministic_and_silent() {
+        for basis in [Basis::Z, Basis::X] {
+            let c = parse(&memory_circuit(&BbCode::bb72(), 3, 0.0, basis)).unwrap();
+            let model = error_model(&c).unwrap();
+            assert!(model.mechanisms.iter().all(|m| m.probability == 0.0) || model.mechanisms.is_empty());
+            let det = sample(&c, 64, 1);
+            assert!((0..64).all(|s| det.fired(s).is_empty() && det.flips(s) == 0));
+        }
+    }
+
+    #[test]
+    fn the_error_model_has_the_reference_columns() {
+        // The source's decoding matrices for [[72,12,6]], six cycles, p = 0.003.
+        let reference = include_str!("../tests/data/bb72_p003_c6.dem");
+        type Column = (f64, Vec<u32>);
+        let mut sections: Vec<(char, Vec<Column>)> = Vec::new();
+        for l in reference.lines() {
+            if let Some(h) = l.strip_prefix("# ") {
+                sections.push((h.chars().next().unwrap(), Vec::new()));
+                continue;
+            }
+            let mut f = l.split_whitespace();
+            let p: f64 = f.next().unwrap().parse().unwrap();
+            let rows: Vec<u32> = f.map(|x| x.parse().unwrap()).collect();
+            // Its matrix drops the logical rows, so columns that differ only
+            // there repeat; their faults add.
+            let cols = &mut sections.last_mut().unwrap().1;
+            match cols.iter_mut().find(|(_, r)| *r == rows) {
+                Some((q, _)) => *q += p,
+                None => cols.push((p, rows)),
+            }
+        }
+        for (kind, columns) in sections {
+            // X faults are seen by Z checks: the Z-basis memory.
+            let basis = if kind == 'X' { Basis::Z } else { Basis::X };
+            let c = parse(&memory_circuit(&BbCode::bb72(), 6, 0.003, basis)).unwrap();
+            let model = error_model(&c).unwrap();
+            let mut ours: Vec<(Vec<u32>, f64)> = Vec::new();
+            for m in model.mechanisms.iter().filter(|m| !m.detectors.is_empty()) {
+                match ours.iter_mut().find(|(d, _)| *d == m.detectors) {
+                    // The same detectors with different observables: the
+                    // source keeps logical rows out of its decoding matrix.
+                    Some((_, q)) => *q = *q + m.probability - 2.0 * *q * m.probability,
+                    None => ours.push((m.detectors.clone(), m.probability)),
+                }
+            }
+            // The source keeps one column of faults that light nothing (its
+            // probabilities summed); here they leave the model. None of them
+            // flips a logical.
+            let columns: Vec<_> = columns.into_iter().filter(|(_, rows)| !rows.is_empty()).collect();
+            assert!(model.mechanisms.iter().all(|m| !m.detectors.is_empty() || m.observables == 0));
+            assert_eq!(ours.len(), columns.len(), "{kind}");
+            for (p, rows) in &columns {
+                let (_, q) = ours.iter().find(|(d, _)| d == rows).unwrap_or_else(|| panic!("{kind}: no column {rows:?}"));
+                assert!((q / p - 1.0).abs() < 0.01, "{kind} {rows:?}: {q} vs {p}");
+            }
+        }
+    }
+
+    #[test]
+    fn bp_osd_corrects_every_single_fault() {
+        let c = parse(&memory_circuit(&BbCode::bb72(), 2, 0.003, Basis::Z)).unwrap();
+        let model = error_model(&c).unwrap();
+        let dec = BpOsd::new(&model, BpOsdConfig::default());
+        for m in model.mechanisms.iter().filter(|m| !m.detectors.is_empty()) {
+            assert_eq!(dec.decode(&m.detectors).observables, m.observables, "{:?}", m.detectors);
+        }
+    }
+
+    #[test]
+    fn the_tally_does_not_depend_on_threads() {
+        let cfg = BpOsdConfig { max_iter: 50, ..BpOsdConfig::default() };
+        let a = memory_experiment(&BbCode::bb72(), 2, 0.004, Basis::Z, 600, 9, cfg, 1).unwrap();
+        let b = memory_experiment(&BbCode::bb72(), 2, 0.004, Basis::Z, 600, 9, cfg, 3).unwrap();
+        assert_eq!(a, b);
+        assert!(a.failures > 0 && a.converged > 0);
+    }
+}

wai/wai-quantum/tests/bit_identity.rs

@@ -20,7 +20,8 @@
     feature = "quantum_tn",
     feature = "quantum_gbs",
     feature = "quantum_sv",
-    feature = "quantum_arch"
+    feature = "quantum_arch",
+    feature = "quantum_ldpc"
 ))]
 
 use wai_quantum::*;
@@ -339,3 +340,27 @@ fn architecture_estimates() {
     }
     assert_eq!(digest(&v), 0xd166_b4c6_c731_7a52);
 }
+
+/// A bivariate-bicycle memory: the circuit's error model (probabilities from
+/// the independent-channel conversion), BP over it with in-crate logarithms,
+/// and OSD on the shots BP leaves unexplained.
+#[test]
+fn qldpc_memory_decoding() {
+    use quantum_frame::{error_model, parse, sample, Basis};
+    use quantum_ldpc::*;
+    let c = parse(&memory_circuit(&BbCode::bb72(), 3, 0.005, Basis::Z)).unwrap();
+    let model = error_model(&c).unwrap();
+    let dec = BpOsd::new(&model, BpOsdConfig { max_iter: 40, ..BpOsdConfig::default() });
+    let det = sample(&c, 192, 11);
+    let mut v: Vec<f64> = model.mechanisms.iter().step_by(37).map(|m| m.probability).collect();
+    let mut osd = 0;
+    for s in 0..192 {
+        let d = dec.decode(&det.fired(s));
+        osd += u32::from(!d.converged);
+        v.push(d.observables as f64);
+        v.push(f64::from(d.iterations));
+        v.push(det.flips(s) as f64);
+    }
+    assert!(osd > 0, "the sample must exercise OSD");
+    assert_eq!(digest(&v), 0x5479_ebe5_f648_b739);
+}

wai/wai-rs/Cargo.lock

@@ -3900,7 +3900,7 @@ dependencies = [
 
 [[package]]
 name = "wai-quantum"
-version = "0.3.29"
+version = "0.3.30"
 dependencies = [
  "blake3",
  "ed25519-dalek",

wai/wai-rs/Cargo.toml

@@ -360,6 +360,8 @@ quantum_sv = ["dep:wai-quantum", "wai-quantum/quantum_sv"]
 # Architecture-level estimation: yoked storage, cultivation, CCZ factories,
 # reaction-limited time, whole machines.
 quantum_arch = ["dep:wai-quantum", "wai-quantum/quantum_arch"]
+# qLDPC memory at circuit level: bivariate-bicycle codes, BP-OSD.
+quantum_ldpc = ["quantum_frame", "wai-quantum/quantum_ldpc"]
 quantum_qir = ["quantum_toolchain", "wai-quantum/quantum_qir"]
 # Variational quantum machine learning (extensions/quantum-sim § VML,
 # `wai.quantum.vml`): the interference-based ML the thermodynamic-computing pivot
@@ -493,6 +495,7 @@ required-features = [
     "quantum_gbs",
     "quantum_sv",
     "quantum_arch",
+    "quantum_ldpc",
 ]
 
 [[bin]]

wai/wai-rs/src/bin/wai_quantum.rs

@@ -165,6 +165,10 @@ fn usage() -> ! {
          wai-quantum gbs --squeezing r1,r2,… [--eta E] [--seed S | --unitary U.json] [--json]\n  \
          wai-quantum gbs … --sample K                 (exact click patterns)\n  \
          wai-quantum gbs … --check SAMPLES.txt        (observed patterns vs exact marginals)\n\n\
+         QLDPC MEMORY (circuit level, BP-OSD):\n  \
+         wai-quantum ldpc [--code 72|90|108|144|288] [--p P] [--cycles N] [--shots S] [--basis z|x|both]\n    \
+         [--max-iter I] [--osd-order K] [--threads T] [--seed S] [--json]\n  \
+         wai-quantum ldpc --circuit [--code …] [--p P] [--cycles N] [--basis z|x]   (the circuit text)\n\n\
          ARCHITECTURE:\n  \
          wai-quantum arch factoring [--bits 1024|1536|2048|3072|4096|6144|8192] [--json]\n  \
          wai-quantum arch storage [--target T] [--fits refit|uniform|si1000] [--json]\n  \
@@ -961,6 +965,66 @@ fn main() {
             }
         }
 
+        // qLDPC memory: a bivariate-bicycle code's depth-8 syndrome cycle,
+        // sampled with circuit-level noise and decoded by BP-OSD.
+        "ldpc" => {
+            use wai::quantum_frame::Basis;
+            use wai::quantum_ldpc::{memory_circuit, memory_experiment, BbCode, BpOsdConfig};
+            let code = match opt(args, "--code", "").as_deref().unwrap_or("72") {
+                "72" => BbCode::bb72(),
+                "90" => BbCode::bb90(),
+                "108" => BbCode::bb108(),
+                "144" => BbCode::bb144(),
+                "288" => BbCode::bb288(),
+                _ => die("--code: 72, 90, 108, 144 or 288"),
+            };
+            let p: f64 = num(args, "--p", "", 0.003);
+            let cycles: u32 = num(args, "--cycles", "", 6);
+            let bases: Vec<(&str, Basis)> = match opt(args, "--basis", "").as_deref().unwrap_or("both") {
+                "z" | "Z" => vec![("z", Basis::Z)],
+                "x" | "X" => vec![("x", Basis::X)],
+                "both" => vec![("z", Basis::Z), ("x", Basis::X)],
+                _ => die("--basis: z, x or both"),
+            };
+            if flag(args, "--circuit") {
+                print!("{}", memory_circuit(&code, cycles, p, bases[0].1));
+                return;
+            }
+            let shots: usize = num(args, "--shots", "", 1000);
+            let seed: u64 = num(args, "--seed", "", 0);
+            let threads: usize = num(args, "--threads", "", std::thread::available_parallelism().map_or(1, |n| n.get()));
+            let cfg = BpOsdConfig { max_iter: num(args, "--max-iter", "", 10_000), osd_order: num(args, "--osd-order", "", 7), ..BpOsdConfig::default() };
+            let mut rows = Vec::new();
+            for (i, (name, basis)) in bases.iter().enumerate() {
+                let r = memory_experiment(&code, cycles, p, *basis, shots, seed ^ ((i as u64 + 1) << 32), cfg, threads).unwrap_or_else(|e| die(&e.to_string()));
+                rows.push((*name, r));
+            }
+            let either = 1.0 - rows.iter().map(|(_, r)| 1.0 - r.per_shot()).product::<f64>();
+            let either_cycle = 1.0 - (1.0 - either).powf(1.0 / f64::from(cycles));
+            if flag(args, "--json") {
+                let items: Vec<String> = rows
+                    .iter()
+                    .map(|(b, r)| format!("{{\"basis\":\"{b}\",\"bp_converged\":{},\"failures\":{},\"per_cycle\":{},\"per_shot\":{},\"shots\":{}}}", r.converged, r.failures, r.per_cycle(), r.per_shot(), r.shots))
+                    .collect();
+                println!("{{\"code\":[{},{}],\"cycles\":{cycles},\"either_per_cycle\":{either_cycle},\"either_per_shot\":{either},\"p\":{p},\"runs\":[{}]}}", code.n(), code.k(), items.join(","));
+            } else {
+                println!("[[{},{}]] bivariate-bicycle memory, p = {p}, {cycles} noisy cycles", code.n(), code.k());
+                for (b, r) in &rows {
+                    println!(
+                        "{b}-basis   {} failures in {} shots   {:.3e} per shot   {:.3e} per cycle   BP alone {:.1}%",
+                        r.failures,
+                        r.shots,
+                        r.per_shot(),
+                        r.per_cycle(),
+                        100.0 * r.converged as f64 / r.shots as f64
+                    );
+                }
+                if rows.len() == 2 {
+                    println!("either    {either:.3e} per shot   {either_cycle:.3e} per cycle (the two halves as independent)");
+                }
+            }
+        }
+
         // The architecture around the code: yoked storage, cultivated magic
         // states, CCZ factories and the clock, priced; and the published
         // factoring machine re-run region by region.

wai/wai-rs/src/lib.rs

@@ -637,6 +637,9 @@ pub use wai_quantum::quantum_sv;
 #[cfg(feature = "quantum_arch")]
 pub use wai_quantum::quantum_arch;
 
+#[cfg(feature = "quantum_ldpc")]
+pub use wai_quantum::quantum_ldpc;
+
 #[cfg(any(feature = "quantum_tn", feature = "quantum_sv"))]
 pub use wai_quantum::gates;
 #[cfg(feature = "quantum_qir")]

wai/wai-web/Cargo.lock

@@ -734,7 +734,7 @@ dependencies = [
 
 [[package]]
 name = "wai-quantum"
-version = "0.3.29"
+version = "0.3.30"
 dependencies = [
  "blake3",
  "ed25519-dalek",