wai/wai-quantum-serve/Cargo.lock

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

wai/wai-quantum/Cargo.lock

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

wai/wai-quantum/Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "wai-quantum"
-version = "0.3.32"
+version = "0.3.33"
 edition = "2024"
 rust-version = "1.97"
 license = "Apache-2.0"
@@ -76,6 +76,9 @@ quantum_ldpc = ["quantum_frame"]
 quantum_qio = []
 # Union-find decoding in near-linear time over the matching graph.
 quantum_uf = ["quantum_frame"]
+# Noisy circuits: quantum trajectories on the state vector, an exact density-
+# matrix referee, and cross-entropy fidelity against the digital error model.
+quantum_noisy = ["quantum_sv"]
 # 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"]
@@ -115,7 +118,7 @@ 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_ldpc", "quantum_qio", "quantum_uf", "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_qio", "quantum_uf", "quantum_noisy", "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",

wai/wai-quantum/README.md

@@ -66,6 +66,7 @@ never creates information.
 | `quantum_resource_receipt` | a signed, energy-accounted receipt over an estimate, re-checkable bit for bit on any machine |
 | `quantum_qio` | **quantum-inspired optimisation**: Ising and Max-Cut problems (G-set reader, spin-glass generator), simulated bifurcation in its ballistic, discrete and heated forms (batched eight trials per pass, bit-identical to one at a time), simulated annealing, single-flip descent, an exact Gray-code referee up to 34 spins, and step-to-solution. Reaches the best-known cut on seven G-set instances |
 | `quantum_uf` | **union-find decoding in near-linear time** over the same matching graph: event-driven cluster growth on integer weights (a radix heap orders the fill events), merges, and peeling. Measured against exact matching on the same shots: within 0.92–1.18× the failures, up to 77× faster at `d = 11` |
