From bf22681fae634d7c5f2f21f676f71fb9a734e172 Mon Sep 17 00:00:00 2001 From: Transaction Science Date: Mon, 19 Oct 2026 06:51:04 -0400 Subject: [PATCH 1/5] wai-quantum: exact amplitudes and samples of a Clifford circuit measured in a product basis, by cut-rank elimination MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Tableau::graph_form writes the state as local Cliffords on a graph state. quantum_cutrank (feature quantum_cutrank, in full) then sums ⟨x|⊗V_q|ψ⟩ = 2^{-n/2} Σ_z ∏ W_q[x_q][z_q] (−1)^{Σ_{a~b} z_a z_b} one qubit at a time along an order, in a table of 2^{cut-rank} entries: each qubit is an in-place pass (a merge where the boundary space loses a dimension, the qubit's signed terms, a doubling or an XOR-shifted sum for its edges ahead). By Parseval the table's norm is the prefix probability, so the same passes sample exactly. Reductions are summed in fixed chunks. Tests: the graph form rebuilds the state (40 random circuits); amplitudes equal the state vector's in probability and relative phase (30 circuits up to 12 qubits); forward, reverse and shuffled orders agree (24–54 qubits, at least 30 comparisons); samples follow the exact law (60k shots) and are seeded; cut-ranks of GHZ and product states. Mutants caught: merge signs, the sampling choice, the graph form without its S, the doubling's sign. examples/cutrank_qasm reads an OpenQASM 2 Clifford bulk with a closing single-qubit layer. On the tracker's 70-qubit random graph-state instance (max cut 32 bits, matching an independent tableau): an exact amplitude in about 7.6 s and an exact sample in about 14 s on the build host, 69 GB peak in f64. Co-Authored-By: Claude Opus 5.5 --- wai/wai-quantum/Cargo.toml | 9 +- wai/wai-quantum/examples/cutrank_qasm.rs | 145 ++++ wai/wai-quantum/src/lib.rs | 2 + wai/wai-quantum/src/quantum_cutrank.rs | 883 +++++++++++++++++++++++ wai/wai-quantum/src/quantum_tableau.rs | 94 +++ 5 files changed, 1132 insertions(+), 1 deletion(-) create mode 100644 wai/wai-quantum/examples/cutrank_qasm.rs create mode 100644 wai/wai-quantum/src/quantum_cutrank.rs diff --git a/wai/wai-quantum/Cargo.toml b/wai/wai-quantum/Cargo.toml index 66f7547ff..40fbf2860 100644 --- a/wai/wai-quantum/Cargo.toml +++ b/wai/wai-quantum/Cargo.toml @@ -77,6 +77,9 @@ quantum_toolchain = ["quantum"] quantum_mps = ["quantum"] quantum_stabilizer = ["quantum"] quantum_tableau = [] +# Exact amplitudes, probabilities and samples of a Clifford circuit measured in a +# product basis, at a cost set by its cut-rank (src/quantum_cutrank.rs). +quantum_cutrank = ["quantum_tableau", "quantum_sv"] quantum_mpsim = ["quantum_sv", "quantum_mps"] quantum_pauli = ["quantum"] # Sparse Pauli dynamics at utility scale: arbitrary-angle Pauli rotations on up @@ -203,7 +206,7 @@ quantum_phasor_meter = ["_signing"] # Everything. full = [ - "quantum", "quantum_toolchain", "quantum_mps", "quantum_stabilizer", "quantum_tableau", "quantum_mpsim", "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_dmrg", "quantum_tdvp", "quantum_chem", "quantum_nqs", "quantum_qir", "quantum_info", "quantum_qudit", "quantum_comm", "quantum_source", "quantum_channel", + "quantum", "quantum_toolchain", "quantum_mps", "quantum_stabilizer", "quantum_tableau", "quantum_cutrank", "quantum_mpsim", "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_dmrg", "quantum_tdvp", "quantum_chem", "quantum_nqs", "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_fabric", "quantum_vml", "quantum_phasor", "quantum_qfhrr", "quantum_kernel", "quantum_qdata", @@ -302,6 +305,10 @@ required-features = ["quantum_fabric"] name = "rv32_listing" required-features = ["quantum_fabric_dataflow"] +[[example]] +name = "cutrank_qasm" +required-features = ["quantum_cutrank"] + [[example]] name = "npu_onnx" required-features = ["quantum_fabric"] diff --git a/wai/wai-quantum/examples/cutrank_qasm.rs b/wai/wai-quantum/examples/cutrank_qasm.rs new file mode 100644 index 000000000..a3e45e88c --- /dev/null +++ b/wai/wai-quantum/examples/cutrank_qasm.rs @@ -0,0 +1,145 @@ +//! A Clifford circuit with one closing layer of single-qubit gates, read +//! from OpenQASM 2, summed exactly by `quantum_cutrank`. +//! +//! ```text +//! cargo run --release --example cutrank_qasm --features quantum_cutrank -- FILE cuts +//! cargo run --release --example cutrank_qasm --features quantum_cutrank -- FILE amplitudes N [SEED] +//! cargo run --release --example cutrank_qasm --features quantum_cutrank -- FILE sample N [SEED] +//! ``` +//! +//! The bulk takes `h s sdg sx x y z cx cz swap`; a qubit's closing layer +//! starts at its first non-Clifford gate (`t tdg rz rx ry p u u3`) and may hold +//! any single-qubit gates. `2^n · p` is printed per bit string: for exact +//! samples of a state whose probabilities follow Porter–Thomas its mean is +//! about 2 (a linear cross-entropy of about 1). + +use std::time::Instant; +use wai_quantum::gates::{C, complex}; +use wai_quantum::quantum_cutrank::{ProductMeasured, Unitary}; +use wai_quantum::quantum_tableau::{Clifford, Tableau}; + +fn mul(a: &Unitary, b: &Unitary) -> Unitary { + let m = |x: C, y: C| complex(x.re * y.re - x.im * y.im, x.re * y.im + x.im * y.re); + let add = |x: C, y: C| complex(x.re + y.re, x.im + y.im); + [add(m(a[0], b[0]), m(a[1], b[2])), add(m(a[0], b[1]), m(a[1], b[3])), add(m(a[2], b[0]), m(a[3], b[2])), add(m(a[2], b[1]), m(a[3], b[3]))] +} + +fn one_qubit(name: &str, p: &[f64]) -> Option { + let c = complex; + let r = std::f64::consts::FRAC_1_SQRT_2; + let e = |x: f64| c(x.cos(), x.sin()); + let (o, one) = (c(0.0, 0.0), c(1.0, 0.0)); + Some(match name { + "h" => [c(r, 0.0), c(r, 0.0), c(r, 0.0), c(-r, 0.0)], + "x" => [o, one, one, o], + "y" => [o, c(0.0, -1.0), c(0.0, 1.0), o], + "z" => [one, o, o, c(-1.0, 0.0)], + "s" => [one, o, o, c(0.0, 1.0)], + "sdg" => [one, o, o, c(0.0, -1.0)], + "t" => [one, o, o, e(std::f64::consts::FRAC_PI_4)], + "tdg" => [one, o, o, e(-std::f64::consts::FRAC_PI_4)], + "sx" => [c(0.5, 0.5), c(0.5, -0.5), c(0.5, -0.5), c(0.5, 0.5)], + "rz" | "p" | "u1" => { + let t = p[0]; + if name == "rz" { [e(-t / 2.0), o, o, e(t / 2.0)] } else { [one, o, o, e(t)] } + } + "rx" => { + let (cs, sn) = ((p[0] / 2.0).cos(), (p[0] / 2.0).sin()); + [c(cs, 0.0), c(0.0, -sn), c(0.0, -sn), c(cs, 0.0)] + } + "ry" => { + let (cs, sn) = ((p[0] / 2.0).cos(), (p[0] / 2.0).sin()); + [c(cs, 0.0), c(-sn, 0.0), c(sn, 0.0), c(cs, 0.0)] + } + _ => return None, + }) +} + +fn main() { + let args: Vec = std::env::args().collect(); + let (path, mode) = (&args[1], args[2].as_str()); + let count: usize = args.get(3).and_then(|s| s.parse().ok()).unwrap_or(4); + let seed: u64 = args.get(4).and_then(|s| s.parse().ok()).unwrap_or(1); + let src = std::fs::read_to_string(path).expect("the file reads"); + let mut n = 0u32; + let mut bulk = Vec::new(); + let mut last: Vec = Vec::new(); + let mut closing: Vec = Vec::new(); + let id = [complex(1.0, 0.0), complex(0.0, 0.0), complex(0.0, 0.0), complex(1.0, 0.0)]; + for line in src.lines().map(str::trim).filter(|l| !l.is_empty()) { + let line = line.trim_end_matches(';'); + if let Some(rest) = line.strip_prefix("qreg q[") { + n = rest.trim_end_matches(']').parse().unwrap(); + last = vec![id; n as usize]; + closing = vec![false; n as usize]; + continue; + } + if line.starts_with("OPENQASM") || line.starts_with("include") || line.starts_with("creg") || line.starts_with("barrier") || line.starts_with("measure") { + continue; + } + let (head, args) = line.split_once(' ').expect("a gate and its qubits"); + let (name, params): (&str, Vec) = match head.split_once('(') { + Some((n, p)) => (n, p.trim_end_matches(')').split(',').map(|x| x.trim().parse().expect("a numeric angle")).collect()), + None => (head, Vec::new()), + }; + let qs: Vec = args.split(',').map(|a| a.trim().trim_start_matches("q[").trim_end_matches(']').parse().unwrap()).collect(); + let a = qs[0]; + let non_clifford = matches!(name, "t" | "tdg" | "rz" | "rx" | "ry" | "p" | "u1"); + if qs.len() == 1 && (closing[a as usize] || non_clifford) { + closing[a as usize] = true; + let u = one_qubit(name, ¶ms).unwrap_or_else(|| panic!("gate {name}")); + last[a as usize] = mul(&u, &last[a as usize]); + continue; + } + assert!(qs.iter().all(|&q| !closing[q as usize]), "a gate on {qs:?} after its closing layer began: {line}"); + bulk.push(match (name, qs.as_slice()) { + ("h", _) => Clifford::H(a), + ("s", _) => Clifford::S(a), + ("sdg", _) => Clifford::Sdg(a), + ("sx", _) => Clifford::SX(a), + ("x", _) => Clifford::X(a), + ("y", _) => Clifford::Y(a), + ("z", _) => Clifford::Z(a), + ("cx", [_, b]) => Clifford::CX(a, *b), + ("cz", [_, b]) => Clifford::CZ(a, *b), + ("swap", [_, b]) => Clifford::Swap(a, *b), + _ => panic!("not a bulk gate: {line}"), + }); + } + let t0 = Instant::now(); + let pm = ProductMeasured::new(&Tableau::run(n, &bulk).unwrap(), &last).unwrap(); + println!("{n} qubits, {} bulk gates; graph form and cuts in {:?}", bulk.len(), t0.elapsed()); + println!("cut-ranks: {:?}", pm.cut_ranks()); + println!("max cut {} bits: a table of {} entries ({} GiB at f64)", pm.max_cut(), 1u64 << pm.max_cut(), (16u64 << pm.max_cut()) >> 30); + let limit = std::env::var("CUTRANK_MAX").ok().and_then(|s| s.parse().ok()).unwrap_or(32); + let scale = (n as f64).exp2(); + match mode { + "cuts" => {} + "amplitudes" => { + let mut s = seed; + for i in 0..count { + let x: Vec = (0..n).map(|_| { + s = s.wrapping_add(0x9e37_79b9_7f4a_7c15); + let mut z = s; + z = (z ^ (z >> 30)).wrapping_mul(0xbf58_476d_1ce4_e5b9); + (z ^ (z >> 27)) & 1 == 1 + }).collect(); + let t = Instant::now(); + let a = pm.amplitude(&x, limit).unwrap(); + println!("amplitude {i}: {a:?}, 2^n·p = {:.6}, in {:?}", scale * (a.re * a.re + a.im * a.im), t.elapsed()); + } + } + "sample" => { + let t = Instant::now(); + let samples = pm.sample(count, seed, limit).unwrap(); + let el = t.elapsed(); + let mean = samples.iter().map(|(_, p)| scale * p).sum::() / samples.len() as f64; + for (x, p) in samples.iter().take(8) { + let bits: String = x.iter().map(|&b| if b { '1' } else { '0' }).collect(); + println!("{bits} 2^n·p = {:.6}", scale * p); + } + println!("{count} samples in {el:?} ({:?} each); mean 2^n·p = {mean:.4}, linear XEB = {:.4}", el / count as u32, mean - 1.0); + } + other => panic!("mode {other}"), + } +} diff --git a/wai/wai-quantum/src/lib.rs b/wai/wai-quantum/src/lib.rs index b62721ac3..876f07274 100644 --- a/wai/wai-quantum/src/lib.rs +++ b/wai/wai-quantum/src/lib.rs @@ -251,6 +251,8 @@ pub mod quantum_receipt; pub mod quantum_stabilizer; #[cfg(feature = "quantum_tableau")] pub mod quantum_tableau; +#[cfg(feature = "quantum_cutrank")] +pub mod quantum_cutrank; #[cfg(feature = "quantum_sv_metal")] pub mod quantum_sv_metal; diff --git a/wai/wai-quantum/src/quantum_cutrank.rs b/wai/wai-quantum/src/quantum_cutrank.rs new file mode 100644 index 000000000..fe640365e --- /dev/null +++ b/wai/wai-quantum/src/quantum_cutrank.rs @@ -0,0 +1,883 @@ +//! Exact amplitudes, probabilities and samples of a stabilizer state measured +//! in a product basis — a Clifford circuit, then one layer of single-qubit +//! gates — at a cost set by the state's cut-rank along an order of its +//! qubits, not by its size. +//! +//! The Clifford part's state is `⊗_q L_q |G⟩` ([`Tableau::graph_form`]): a +//! graph state `|G⟩ = 2^{-n/2} Σ_z (−1)^{Σ_{a~b} z_a z_b} |z⟩` and local +//! Cliffords. With `W_q = V_q L_q` (`V_q` the qubit's last gate), +//! +//! ```text +//! ⟨x| ⊗ V_q |ψ⟩ = 2^{-n/2} Σ_z ∏_q W_q[x_q][z_q] · (−1)^{Σ_{a~b} z_a z_b}. +//! ``` +//! +//! The sum is taken one qubit at a time along an order. After the first k +//! qubits, the part summed depends on the rest only through the image of +//! the summed bits under the graph's edges between the two sides: a space +//! over GF(2) whose dimension is the cut-rank there, the entanglement (in +//! bits) of the state across that cut. A table of `2^{cut-rank}` complex +//! numbers holds it, and each qubit is one pass over the table: a merge +//! where the space loses a dimension, the qubit's two terms with their +//! signs, and a doubling (or an XOR-shifted sum) for its edges ahead. +//! +//! By Parseval, the probability of a prefix `x_1 … x_k` is `Σ |T|²` over +//! the table after it (the gates still to come are unitary and the graph's +//! signs have modulus 1), so the same passes sample exactly: each qubit's +//! conditional probabilities from one read of the table, then the pass for +//! the outcome drawn. +//! +//! Every reduction is summed in fixed chunks, in order, so a result is the +//! same bits at any thread count. + +use crate::gates::{C, complex}; +use crate::quantum_tableau::{GraphForm, Tableau}; + +/// A single-qubit gate, row-major: `[m00, m01, m10, m11]`. +pub type Unitary = [C; 4]; + +/// Why a computation was refused. +#[derive(Clone, Debug, PartialEq, Eq)] +pub enum CutRankError { + /// The gates given are not one per qubit. + Gates { qubits: usize, gates: usize }, + /// The order is not a permutation of the qubits. + Order, + /// A cut along the order is wider than the table may be. + TooWide { cut: usize, limit: usize }, + /// A bit string is not one bit per qubit. + Bits, +} + +impl std::fmt::Display for CutRankError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + match self { + CutRankError::Gates { qubits, gates } => write!(f, "{gates} final gates for {qubits} qubits"), + CutRankError::Order => write!(f, "the order is not a permutation of the qubits"), + CutRankError::TooWide { cut, limit } => write!(f, "a cut of {cut} bits along the order; at most {limit}"), + CutRankError::Bits => write!(f, "a bit string is not one bit per qubit"), + } + } +} + +impl std::error::Error for CutRankError {} + +fn mul(a: &Unitary, b: &Unitary) -> Unitary { + let m = |x: C, y: C| complex(x.re * y.re - x.im * y.im, x.re * y.im + x.im * y.re); + let add = |x: C, y: C| complex(x.re + y.re, x.im + y.im); + [add(m(a[0], b[0]), m(a[1], b[2])), add(m(a[0], b[1]), m(a[1], b[3])), add(m(a[2], b[0]), m(a[3], b[2])), add(m(a[2], b[1]), m(a[3], b[3]))] +} + +fn get(v: &[u64], i: usize) -> bool { + v[i / 64] >> (i % 64) & 1 == 1 +} + +/// A stabilizer state measured in a product basis, ready to be summed. +#[derive(Clone, Debug)] +pub struct ProductMeasured { + n: usize, + adjacency: Vec>, + /// `W_q / √2`. + w: Vec, + order: Vec, + /// The cut-rank after each qubit of the order. + cuts: Vec, +} + +impl ProductMeasured { + /// The state `tableau` holds with `last[q]` applied to each qubit q, to be + /// summed in the qubits' own order. + pub fn new(tableau: &Tableau, last: &[Unitary]) -> Result { + let form = tableau.graph_form(); + let n = form.n as usize; + if last.len() != n { + return Err(CutRankError::Gates { qubits: n, gates: last.len() }); + } + Ok(ProductMeasured::from_form(&form, last)) + } + + /// From a graph form already taken. + pub fn from_form(form: &GraphForm, last: &[Unitary]) -> ProductMeasured { + let n = form.n as usize; + let r = std::f64::consts::FRAC_1_SQRT_2; + let (o, one) = (complex(0.0, 0.0), complex(1.0, 0.0)); + let h = [complex(r, 0.0), complex(r, 0.0), complex(r, 0.0), complex(-r, 0.0)]; + let s = [one, o, o, complex(0.0, 1.0)]; + let z = [one, o, o, complex(-1.0, 0.0)]; + let w = (0..n) + .map(|q| { + let l = form.local[q]; + let mut m = [one, o, o, one]; + if l.h { + m = mul(&m, &h); + } + if l.s { + m = mul(&m, &s); + } + if l.z { + m = mul(&m, &z); + } + let wq = mul(&last[q], &m); + wq.map(|c| complex(c.re * r, c.im * r)) + }) + .collect(); + let mut pm = ProductMeasured { n, adjacency: form.adjacency.clone(), w, order: (0..n).collect(), cuts: Vec::new() }; + pm.cuts = pm.cut_ranks_of(&pm.order); + pm + } + + /// The same state, summed in `order`. + pub fn with_order(mut self, order: &[usize]) -> Result { + let mut seen = vec![false; self.n]; + if order.len() != self.n || order.iter().any(|&q| q >= self.n || std::mem::replace(&mut seen[q], true)) { + return Err(CutRankError::Order); + } + self.order = order.to_vec(); + self.cuts = self.cut_ranks_of(order); + Ok(self) + } + + pub fn qubits(&self) -> usize { + self.n + } + + /// The cut-rank after each qubit of the order: the table's size, in bits, + /// between that qubit and the next. + pub fn cut_ranks(&self) -> &[usize] { + &self.cuts + } + + /// The widest cut: `2^max_cut` entries is the largest table. + pub fn max_cut(&self) -> usize { + self.cuts.iter().copied().max().unwrap_or(0) + } + + fn cut_ranks_of(&self, order: &[usize]) -> Vec { + let words = self.n.div_ceil(64); + let mut ahead = vec![u64::MAX; words]; + if self.n % 64 != 0 { + ahead[words - 1] = (1u64 << (self.n % 64)) - 1; + } + let mut cuts = Vec::with_capacity(self.n); + let mut rows: Vec> = Vec::new(); + for &v in order { + ahead[v / 64] &= !(1u64 << (v % 64)); + rows.push(self.adjacency[v].clone()); + // The rank of the edges from the summed qubits to the rest. + let mut basis: Vec> = Vec::new(); + for row in &rows { + let mut x: Vec = row.iter().zip(&ahead).map(|(a, b)| a & b).collect(); + for b in &basis { + let p = pivot(b); + if get(&x, p) { + x.iter_mut().zip(b).for_each(|(a, b)| *a ^= b); + } + } + if x.iter().any(|&w| w != 0) { + basis.push(x); + } + } + cuts.push(basis.len()); + } + cuts + } + + /// `⟨x|ψ⟩` up to the state's global phase: `x[q]` is qubit q's bit. + pub fn amplitude(&self, x: &[bool], max_cut: usize) -> Result { + if x.len() != self.n { + return Err(CutRankError::Bits); + } + self.check(max_cut)?; + let mut walk = Walk::new(self); + for &v in &self.order { + walk.step(v, x[v], None); + } + Ok(walk.finish()) + } + + /// `|⟨x|ψ⟩|²`. + pub fn probability(&self, x: &[bool], max_cut: usize) -> Result { + let a = self.amplitude(x, max_cut)?; + Ok(a.re * a.re + a.im * a.im) + } + + /// `shots` exact samples, each with its probability, from a seeded stream: + /// qubit by qubit along the order, each outcome drawn from its + /// conditional probability. + pub fn sample(&self, shots: usize, seed: u64, max_cut: usize) -> Result, f64)>, CutRankError> { + self.check(max_cut)?; + let mut rng = seed; + let mut out = Vec::with_capacity(shots); + for _ in 0..shots { + let mut walk = Walk::new(self); + let mut x = vec![false; self.n]; + for &v in &self.order { + x[v] = walk.step(v, false, Some(&mut rng)); + } + let a = walk.finish(); + out.push((x, a.re * a.re + a.im * a.im)); + } + Ok(out) + } + + fn check(&self, max_cut: usize) -> Result<(), CutRankError> { + let cut = self.max_cut(); + if cut > max_cut { + return Err(CutRankError::TooWide { cut, limit: max_cut }); + } + Ok(()) + } +} + +fn pivot(v: &[u64]) -> usize { + let w = v.iter().position(|&x| x != 0).expect("a nonzero vector"); + w * 64 + v[w].trailing_zeros() as usize +} + +fn splitmix(s: &mut u64) -> u64 { + *s = s.wrapping_add(0x9e37_79b9_7f4a_7c15); + let mut z = *s; + z = (z ^ (z >> 30)).wrapping_mul(0xbf58_476d_1ce4_e5b9); + z = (z ^ (z >> 27)).wrapping_mul(0x94d0_49bb_1331_11eb); + z ^ (z >> 31) +} + +/// `(−1)^{|c ∧ ℓ|}`. +#[inline(always)] +fn sign(c: usize, l: usize) -> f64 { + if (c & l).count_ones() % 2 == 0 { 1.0 } else { -1.0 } +} + +/// Insert a zero bit at position `p` of `i`. +#[inline(always)] +fn insert_zero(i: usize, p: u32) -> usize { + let low = i & ((1usize << p) - 1); + ((i >> p) << (p + 1)) | low +} + +/// The table's two arrays, shared by the threads of a pass: each thread +/// writes only the indices of the groups it was given, and no two groups +/// share an index. +#[derive(Clone, Copy)] +struct Shared { + re: *mut f64, + im: *mut f64, +} + +// SAFETY: access is partitioned by group, as above. +unsafe impl Send for Shared {} +unsafe impl Sync for Shared {} + +impl Shared { + #[inline(always)] + fn get(self, i: usize) -> (f64, f64) { + // SAFETY: callers pass indices inside the table, of their own groups. + unsafe { (*self.re.add(i), *self.im.add(i)) } + } + #[inline(always)] + fn set(self, i: usize, v: (f64, f64)) { + // SAFETY: as for `get`. + unsafe { + *self.re.add(i) = v.0; + *self.im.add(i) = v.1; + } + } +} + +fn threads() -> usize { + #[cfg(not(target_arch = "wasm32"))] + { + std::thread::available_parallelism().map_or(1, |t| t.get()) + } + #[cfg(target_arch = "wasm32")] + { + 1 + } +} + +/// Run `f(lo, hi)` over `0..len` in chunks, on the host's threads. +fn par(len: usize, f: &(dyn Fn(usize, usize) + Sync)) { + let t = threads(); + if t == 1 || len < 1 << 15 || cfg!(target_arch = "wasm32") { + f(0, len); + return; + } + let chunk = len.div_ceil(t * 8).max(1 << 12); + let next = std::sync::atomic::AtomicUsize::new(0); + std::thread::scope(|sc| { + for _ in 0..t { + sc.spawn(|| loop { + let lo = next.fetch_add(chunk, std::sync::atomic::Ordering::Relaxed); + if lo >= len { + break; + } + f(lo, (lo + chunk).min(len)); + }); + } + }); +} + +/// Fixed chunks a reduction is summed in, so its bits do not depend on the +/// threads. +const CHUNKS: usize = 256; + +/// `Σ_{i [f64; 4] + Sync)) -> [f64; 4] { + let chunks = len.clamp(1, CHUNKS); + let per = len.div_ceil(chunks); + let slots: Vec> = (0..chunks).map(|_| std::sync::Mutex::new([0.0; 4])).collect(); + // `par` hands out chunk indices; each chunk is summed in order. + let run = |lo: usize, hi: usize| { + for c in lo..hi { + let mut acc = [0.0f64; 4]; + for i in c * per..((c + 1) * per).min(len) { + let t = term(i); + for k in 0..4 { + acc[k] += t[k]; + } + } + *slots[c].lock().unwrap() = acc; + } + }; + if len >= 1 << 15 { + let t = threads().min(chunks); + let next = std::sync::atomic::AtomicUsize::new(0); + std::thread::scope(|sc| { + for _ in 0..t { + sc.spawn(|| loop { + let c = next.fetch_add(1, std::sync::atomic::Ordering::Relaxed); + if c >= chunks { + break; + } + run(c, c + 1); + }); + } + }); + } else { + run(0, chunks); + } + slots.into_iter().map(|s| s.into_inner().unwrap()).fold([0.0; 4], |a, p| [a[0] + p[0], a[1] + p[1], a[2] + p[2], a[3] + p[3]]) +} + +struct Walk<'a> { + pm: &'a ProductMeasured, + /// The table's space: one vector a coordinate, over the qubits not yet summed. + basis: Vec>, + re: Vec, + im: Vec, + /// The qubits not yet summed. + ahead: Vec, +} + +impl<'a> Walk<'a> { + fn new(pm: &'a ProductMeasured) -> Walk<'a> { + let words = pm.n.div_ceil(64); + let mut ahead = vec![u64::MAX; words]; + if pm.n % 64 != 0 { + ahead[words - 1] = (1u64 << (pm.n % 64)) - 1; + } + let cap = 1usize << pm.max_cut(); + let (mut re, mut im) = (Vec::with_capacity(cap), Vec::with_capacity(cap)); + re.push(1.0); + im.push(0.0); + Walk { pm, basis: Vec::new(), re, im, ahead } + } + + fn shared(&mut self) -> Shared { + Shared { re: self.re.as_mut_ptr(), im: self.im.as_mut_ptr() } + } + + fn finish(self) -> C { + debug_assert_eq!(self.re.len(), 1); + complex(self.re[0], self.im[0]) + } + + /// Swap coordinates `a` and `b` (table, basis). + fn swap_coords(&mut self, a: usize, b: usize) { + if a == b { + return; + } + let (lo, hi) = (a.min(b) as u32, a.max(b) as u32); + self.basis.swap(a, b); + let len = self.re.len() >> 2; + let t = self.shared(); + par(len, &|s, e| { + for i in s..e { + // Bit lo = 1, bit hi = 0, against bit lo = 0, bit hi = 1. + let base = insert_zero(insert_zero(i, lo), hi); + let (p, q) = (base | 1 << lo, base | 1 << hi); + let (x, y) = (t.get(p), t.get(q)); + t.set(p, y); + t.set(q, x); + } + }); + } + + /// Sum qubit `v` with outcome `bit`; or, given `rng`, draw its outcome + /// from its conditional probability and return it. + fn step(&mut self, v: usize, bit: bool, rng: Option<&mut u64>) -> bool { + let w = self.pm.w[v]; + let mut ell = 0usize; + for (j, b) in self.basis.iter().enumerate() { + if get(b, v) { + ell |= 1 << j; + } + } + self.ahead[v / 64] &= !(1u64 << (v % 64)); + for b in &mut self.basis { + b[v / 64] &= !(1u64 << (v % 64)); + } + // A combination of the projected basis that vanishes (at most one). + let kernel = dependency(&self.basis); + let merged = if let Some(mut kappa) = kernel { + let top = self.basis.len() - 1; + let j0 = usize::BITS as usize - 1 - kappa.leading_zeros() as usize; + if j0 != top { + self.swap_coords(j0, top); + let swap = |m: usize| { + let (x, y) = (m >> j0 & 1, m >> top & 1); + (m & !(1 << j0) & !(1 << top)) | (y << j0) | (x << top) + }; + ell = swap(ell); + kappa = swap(kappa); + } + self.merge(kappa & !(1 << top), ell); + self.basis.pop(); + true + } else { + false + }; + // The qubit's edges ahead, in the table's coordinates. + let a: Vec = self.pm.adjacency[v].iter().zip(&self.ahead).map(|(x, y)| x & y).collect(); + let alpha = coordinates(&self.basis, &a); + // The outcome: given, or drawn from the conditional probabilities. + let chosen = match rng { + None => bit, + Some(rng) => { + let (sa, sb, x) = self.overlap(merged, alpha, ell); + // Σ|w0·A + w1·B|² = |w0|²Σ|A|² + |w1|²Σ|B|² + 2·Re(w0·conj(w1)·ΣA·conj(B)). + let norm = |b: bool| { + let (w0, w1) = (w[2 * usize::from(b)], w[2 * usize::from(b) + 1]); + let (cr, ci) = (w0.re * w1.re + w0.im * w1.im, w0.im * w1.re - w0.re * w1.im); + (w0.re * w0.re + w0.im * w0.im) * sa + (w1.re * w1.re + w1.im * w1.im) * sb + 2.0 * (cr * x.0 - ci * x.1) + }; + let (p0, p1) = (norm(false).max(0.0), norm(true).max(0.0)); + let u = ((splitmix(rng) >> 11) as f64 + 0.5) / 9_007_199_254_740_992.0; + u * (p0 + p1) >= p0 + } + }; + let (w0, w1) = (w[2 * usize::from(chosen)], w[2 * usize::from(chosen) + 1]); + self.apply(merged, alpha, ell, w0, w1); + if alpha.is_none() { + self.basis.push(a); + } + chosen + } + + /// Merge along the kernel: the table's top coordinate goes, and its two + /// halves become the qubit's z = 0 and z = 1 sums (the latter signed by + /// `ell`), each indexed by the remaining coordinates. + fn merge(&mut self, kappa_low: usize, ell: usize) { + let top = self.re.len() >> 1; + let t = self.shared(); + if kappa_low == 0 { + par(top, &|s, e| { + for c in s..e { + let (a, b) = (t.get(c), t.get(c | top)); + let (sa, sb) = (sign(c, ell), sign(c | top, ell)); + t.set(c, (a.0 + b.0, a.1 + b.1)); + t.set(c | top, (sa * a.0 + sb * b.0, sa * a.1 + sb * b.1)); + } + }); + } else { + let hk = usize::BITS - 1 - kappa_low.leading_zeros(); + par(top >> 1, &|s, e| { + for i in s..e { + let c = insert_zero(i, hk); + let c2 = c ^ kappa_low; + let (a, b, aa, bb) = (t.get(c), t.get(c2), t.get(c | top), t.get(c2 | top)); + // c pairs with c2 | top, and c2 with c | top. + t.set(c, (a.0 + bb.0, a.1 + bb.1)); + t.set(c2, (b.0 + aa.0, b.1 + aa.1)); + let (s1, s2, s3, s4) = (sign(c, ell), sign(c2 | top, ell), sign(c2, ell), sign(c | top, ell)); + t.set(c | top, (s1 * a.0 + s2 * bb.0, s1 * a.1 + s2 * bb.1)); + t.set(c2 | top, (s3 * b.0 + s4 * aa.0, s3 * b.1 + s4 * aa.1)); + } + }); + } + } + + /// `(Σ|A|², Σ|B|², Σ A·conj(B))` over the outputs, where an output is + /// `w0·A + w1·B`: what the conditional probabilities need. + fn overlap(&self, merged: bool, alpha: Option, ell: usize) -> (f64, f64, (f64, f64)) { + let len = self.re.len(); + let (re, im) = (&self.re, &self.im); + let m2 = |i: usize| re[i] * re[i] + im[i] * im[i]; + // conj-product A·conj(B) of entries i (A) and j (B), times s. + let cp = |i: usize, j: usize, s: f64| [s * (re[i] * re[j] + im[i] * im[j]), s * (im[i] * re[j] - re[i] * im[j])]; + match (merged, alpha) { + // A and B on separate outputs. + (true, None) => { + let half = len >> 1; + let [a, b, _, _] = reduce4(half, &|c| [m2(c), m2(c | half), 0.0, 0.0]); + (a, b, (0.0, 0.0)) + } + (false, None) => { + let [s, _, _, _] = reduce4(len, &|c| [m2(c), 0.0, 0.0, 0.0]); + (s, s, (0.0, 0.0)) + } + // A = lower[c], B = upper[c ⊕ α]. + (true, Some(al)) => { + let half = len >> 1; + let [a, b, xr, xi] = reduce4(half, &|c| { + let d = (c ^ al) | half; + let [xr, xi] = cp(c, d, 1.0); + [m2(c), m2(d), xr, xi] + }); + (a, b, (xr, xi)) + } + // A = T[c], B = sign(c ⊕ α)·T[c ⊕ α]. + (false, Some(al)) => { + let [s, xr, xi, _] = reduce4(len, &|c| { + let d = c ^ al; + let [xr, xi] = cp(c, d, sign(d, ell)); + [m2(c), xr, xi, 0.0] + }); + (s, s, (xr, xi)) + } + } + } + + /// The qubit's pass: outputs `w0·A + w1·B` (see [`Walk::overlap`]). + fn apply(&mut self, merged: bool, alpha: Option, ell: usize, w0: C, w1: C) { + let len = self.re.len(); + let cm = |w: C, v: (f64, f64)| (w.re * v.0 - w.im * v.1, w.re * v.1 + w.im * v.0); + let add = |a: (f64, f64), b: (f64, f64)| (a.0 + b.0, a.1 + b.1); + match (merged, alpha) { + (true, None) => { + let half = len >> 1; + let t = self.shared(); + par(half, &|s, e| { + for c in s..e { + t.set(c, cm(w0, t.get(c))); + t.set(c | half, cm(w1, t.get(c | half))); + } + }); + } + (true, Some(al)) => { + let half = len >> 1; + let t = self.shared(); + par(half, &|s, e| { + for c in s..e { + t.set(c, add(cm(w0, t.get(c)), cm(w1, t.get((c ^ al) | half)))); + } + }); + self.re.truncate(half); + self.im.truncate(half); + } + (false, None) => { + self.re.resize(2 * len, 0.0); + self.im.resize(2 * len, 0.0); + let t = self.shared(); + par(len, &|s, e| { + for c in s..e { + let v = t.get(c); + let sv = sign(c, ell); + t.set(c | len, cm(w1, (sv * v.0, sv * v.1))); + t.set(c, cm(w0, v)); + } + }); + } + (false, Some(0)) => { + let t = self.shared(); + par(len, &|s, e| { + for c in s..e { + let v = t.get(c); + let sv = sign(c, ell); + t.set(c, add(cm(w0, v), cm(w1, (sv * v.0, sv * v.1)))); + } + }); + } + (false, Some(al)) => { + let hb = usize::BITS - 1 - al.leading_zeros(); + let t = self.shared(); + par(len >> 1, &|s, e| { + for i in s..e { + let c = insert_zero(i, hb); + let d = c ^ al; + let (x, y) = (t.get(c), t.get(d)); + let (sx, sy) = (sign(c, ell), sign(d, ell)); + t.set(c, add(cm(w0, x), cm(w1, (sy * y.0, sy * y.1)))); + t.set(d, add(cm(w0, y), cm(w1, (sx * x.0, sx * x.1)))); + } + }); + } + } + } +} + +/// A nonzero combination (as a coordinate mask) of `vectors` that sums to +/// zero, if there is one. +fn dependency(vectors: &[Vec]) -> Option { + let mut rows: Vec<(Vec, usize)> = Vec::new(); + for (j, v) in vectors.iter().enumerate() { + let mut x = v.clone(); + let mut combo = 1usize << j; + for (b, bc) in &rows { + if get(&x, pivot(b)) { + x.iter_mut().zip(b).for_each(|(a, b)| *a ^= b); + combo ^= bc; + } + } + if x.iter().all(|&w| w == 0) { + return Some(combo); + } + rows.push((x, combo)); + } + None +} + +/// `a`'s coordinates in `basis` (independent vectors), or `None` if `a` is +/// outside their span. +fn coordinates(basis: &[Vec], a: &[u64]) -> Option { + let mut rows: Vec<(Vec, usize)> = Vec::new(); + for (j, v) in basis.iter().enumerate() { + let mut x = v.clone(); + let mut combo = 1usize << j; + for (b, bc) in &rows { + if get(&x, pivot(b)) { + x.iter_mut().zip(b).for_each(|(a, b)| *a ^= b); + combo ^= bc; + } + } + debug_assert!(x.iter().any(|&w| w != 0), "the basis is independent"); + rows.push((x, combo)); + } + let mut x = a.to_vec(); + let mut combo = 0usize; + for (b, bc) in &rows { + if get(&x, pivot(b)) { + x.iter_mut().zip(b).for_each(|(a, b)| *a ^= b); + combo ^= bc; + } + } + x.iter().all(|&w| w == 0).then_some(combo) +} + +#[cfg(all(test, feature = "quantum_sv"))] +mod tests { + use super::*; + use crate::gates::Gate; + use crate::quantum_sv::StateVector; + use crate::quantum_tableau::Clifford; + + fn next(s: &mut u64) -> u64 { + splitmix(s) + } + + fn random_clifford(n: u32, ops: usize, seed: u64) -> Vec { + let mut s = seed; + (0..ops) + .map(|_| { + let a = (next(&mut s) % n as u64) as u32; + let mut b = (next(&mut s) % n as u64) as u32; + if b == a { + b = (a + 1) % n; + } + match next(&mut s) % 9 { + 0 => Clifford::H(a), + 1 => Clifford::S(a), + 2 => Clifford::Sdg(a), + 3 => Clifford::X(a), + 4 => Clifford::Y(a), + 5 => Clifford::SX(a), + 6 => Clifford::CX(a, b), + 7 => Clifford::CZ(a, b), + _ => Clifford::Z(a), + } + }) + .collect() + } + + fn matrix(op: Clifford) -> Gate { + let c = complex; + match op { + Clifford::H(a) => Gate::h(a), + Clifford::S(a) => Gate::s(a), + Clifford::Sdg(a) => Gate::sdg(a), + Clifford::X(a) => Gate::x(a), + Clifford::Y(a) => Gate::y(a), + Clifford::Z(a) => Gate::z(a), + Clifford::SX(a) => Gate { qubits: vec![a], matrix: vec![c(0.5, 0.5), c(0.5, -0.5), c(0.5, -0.5), c(0.5, 0.5)] }, + Clifford::CX(a, b) => Gate::cx(a, b), + Clifford::CZ(a, b) => Gate::cz(a, b), + other => unreachable!("{other:?}"), + } + } + + /// A random single-qubit unitary, `U(θ, φ, λ)`. + fn random_unitary(s: &mut u64) -> Unitary { + let mut ang = || (next(s) >> 11) as f64 / (1u64 << 53) as f64 * std::f64::consts::TAU; + let (t, p, l) = (ang(), ang(), ang()); + let (c, sn) = ((t / 2.0).cos(), (t / 2.0).sin()); + let e = |x: f64| complex(x.cos(), x.sin()); + let sc = |z: C, k: f64| complex(z.re * k, z.im * k); + [complex(c, 0.0), sc(e(l), -sn), sc(e(p), sn), sc(e(p + l), c)] + } + + /// The circuit and the last gates on the dense state vector. + fn dense(n: u32, ops: &[Clifford], last: &[Unitary]) -> StateVector { + let mut sv = StateVector::zero(n).unwrap(); + for &o in ops { + sv.apply(&matrix(o)).unwrap(); + } + for (q, u) in last.iter().enumerate() { + sv.apply(&Gate { qubits: vec![q as u32], matrix: u.to_vec() }).unwrap(); + } + sv + } + + fn bits(x: usize, n: usize) -> Vec { + // The state vector's index has qubit 0 as its most significant bit. + (0..n).map(|q| x >> (n - 1 - q) & 1 == 1).collect() + } + + #[test] + fn the_graph_form_rebuilds_the_state() { + for seed in 0..40u64 { + let n = 2 + (seed % 9) as u32; + let ops = random_clifford(n, 40, seed); + let t = Tableau::run(n, &ops).unwrap(); + let f = t.graph_form(); + // |G⟩, then each L_q = H^h S^s Z^z (Z first). + let mut g = StateVector::zero(n).unwrap(); + for q in 0..n { + g.apply(&Gate::h(q)).unwrap(); + } + for a in 0..n as usize { + for b in a + 1..n as usize { + if get(&f.adjacency[a], b) { + g.apply(&Gate::cz(a as u32, b as u32)).unwrap(); + } + } + } + for (q, l) in f.local.iter().enumerate() { + let q = q as u32; + if l.z { + g.apply(&Gate::z(q)).unwrap(); + } + if l.s { + g.apply(&Gate::s(q)).unwrap(); + } + if l.h { + g.apply(&Gate::h(q)).unwrap(); + } + } + let sv = dense(n, &ops, &[]); + let (mut re, mut im) = (0.0, 0.0); + for i in 0..sv.re.len() { + re += g.re[i] * sv.re[i] + g.im[i] * sv.im[i]; + im += g.re[i] * sv.im[i] - g.im[i] * sv.re[i]; + } + assert!(((re * re + im * im).sqrt() - 1.0).abs() < 1e-9, "seed {seed}: |⟨G|ψ⟩| = {}", (re * re + im * im).sqrt()); + } + } + + #[test] + fn amplitudes_are_the_state_vectors() { + for seed in 0..30u64 { + let n = 1 + (seed % 12) as usize; + let mut s = seed ^ 0xa5a5; + let ops = if n > 1 { random_clifford(n as u32, 8 * n, seed) } else { vec![Clifford::H(0)] }; + let last: Vec = (0..n).map(|_| random_unitary(&mut s)).collect(); + let pm = ProductMeasured::new(&Tableau::run(n as u32, &ops).unwrap(), &last).unwrap(); + let sv = dense(n as u32, &ops, &last); + // One global phase relates the two: fix it on the largest amplitude. + let x0 = (0..sv.re.len()).max_by(|&a, &b| sv.amplitude(a).re.hypot(sv.amplitude(a).im).total_cmp(&sv.amplitude(b).re.hypot(sv.amplitude(b).im))).unwrap(); + let (a0, b0) = (pm.amplitude(&bits(x0, n), 40).unwrap(), sv.amplitude(x0)); + for x in 0..sv.re.len() { + let (a, b) = (pm.amplitude(&bits(x, n), 40).unwrap(), sv.amplitude(x)); + assert!(((a.re * a.re + a.im * a.im) - (b.re * b.re + b.im * b.im)).abs() < 1e-12, "seed {seed} x {x}: probability"); + // a/a0 = b/b0: a·b0 = b·a0. + let l = (a.re * b0.re - a.im * b0.im, a.re * b0.im + a.im * b0.re); + let r = (b.re * a0.re - b.im * a0.im, b.re * a0.im + b.im * a0.re); + assert!((l.0 - r.0).abs() < 1e-12 && (l.1 - r.1).abs() < 1e-12, "seed {seed} x {x}: relative phase"); + } + } + } + + #[test] + fn every_order_gives_the_same_amplitude() { + let mut compared = 0; + for seed in 0..6u64 { + let n = 24 + 6 * seed as usize; + let mut s = seed ^ 0x77; + let ops = random_clifford(n as u32, 3 * n, seed); + let last: Vec = (0..n).map(|_| random_unitary(&mut s)).collect(); + let t = Tableau::run(n as u32, &ops).unwrap(); + let fwd = ProductMeasured::new(&t, &last).unwrap(); + let rev = fwd.clone().with_order(&(0..n).rev().collect::>()).unwrap(); + let mut perm: Vec = (0..n).collect(); + for i in (1..n).rev() { + perm.swap(i, (next(&mut s) % (i as u64 + 1)) as usize); + } + let shuffled = fwd.clone().with_order(&perm).unwrap(); + for _ in 0..4 { + let x: Vec = (0..n).map(|_| next(&mut s) & 1 == 1).collect(); + let a = fwd.amplitude(&x, 24).unwrap(); + for other in [&rev, &shuffled] { + if other.max_cut() > 22 { + continue; + } + let b = other.amplitude(&x, 24).unwrap(); + let scale = a.re.hypot(a.im).max(1e-300); + assert!((a.re - b.re).abs() / scale < 1e-9 && (a.im - b.im).abs() / scale < 1e-9, "seed {seed}: {a:?} vs {b:?}"); + compared += 1; + } + } + } + assert!(compared >= 30, "only {compared} comparisons ran"); + } + + #[test] + fn samples_follow_the_probabilities() { + let n = 6usize; + let mut s = 99u64; + let ops = random_clifford(n as u32, 40, 7); + let last: Vec = (0..n).map(|_| random_unitary(&mut s)).collect(); + let pm = ProductMeasured::new(&Tableau::run(n as u32, &ops).unwrap(), &last).unwrap(); + let sv = dense(n as u32, &ops, &last); + let shots = 60_000; + let samples = pm.sample(shots, 5, 10).unwrap(); + let mut counts = vec![0usize; 1 << n]; + for (x, p) in &samples { + let idx = x.iter().enumerate().fold(0usize, |i, (q, &b)| i | (usize::from(b) << (n - 1 - q))); + counts[idx] += 1; + let want = sv.amplitude(idx); + assert!((p - (want.re * want.re + want.im * want.im)).abs() < 1e-12, "a sample's probability is its own"); + } + for (i, &c) in counts.iter().enumerate() { + let a = sv.amplitude(i); + let p = a.re * a.re + a.im * a.im; + let sd = (shots as f64 * p * (1.0 - p)).sqrt().max(1.0); + assert!((c as f64 - shots as f64 * p).abs() < 5.0 * sd, "{i:06b}: {c} vs {:.1}", shots as f64 * p); + } + // Deterministic in the seed. + assert_eq!(pm.sample(50, 5, 10).unwrap(), pm.sample(50, 5, 10).unwrap()); + } + + #[test] + fn cut_ranks_are_the_entanglement() { + // A GHZ chain: one bit across every cut. A product state: none. + let n = 9u32; + let mut ghz = vec![Clifford::H(0)]; + ghz.extend((0..n - 1).map(|q| Clifford::CX(q, q + 1))); + let id = [complex(1.0, 0.0), complex(0.0, 0.0), complex(0.0, 0.0), complex(1.0, 0.0)]; + let pm = ProductMeasured::new(&Tableau::run(n, &ghz).unwrap(), &vec![id; n as usize]).unwrap(); + assert!(pm.cut_ranks()[..n as usize - 1].iter().all(|&c| c == 1), "{:?}", pm.cut_ranks()); + let prod = ProductMeasured::new(&Tableau::run(n, &[Clifford::H(3)]).unwrap(), &vec![id; n as usize]).unwrap(); + assert_eq!(prod.max_cut(), 0); + let p = prod.probability(&(0..n).map(|q| q == 3).collect::>(), 4).unwrap(); + assert!((p - 0.5).abs() < 1e-12); + assert!(matches!(pm.amplitude(&[false; 9], 0), Err(CutRankError::TooWide { cut: 1, limit: 0 }))); + } +} diff --git a/wai/wai-quantum/src/quantum_tableau.rs b/wai/wai-quantum/src/quantum_tableau.rs index 66e9a7bab..34769bd1a 100644 --- a/wai/wai-quantum/src/quantum_tableau.rs +++ b/wai/wai-quantum/src/quantum_tableau.rs @@ -118,6 +118,29 @@ fn set(v: &mut [u64], i: usize) { v[i / 64] |= 1 << (i % 64); } +fn put(v: &mut [u64], i: usize, bit: bool) { + v[i / 64] = (v[i / 64] & !(1u64 << (i % 64))) | (u64::from(bit) << (i % 64)); +} + +/// One qubit's local Clifford in a [`GraphForm`]: `H^h · S^s · Z^z`. +#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)] +pub struct LocalClifford { + pub h: bool, + pub s: bool, + pub z: bool, +} + +/// A stabilizer state as `⊗_q L_q |G⟩` up to a global phase, where `|G⟩ = +/// ∏_{a~b} CZ_ab |+⟩^n` is the graph state of `adjacency` and `L_q` is +/// qubit q's [`LocalClifford`]. +#[derive(Clone, Debug, PartialEq, Eq)] +pub struct GraphForm { + pub n: u32, + /// Row q: q's neighbours, packed (qubit j at bit `j % 64` of word `j / 64`). + pub adjacency: Vec>, + pub local: Vec, +} + /// The first set bit of `v` in `from..to`. fn first_set(v: &[u64], from: usize, to: usize) -> Option { let mut w = from / 64; @@ -355,6 +378,77 @@ impl Tableau { Ok(()) } + /// The state as a graph state and local Cliffords ([`GraphForm`]). + /// + /// The stabilizer generators are row-reduced on their X part; the Z-only + /// rows left over name, by their pivots, the qubits a Hadamard turns into + /// X pivots; the X part is then reduced to the identity, so row q is + /// `±X_q` (or `±Y_q`) times Z on q's neighbours. A Y on its own qubit + /// takes an S, a minus sign a Z. + pub fn graph_form(&self) -> GraphForm { + let n = self.n; + let mut s = self.rows(n, 2 * n); + let wn = s.wn; + // The X part, reduced. + let mut k = 0; + for q in 0..n { + let Some(p) = (k..n).find(|&i| get(s.xr(i), q)) else { continue }; + s.swap(k, p); + for i in 0..n { + if i != k && get(s.xr(i), q) { + s.mul_into(i, k); + } + } + k += 1; + } + // The Z-only rows, reduced: their pivot qubits take a Hadamard. + let mut h = vec![false; n]; + let mut k2 = k; + for q in 0..n { + let Some(p) = (k2..n).find(|&i| get(s.zr(i), q)) else { continue }; + s.swap(k2, p); + for i in k..n { + if i != k2 && get(s.zr(i), q) { + s.mul_into(i, k2); + } + } + h[q] = true; + k2 += 1; + } + debug_assert_eq!(k2, n, "n independent generators"); + // H on those qubits, in every row: X ↔ Z, and Y → −Y. + for q in (0..n).filter(|&q| h[q]) { + for i in 0..n { + let (x, z) = (get(s.xr(i), q), get(s.zr(i), q)); + s.r[i] ^= x && z; + put(&mut s.x[i * wn..(i + 1) * wn], q, z); + put(&mut s.z[i * wn..(i + 1) * wn], q, x); + } + } + // The X part to the identity. + for q in 0..n { + let p = (q..n).find(|&i| get(s.xr(i), q)).expect("after the Hadamards the X part has full rank"); + s.swap(q, p); + for i in 0..n { + if i != q && get(s.xr(i), q) { + s.mul_into(i, q); + } + } + } + // Row q holds X_q or Y_q on its own qubit, and Z on the others: S† + // turns a Y into an X (only row q has an X part on q, and S† leaves + // Z alone), and Z on q flips row q's sign alone. + let mut local = vec![LocalClifford::default(); n]; + let mut adjacency = vec![vec![0u64; wn]; n]; + for q in 0..n { + local[q] = LocalClifford { h: h[q], s: get(s.zr(q), q), z: s.r[q] }; + adjacency[q].copy_from_slice(s.zr(q)); + put(&mut adjacency[q], q, false); + } + debug_assert!((0..n).all(|a| (0..n).all(|b| get(&adjacency[a], b) == get(&adjacency[b], a))), "the graph is undirected"); + GraphForm { n: n as u32, adjacency, local } + } + /// Whether a measurement of qubit `a` is random (each outcome ½) rather /// than determined: some stabilizer anticommutes with `Z_a`. One scan of /// a column, so far cheaper than reading the determined outcome. -- 2.54.0 (Apple Git-157) From 70af6b4dcd7b7b812b09e3a8fc85d2eb3955dd30 Mon Sep 17 00:00:00 2001 From: Transaction Science Date: Mon, 19 Oct 2026 07:11:13 -0400 Subject: [PATCH 2/5] wai-quantum cutrank: the steps planned once, tables stored at a chosen precision, and orders annealed for less work MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - The table operations (coordinate swaps, kernel merges, doublings, XOR sums, the sign masks) depend on the graph and the order alone, so they are planned once per state; an amplitude or a sample only runs passes. - Precision::Single stores the table in f32 (arithmetic stays f64): half the memory and traffic, amplitudes within 1e-5 of the double ones (test). - A merge followed by a doubling is one pass; a doubling no longer zero-fills. Every pass is cut into fixed chunks, so reductions do not depend on the threads. - work() = Σ 2^cut, passes() tallies the passes by kind, and anneal_order(moves, seed) searches orders by simulated annealing, never returning a worse one (test: lower work, same amplitudes). On the tracker's 70-qubit random graph-state instance (build host): natural order, max cut 32: amplitude 6.2 s (f64) / 4.7 s (f32), sample 7.5 s (f32). Annealed order (seed 3, 20 000 moves): max cut 29, work 6.5× less, largest table 4 GiB in f32: amplitude 0.53 s, exact sample 1.49 s, amplitudes agreeing with the natural order's double ones to 1e-7. Co-Authored-By: Claude Opus 5.5 --- wai/wai-quantum/examples/cutrank_qasm.rs | 36 +- wai/wai-quantum/src/quantum_cutrank.rs | 826 +++++++++++++++-------- 2 files changed, 575 insertions(+), 287 deletions(-) diff --git a/wai/wai-quantum/examples/cutrank_qasm.rs b/wai/wai-quantum/examples/cutrank_qasm.rs index a3e45e88c..7370ce9d0 100644 --- a/wai/wai-quantum/examples/cutrank_qasm.rs +++ b/wai/wai-quantum/examples/cutrank_qasm.rs @@ -15,7 +15,7 @@ use std::time::Instant; use wai_quantum::gates::{C, complex}; -use wai_quantum::quantum_cutrank::{ProductMeasured, Unitary}; +use wai_quantum::quantum_cutrank::{Precision, ProductMeasured, Unitary}; use wai_quantum::quantum_tableau::{Clifford, Tableau}; fn mul(a: &Unitary, b: &Unitary) -> Unitary { @@ -107,14 +107,37 @@ fn main() { }); } let t0 = Instant::now(); - let pm = ProductMeasured::new(&Tableau::run(n, &bulk).unwrap(), &last).unwrap(); + let mut pm = ProductMeasured::new(&Tableau::run(n, &bulk).unwrap(), &last).unwrap(); + // CUTRANK_ORDER=FILE: an order (as `order` prints it) to sum in. + if let Ok(path) = std::env::var("CUTRANK_ORDER") { + let text = std::fs::read_to_string(&path).expect("the order file reads"); + let line = text.lines().find(|l| l.starts_with("order: ")).expect("an `order: [...]` line"); + let order: Vec = line["order: ".len()..].trim_matches(|c| c == '[' || c == ']').split(',').map(|x| x.trim().parse().unwrap()).collect(); + pm = pm.with_order(&order).expect("a permutation of the qubits"); + } println!("{n} qubits, {} bulk gates; graph form and cuts in {:?}", bulk.len(), t0.elapsed()); println!("cut-ranks: {:?}", pm.cut_ranks()); println!("max cut {} bits: a table of {} entries ({} GiB at f64)", pm.max_cut(), 1u64 << pm.max_cut(), (16u64 << pm.max_cut()) >> 30); let limit = std::env::var("CUTRANK_MAX").ok().and_then(|s| s.parse().ok()).unwrap_or(32); + let precision = if std::env::var("CUTRANK_PRECISION").as_deref() == Ok("single") { Precision::Single } else { Precision::Double }; + println!("tables at {precision:?}: the largest {} GiB", pm.table_bytes(precision) >> 30); + for (kind, passes, entries) in pm.passes() { + println!(" {kind:8} {passes:3} passes, {:.2}e9 entries", entries as f64 / 1e9); + } let scale = (n as f64).exp2(); + // `order N`: anneal the order for the least work, Σ 2^cut. + let cost = |c: &[usize]| c.iter().map(|&k| (k as f64).exp2()).sum::(); match mode { "cuts" => {} + "order" => { + let t0 = Instant::now(); + let best_order = pm.anneal_order(count, seed); + let best = cost(pm.clone().with_order(&best_order).unwrap().cut_ranks()); + let b = pm.clone().with_order(&best_order).unwrap(); + println!("annealed {count} moves in {:?}: work Σ2^cut {best:.3e} (natural {:.3e}), max cut {}", t0.elapsed(), cost(pm.cut_ranks()), b.max_cut()); + println!("cuts: {:?}", b.cut_ranks()); + println!("order: {best_order:?}"); + } "amplitudes" => { let mut s = seed; for i in 0..count { @@ -125,15 +148,20 @@ fn main() { (z ^ (z >> 27)) & 1 == 1 }).collect(); let t = Instant::now(); - let a = pm.amplitude(&x, limit).unwrap(); + let a = pm.amplitude_in(&x, limit, precision).unwrap(); println!("amplitude {i}: {a:?}, 2^n·p = {:.6}, in {:?}", scale * (a.re * a.re + a.im * a.im), t.elapsed()); } } "sample" => { let t = Instant::now(); - let samples = pm.sample(count, seed, limit).unwrap(); + let samples = pm.sample_in(count, seed, limit, precision).unwrap(); let el = t.elapsed(); let mean = samples.iter().map(|(_, p)| scale * p).sum::() / samples.len() as f64; + // CUTRANK_OUT=FILE: every sample, `bits 2^n·p` a line (qubit 0 first). + if let Ok(out) = std::env::var("CUTRANK_OUT") { + let text: String = samples.iter().map(|(x, p)| format!("{} {:.9}\n", x.iter().map(|&b| if b { '1' } else { '0' }).collect::(), scale * p)).collect(); + std::fs::write(&out, text).expect("the samples write"); + } for (x, p) in samples.iter().take(8) { let bits: String = x.iter().map(|&b| if b { '1' } else { '0' }).collect(); println!("{bits} 2^n·p = {:.6}", scale * p); diff --git a/wai/wai-quantum/src/quantum_cutrank.rs b/wai/wai-quantum/src/quantum_cutrank.rs index fe640365e..6462ff4c5 100644 --- a/wai/wai-quantum/src/quantum_cutrank.rs +++ b/wai/wai-quantum/src/quantum_cutrank.rs @@ -81,6 +81,8 @@ pub struct ProductMeasured { order: Vec, /// The cut-rank after each qubit of the order. cuts: Vec, + /// The steps along the order, the same for every bit string. + steps: Vec, } impl ProductMeasured { @@ -120,8 +122,9 @@ impl ProductMeasured { wq.map(|c| complex(c.re * r, c.im * r)) }) .collect(); - let mut pm = ProductMeasured { n, adjacency: form.adjacency.clone(), w, order: (0..n).collect(), cuts: Vec::new() }; + let mut pm = ProductMeasured { n, adjacency: form.adjacency.clone(), w, order: (0..n).collect(), cuts: Vec::new(), steps: Vec::new() }; pm.cuts = pm.cut_ranks_of(&pm.order); + pm.steps = plan(&pm.adjacency, n, &pm.order); pm } @@ -133,6 +136,7 @@ impl ProductMeasured { } self.order = order.to_vec(); self.cuts = self.cut_ranks_of(order); + self.steps = plan(&self.adjacency, self.n, order); Ok(self) } @@ -181,17 +185,106 @@ impl ProductMeasured { cuts } + /// The work of an order: `Σ_k 2^{cut_k}`, the table entries its passes + /// touch, up to a small factor. + pub fn work(&self) -> f64 { + self.cuts.iter().map(|&k| (k as f64).exp2()).sum() + } + + /// An order with less work, by simulated annealing from the current one: + /// `moves` random swaps and moves of one qubit, seeded, accepted by the + /// Metropolis rule on the work's logarithm as the temperature falls. The + /// cheapest order seen is returned; it is never worse than the current. + pub fn anneal_order(&self, moves: usize, seed: u64) -> Vec { + let n = self.n; + let mut s = seed; + let cost = |o: &[usize]| self.cut_ranks_of(o).iter().map(|&k| (k as f64).exp2()).sum::(); + let mut order = self.order.clone(); + let mut cur = self.work(); + let (mut best, mut best_order) = (cur, order.clone()); + if n < 2 { + return best_order; + } + for it in 0..moves { + let temp = 0.3 * (1.0 - it as f64 / moves as f64) + 1e-3; + let (i, j) = ((splitmix(&mut s) % n as u64) as usize, (splitmix(&mut s) % n as u64) as usize); + let mut cand = order.clone(); + if splitmix(&mut s) & 1 == 0 { + cand.swap(i, j); + } else { + let v = cand.remove(i); + cand.insert(j, v); + } + let c = cost(&cand); + // The acceptance draw only for a worse order. + if c <= cur || ((cur.log2() - c.log2()) / temp).exp() > (splitmix(&mut s) >> 11) as f64 / (1u64 << 53) as f64 { + order = cand; + cur = c; + if c < best { + best = c; + best_order = order.clone(); + } + } + } + best_order + } + + /// The passes one amplitude makes over the table, by kind: + /// `(kind, passes, entries visited)`. + pub fn passes(&self) -> Vec<(&'static str, usize, u64)> { + let mut out = vec![("swap", 0usize, 0u64), ("merge", 0, 0), ("fold", 0, 0), ("double", 0, 0), ("xor-sum", 0, 0)]; + let mut r = 0usize; + for st in &self.steps { + let len = 1u64 << r; + if st.swap.is_some() { + out[0].1 += 1; + out[0].2 += len; + } + if st.kernel.is_some() { + out[1].1 += 1; + out[1].2 += len; + r -= 1; + } + match st.alpha { + None if st.kernel.is_some() => r += 1, + None => { + out[3].1 += 1; + out[3].2 += 2 * len; + r += 1; + } + Some(_) if st.kernel.is_some() => { + out[2].1 += 1; + out[2].2 += len; + } + Some(_) => { + out[4].1 += 1; + out[4].2 += len; + } + } + } + out + } + + /// The largest table, in bytes, at `precision`. + pub fn table_bytes(&self, precision: Precision) -> u64 { + precision.entry_bytes() << self.max_cut() + } + /// `⟨x|ψ⟩` up to the state's global phase: `x[q]` is qubit q's bit. pub fn amplitude(&self, x: &[bool], max_cut: usize) -> Result { + self.amplitude_in(x, max_cut, Precision::Double) + } + + /// [`ProductMeasured::amplitude`], the table stored at `precision`. + pub fn amplitude_in(&self, x: &[bool], max_cut: usize, precision: Precision) -> Result { if x.len() != self.n { return Err(CutRankError::Bits); } self.check(max_cut)?; - let mut walk = Walk::new(self); - for &v in &self.order { - walk.step(v, x[v], None); - } - Ok(walk.finish()) + Ok(match precision { + Precision::Double => self.run::(Some(x), &mut 0).1, + Precision::Single => self.run::(Some(x), &mut 0).1, + }) } /// `|⟨x|ψ⟩|²`. @@ -204,19 +297,35 @@ impl ProductMeasured { /// qubit by qubit along the order, each outcome drawn from its /// conditional probability. pub fn sample(&self, shots: usize, seed: u64, max_cut: usize) -> Result, f64)>, CutRankError> { + self.sample_in(shots, seed, max_cut, Precision::Double) + } + + /// [`ProductMeasured::sample`], the table stored at `precision`. + pub fn sample_in(&self, shots: usize, seed: u64, max_cut: usize, precision: Precision) -> Result, f64)>, CutRankError> { self.check(max_cut)?; let mut rng = seed; - let mut out = Vec::with_capacity(shots); - for _ in 0..shots { - let mut walk = Walk::new(self); - let mut x = vec![false; self.n]; - for &v in &self.order { - x[v] = walk.step(v, false, Some(&mut rng)); - } - let a = walk.finish(); - out.push((x, a.re * a.re + a.im * a.im)); + Ok((0..shots) + .map(|_| { + let (x, a) = match precision { + Precision::Double => self.run::(None, &mut rng), + Precision::Single => self.run::(None, &mut rng), + }; + (x, a.re * a.re + a.im * a.im) + }) + .collect()) + } + + /// One walk: with `x`, its amplitude; without, a sample drawn from `rng`. + fn run(&self, x: Option<&[bool]>, rng: &mut u64) -> (Vec, C) { + let mut walk = Walk::::new(self.max_cut()); + let mut out = vec![false; self.n]; + for st in &self.steps { + out[st.v] = match x { + Some(x) => walk.step(st, &self.w[st.v], x[st.v], None), + None => walk.step(st, &self.w[st.v], false, Some(rng)), + }; } - Ok(out) + (out, walk.finish()) } fn check(&self, max_cut: usize) -> Result<(), CutRankError> { @@ -254,35 +363,6 @@ fn insert_zero(i: usize, p: u32) -> usize { ((i >> p) << (p + 1)) | low } -/// The table's two arrays, shared by the threads of a pass: each thread -/// writes only the indices of the groups it was given, and no two groups -/// share an index. -#[derive(Clone, Copy)] -struct Shared { - re: *mut f64, - im: *mut f64, -} - -// SAFETY: access is partitioned by group, as above. -unsafe impl Send for Shared {} -unsafe impl Sync for Shared {} - -impl Shared { - #[inline(always)] - fn get(self, i: usize) -> (f64, f64) { - // SAFETY: callers pass indices inside the table, of their own groups. - unsafe { (*self.re.add(i), *self.im.add(i)) } - } - #[inline(always)] - fn set(self, i: usize, v: (f64, f64)) { - // SAFETY: as for `get`. - unsafe { - *self.re.add(i) = v.0; - *self.im.add(i) = v.1; - } - } -} - fn threads() -> usize { #[cfg(not(target_arch = "wasm32"))] { @@ -294,52 +374,24 @@ fn threads() -> usize { } } -/// Run `f(lo, hi)` over `0..len` in chunks, on the host's threads. -fn par(len: usize, f: &(dyn Fn(usize, usize) + Sync)) { - let t = threads(); - if t == 1 || len < 1 << 15 || cfg!(target_arch = "wasm32") { - f(0, len); - return; - } - let chunk = len.div_ceil(t * 8).max(1 << 12); - let next = std::sync::atomic::AtomicUsize::new(0); - std::thread::scope(|sc| { - for _ in 0..t { - sc.spawn(|| loop { - let lo = next.fetch_add(chunk, std::sync::atomic::Ordering::Relaxed); - if lo >= len { - break; - } - f(lo, (lo + chunk).min(len)); - }); - } - }); -} - -/// Fixed chunks a reduction is summed in, so its bits do not depend on the -/// threads. +/// Fixed chunks a pass is cut into, so a reduction's bits do not depend on +/// the threads. const CHUNKS: usize = 256; -/// `Σ_{i [f64; 4] + Sync)) -> [f64; 4] { +/// Run `f(lo, hi) -> partial sums` over `0..len` in [`CHUNKS`] fixed chunks, +/// on the host's threads, and add the chunks' sums in order. +fn chunked(len: usize, f: &(dyn Fn(usize, usize) -> [f64; 4] + Sync)) -> [f64; 4] { let chunks = len.clamp(1, CHUNKS); let per = len.div_ceil(chunks); let slots: Vec> = (0..chunks).map(|_| std::sync::Mutex::new([0.0; 4])).collect(); - // `par` hands out chunk indices; each chunk is summed in order. - let run = |lo: usize, hi: usize| { - for c in lo..hi { - let mut acc = [0.0f64; 4]; - for i in c * per..((c + 1) * per).min(len) { - let t = term(i); - for k in 0..4 { - acc[k] += t[k]; - } - } - *slots[c].lock().unwrap() = acc; + let run = |c: usize| { + let (lo, hi) = (c * per, ((c + 1) * per).min(len)); + if lo < hi { + *slots[c].lock().unwrap() = f(lo, hi); } }; - if len >= 1 << 15 { - let t = threads().min(chunks); + let t = threads().min(chunks); + if t > 1 && len >= 1 << 15 && !cfg!(target_arch = "wasm32") { let next = std::sync::atomic::AtomicUsize::new(0); std::thread::scope(|sc| { for _ in 0..t { @@ -348,59 +400,189 @@ fn reduce4(len: usize, term: &(dyn Fn(usize) -> [f64; 4] + Sync)) -> [f64; 4] { if c >= chunks { break; } - run(c, c + 1); + run(c); }); } }); } else { - run(0, chunks); + (0..chunks).for_each(run); } slots.into_iter().map(|s| s.into_inner().unwrap()).fold([0.0; 4], |a, p| [a[0] + p[0], a[1] + p[1], a[2] + p[2], a[3] + p[3]]) } -struct Walk<'a> { - pm: &'a ProductMeasured, - /// The table's space: one vector a coordinate, over the qubits not yet summed. - basis: Vec>, - re: Vec, - im: Vec, - /// The qubits not yet summed. - ahead: Vec, +/// How the table's entries are stored. Arithmetic is in double precision +/// either way; `Single` halves the memory and the traffic. +#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)] +pub enum Precision { + #[default] + Double, + Single, } -impl<'a> Walk<'a> { - fn new(pm: &'a ProductMeasured) -> Walk<'a> { - let words = pm.n.div_ceil(64); - let mut ahead = vec![u64::MAX; words]; - if pm.n % 64 != 0 { - ahead[words - 1] = (1u64 << (pm.n % 64)) - 1; +impl Precision { + /// Bytes of one table entry (its real and imaginary parts). + pub fn entry_bytes(self) -> u64 { + match self { + Precision::Double => 16, + Precision::Single => 8, } - let cap = 1usize << pm.max_cut(); + } +} + +/// A table entry's storage. +trait Elem: Copy + Send + Sync + 'static { + fn ld(self) -> f64; + fn st(x: f64) -> Self; +} + +impl Elem for f64 { + #[inline(always)] + fn ld(self) -> f64 { + self + } + #[inline(always)] + fn st(x: f64) -> f64 { + x + } +} + +impl Elem for f32 { + #[inline(always)] + fn ld(self) -> f64 { + self as f64 + } + #[inline(always)] + fn st(x: f64) -> f32 { + x as f32 + } +} + +/// The table's two arrays, shared by the threads of a pass: each thread +/// writes only the indices of the groups it was given, and no two groups +/// share an index. +struct Shared { + re: *mut T, + im: *mut T, +} + +impl Clone for Shared { + fn clone(&self) -> Self { + *self + } +} +impl Copy for Shared {} + +// SAFETY: access is partitioned by group, as above. +unsafe impl Send for Shared {} +unsafe impl Sync for Shared {} + +impl Shared { + #[inline(always)] + fn get(self, i: usize) -> (f64, f64) { + // SAFETY: callers pass indices inside the table, of their own groups. + unsafe { ((*self.re.add(i)).ld(), (*self.im.add(i)).ld()) } + } + #[inline(always)] + fn set(self, i: usize, v: (f64, f64)) { + // SAFETY: as for `get`. + unsafe { + *self.re.add(i) = T::st(v.0); + *self.im.add(i) = T::st(v.1); + } + } +} + +/// One qubit's step, the same for every bit string: fixed by the graph and +/// the order alone. +#[derive(Clone, Debug, PartialEq, Eq)] +struct Step { + v: usize, + /// Coordinates swapped first, so the kernel's last coordinate is the top. + swap: Option<(usize, usize)>, + /// The signs `(−1)^{z·ℓ(c)}`, in the coordinates after the swap. + ell: usize, + /// Merge along the kernel (its coordinates below the top), dropping the top. + kernel: Option, + /// The edges ahead: `None` a new coordinate (the table doubles), else + /// their coordinates. + alpha: Option, +} + +/// The steps along `order` (see [`Walk`]). +fn plan(adjacency: &[Vec], n: usize, order: &[usize]) -> Vec { + let words = n.div_ceil(64); + let mut ahead = vec![u64::MAX; words]; + if n % 64 != 0 { + ahead[words - 1] = (1u64 << (n % 64)) - 1; + } + let mut basis: Vec> = Vec::new(); + let mut steps = Vec::with_capacity(n); + for &v in order { + let mut ell = 0usize; + for (j, b) in basis.iter().enumerate() { + if get(b, v) { + ell |= 1 << j; + } + } + ahead[v / 64] &= !(1u64 << (v % 64)); + for b in &mut basis { + b[v / 64] &= !(1u64 << (v % 64)); + } + let (mut swap, mut kernel) = (None, None); + if let Some(mut kappa) = dependency(&basis) { + let top = basis.len() - 1; + let j0 = usize::BITS as usize - 1 - kappa.leading_zeros() as usize; + if j0 != top { + basis.swap(j0, top); + let sw = |m: usize| (m & !(1 << j0) & !(1 << top)) | ((m >> top & 1) << j0) | ((m >> j0 & 1) << top); + ell = sw(ell); + kappa = sw(kappa); + swap = Some((j0, top)); + } + kernel = Some(kappa & !(1 << top)); + basis.pop(); + } + let a: Vec = adjacency[v].iter().zip(&ahead).map(|(x, y)| x & y).collect(); + let alpha = coordinates(&basis, &a); + if alpha.is_none() { + basis.push(a); + } + steps.push(Step { v, swap, ell, kernel, alpha }); + } + debug_assert!(basis.is_empty(), "every qubit summed"); + steps +} + +/// A sum in progress: the table over the current boundary space. +struct Walk { + re: Vec, + im: Vec, +} + +impl Walk { + fn new(max_cut: usize) -> Walk { + let cap = 1usize << max_cut; let (mut re, mut im) = (Vec::with_capacity(cap), Vec::with_capacity(cap)); - re.push(1.0); - im.push(0.0); - Walk { pm, basis: Vec::new(), re, im, ahead } + re.push(T::st(1.0)); + im.push(T::st(0.0)); + Walk { re, im } } - fn shared(&mut self) -> Shared { + fn shared(&mut self) -> Shared { Shared { re: self.re.as_mut_ptr(), im: self.im.as_mut_ptr() } } fn finish(self) -> C { debug_assert_eq!(self.re.len(), 1); - complex(self.re[0], self.im[0]) + complex(self.re[0].ld(), self.im[0].ld()) } - /// Swap coordinates `a` and `b` (table, basis). + /// Swap coordinates `a` and `b`. fn swap_coords(&mut self, a: usize, b: usize) { - if a == b { - return; - } let (lo, hi) = (a.min(b) as u32, a.max(b) as u32); - self.basis.swap(a, b); let len = self.re.len() >> 2; let t = self.shared(); - par(len, &|s, e| { + chunked(len, &|s, e| { for i in s..e { // Bit lo = 1, bit hi = 0, against bit lo = 0, bit hi = 1. let base = insert_zero(insert_zero(i, lo), hi); @@ -409,207 +591,248 @@ impl<'a> Walk<'a> { t.set(p, y); t.set(q, x); } + [0.0; 4] }); } - /// Sum qubit `v` with outcome `bit`; or, given `rng`, draw its outcome - /// from its conditional probability and return it. - fn step(&mut self, v: usize, bit: bool, rng: Option<&mut u64>) -> bool { - let w = self.pm.w[v]; - let mut ell = 0usize; - for (j, b) in self.basis.iter().enumerate() { - if get(b, v) { - ell |= 1 << j; - } - } - self.ahead[v / 64] &= !(1u64 << (v % 64)); - for b in &mut self.basis { - b[v / 64] &= !(1u64 << (v % 64)); - } - // A combination of the projected basis that vanishes (at most one). - let kernel = dependency(&self.basis); - let merged = if let Some(mut kappa) = kernel { - let top = self.basis.len() - 1; - let j0 = usize::BITS as usize - 1 - kappa.leading_zeros() as usize; - if j0 != top { - self.swap_coords(j0, top); - let swap = |m: usize| { - let (x, y) = (m >> j0 & 1, m >> top & 1); - (m & !(1 << j0) & !(1 << top)) | (y << j0) | (x << top) - }; - ell = swap(ell); - kappa = swap(kappa); - } - self.merge(kappa & !(1 << top), ell); - self.basis.pop(); - true - } else { - false - }; - // The qubit's edges ahead, in the table's coordinates. - let a: Vec = self.pm.adjacency[v].iter().zip(&self.ahead).map(|(x, y)| x & y).collect(); - let alpha = coordinates(&self.basis, &a); - // The outcome: given, or drawn from the conditional probabilities. - let chosen = match rng { - None => bit, - Some(rng) => { - let (sa, sb, x) = self.overlap(merged, alpha, ell); - // Σ|w0·A + w1·B|² = |w0|²Σ|A|² + |w1|²Σ|B|² + 2·Re(w0·conj(w1)·ΣA·conj(B)). - let norm = |b: bool| { - let (w0, w1) = (w[2 * usize::from(b)], w[2 * usize::from(b) + 1]); - let (cr, ci) = (w0.re * w1.re + w0.im * w1.im, w0.im * w1.re - w0.re * w1.im); - (w0.re * w0.re + w0.im * w0.im) * sa + (w1.re * w1.re + w1.im * w1.im) * sb + 2.0 * (cr * x.0 - ci * x.1) - }; - let (p0, p1) = (norm(false).max(0.0), norm(true).max(0.0)); - let u = ((splitmix(rng) >> 11) as f64 + 0.5) / 9_007_199_254_740_992.0; - u * (p0 + p1) >= p0 - } - }; - let (w0, w1) = (w[2 * usize::from(chosen)], w[2 * usize::from(chosen) + 1]); - self.apply(merged, alpha, ell, w0, w1); - if alpha.is_none() { - self.basis.push(a); - } - chosen - } - - /// Merge along the kernel: the table's top coordinate goes, and its two - /// halves become the qubit's z = 0 and z = 1 sums (the latter signed by - /// `ell`), each indexed by the remaining coordinates. - fn merge(&mut self, kappa_low: usize, ell: usize) { + /// Merge along the kernel: the top coordinate goes; the lower half + /// becomes the qubit's z = 0 sum times `w0`, the upper its z = 1 sum + /// (signed by `ell`) times `w1`, each indexed by the remaining + /// coordinates. Returns `[Σ|lower|², Σ|upper|²]` before the weights. + fn merge(&mut self, kappa_low: usize, ell: usize, w0: C, w1: C) -> [f64; 2] { let top = self.re.len() >> 1; let t = self.shared(); - if kappa_low == 0 { - par(top, &|s, e| { + let cm = |w: C, v: (f64, f64)| (w.re * v.0 - w.im * v.1, w.re * v.1 + w.im * v.0); + let m2 = |v: (f64, f64)| v.0 * v.0 + v.1 * v.1; + let [a, b, _, _] = if kappa_low == 0 { + chunked(top, &|s, e| { + let mut acc = [0.0; 4]; for c in s..e { - let (a, b) = (t.get(c), t.get(c | top)); - let (sa, sb) = (sign(c, ell), sign(c | top, ell)); - t.set(c, (a.0 + b.0, a.1 + b.1)); - t.set(c | top, (sa * a.0 + sb * b.0, sa * a.1 + sb * b.1)); + let (x, y) = (t.get(c), t.get(c | top)); + let (sx, sy) = (sign(c, ell), sign(c | top, ell)); + let m0 = (x.0 + y.0, x.1 + y.1); + let m1 = (sx * x.0 + sy * y.0, sx * x.1 + sy * y.1); + acc[0] += m2(m0); + acc[1] += m2(m1); + t.set(c, cm(w0, m0)); + t.set(c | top, cm(w1, m1)); } - }); + acc + }) } else { let hk = usize::BITS - 1 - kappa_low.leading_zeros(); - par(top >> 1, &|s, e| { + chunked(top >> 1, &|s, e| { + let mut acc = [0.0; 4]; for i in s..e { let c = insert_zero(i, hk); let c2 = c ^ kappa_low; - let (a, b, aa, bb) = (t.get(c), t.get(c2), t.get(c | top), t.get(c2 | top)); + let (x, y, xx, yy) = (t.get(c), t.get(c2), t.get(c | top), t.get(c2 | top)); // c pairs with c2 | top, and c2 with c | top. - t.set(c, (a.0 + bb.0, a.1 + bb.1)); - t.set(c2, (b.0 + aa.0, b.1 + aa.1)); + let m0c = (x.0 + yy.0, x.1 + yy.1); + let m0c2 = (y.0 + xx.0, y.1 + xx.1); let (s1, s2, s3, s4) = (sign(c, ell), sign(c2 | top, ell), sign(c2, ell), sign(c | top, ell)); - t.set(c | top, (s1 * a.0 + s2 * bb.0, s1 * a.1 + s2 * bb.1)); - t.set(c2 | top, (s3 * b.0 + s4 * aa.0, s3 * b.1 + s4 * aa.1)); + let m1c = (s1 * x.0 + s2 * yy.0, s1 * x.1 + s2 * yy.1); + let m1c2 = (s3 * y.0 + s4 * xx.0, s3 * y.1 + s4 * xx.1); + acc[0] += m2(m0c) + m2(m0c2); + acc[1] += m2(m1c) + m2(m1c2); + t.set(c, cm(w0, m0c)); + t.set(c2, cm(w0, m0c2)); + t.set(c | top, cm(w1, m1c)); + t.set(c2 | top, cm(w1, m1c2)); } - }); - } + acc + }) + }; + [a, b] } - /// `(Σ|A|², Σ|B|², Σ A·conj(B))` over the outputs, where an output is - /// `w0·A + w1·B`: what the conditional probabilities need. - fn overlap(&self, merged: bool, alpha: Option, ell: usize) -> (f64, f64, (f64, f64)) { - let len = self.re.len(); - let (re, im) = (&self.re, &self.im); - let m2 = |i: usize| re[i] * re[i] + im[i] * im[i]; - // conj-product A·conj(B) of entries i (A) and j (B), times s. - let cp = |i: usize, j: usize, s: f64| [s * (re[i] * re[j] + im[i] * im[j]), s * (im[i] * re[j] - re[i] * im[j])]; - match (merged, alpha) { - // A and B on separate outputs. - (true, None) => { - let half = len >> 1; - let [a, b, _, _] = reduce4(half, &|c| [m2(c), m2(c | half), 0.0, 0.0]); - (a, b, (0.0, 0.0)) + /// `[Re, Im]` of `Σ_c lower[c] · conj(upper[c ⊕ α])`. + fn cross_halves(&mut self, al: usize) -> (f64, f64) { + let half = self.re.len() >> 1; + let t = self.shared(); + let [r, i, _, _] = chunked(half, &|s, e| { + let mut acc = [0.0; 4]; + for c in s..e { + let (x, y) = (t.get(c), t.get((c ^ al) | half)); + acc[0] += x.0 * y.0 + x.1 * y.1; + acc[1] += x.1 * y.0 - x.0 * y.1; } - (false, None) => { - let [s, _, _, _] = reduce4(len, &|c| [m2(c), 0.0, 0.0, 0.0]); - (s, s, (0.0, 0.0)) - } - // A = lower[c], B = upper[c ⊕ α]. - (true, Some(al)) => { - let half = len >> 1; - let [a, b, xr, xi] = reduce4(half, &|c| { - let d = (c ^ al) | half; - let [xr, xi] = cp(c, d, 1.0); - [m2(c), m2(d), xr, xi] - }); - (a, b, (xr, xi)) - } - // A = T[c], B = sign(c ⊕ α)·T[c ⊕ α]. - (false, Some(al)) => { - let [s, xr, xi, _] = reduce4(len, &|c| { - let d = c ^ al; - let [xr, xi] = cp(c, d, sign(d, ell)); - [m2(c), xr, xi, 0.0] - }); - (s, s, (xr, xi)) - } - } + acc + }); + (r, i) } - /// The qubit's pass: outputs `w0·A + w1·B` (see [`Walk::overlap`]). - fn apply(&mut self, merged: bool, alpha: Option, ell: usize, w0: C, w1: C) { + /// After a merge, the edges ahead in span: `lower[c] ← w0·lower[c] + + /// w1·upper[c ⊕ α]`, and the upper half goes. + fn fold(&mut self, al: usize, w0: C, w1: C) { + let half = self.re.len() >> 1; + let t = self.shared(); + let cm = |w: C, v: (f64, f64)| (w.re * v.0 - w.im * v.1, w.re * v.1 + w.im * v.0); + chunked(half, &|s, e| { + for c in s..e { + let (a, b) = (cm(w0, t.get(c)), cm(w1, t.get((c ^ al) | half))); + t.set(c, (a.0 + b.0, a.1 + b.1)); + } + [0.0; 4] + }); + self.re.truncate(half); + self.im.truncate(half); + } + + /// No merge, the edges ahead new: the table doubles, `upper[c] = w1·s(c)·T[c]`, `lower[c] = w0·T[c]`. + fn double(&mut self, ell: usize, w0: C, w1: C) { let len = self.re.len(); + debug_assert!(self.re.capacity() >= 2 * len); + // SAFETY: the capacity holds 2·len (reserved for the widest cut), and + // every new entry is written by the pass below before anything reads it. + unsafe { + self.re.set_len(2 * len); + self.im.set_len(2 * len); + } + let t = self.shared(); + let cm = |w: C, v: (f64, f64)| (w.re * v.0 - w.im * v.1, w.re * v.1 + w.im * v.0); + chunked(len, &|s, e| { + for c in s..e { + let v = t.get(c); + let sv = sign(c, ell); + t.set(c | len, cm(w1, (sv * v.0, sv * v.1))); + t.set(c, cm(w0, v)); + } + [0.0; 4] + }); + } + + /// No merge, the edges ahead in span: `T[c] ← w0·T[c] + w1·s(c ⊕ α)·T[c ⊕ α]`. + fn xor_sum(&mut self, al: usize, ell: usize, w0: C, w1: C) { + let len = self.re.len(); + let t = self.shared(); let cm = |w: C, v: (f64, f64)| (w.re * v.0 - w.im * v.1, w.re * v.1 + w.im * v.0); let add = |a: (f64, f64), b: (f64, f64)| (a.0 + b.0, a.1 + b.1); - match (merged, alpha) { - (true, None) => { - let half = len >> 1; - let t = self.shared(); - par(half, &|s, e| { - for c in s..e { - t.set(c, cm(w0, t.get(c))); - t.set(c | half, cm(w1, t.get(c | half))); - } - }); + if al == 0 { + chunked(len, &|s, e| { + for c in s..e { + let v = t.get(c); + let sv = sign(c, ell); + t.set(c, add(cm(w0, v), cm(w1, (sv * v.0, sv * v.1)))); + } + [0.0; 4] + }); + return; + } + let hb = usize::BITS - 1 - al.leading_zeros(); + chunked(len >> 1, &|s, e| { + for i in s..e { + let c = insert_zero(i, hb); + let d = c ^ al; + let (x, y) = (t.get(c), t.get(d)); + let (sx, sy) = (sign(c, ell), sign(d, ell)); + t.set(c, add(cm(w0, x), cm(w1, (sy * y.0, sy * y.1)))); + t.set(d, add(cm(w0, y), cm(w1, (sx * x.0, sx * x.1)))); } - (true, Some(al)) => { - let half = len >> 1; - let t = self.shared(); - par(half, &|s, e| { - for c in s..e { - t.set(c, add(cm(w0, t.get(c)), cm(w1, t.get((c ^ al) | half)))); - } - }); - self.re.truncate(half); - self.im.truncate(half); + [0.0; 4] + }); + } + + /// Without a merge: `[Σ|T|², Re, Im]` of `Σ_c T[c]·conj(s(c⊕α)·T[c⊕α])`. + fn norms(&mut self, al: Option, ell: usize) -> (f64, (f64, f64)) { + let len = self.re.len(); + let t = self.shared(); + let [s, r, i, _] = chunked(len, &|lo, hi| { + let mut acc = [0.0; 4]; + for c in lo..hi { + let x = t.get(c); + acc[0] += x.0 * x.0 + x.1 * x.1; + if let Some(al) = al { + let d = c ^ al; + let (y, sy) = (t.get(d), sign(d, ell)); + acc[1] += sy * (x.0 * y.0 + x.1 * y.1); + acc[2] += sy * (x.1 * y.0 - x.0 * y.1); + } } - (false, None) => { - self.re.resize(2 * len, 0.0); - self.im.resize(2 * len, 0.0); - let t = self.shared(); - par(len, &|s, e| { - for c in s..e { - let v = t.get(c); - let sv = sign(c, ell); - t.set(c | len, cm(w1, (sv * v.0, sv * v.1))); - t.set(c, cm(w0, v)); + acc + }); + (s, (r, i)) + } + + /// One step with outcome `bit`, or, given `rng`, an outcome drawn from + /// its conditional probability; returns the outcome. + fn step(&mut self, st: &Step, w: &Unitary, bit: bool, rng: Option<&mut u64>) -> bool { + if let Some((a, b)) = st.swap { + self.swap_coords(a, b); + } + let ws = |b: bool| (w[2 * usize::from(b)], w[2 * usize::from(b) + 1]); + // Σ|w0·A + w1·B|² = |w0|²Σ|A|² + |w1|²Σ|B|² + 2·Re(w0·conj(w1)·ΣA·conj(B)). + let norm = |b: bool, sa: f64, sb: f64, x: (f64, f64)| { + let (w0, w1) = ws(b); + let (cr, ci) = (w0.re * w1.re + w0.im * w1.im, w0.im * w1.re - w0.re * w1.im); + ((w0.re * w0.re + w0.im * w0.im) * sa + (w1.re * w1.re + w1.im * w1.im) * sb + 2.0 * (cr * x.0 - ci * x.1)).max(0.0) + }; + let draw = |rng: &mut u64, p0: f64, p1: f64| { + let u = ((splitmix(rng) >> 11) as f64 + 0.5) / 9_007_199_254_740_992.0; + u * (p0 + p1) >= p0 + }; + let one = complex(1.0, 0.0); + match (st.kernel, st.alpha) { + (Some(k), None) => match rng { + None => { + let (w0, w1) = ws(bit); + self.merge(k, st.ell, w0, w1); + bit + } + Some(rng) => { + let [sa, sb] = self.merge(k, st.ell, one, one); + let b = draw(rng, norm(false, sa, sb, (0.0, 0.0)), norm(true, sa, sb, (0.0, 0.0))); + let (w0, w1) = ws(b); + // Scale the halves now the outcome is known. + let half = self.re.len() >> 1; + let t = self.shared(); + let cm = |w: C, v: (f64, f64)| (w.re * v.0 - w.im * v.1, w.re * v.1 + w.im * v.0); + chunked(half, &|s, e| { + for c in s..e { + t.set(c, cm(w0, t.get(c))); + t.set(c | half, cm(w1, t.get(c | half))); + } + [0.0; 4] + }); + b + } + }, + (Some(k), Some(al)) => { + let [sa, sb] = self.merge(k, st.ell, one, one); + let b = match rng { + None => bit, + Some(rng) => { + let x = self.cross_halves(al); + draw(rng, norm(false, sa, sb, x), norm(true, sa, sb, x)) } - }); + }; + let (w0, w1) = ws(b); + self.fold(al, w0, w1); + b } - (false, Some(0)) => { - let t = self.shared(); - par(len, &|s, e| { - for c in s..e { - let v = t.get(c); - let sv = sign(c, ell); - t.set(c, add(cm(w0, v), cm(w1, (sv * v.0, sv * v.1)))); + (None, None) => { + let b = match rng { + None => bit, + Some(rng) => { + let (s, _) = self.norms(None, st.ell); + draw(rng, norm(false, s, s, (0.0, 0.0)), norm(true, s, s, (0.0, 0.0))) } - }); + }; + let (w0, w1) = ws(b); + self.double(st.ell, w0, w1); + b } - (false, Some(al)) => { - let hb = usize::BITS - 1 - al.leading_zeros(); - let t = self.shared(); - par(len >> 1, &|s, e| { - for i in s..e { - let c = insert_zero(i, hb); - let d = c ^ al; - let (x, y) = (t.get(c), t.get(d)); - let (sx, sy) = (sign(c, ell), sign(d, ell)); - t.set(c, add(cm(w0, x), cm(w1, (sy * y.0, sy * y.1)))); - t.set(d, add(cm(w0, y), cm(w1, (sx * x.0, sx * x.1)))); + (None, Some(al)) => { + let b = match rng { + None => bit, + Some(rng) => { + let (s, x) = self.norms(Some(al), st.ell); + draw(rng, norm(false, s, s, x), norm(true, s, s, x)) } - }); + }; + let (w0, w1) = ws(b); + self.xor_sum(al, st.ell, w0, w1); + b } } } @@ -838,6 +1061,27 @@ mod tests { assert!(compared >= 30, "only {compared} comparisons ran"); } + /// Single-precision storage, double-precision arithmetic: amplitudes + /// within single-precision rounding of the double ones, and samples from + /// the same law. + #[test] + fn single_precision_tables_agree() { + let mut s = 31u64; + for seed in 0..8u64 { + let n = 10 + 3 * seed as usize; + let ops = random_clifford(n as u32, 4 * n, seed); + let last: Vec = (0..n).map(|_| random_unitary(&mut s)).collect(); + let pm = ProductMeasured::new(&Tableau::run(n as u32, &ops).unwrap(), &last).unwrap(); + for _ in 0..3 { + let x: Vec = (0..n).map(|_| next(&mut s) & 1 == 1).collect(); + let (a, b) = (pm.amplitude(&x, 30).unwrap(), pm.amplitude_in(&x, 30, Precision::Single).unwrap()); + let scale = a.re.hypot(a.im).max(1e-300); + assert!((a.re - b.re).hypot(a.im - b.im) / scale < 1e-5, "seed {seed}: {a:?} vs {b:?}"); + } + } + assert_eq!(Precision::Single.entry_bytes() * 2, Precision::Double.entry_bytes()); + } + #[test] fn samples_follow_the_probabilities() { let n = 6usize; @@ -865,6 +1109,22 @@ mod tests { assert_eq!(pm.sample(50, 5, 10).unwrap(), pm.sample(50, 5, 10).unwrap()); } + #[test] + fn annealing_lowers_the_work_and_keeps_the_amplitudes() { + let n = 28usize; + let mut s = 5u64; + let ops = random_clifford(n as u32, 5 * n, 17); + let last: Vec = (0..n).map(|_| random_unitary(&mut s)).collect(); + let pm = ProductMeasured::new(&Tableau::run(n as u32, &ops).unwrap(), &last).unwrap(); + let order = pm.anneal_order(400, 3); + assert_eq!(order, pm.anneal_order(400, 3), "seeded"); + let better = pm.clone().with_order(&order).unwrap(); + assert!(better.work() <= pm.work(), "{} > {}", better.work(), pm.work()); + let x: Vec = (0..n).map(|_| next(&mut s) & 1 == 1).collect(); + let (a, b) = (pm.amplitude(&x, 28).unwrap(), better.amplitude(&x, 28).unwrap()); + assert!((a.re - b.re).hypot(a.im - b.im) < 1e-9 * a.re.hypot(a.im).max(1e-300), "{a:?} vs {b:?}"); + } + #[test] fn cut_ranks_are_the_entanglement() { // A GHZ chain: one bit across every cut. A product state: none. -- 2.54.0 (Apple Git-157) From bdd46363b505d61816446a30522d3d91a7a83d9f Mon Sep 17 00:00:00 2001 From: Transaction Science Date: Mon, 19 Oct 2026 07:31:33 -0400 Subject: [PATCH 3/5] wai-quantum examples: patched cross-entropy benchmarking from shots (patch_xeb) and one amplitude's contraction cost (tn_cost) patch_xeb simulates each patch of a patched circuit exactly on the dense state vector and scores the measured shots (normalised linear XEB per patch, the product, its error). On a published 61-qubit experiment's release it reproduces the released per-circuit fidelities to 5e-7 (their JSON's rounding). tn_cost prints log10 of one amplitude's contraction cost as find_path sees it, with the largest intermediate and the slices. Co-Authored-By: Claude Opus 5.5 --- wai/wai-quantum/Cargo.toml | 8 +++ wai/wai-quantum/examples/patch_xeb.rs | 89 +++++++++++++++++++++++++++ wai/wai-quantum/examples/tn_cost.rs | 65 +++++++++++++++++++ 3 files changed, 162 insertions(+) create mode 100644 wai/wai-quantum/examples/patch_xeb.rs create mode 100644 wai/wai-quantum/examples/tn_cost.rs diff --git a/wai/wai-quantum/Cargo.toml b/wai/wai-quantum/Cargo.toml index 40fbf2860..fa1cd8ecb 100644 --- a/wai/wai-quantum/Cargo.toml +++ b/wai/wai-quantum/Cargo.toml @@ -305,6 +305,14 @@ required-features = ["quantum_fabric"] name = "rv32_listing" required-features = ["quantum_fabric_dataflow"] +[[example]] +name = "tn_cost" +required-features = ["quantum_tn"] + +[[example]] +name = "patch_xeb" +required-features = ["quantum_sv"] + [[example]] name = "cutrank_qasm" required-features = ["quantum_cutrank"] diff --git a/wai/wai-quantum/examples/patch_xeb.rs b/wai/wai-quantum/examples/patch_xeb.rs new file mode 100644 index 000000000..9554b4de3 --- /dev/null +++ b/wai/wai-quantum/examples/patch_xeb.rs @@ -0,0 +1,89 @@ +//! Patched cross-entropy benchmarking from a circuit, its patches and the +//! measured shots: each patch simulated exactly on the dense state vector. +//! +//! ```text +//! cargo run --release --example patch_xeb --features quantum_sv -- CIRCUIT.qasm SHOTS.u64 PATCHES +//! ``` +//! +//! `CIRCUIT.qasm`: OpenQASM 3 with `rz`, `rx`, `ry`, `h`, `x`, `sx`, `cz`, +//! `cx` (barriers and measurements are skipped). Its two-qubit gates must stay +//! inside the patches: a patched circuit has its boundary couplers removed. +//! `SHOTS.u64`: one little-endian `u64` a shot, qubit q at bit q. +//! `PATCHES`: the patches' qubits, `0,1,2;3,4,5`. +//! +//! Prints, per patch: its qubits, the mean of `2^n·p(x)` over the shots, its +//! ideal value `2^n·Σp² − 1` (the XEB of the patch's own exact samples), and +//! the normalised fidelity `(2^n·mean p − 1) / (2^n·Σp² − 1)`; then the +//! product of the fidelities and its standard error (the patches' errors +//! added in quadrature, relative). + +use wai_quantum::gates::Gate; +use wai_quantum::quantum_sv::StateVector; + +fn main() { + let args: Vec = std::env::args().collect(); + let src = std::fs::read_to_string(&args[1]).expect("the circuit reads"); + let raw = std::fs::read(&args[2]).expect("the shots read"); + let shots: Vec = raw.chunks_exact(8).map(|b| u64::from_le_bytes(b.try_into().unwrap())).collect(); + let patches: Vec> = args[3].split(';').map(|p| p.split(',').map(|q| q.trim().parse().unwrap()).collect()).collect(); + // The gates, with their qubits. + let mut gates: Vec<(String, f64, Vec)> = Vec::new(); + for line in src.lines().map(str::trim) { + let line = line.trim_end_matches(';'); + if line.is_empty() || ["OPENQASM", "include", "bit", "qubit", "barrier", "c["].iter().any(|p| line.starts_with(p)) || line.contains("measure") { + continue; + } + let (head, args) = line.split_once(' ').expect("a gate and its qubits"); + let (name, theta) = match head.split_once('(') { + Some((n, p)) => (n.to_string(), p.trim_end_matches(')').parse::().expect("one angle")), + None => (head.to_string(), 0.0), + }; + let qs: Vec = args.split(',').map(|a| a.trim().trim_start_matches("q[").trim_end_matches(']').parse().unwrap()).collect(); + gates.push((name, theta, qs)); + } + let mut product = 1.0; + let mut rel_var = 0.0; + for patch in &patches { + let n = patch.len() as u32; + let local = |q: u32| patch.iter().position(|&p| p == q).map(|i| i as u32); + let mut sv = StateVector::zero(n).unwrap(); + for (name, theta, qs) in &gates { + let l: Vec> = qs.iter().map(|&q| local(q)).collect(); + if l.iter().all(Option::is_none) { + continue; + } + let l: Vec = l.into_iter().map(|x| x.expect("a gate crossing the patch boundary")).collect(); + let g = match (name.as_str(), l.as_slice()) { + ("rz", [a]) => Gate::rz(*a, *theta), + ("rx", [a]) => Gate::rx(*a, *theta), + ("ry", [a]) => Gate::ry(*a, *theta), + ("h", [a]) => Gate::h(*a), + ("x", [a]) => Gate::x(*a), + ("sx", [a]) => Gate::sqrt_x(*a), + ("cz", [a, b]) => Gate::cz(*a, *b), + ("cx", [a, b]) => Gate::cx(*a, *b), + _ => panic!("gate {name} on {qs:?}"), + }; + sv.apply(&g).unwrap(); + } + let dim = (n as f64).exp2(); + let prob = |i: usize| sv.re[i] * sv.re[i] + sv.im[i] * sv.im[i]; + let ideal = dim * (0..sv.re.len()).map(|i| prob(i) * prob(i)).sum::() - 1.0; + // Each shot's patch bits, as the state vector's index (local qubit 0 most significant). + let vals: Vec = shots + .iter() + .map(|&x| { + let idx = patch.iter().enumerate().fold(0usize, |i, (l, &q)| i | (((x >> q) & 1) as usize) << (n as usize - 1 - l)); + dim * prob(idx) + }) + .collect(); + let mean = vals.iter().sum::() / vals.len() as f64; + let var = vals.iter().map(|v| (v - mean).powi(2)).sum::() / (vals.len() as f64 - 1.0); + let f = (mean - 1.0) / ideal; + let se = (var / vals.len() as f64).sqrt() / ideal; + product *= f; + rel_var += (se / f).powi(2); + println!("patch {n} qubits: mean 2^n·p {mean:.6}, ideal xeb {ideal:.6}, fidelity {f:.6} ± {se:.6}"); + } + println!("fidelity {product:.6e} ± {:.6e} ({} shots)", product.abs() * rel_var.sqrt(), shots.len()); +} diff --git a/wai/wai-quantum/examples/tn_cost.rs b/wai/wai-quantum/examples/tn_cost.rs new file mode 100644 index 000000000..595bcae35 --- /dev/null +++ b/wai/wai-quantum/examples/tn_cost.rs @@ -0,0 +1,65 @@ +//! The cost of one amplitude of a circuit by tensor-network contraction, as +//! the hyper-optimised path search finds it. +//! +//! ```text +//! cargo run --release --example tn_cost --features quantum_tn -- CIRCUIT.qasm TRIALS SEED [LOG2_MAX_SIZE] +//! ``` +//! +//! `CIRCUIT.qasm`: OpenQASM 3 with `rz`, `rx`, `ry`, `h`, `x`, `sx`, `cz`, +//! `cx` (barriers and measurements skipped). Prints log10 of the complex +//! multiply-adds for one amplitude (every slice), log2 of the largest +//! intermediate, and the number of slices. + +use std::time::Instant; +use wai_quantum::gates::Gate; +use wai_quantum::quantum_tn::{Network, PathOptions, find_path}; + +fn main() { + let args: Vec = std::env::args().collect(); + let src = std::fs::read_to_string(&args[1]).expect("the circuit reads"); + let trials: usize = args[2].parse().unwrap(); + let seed: u64 = args[3].parse().unwrap(); + let max_size = args.get(4).map(|s| s.parse::().unwrap().exp2()); + let mut n = 0u32; + let mut gates = Vec::new(); + for line in src.lines().map(str::trim) { + let line = line.trim_end_matches(';'); + if let Some(rest) = line.strip_prefix("qubit[") { + n = rest.split(']').next().unwrap().parse().unwrap(); + continue; + } + if line.is_empty() || ["OPENQASM", "include", "bit", "barrier", "c["].iter().any(|p| line.starts_with(p)) || line.contains("measure") { + continue; + } + let (head, args) = line.split_once(' ').expect("a gate and its qubits"); + let (name, theta) = match head.split_once('(') { + Some((n, p)) => (n, p.trim_end_matches(')').parse::().expect("one angle")), + None => (head, 0.0), + }; + let qs: Vec = args.split(',').map(|a| a.trim().trim_start_matches("q[").trim_end_matches(']').parse().unwrap()).collect(); + gates.push(match (name, qs.as_slice()) { + ("rz", [a]) => Gate::rz(*a, theta), + ("rx", [a]) => Gate::rx(*a, theta), + ("ry", [a]) => Gate::ry(*a, theta), + ("h", [a]) => Gate::h(*a), + ("x", [a]) => Gate::x(*a), + ("sx", [a]) => Gate::sqrt_x(*a), + ("cz", [a, b]) => Gate::cz(*a, *b), + ("cx", [a, b]) => Gate::cx(*a, *b), + _ => panic!("gate {name} on {qs:?}"), + }); + } + let bits = vec![0u8; n as usize]; + let mut net = Network::amplitude(n, &gates, &bits, &[]).expect("the network builds"); + net.simplify(); + let t = Instant::now(); + let plan = find_path(&net.shape(), &net.dims, &net.open, &PathOptions { trials, seed, max_size, ..PathOptions::default() }); + println!( + "{n} qubits, {} gates: log10 flops {:.3}, log2 largest {:.1}, slices {:.3e}, searched in {:?}", + gates.len(), + plan.total_flops().log10(), + plan.per_slice.max_size.log2(), + plan.slices, + t.elapsed() + ); +} -- 2.54.0 (Apple Git-157) From f599f0a6fce2a5cd43f75d53f6a724696e01c4b0 Mon Sep 17 00:00:00 2001 From: Transaction Science Date: Mon, 19 Oct 2026 11:23:35 -0400 Subject: [PATCH 4/5] wai-quantum 0.6.6: the cut-rank engine in the README and changelog A feature row and a section for quantum_cutrank; the changelog entry; the version in every lockfile that pins it. Co-Authored-By: Claude Opus 5.5 --- wai/CHANGELOG.md | 41 +++++++++++++++++++++++++ wai/tools/quantum-meter/Cargo.lock | 2 +- wai/tools/quantum-wgsl-check/Cargo.lock | 2 +- wai/wai-pack/Cargo.lock | 2 +- wai/wai-quantum-qpu/Cargo.lock | 2 +- wai/wai-quantum-serve/Cargo.lock | 2 +- wai/wai-quantum-service/Cargo.lock | 2 +- wai/wai-quantum/Cargo.lock | 2 +- wai/wai-quantum/Cargo.toml | 2 +- wai/wai-quantum/README.md | 29 +++++++++++++++++ wai/wai-rs/Cargo.lock | 2 +- wai/wai-web/Cargo.lock | 2 +- 12 files changed, 80 insertions(+), 10 deletions(-) diff --git a/wai/CHANGELOG.md b/wai/CHANGELOG.md index edeb744a4..87e9a5ed6 100644 --- a/wai/CHANGELOG.md +++ b/wai/CHANGELOG.md @@ -6,6 +6,47 @@ wire format may change and such changes are recorded here. ## v0.4.0 — current +### wai-quantum 0.6.6: Clifford circuits measured in a product basis, exactly, by cut-rank + +Additive. Every pinned digest is unchanged. + +- **`Tableau::graph_form`** (`quantum_tableau`) writes a stabilizer state as local Cliffords + (`H^h S^s Z^z` per qubit) on a graph state. Tested by rebuilding the state from the form + over 40 random circuits. +- **`quantum_cutrank`** (a new feature, in `full`) sums a Clifford circuit's state measured + after one closing layer of single-qubit gates at any angles: + `⟨x| ⊗ V_q |ψ⟩ = 2^{-n/2} Σ_z ∏_q W_q[x_q][z_q] (−1)^{Σ_{a~b} z_a z_b}`. + - The sum goes one qubit at a time along an order, in a table of `2^{cut-rank}` entries. + Each qubit is an in-place pass: a merge where the boundary space loses a dimension, the + qubit's signed terms, and a doubling or an XOR-shifted sum for its edges ahead. + - The passes depend on the graph and the order alone, so they are planned once per state. + - By Parseval, the table's norm is a prefix's probability, so the same passes sample + exactly: `amplitude`, `probability`, `sample`. + - `Precision::Single` stores the table in `f32` with `f64` arithmetic: half the memory and + traffic. `anneal_order` searches orders for less work (`work` = Σ 2^cut) and never + returns a worse one. + - Every reduction is summed in fixed chunks, so a result is the same bits at any thread + count. + - Tests: + - amplitudes equal the state vector's in probability and relative phase (30 circuits up + to 12 qubits); + - forward, reverse and shuffled orders agree (24 to 54 qubits); + - single-precision tables agree with double; + - samples follow the exact law and are seeded; + - annealing lowers the work and keeps the amplitudes; + - cut-ranks are the entanglement of GHZ and product states. + - Five mutants are caught: the merge's signs, the sampling choice, the doubling's sign, + the fold's coefficient, and the graph form without its S. +- **Examples:** + - `cutrank_qasm` reads an OpenQASM 2 Clifford bulk with a closing single-qubit layer, and + prints cut-ranks, annealed orders, amplitudes and exact samples. On a 70-qubit + graph-state sampling circuit, an annealed order keeps every cut at 29 bits or fewer: one + amplitude in about half a second, one exact sample in 1.1 to 1.5 s, with a largest + table of 4 GiB. + - `patch_xeb` scores measured shots of a patched circuit by patched cross-entropy, each + patch simulated exactly on the dense state vector. + - `tn_cost` prints the contraction cost of one amplitude as `find_path` finds it. + ### wai-quantum 0.6.5 and the service: programs that measure in the middle Additive. Every pinned digest is unchanged, and a program that measures only at the end diff --git a/wai/tools/quantum-meter/Cargo.lock b/wai/tools/quantum-meter/Cargo.lock index 904f1b44e..ce9265c79 100644 --- a/wai/tools/quantum-meter/Cargo.lock +++ b/wai/tools/quantum-meter/Cargo.lock @@ -1017,7 +1017,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a" [[package]] name = "wai-quantum" -version = "0.6.5" +version = "0.6.6" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/tools/quantum-wgsl-check/Cargo.lock b/wai/tools/quantum-wgsl-check/Cargo.lock index e9f467917..80a09f5b7 100644 --- a/wai/tools/quantum-wgsl-check/Cargo.lock +++ b/wai/tools/quantum-wgsl-check/Cargo.lock @@ -599,7 +599,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a" [[package]] name = "wai-quantum" -version = "0.6.5" +version = "0.6.6" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-pack/Cargo.lock b/wai/wai-pack/Cargo.lock index a0356395f..bb1920836 100644 --- a/wai/wai-pack/Cargo.lock +++ b/wai/wai-pack/Cargo.lock @@ -713,7 +713,7 @@ dependencies = [ [[package]] name = "wai-quantum" -version = "0.6.5" +version = "0.6.6" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-quantum-qpu/Cargo.lock b/wai/wai-quantum-qpu/Cargo.lock index 143d92e94..c76c61b77 100644 --- a/wai/wai-quantum-qpu/Cargo.lock +++ b/wai/wai-quantum-qpu/Cargo.lock @@ -422,7 +422,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a" [[package]] name = "wai-quantum" -version = "0.6.5" +version = "0.6.6" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-quantum-serve/Cargo.lock b/wai/wai-quantum-serve/Cargo.lock index 2ac649c7d..b3bc9297d 100644 --- a/wai/wai-quantum-serve/Cargo.lock +++ b/wai/wai-quantum-serve/Cargo.lock @@ -599,7 +599,7 @@ dependencies = [ [[package]] name = "wai-quantum" -version = "0.6.5" +version = "0.6.6" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-quantum-service/Cargo.lock b/wai/wai-quantum-service/Cargo.lock index 55d4fabf3..1c893cea6 100644 --- a/wai/wai-quantum-service/Cargo.lock +++ b/wai/wai-quantum-service/Cargo.lock @@ -857,7 +857,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a" [[package]] name = "wai-quantum" -version = "0.6.5" +version = "0.6.6" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-quantum/Cargo.lock b/wai/wai-quantum/Cargo.lock index 3294d0fa5..04e9e8a43 100644 --- a/wai/wai-quantum/Cargo.lock +++ b/wai/wai-quantum/Cargo.lock @@ -1079,7 +1079,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a" [[package]] name = "wai-quantum" -version = "0.6.5" +version = "0.6.6" dependencies = [ "blake3", "block2", diff --git a/wai/wai-quantum/Cargo.toml b/wai/wai-quantum/Cargo.toml index fa1cd8ecb..84d94b9c8 100644 --- a/wai/wai-quantum/Cargo.toml +++ b/wai/wai-quantum/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "wai-quantum" -version = "0.6.5" +version = "0.6.6" edition = "2024" rust-version = "1.97" license = "Apache-2.0" diff --git a/wai/wai-quantum/README.md b/wai/wai-quantum/README.md index 7f2b776da..efdf53b6e 100644 --- a/wai/wai-quantum/README.md +++ b/wai/wai-quantum/README.md @@ -91,6 +91,7 @@ never creates information. | `quantum_qec` | CSS codes (distance verified by exhaustive search, not asserted) — **rotated surface** `[[d²,1,d]]`, toric, bivariate-bicycle, and the **`[[7,1,3]]` colour code**, whose self-duality buys a *transversal* logical Hadamard — one round, no ancillas, no surgery — verified on the simulator, with the surface code shown failing the same move — with **two decoders**: union-find (weighted matching, single-shot, space-time with faulty measurement, **and circuit-level with hook errors**; **imports and exports** Stim-syntax detector error models, so other tools' noise models decode here and ours decode there) and Relay-BP (for qLDPC) | | `quantum_frame` | **noisy Clifford circuits at QEC scale**: reads the ecosystem's circuit text; samples detection events by bit-packed Pauli-frame simulation (64 shots per word); derives the detector error model by backward sensitivity, refusing non-deterministic detectors; generates the rotated surface-code memory experiment byte for byte as the reference layout; and decodes by **exact minimum-weight perfect matching** (Edmonds' blossom on integer weights). Measures logical error rates and fits the code model resource estimates assume | | `quantum_tn` | **tensor-network contraction**: circuits (any unitary on any qubits) become networks; a seeded hyper-optimiser finds contraction trees from randomised-greedy and multilevel-partition families, reconfigures them optimally in subtrees, and slices them to a memory cap with a second search around the slices; slices contract in a fixed order, threaded, to the same bits everywhere. Single amplitudes or open-output batches | +| `quantum_cutrank` | **Clifford circuits measured in a product basis, exactly**: the stabilizer state's graph form (`Tableau::graph_form`) and one closing layer of single-qubit gates at any angles, summed one qubit at a time along an order in a table of `2^{cut-rank}` entries. Exact amplitudes, probabilities and samples (by Parseval), at a cost set by the entanglement across the order's cuts, not by the number of qubits; orders annealed for less work; tables in double or single precision; the same bits at any thread count | | `quantum_gbs` | **Gaussian boson sampling, exactly**: hafnians (power traces, Hessenberg + characteristic polynomial, `O(n³ 2^{n/2})`), torontonians, lossy Gaussian states from squeezers and an interferometer, exact click and photon-number probabilities, exact samples drawn mode by mode, and a check of observed samples against the exact one- and two-mode click statistics | | `quantum_sv` | **the dense state vector at any angle**: one-, two- and `k`-qubit kernels with a diagonal fast path, gate fusion around multi-qubit gates, threads over independent blocks; every kernel path computes each amplitude in the same order, so threading never changes a bit. Agrees with an independent simulator to 2·10⁻¹⁶ | | `quantum_compile` | Clifford routing + stabilizer-tableau equivalence proof | @@ -432,6 +433,34 @@ agrees with an independent contraction to 12 significant figures. The module documentation states what slicing costs on deep circuits and what the scalar executor can and cannot contract. +## Clifford circuits in a product basis + +A Clifford circuit followed by one layer of single-qubit gates, any gates at any +angles, is out of the stabilizer tableau's reach. Its amplitudes are still a +sum over a graph state, and `quantum_cutrank` takes that sum one qubit at a +time. The table it carries has `2^r` entries, where `r` is the state's cut-rank +(its entanglement, in bits) across the cut the sum has reached. So the cost is +set by the widest cut along the order, not by the number of qubits. +`anneal_order` searches for an order with narrower cuts, and `work` reports the +entries the passes touch. + +```sh +cargo run --release --example cutrank_qasm --features quantum_cutrank -- circuit.qasm cuts +cargo run --release --example cutrank_qasm --features quantum_cutrank -- circuit.qasm order 20000 3 > order.txt +CUTRANK_ORDER=order.txt CUTRANK_PRECISION=single \ + cargo run --release --example cutrank_qasm --features quantum_cutrank -- circuit.qasm sample 1000 7 +``` + +`order` anneals for the least work and prints the order it found; `sample` +prints each exact sample with `2^n · p`. + +A 70-qubit graph-state sampling circuit has a nearest-neighbour Clifford bulk, +then a non-Clifford basis on every qubit. With an annealed order, every cut is +29 bits or fewer. On one workstation, one amplitude takes about half a second +and one exact sample 1.1 to 1.5 s, with a largest table of 4 GiB in single +precision. Amplitudes agree between orders and between precisions to about +10⁻⁷ relative. + ## Gaussian boson sampling ```sh diff --git a/wai/wai-rs/Cargo.lock b/wai/wai-rs/Cargo.lock index d27c01998..358547e67 100644 --- a/wai/wai-rs/Cargo.lock +++ b/wai/wai-rs/Cargo.lock @@ -4310,7 +4310,7 @@ version = "0.1.0" [[package]] name = "wai-quantum" -version = "0.6.5" +version = "0.6.6" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-web/Cargo.lock b/wai/wai-web/Cargo.lock index a18a919fd..38c97d7f3 100644 --- a/wai/wai-web/Cargo.lock +++ b/wai/wai-web/Cargo.lock @@ -1241,7 +1241,7 @@ version = "0.1.0" [[package]] name = "wai-quantum" -version = "0.6.5" +version = "0.6.6" dependencies = [ "blake3", "ed25519-dalek", -- 2.54.0 (Apple Git-157) From c1c325002f27ebca46f3cb35ab41d9a29a040dfb Mon Sep 17 00:00:00 2001 From: Transaction Science Date: Mon, 19 Oct 2026 11:29:49 -0400 Subject: [PATCH 5/5] wai-quantum cutrank: clippy's is_multiple_of and an enumerated loop Co-Authored-By: Claude Opus 5.5 --- wai/wai-quantum/src/quantum_cutrank.rs | 6 +++--- wai/wai-quantum/src/quantum_tableau.rs | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/wai/wai-quantum/src/quantum_cutrank.rs b/wai/wai-quantum/src/quantum_cutrank.rs index 6462ff4c5..7505df12f 100644 --- a/wai/wai-quantum/src/quantum_cutrank.rs +++ b/wai/wai-quantum/src/quantum_cutrank.rs @@ -158,7 +158,7 @@ impl ProductMeasured { fn cut_ranks_of(&self, order: &[usize]) -> Vec { let words = self.n.div_ceil(64); let mut ahead = vec![u64::MAX; words]; - if self.n % 64 != 0 { + if !self.n.is_multiple_of(64) { ahead[words - 1] = (1u64 << (self.n % 64)) - 1; } let mut cuts = Vec::with_capacity(self.n); @@ -353,7 +353,7 @@ fn splitmix(s: &mut u64) -> u64 { /// `(−1)^{|c ∧ ℓ|}`. #[inline(always)] fn sign(c: usize, l: usize) -> f64 { - if (c & l).count_ones() % 2 == 0 { 1.0 } else { -1.0 } + if (c & l).count_ones().is_multiple_of(2) { 1.0 } else { -1.0 } } /// Insert a zero bit at position `p` of `i`. @@ -512,7 +512,7 @@ struct Step { fn plan(adjacency: &[Vec], n: usize, order: &[usize]) -> Vec { let words = n.div_ceil(64); let mut ahead = vec![u64::MAX; words]; - if n % 64 != 0 { + if !n.is_multiple_of(64) { ahead[words - 1] = (1u64 << (n % 64)) - 1; } let mut basis: Vec> = Vec::new(); diff --git a/wai/wai-quantum/src/quantum_tableau.rs b/wai/wai-quantum/src/quantum_tableau.rs index 34769bd1a..3a8f8c4ac 100644 --- a/wai/wai-quantum/src/quantum_tableau.rs +++ b/wai/wai-quantum/src/quantum_tableau.rs @@ -404,7 +404,7 @@ impl Tableau { // The Z-only rows, reduced: their pivot qubits take a Hadamard. let mut h = vec![false; n]; let mut k2 = k; - for q in 0..n { + for (q, hq) in h.iter_mut().enumerate() { let Some(p) = (k2..n).find(|&i| get(s.zr(i), q)) else { continue }; s.swap(k2, p); for i in k..n { @@ -412,7 +412,7 @@ impl Tableau { s.mul_into(i, k2); } } - h[q] = true; + *hq = true; k2 += 1; } debug_assert_eq!(k2, n, "n independent generators"); -- 2.54.0 (Apple Git-157)