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

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",

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",

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",

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",

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",

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",

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",

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"
@@ -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,18 @@ 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"]
+
 [[example]]
 name = "npu_onnx"
 required-features = ["quantum_fabric"]

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

wai/wai-quantum/examples/cutrank_qasm.rs

@@ -0,0 +1,173 @@
+//! 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::{Precision, 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<Unitary> {
+    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<String> = 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<Unitary> = Vec::new();
+    let mut closing: Vec<bool> = 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<f64>) = 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<u32> = 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, &params).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 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<usize> = 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::<f64>();
+    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 {
+                let x: Vec<bool> = (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_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_in(count, seed, limit, precision).unwrap();
+            let el = t.elapsed();
+            let mean = samples.iter().map(|(_, p)| scale * p).sum::<f64>() / 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::<String>(), 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);
+            }
+            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}"),
+    }
+}

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<String> = 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<u64> = raw.chunks_exact(8).map(|b| u64::from_le_bytes(b.try_into().unwrap())).collect();
+    let patches: Vec<Vec<u32>> = 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<u32>)> = 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::<f64>().expect("one angle")),
+            None => (head.to_string(), 0.0),
+        };
+        let qs: Vec<u32> = 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<Option<u32>> = qs.iter().map(|&q| local(q)).collect();
+            if l.iter().all(Option::is_none) {
+                continue;
+            }
+            let l: Vec<u32> = 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::<f64>() - 1.0;
+        // Each shot's patch bits, as the state vector's index (local qubit 0 most significant).
+        let vals: Vec<f64> = 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::<f64>() / vals.len() as f64;
+        let var = vals.iter().map(|v| (v - mean).powi(2)).sum::<f64>() / (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());
+}

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<String> = 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::<f64>().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::<f64>().expect("one angle")),
+            None => (head, 0.0),
+        };
+        let qs: Vec<u32> = 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()
+    );
+}

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;

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<Vec<u64>>,
+    pub local: Vec<LocalClifford>,
+}
+
 /// The first set bit of `v` in `from..to`.
 fn first_set(v: &[u64], from: usize, to: usize) -> Option<usize> {
     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, 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 {
+                if i != k2 && get(s.zr(i), q) {
+                    s.mul_into(i, k2);
+                }
+            }
+            *hq = 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.

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",

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",