+| `quantum_noisy` | **noisy circuits**: quantum trajectories on the state vector (depolarizing, dephasing, amplitude damping, readout), an exact density-matrix referee up to 12 qubits that reuses the same kernels, linear cross-entropy fidelity (sampled, exact and the circuit's own), and the digital error model. Measures where that model is off: below the true fidelity by a margin that grows with the error rate, not with depth |
 | `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 |
@@ -177,6 +178,27 @@ in one pass. On the G-set Max-Cut benchmark the best-known cut is reached on
 G1, G2, G3, G11, G14, G22 and G43. On G55 the best found is 10,292 against
 10,299. The module documents each setting and success rate, and the bound
 that sparse graphs put on the time step.
+## Noisy circuits
+
+```sh
+wai-quantum noisy --file circuit.jsonl --depol2 0.01 --depol1 0.001 --readout 0.01 --trajectories 4000 --xeb
+wai-quantum noisy --file circuit.jsonl --damping 0.02 --expect "Z0; X3 X4" --exact
+```
+
+Every quantum-advantage experiment so far runs on noisy hardware and reports
+a fidelity. For random-circuit sampling that fidelity is a cross-entropy
+estimate, compared against a digital error model that multiplies the success
+probabilities of every gate and readout. `quantum_noisy` checks that
+arithmetic three ways:
+- quantum trajectories on the state vector;
+- the exact density matrix, for registers up to 12 qubits;
+- the model itself.
+
+On random brickwork circuits under depolarizing noise, the exact fidelity sits
+above the model's prediction by a few percent. The margin grows with the error
+rate, not with depth, and it is wider for circuits that do not scramble. The
+module tabulates the margin.
+
 ## Fast decoding
 
 ```sh

wai/wai-quantum/src/lib.rs

@@ -116,6 +116,9 @@ pub mod quantum_qio;
 #[cfg(feature = "quantum_uf")]
 pub mod quantum_uf;
 
+#[cfg(feature = "quantum_noisy")]
+pub mod quantum_noisy;
+
 #[cfg(feature = "quantum_frame")]
 pub mod quantum_match;
 

wai/wai-quantum/src/quantum_noisy.rs

@@ -0,0 +1,664 @@
+//! Noisy circuits — `wai.quantum.noisy`.
+//!
+//! Every quantum-advantage experiment to date runs on noisy hardware, and its
+//! headline is a fidelity. For random-circuit sampling, the estimate is the
+//! linear cross-entropy benchmark,
+//!
+//! ```text
+//! F = 2ⁿ · ⟨p_ideal(x)⟩ − 1
+//! ```
+//!
+//! averaged over the measured bitstrings. Behind it sits a simple error model:
+//! a circuit's fidelity is the product of its gates' and readouts' success
+//! probabilities. Whether that model holds for a given circuit and noise is a
+//! question a classical computer can answer. This module answers it three
+//! ways:
+//!
+//! - **Trajectories.** [`trajectories`] runs the circuit on a pure state many
+//!   times. After each gate it draws one branch of each noise channel on the
+//!   qubits that gate touched.
+//!   - Pauli channels (depolarizing noise, dephasing) pick a Pauli and stay
+//!     unitary.
+//!   - Amplitude damping jumps with the state's own probability `γ·P(1)` and
+//!     is renormalised.
+//!
+//!   Averages over trajectories converge to the noisy expectation, with a
+//!   standard error.
+//! - **The exact density matrix.** [`DensityMatrix`] holds `ρ` as a state on
+//!   `2n` qubits: rows are qubits `0..n`, columns `n..2n`. A gate `U` acts as
+//!   `U ⊗ U*`, and a channel as its Kraus sum. It reuses the state vector's
+//!   kernels, and it is the referee for small registers.
+//! - **The digital error model.** [`digital_fidelity`] is the product of
+//!   success probabilities. [`linear_xeb`] scores samples against the ideal
+//!   distribution. [`exact_xeb`] gives the value an infinite sample would
+//!   reach.
+//!
+//! Each trajectory is seeded on its own and results are reduced in trajectory
+//! order, so a run is the same bits at any thread count and on every machine.
+//!
+//! # Checked
+//!
+//! - **Trajectories against the exact density matrix.** On 5-qubit random
+//!   circuits with every noise channel, 4,000 trajectories land within four
+//!   standard errors of the exact expectations.
+//! - **Sampling.** Sampled bitstrings, with readout flips, follow the exact
+//!   distribution, and the sampled XEB converges to [`exact_xeb`].
+//!
+//! # What the digital error model gets wrong
+//!
+//! The circuits here are 8-qubit brickwork circuits: random rotations, then CZ
+//! on alternating pairs. Under depolarizing noise, the exact XEB (normalised
+//! by the noiseless circuit's own) sits *above* the digital model:
+//!
+//! | layers | two-qubit error 0.005 | 0.01  |
+//! |--------|-----------------------|-------|
+//! | 4      | 1.035                 | 1.073 |
+//! | 8      | 1.045                 | 1.095 |
+//! | 16     | 1.034                 | 1.073 |
+//! | 32     | 1.036                 | 1.077 |
+//! | 48     | 1.041                 | 1.089 |
+//!
+//! The margin grows with the error rate but not with depth. A per-gate bias
+//! would compound with every layer; this one does not. It behaves like a
+//! fixed set of errors near the measurement that the model wrongly counts as
+//! destroying the signal. A Z error just before a computational-basis
+//! measurement, for one, changes no outcome at all.
+//!
+//! The model therefore understates the fidelity of a shallow-noise circuit by
+//! a few percent, independent of size. An experiment that reports fidelity
+//! from that model inherits the same offset.
+//!
+//! # Honest boundaries
+//!
+//! - **Markovian, gate-local noise.** Errors follow the gate that caused them,
+//!   on that gate's qubits. Crosstalk, leakage, coherent over-rotation and
+//!   correlated noise in time are not modelled.
+//! - **The exact referee is small.** The density matrix is `4ⁿ` amplitudes,
+//!   so it stops at 12 qubits.
+
+use crate::gates::{complex, Gate, C};
+use crate::quantum_sv::{StateVector, SvError};
+
+/// Per-gate noise, applied after each gate on the qubits it touched, plus
+/// readout errors.
+#[derive(Clone, Copy, Debug, PartialEq, Default)]
+pub struct NoiseModel {
+    /// Depolarizing probability after a one-qubit gate (X, Y, Z each `p/3`);
+    /// also applied to each qubit of a gate on three or more qubits.
+    pub depol1: f64,
+    /// Depolarizing probability after a two-qubit gate (each of the 15
+    /// non-identity two-qubit Paulis `p/15`).
+    pub depol2: f64,
+    /// Amplitude-damping strength `γ` per touched qubit per gate.
+    pub damping: f64,
+    /// Phase-flip probability per touched qubit per gate.
+    pub dephasing: f64,
+    /// Bit-flip probability of each measured bit.
+    pub readout: f64,
+}
+
+/// The fidelity the digital error model predicts: the product of `1 − e`
+/// over every error location. The locations are:
+/// - each gate's Pauli error;
+/// - each touched qubit's dephasing, and its damping at the average jump
+///   probability `γ/2`;
+/// - each readout.
+///
+/// Each error is taken to destroy the signal entirely; see the module notes
+/// for what that misses.
+pub fn digital_fidelity(gates: &[Gate], n: u32, noise: &NoiseModel) -> f64 {
+    let mut f = 1.0;
+    for g in gates {
+        match g.qubits.len() {
+            1 => f *= 1.0 - noise.depol1,
+            2 => f *= 1.0 - noise.depol2,
+            k => {
+                for _ in 0..k {
+                    f *= 1.0 - noise.depol1;
+                }
+            }
+        }
+        for _ in &g.qubits {
+            f *= (1.0 - noise.dephasing) * (1.0 - noise.damping / 2.0);
+        }
+    }
+    for _ in 0..n {
+        f *= 1.0 - noise.readout;
+    }
+    f
+}
+
+// ---------------------------------------------------------------------------
+// Randomness
+// ---------------------------------------------------------------------------
+
+#[derive(Clone)]
+struct Rng(u64);
+
+impl Rng {
+    fn next(&mut self) -> u64 {
+        self.0 = self.0.wrapping_add(0x9e37_79b9_7f4a_7c15);
+        let mut z = self.0;
+        z = (z ^ (z >> 30)).wrapping_mul(0xbf58_476d_1ce4_e5b9);
+        z = (z ^ (z >> 27)).wrapping_mul(0x94d0_49bb_1331_11eb);
+        z ^ (z >> 31)
+    }
+    /// Uniform on `[0, 1)`.
+    fn unit(&mut self) -> f64 {
+        (self.next() >> 11) as f64 * (1.0 / 9_007_199_254_740_992.0)
+    }
+}
+
+fn trajectory_rng(seed: u64, t: usize) -> Rng {
+    let mut r = Rng(seed ^ (t as u64).wrapping_mul(0xd1b5_4a32_d192_ed03));
+    r.next();
+    r
+}
+
+const PAULI_I: [C; 4] = [C { re: 1.0, im: 0.0 }, C { re: 0.0, im: 0.0 }, C { re: 0.0, im: 0.0 }, C { re: 1.0, im: 0.0 }];
+
+fn pauli(q: u32, p: u8) -> Gate {
+    match p {
+        1 => Gate::x(q),
+        2 => Gate::y(q),
+        3 => Gate::z(q),
+        _ => Gate { qubits: vec![q], matrix: PAULI_I.to_vec() },
+    }
+}
+
+/// The two Kraus operators of amplitude damping.
+fn damping_kraus(q: u32, gamma: f64) -> [Gate; 2] {
+    let (z, o) = (complex(0.0, 0.0), complex(1.0, 0.0));
+    [
+        Gate { qubits: vec![q], matrix: vec![o, z, z, complex((1.0 - gamma).sqrt(), 0.0)] },
+        Gate { qubits: vec![q], matrix: vec![z, complex(gamma.sqrt(), 0.0), z, z] },
+    ]
+}
+
+// ---------------------------------------------------------------------------
+// Trajectories
+// ---------------------------------------------------------------------------
+
+/// Probability that qubit `q` reads 1.
+fn prob_one(sv: &StateVector, q: u32) -> f64 {
+    let bit = 1usize << (sv.n - 1 - q);
+    let mut p = 0.0;
+    for (x, (r, i)) in sv.re.iter().zip(&sv.im).enumerate() {
+        if x & bit != 0 {
+            p += r * r + i * i;
+        }
+    }
+    p
+}
+
+fn scale(sv: &mut StateVector, s: f64) {
+    sv.re.iter_mut().for_each(|v| *v *= s);
+    sv.im.iter_mut().for_each(|v| *v *= s);
+}
+
+/// Run one trajectory of `gates` on `n` qubits under `noise`.
+fn run_trajectory(gates: &[Gate], n: u32, noise: &NoiseModel, rng: &mut Rng) -> Result<StateVector, SvError> {
+    let mut sv = StateVector::zero(n)?;
+    for g in gates {
+        sv.apply(g)?;
+        let k = g.qubits.len();
+        if k == 2 && noise.depol2 > 0.0 && rng.unit() < noise.depol2 {
+            let which = 1 + (rng.next() % 15) as u8;
+            sv.apply(&pauli(g.qubits[0], which >> 2))?;
+            sv.apply(&pauli(g.qubits[1], which & 3))?;
+        }
+        for &q in &g.qubits {
+            if k != 2 && noise.depol1 > 0.0 && rng.unit() < noise.depol1 {
+                sv.apply(&pauli(q, 1 + (rng.next() % 3) as u8))?;
+            }
+            if noise.dephasing > 0.0 && rng.unit() < noise.dephasing {
+                sv.apply(&Gate::z(q))?;
+            }
+            if noise.damping > 0.0 {
+                let p1 = prob_one(&sv, q);
+                let [k0, k1] = damping_kraus(q, noise.damping);
+                let jump = noise.damping * p1;
+                if rng.unit() < jump {
+                    sv.apply(&k1)?;
+                    scale(&mut sv, 1.0 / jump.sqrt());
+                } else {
+                    sv.apply(&k0)?;
+                    scale(&mut sv, 1.0 / (1.0 - jump).sqrt());
+                }
+            }
+        }
+    }
+    Ok(sv)
+}
+
+/// A Pauli observable: `Σ c·P`, each string a list of `(qubit, 'X'|'Y'|'Z')`.
+pub type Observable = Vec<(f64, Vec<(u32, char)>)>;
+
+/// One trajectory's observable values and measured bitstring.
+type Outcome = Result<(Vec<f64>, u64), SvError>;
+
+/// A mean over trajectories and its standard error.
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub struct Estimate {
+    pub mean: f64,
+    pub stderr: f64,
+}
+
+fn estimate(values: &[f64]) -> Estimate {
+    let n = values.len() as f64;
+    let mean = values.iter().sum::<f64>() / n;
+    let var = if values.len() > 1 { values.iter().map(|v| (v - mean) * (v - mean)).sum::<f64>() / (n - 1.0) } else { 0.0 };
+    Estimate { mean, stderr: (var / n).sqrt() }
+}
+
+/// The outcome of a trajectory run.
+#[derive(Clone, Debug, PartialEq)]
+pub struct TrajectoryRun {
+    /// Each observable's mean over trajectories.
+    pub expectations: Vec<Estimate>,
+    /// One measured bitstring per trajectory, readout noise included.
+    pub samples: Vec<u64>,
+}
+
+/// Run `count` noisy trajectories of `gates` on `n` qubits. Each yields the
+/// observables' values and one measured bitstring (with readout flips).
+/// Trajectory `t` is seeded from `(seed, t)`; results come back in trajectory
+/// order whatever `threads` is (wasm32 runs one).
+pub fn trajectories(gates: &[Gate], n: u32, noise: &NoiseModel, observables: &[Observable], count: usize, seed: u64, threads: usize) -> Result<TrajectoryRun, SvError> {
+    let run = |t: usize| -> Outcome {
+        let mut rng = trajectory_rng(seed, t);
+        let sv = run_trajectory(gates, n, noise, &mut rng)?;
+        let values: Vec<f64> = observables.iter().map(|o| sv.expectation(o)).collect();
+        let mut x = sv.sample(1, rng.next())[0];
+        if noise.readout > 0.0 {
+            for q in 0..n {
+                if rng.unit() < noise.readout {
+                    x ^= 1 << (n - 1 - q);
+                }
+            }
+        }
+        Ok((values, x))
+    };
+    let threads = if cfg!(target_arch = "wasm32") { 1 } else { threads.max(1).min(count.max(1)) };
+    let results: Vec<Outcome> = if threads == 1 {
+        (0..count).map(run).collect()
+    } else {
+        let next = std::sync::atomic::AtomicUsize::new(0);
+        let slots: Vec<std::sync::Mutex<Option<Outcome>>> = (0..count).map(|_| std::sync::Mutex::new(None)).collect();
+        std::thread::scope(|scope| {
+            for _ in 0..threads {
+                scope.spawn(|| loop {
+                    let t = next.fetch_add(1, std::sync::atomic::Ordering::Relaxed);
+                    if t >= count {
+                        break;
+                    }
+                    *slots[t].lock().unwrap() = Some(run(t));
+                });
+            }
+        });
+        slots.into_iter().map(|m| m.into_inner().unwrap().expect("every trajectory ran")).collect()
+    };
+    let mut per_obs: Vec<Vec<f64>> = vec![Vec::with_capacity(count); observables.len()];
+    let mut samples = Vec::with_capacity(count);
+    for r in results {
+        let (values, x) = r?;
+        for (k, v) in values.into_iter().enumerate() {
+            per_obs[k].push(v);
+        }
+        samples.push(x);
+    }
+    Ok(TrajectoryRun { expectations: per_obs.iter().map(|v| estimate(v)).collect(), samples })
+}
+
+/// Linear cross-entropy fidelity of `samples` against the ideal state:
+/// `2ⁿ · mean(p_ideal(x)) − 1`.
+pub fn linear_xeb(ideal: &StateVector, samples: &[u64]) -> Estimate {
+    let dim = (1u64 << ideal.n) as f64;
+    let v: Vec<f64> = samples
+        .iter()
+        .map(|&x| {
+            let a = ideal.amplitude(x as usize);
+            dim * (a.re * a.re + a.im * a.im) - 1.0
+        })
+        .collect();
+    estimate(&v)
+}
+
+// ---------------------------------------------------------------------------
+// The exact density matrix
+// ---------------------------------------------------------------------------
+
+/// Largest register the exact density matrix accepts.
+pub const MAX_DENSITY_QUBITS: u32 = 12;
+
+/// `ρ` on `n` qubits, stored as a state on `2n` qubits: entry `(r, c)` at
+/// index `r·2ⁿ + c`.
+#[derive(Clone, Debug, PartialEq)]
+pub struct DensityMatrix {
+    pub n: u32,
+    vec: StateVector,
+}
+
+fn conj(g: &Gate, shift: u32) -> Gate {
+    Gate { qubits: g.qubits.iter().map(|&q| q + shift).collect(), matrix: g.matrix.iter().map(|c| complex(c.re, -c.im)).collect() }
+}
+
+impl DensityMatrix {
+    /// `|0…0⟩⟨0…0|`.
+    pub fn zero(n: u32) -> Result<DensityMatrix, SvError> {
+        if n > MAX_DENSITY_QUBITS {
+            return Err(SvError::TooLarge);
+        }
+        Ok(DensityMatrix { n, vec: StateVector::zero(2 * n)? })
+    }
+
+    /// `ρ → K ρ K†`.
+    fn sandwich(&mut self, k: &Gate) -> Result<(), SvError> {
+        self.vec.apply(k)?;
+        self.vec.apply(&conj(k, self.n))
+    }
+
+    /// `ρ → Σ_i w_i K_i ρ K_i†`.
+    fn channel(&mut self, terms: &[(f64, Gate)]) -> Result<(), SvError> {
+        let mut acc = StateVector { n: self.vec.n, re: vec![0.0; self.vec.re.len()], im: vec![0.0; self.vec.im.len()] };
+        for (w, k) in terms {
+            let mut t = DensityMatrix { n: self.n, vec: self.vec.clone() };
+            t.sandwich(k)?;
+            for (a, b) in acc.re.iter_mut().zip(&t.vec.re) {
+                *a += w * b;
+            }
+            for (a, b) in acc.im.iter_mut().zip(&t.vec.im) {
+                *a += w * b;
+            }
+        }
+        self.vec = acc;
+        Ok(())
+    }
+
+    fn depolarize1(&mut self, q: u32, p: f64) -> Result<(), SvError> {
+        let mut terms = vec![(1.0 - p, pauli(q, 0))];
+        for k in 1..4 {
+            terms.push((p / 3.0, pauli(q, k)));
+        }
+        self.channel(&terms)
+    }
+
+    fn depolarize2(&mut self, a: u32, b: u32, p: f64) -> Result<(), SvError> {
+        let mut terms = Vec::with_capacity(16);
+        for k in 0..16u8 {
+            let w = if k == 0 { 1.0 - p } else { p / 15.0 };
+            let (pa, pb) = (pauli(a, k >> 2), pauli(b, k & 3));
+            terms.push((w, two_qubit(&pa, &pb)));
+        }
+        self.channel(&terms)
+    }
+
+    /// Run `gates` under `noise` from `|0…0⟩⟨0…0|`, exactly.
+    pub fn run(gates: &[Gate], n: u32, noise: &NoiseModel) -> Result<DensityMatrix, SvError> {
+        let mut rho = DensityMatrix::zero(n)?;
+        for g in gates {
+            rho.sandwich(g)?;
+            let k = g.qubits.len();
+            if k == 2 && noise.depol2 > 0.0 {
+                rho.depolarize2(g.qubits[0], g.qubits[1], noise.depol2)?;
+            }
+            for &q in &g.qubits {
+                if k != 2 && noise.depol1 > 0.0 {
+                    rho.depolarize1(q, noise.depol1)?;
+                }
+                if noise.dephasing > 0.0 {
+                    rho.channel(&[(1.0 - noise.dephasing, pauli(q, 0)), (noise.dephasing, pauli(q, 3))])?;
+                }
+                if noise.damping > 0.0 {
+                    let [k0, k1] = damping_kraus(q, noise.damping);
+                    rho.channel(&[(1.0, k0), (1.0, k1)])?;
+                }
+            }
+        }
+        Ok(rho)
+    }
+
+    /// Entry `ρ_{r,c}`.
+    pub fn entry(&self, r: usize, c: usize) -> C {
+        self.vec.amplitude((r << self.n) | c)
+    }
+
+    /// `Tr ρ`.
+    pub fn trace(&self) -> f64 {
+        (0..1usize << self.n).map(|x| self.entry(x, x).re).sum()
+    }
+
+    /// `Tr(ρ Σ c·P)`.
+    pub fn expectation(&self, observable: &[(f64, Vec<(u32, char)>)]) -> f64 {
+        let n = self.n;
+        let mut total = 0.0;
+        for (coef, string) in observable {
+            let (mut flip, mut zmask, mut ys) = (0usize, 0usize, 0u32);
+            for &(q, p) in string {
+                let bit = 1usize << (n - 1 - q);
+                match p {
+                    'X' => flip |= bit,
+                    'Y' => {
+                        flip |= bit;
+                        zmask |= bit;
+                        ys += 1;
+                    }
+                    'Z' => zmask |= bit,
+                    _ => {}
+                }
+            }
+            // Tr(Pρ) = Σ_y ⟨y⊕f| P |y⟩ ρ_{y, y⊕f} = Σ_y i^{ys} (−1)^{|y∧z|} ρ_{y, y⊕f}.
+            let mut v = 0.0;
+            for y in 0..1usize << n {
+                let e = self.entry(y, y ^ flip);
+                let sign = if (y & zmask).count_ones().is_multiple_of(2) { 1.0 } else { -1.0 };
+                v += sign
+                    * match ys % 4 {
+                        0 => e.re,
+                        1 => -e.im,
+                        2 => -e.re,
+                        _ => e.im,
+                    };
+            }
+            total += coef * v;
+        }
+        total
+    }
+
+    /// The measured distribution, readout flips included.
+    pub fn probabilities(&self, readout: f64) -> Vec<f64> {
+        let n = self.n;
+        let mut p: Vec<f64> = (0..1usize << n).map(|x| self.entry(x, x).re).collect();
+        if readout > 0.0 {
+            for q in 0..n {
+                let bit = 1usize << (n - 1 - q);
+                for x in 0..p.len() {
+                    if x & bit == 0 {
+                        let (a, b) = (p[x], p[x | bit]);
+                        p[x] = (1.0 - readout) * a + readout * b;
+                        p[x | bit] = readout * a + (1.0 - readout) * b;
+                    }
+                }
+            }
+        }
+        p
+    }
+}
+
+/// The tensor product `a ⊗ b` of two one-qubit gates on different qubits.
+fn two_qubit(a: &Gate, b: &Gate) -> Gate {
+    let mut m = vec![complex(0.0, 0.0); 16];
+    for r in 0..4 {
+        for c in 0..4 {
+            let (x, y) = (a.matrix[(r >> 1) * 2 + (c >> 1)], b.matrix[(r & 1) * 2 + (c & 1)]);
+            m[r * 4 + c] = complex(x.re * y.re - x.im * y.im, x.re * y.im + x.im * y.re);
+        }
+    }
+    Gate { qubits: vec![a.qubits[0], b.qubits[0]], matrix: m }
+}
+
+/// The noiseless circuit's own cross-entropy, `2ⁿ Σ_x p_ideal(x)² − 1`: 1 for
+/// a fully scrambling circuit, larger for one that is not. Noisy values are
+/// compared with the digital model after dividing by it.
+pub fn ideal_xeb(ideal: &StateVector) -> f64 {
+    let dim = (1u64 << ideal.n) as f64;
+    let mut s = 0.0;
+    for (r, i) in ideal.re.iter().zip(&ideal.im) {
+        let q = r * r + i * i;
+        s += q * q;
+    }
+    dim * s - 1.0
+}
+
+/// The cross-entropy fidelity an infinite sample would give:
+/// `2ⁿ Σ_x p_noisy(x) p_ideal(x) − 1`.
+pub fn exact_xeb(ideal: &StateVector, noisy: &DensityMatrix, readout: f64) -> f64 {
+    let dim = (1u64 << ideal.n) as f64;
+    let p = noisy.probabilities(readout);
+    let mut s = 0.0;
+    for (x, px) in p.iter().enumerate() {
+        let a = ideal.amplitude(x);
+        s += px * (a.re * a.re + a.im * a.im);
+    }
+    dim * s - 1.0
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    /// A brickwork random circuit: random one-qubit rotations, then CZ on
+    /// alternating pairs, `layers` times.
+    fn random_circuit(n: u32, layers: usize, seed: u64) -> Vec<Gate> {
+        let mut rng = Rng(seed);
+        let mut g = Vec::new();
+        for l in 0..layers {
+            for q in 0..n {
+                let a = core::f64::consts::TAU * rng.unit();
+                g.push(match rng.next() % 3 {
+                    0 => Gate::rx(q, a),
+                    1 => Gate::ry(q, a),
+                    _ => Gate::sqrt_w(q),
+                });
+            }
+            let mut q = (l % 2) as u32;
+            while q + 1 < n {
+                g.push(Gate::cz(q, q + 1));
+                q += 2;
+            }
+        }
+        g
+    }
+
+    fn observables(n: u32) -> Vec<Observable> {
+        vec![vec![(1.0, vec![(0, 'Z')])], vec![(1.0, vec![(1, 'X'), (2, 'X')])], vec![(1.0, vec![(n - 1, 'Y'), (0, 'Z')])], vec![(0.5, vec![(2, 'Z')]), (0.5, vec![(3, 'Y')])]]
+    }
+
+    #[test]
+    fn noiseless_runs_match_the_pure_state() {
+        let g = random_circuit(5, 4, 1);
+        let mut sv = StateVector::zero(5).unwrap();
+        sv.run(&g, 0).unwrap();
+        let rho = DensityMatrix::run(&g, 5, &NoiseModel::default()).unwrap();
+        let tr = trajectories(&g, 5, &NoiseModel::default(), &observables(5), 3, 2, 2).unwrap();
+        for (k, o) in observables(5).iter().enumerate() {
+            let pure = sv.expectation(o);
+            assert!((rho.expectation(o) - pure).abs() < 1e-12);
+            assert!((tr.expectations[k].mean - pure).abs() < 1e-12);
+            assert!(tr.expectations[k].stderr < 1e-12);
+        }
+        assert!((rho.trace() - 1.0).abs() < 1e-12);
+    }
+
+    #[test]
+    fn trajectories_converge_to_the_density_matrix() {
+        let n = 5;
+        let g = random_circuit(n, 5, 7);
+        let noise = NoiseModel { depol1: 0.01, depol2: 0.03, damping: 0.02, dephasing: 0.01, readout: 0.0 };
+        let rho = DensityMatrix::run(&g, n, &noise).unwrap();
+        assert!((rho.trace() - 1.0).abs() < 1e-12);
+        let tr = trajectories(&g, n, &noise, &observables(n), 4000, 11, 4).unwrap();
+        for (k, o) in observables(n).iter().enumerate() {
+            let exact = rho.expectation(o);
+            let e = tr.expectations[k];
+            assert!((e.mean - exact).abs() < 4.0 * e.stderr + 1e-9, "observable {k}: {} ± {} vs {exact}", e.mean, e.stderr);
+        }
+    }
+
+    #[test]
+    fn readout_noise_and_sampling_match_the_exact_distribution() {
+        let n = 4;
+        let g = random_circuit(n, 3, 3);
+        let noise = NoiseModel { depol2: 0.02, readout: 0.05, ..NoiseModel::default() };
+        let p = DensityMatrix::run(&g, n, &noise).unwrap().probabilities(noise.readout);
+        assert!((p.iter().sum::<f64>() - 1.0).abs() < 1e-12);
+        let shots = 20_000;
+        let tr = trajectories(&g, n, &noise, &[], shots, 5, 4).unwrap();
+        let mut counts = vec![0usize; 1 << n];
+        for &x in &tr.samples {
+            counts[x as usize] += 1;
+        }
+        for (x, &c) in counts.iter().enumerate() {
+            let expected = p[x] * shots as f64;
+            let sd = (expected * (1.0 - p[x])).sqrt().max(1.0);
+            assert!((c as f64 - expected).abs() < 5.0 * sd, "x={x}: {c} vs {expected:.1}");
+        }
+    }
+
+    #[test]
+    fn runs_do_not_depend_on_threads() {
+        let g = random_circuit(6, 4, 9);
+        let noise = NoiseModel { depol1: 0.02, depol2: 0.05, damping: 0.01, dephasing: 0.01, readout: 0.02 };
+        let a = trajectories(&g, 6, &noise, &observables(6), 50, 3, 1).unwrap();
+        let b = trajectories(&g, 6, &noise, &observables(6), 50, 3, 5).unwrap();
+        assert_eq!(a, b);
+    }
+
+    fn self_xeb(ideal: &StateVector) -> f64 {
+        ideal_xeb(ideal)
+    }
+
+    /// The exact XEB, normalised by the noiseless circuit's own, over the
+    /// digital model's prediction.
+    fn model_ratio(n: u32, layers: usize, p2: f64) -> f64 {
+        let g = random_circuit(n, layers, 21);
+        let mut ideal = StateVector::zero(n).unwrap();
+        ideal.run(&g, 0).unwrap();
+        let noise = NoiseModel { depol1: p2 / 10.0, depol2: p2, ..NoiseModel::default() };
+        let exact = exact_xeb(&ideal, &DensityMatrix::run(&g, n, &noise).unwrap(), 0.0);
+        exact / self_xeb(&ideal) / digital_fidelity(&g, n, &noise)
+    }
+
+    #[test]
+    fn the_digital_model_is_low_by_a_depth_independent_margin() {
+        // Under Pauli noise the exact (normalised) XEB sits above the digital
+        // model by a margin that grows with the error rate but not with depth:
+        // a fixed number of errors near the end, not a per-gate bias.
+        for p2 in [0.005, 0.01] {
+            let (shallow, deep) = (model_ratio(7, 8, p2), model_ratio(7, 24, p2));
+            eprintln!("p2={p2}: {shallow:.4} {deep:.4}");
+            assert!(shallow > 1.0 && deep > 1.0, "p2={p2}: {shallow} {deep}");
+            assert!((shallow - deep).abs() < 0.03, "p2={p2}: {shallow} {deep}");
+            assert!(deep < 1.0 + 12.0 * p2, "p2={p2}: {deep}");
+        }
+        // Noiseless: the exact XEB is the ideal circuit's own, 2ⁿ Σp² − 1.
+        let g = random_circuit(7, 8, 21);
+        let mut ideal = StateVector::zero(7).unwrap();
+        ideal.run(&g, 0).unwrap();
+        let noiseless = DensityMatrix::run(&g, 7, &NoiseModel::default()).unwrap();
+        assert!((exact_xeb(&ideal, &noiseless, 0.0) - self_xeb(&ideal)).abs() < 1e-9);
+    }
+
+    #[test]
+    fn sampled_xeb_converges_to_the_exact_value() {
+        let n = 7;
+        let g = random_circuit(n, 7, 4);
+        let mut ideal = StateVector::zero(n).unwrap();
+        ideal.run(&g, 0).unwrap();
+        let noise = NoiseModel { depol1: 0.002, depol2: 0.015, readout: 0.01, ..NoiseModel::default() };
+        let exact = exact_xeb(&ideal, &DensityMatrix::run(&g, n, &noise).unwrap(), noise.readout);
+        let tr = trajectories(&g, n, &noise, &[], 20_000, 8, 4).unwrap();
+        let est = linear_xeb(&ideal, &tr.samples);
+        assert!((est.mean - exact).abs() < 4.0 * est.stderr, "{} ± {} vs {exact}", est.mean, est.stderr);
+    }
+}

wai/wai-quantum/src/quantum_sv.rs

@@ -1,6 +1,6 @@
 //! The dense state vector, any gate at any angle — `wai.quantum.sv`.
 //!
-//! The crate's byte-exact simulator ([`crate::quantum`]) is fixed-point and
+//! The crate's byte-exact simulator (`quantum`) is fixed-point and
 //! limited to the dyadic gate family. That limit is deliberate: it is what
 //! makes its amplitudes a portable hash. Most circuits are not dyadic.
 //! Rotations at arbitrary angles, √W, fSim and random two-qubit unitaries need

wai/wai-quantum/tests/bit_identity.rs

@@ -23,7 +23,8 @@
     feature = "quantum_arch",
     feature = "quantum_ldpc",
     feature = "quantum_qio",
-    feature = "quantum_uf"
+    feature = "quantum_uf",
+    feature = "quantum_noisy"
 ))]
 
 use wai_quantum::*;
@@ -409,3 +410,42 @@ fn union_find_decoding() {
     assert!(v.iter().any(|&x| x != 0.0), "the sample must include flips");
     assert_eq!(digest(&v), 0xe43c_2ad4_9244_3118);
 }
+
+/// Noisy circuits: trajectories under every channel (Pauli branches, damping
+/// jumps with renormalisation, readout flips), the exact density matrix and
+/// both cross-entropy figures.
+#[test]
+fn noisy_circuits() {
+    use gates::Gate;
+    use quantum_noisy::*;
+    let n = 6u32;
+    let mut g = Vec::new();
+    let mut a = 0.3f64;
+    for layer in 0..5u32 {
+        for q in 0..n {
+            a += 0.917;
+            g.push(if (q + layer) % 2 == 0 { Gate::rx(q, a) } else { Gate::ry(q, a) });
+        }
+        let mut q = layer % 2;
+        while q + 1 < n {
+            g.push(Gate::cz(q, q + 1));
+            q += 2;
+        }
+    }
+    let noise = NoiseModel { depol1: 0.01, depol2: 0.03, damping: 0.02, dephasing: 0.01, readout: 0.02 };
+    let obs = vec![vec![(1.0, vec![(0, 'Z')])], vec![(1.0, vec![(1, 'X'), (4, 'Y')])]];
+    let run = trajectories(&g, n, &noise, &obs, 40, 6, 3).unwrap();
+    let mut v: Vec<f64> = run.expectations.iter().flat_map(|e| [e.mean, e.stderr]).collect();
+    v.extend(run.samples.iter().map(|&x| x as f64));
+    let rho = DensityMatrix::run(&g, n, &noise).unwrap();
+    v.extend(obs.iter().map(|o| rho.expectation(o)));
+    v.push(rho.trace());
+    let mut ideal = quantum_sv::StateVector::zero(n).unwrap();
+    ideal.run(&g, 0).unwrap();
+    v.push(exact_xeb(&ideal, &rho, noise.readout));
+    let x = linear_xeb(&ideal, &run.samples);
+    v.push(x.mean);
+    v.push(x.stderr);
+    v.push(digital_fidelity(&g, n, &noise));
+    assert_eq!(digest(&v), 0x78b2_ddfb_e7c3_6bea);
+}

wai/wai-rs/Cargo.lock

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

wai/wai-rs/Cargo.toml

@@ -361,7 +361,7 @@ quantum_tn = ["dep:wai-quantum", "wai-quantum/quantum_tn"]
 # Gaussian boson sampling: hafnians, torontonians, exact probabilities, samples, checks.
 quantum_gbs = ["dep:wai-quantum", "wai-quantum/quantum_gbs"]
 # The dense state vector at any angle: fusion, threads, the same bits everywhere.
-quantum_sv = ["dep:wai-quantum", "wai-quantum/quantum_sv"]
+quantum_sv = ["dep:wai-quantum", "wai-quantum/quantum_sv", "wai-quantum/quantum_noisy"]
 # Architecture-level estimation: yoked storage, cultivation, CCZ factories,
 # reaction-limited time, whole machines.
 quantum_arch = ["dep:wai-quantum", "wai-quantum/quantum_arch"]

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

@@ -177,6 +177,9 @@ fn usage() -> ! {
          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  \
          wai-quantum arch cultivate [--p 5e-4|1e-3|2e-3] [--target T] [--json]\n\n\
+         NOISY CIRCUITS (trajectories, exact density matrix, cross-entropy):\n  \
+         wai-quantum noisy --file CIRCUIT.jsonl [-n N] [--depol1 P] [--depol2 P] [--damping G] [--dephasing P]\n    \
+         [--readout P] [--trajectories T] [--seed S] [--expect \"Z0 Z1; X3\"] [--xeb] [--exact] [--json]\n\n\
          STATE VECTOR (any gate, any angle):\n  \
          wai-quantum sv --file CIRCUIT.jsonl [-n N] [--fuse W] [--top K | --sample S [--seed X] | --expect \"Z0 Z1; X3\"] [--json]\n\n\
          TENSOR NETWORKS:\n  \
@@ -1188,6 +1191,98 @@ fn main() {
             }
         }
 
+        // Noisy circuits: trajectories on the state vector, optionally the
+        // exact density matrix (up to 12 qubits) and cross-entropy fidelity
+        // against the ideal circuit and the digital error model.
+        "noisy" => {
+            use wai::quantum_noisy::*;
+            let path = opt(args, "--file", "").unwrap_or_else(|| die("--file CIRCUIT.jsonl is required"));
+            let (gates, n) = read_gate_file(&path);
+            let n: u32 = num(args, "--qubits", "-n", n);
+            let noise = NoiseModel {
+                depol1: num(args, "--depol1", "", 0.0),
+                depol2: num(args, "--depol2", "", 0.0),
+                damping: num(args, "--damping", "", 0.0),
+                dephasing: num(args, "--dephasing", "", 0.0),
+                readout: num(args, "--readout", "", 0.0),
+            };
+            let count: usize = num(args, "--trajectories", "", 1000);
+            let seed: u64 = num(args, "--seed", "", 1);
+            let threads: usize = num(args, "--threads", "", std::thread::available_parallelism().map_or(1, |n| n.get()));
+            let mut names = Vec::new();
+            let mut observables = Vec::new();
+            if let Some(spec) = opt(args, "--expect", "") {
+                for term in spec.split(';').map(str::trim).filter(|t| !t.is_empty()) {
+                    let string: Vec<(u32, char)> = term
+                        .split_whitespace()
+                        .map(|f| {
+                            let (p, q) = f.split_at(1);
+                            let p = p.chars().next().unwrap().to_ascii_uppercase();
+                            if !"XYZ".contains(p) {
+                                die("Pauli factors are X, Y or Z followed by a qubit, as in Z0");
+                            }
+                            (q.parse().unwrap_or_else(|_| die("bad qubit in --expect")), p)
+                        })
+                        .collect();
+                    if string.iter().any(|&(q, _)| q >= n) {
+                        die("--expect names a qubit outside the register");
+                    }
+                    names.push(term.to_string());
+                    observables.push(vec![(1.0, string)]);
+                }
+            }
+            let run = trajectories(&gates, n, &noise, &observables, count, seed, threads).unwrap_or_else(|e| die(&format!("{e:?}")));
+            let exact = if flag(args, "--exact") || (flag(args, "--xeb") && n <= MAX_DENSITY_QUBITS) {
+                Some(DensityMatrix::run(&gates, n, &noise).unwrap_or_else(|e| die(&format!("{e:?} (the exact density matrix stops at {MAX_DENSITY_QUBITS} qubits)"))))
+            } else {
+                None
+            };
+            let xeb = if flag(args, "--xeb") {
+                let mut ideal = wai::quantum_sv::StateVector::zero(n).unwrap_or_else(|e| die(&format!("{e:?}")));
+                ideal.run(&gates, 2).unwrap_or_else(|e| die(&format!("{e:?}")));
+                let sampled = linear_xeb(&ideal, &run.samples);
+                let exact_x = exact.as_ref().map(|rho| exact_xeb(&ideal, rho, noise.readout));
+                Some((sampled, exact_x, digital_fidelity(&gates, n, &noise), ideal_xeb(&ideal)))
+            } else {
+                None
+            };
+            if flag(args, "--json") {
+                let mut parts = vec![format!("\"qubits\":{n}"), format!("\"trajectories\":{count}")];
+                let ex: Vec<String> = names
+                    .iter()
+                    .enumerate()
+                    .map(|(k, t)| {
+                        let e = run.expectations[k];
+                        let exact_v = exact.as_ref().map_or("null".into(), |rho| rho.expectation(&observables[k]).to_string());
+                        format!("{{\"exact\":{exact_v},\"mean\":{},\"stderr\":{},\"term\":{}}}", e.mean, e.stderr, serde_json::to_string(t).unwrap())
+                    })
+                    .collect();
+                parts.push(format!("\"expectations\":[{}]", ex.join(",")));
+                if let Some((s, e, d, i)) = xeb {
+                    parts.push(format!("\"xeb\":{{\"digital_model\":{d},\"exact\":{},\"ideal\":{i},\"sampled\":{},\"stderr\":{}}}", e.map_or("null".into(), |v| v.to_string()), s.mean, s.stderr));
+                }
+                parts.sort();
+                println!("{{{}}}", parts.join(","));
+            } else {
+                println!("{n} qubits, {} gates, {count} trajectories", gates.len());
+                for (k, t) in names.iter().enumerate() {
+                    let e = run.expectations[k];
+                    match &exact {
+                        Some(rho) => println!("<{t}> = {:+.6} ± {:.6}   exact {:+.6}", e.mean, e.stderr, rho.expectation(&observables[k])),
+                        None => println!("<{t}> = {:+.6} ± {:.6}", e.mean, e.stderr),
+                    }
+                }
+                if let Some((s, e, d, i)) = xeb {
+                    println!("XEB noiseless {i:.4}   (the circuit's own; 1 when fully scrambling)");
+                    println!("XEB sampled   {:.4} ± {:.4}   normalised {:.4}", s.mean, s.stderr, s.mean / i);
+                    if let Some(v) = e {
+                        println!("XEB exact     {v:.4}   normalised {:.4}", v / i);
+                    }
+                    println!("digital model {d:.4}");
+                }
+            }
+        }
+
         // The dense state vector, any gate at any angle: run a circuit file
         // (fused to --fuse qubits), then show the likeliest outcomes, draw
         // samples, or evaluate Pauli observables.

wai/wai-rs/src/lib.rs

@@ -648,6 +648,9 @@ pub use wai_quantum::quantum_gbs;
 #[cfg(feature = "quantum_sv")]
 pub use wai_quantum::quantum_sv;
 
+#[cfg(feature = "quantum_sv")]
+pub use wai_quantum::quantum_noisy;
+
 #[cfg(feature = "quantum_arch")]
 pub use wai_quantum::quantum_arch;
 

wai/wai-web/Cargo.lock

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