wai-quantum 0.6.6: Clifford circuits measured in a product basis, exactly, by cut-rank #340

Merged
dcharlot merged 5 commits from wai/quantum-cutrank-0.6.6 into main 2026-10-19 20:26:27 -04:00
18 changed files with 1662 additions and 11 deletions

View file

@ -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

View file

@ -1017,7 +1017,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a"
[[package]]
name = "wai-quantum"
version = "0.6.5"
version = "0.6.6"
dependencies = [
"blake3",
"ed25519-dalek",

View file

@ -599,7 +599,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a"
[[package]]
name = "wai-quantum"
version = "0.6.5"
version = "0.6.6"
dependencies = [
"blake3",
"ed25519-dalek",

View file

@ -713,7 +713,7 @@ dependencies = [
[[package]]
name = "wai-quantum"
version = "0.6.5"
version = "0.6.6"
dependencies = [
"blake3",
"ed25519-dalek",

View file

@ -422,7 +422,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a"
[[package]]
name = "wai-quantum"
version = "0.6.5"
version = "0.6.6"
dependencies = [
"blake3",
"ed25519-dalek",

View file

@ -599,7 +599,7 @@ dependencies = [
[[package]]
name = "wai-quantum"
version = "0.6.5"
version = "0.6.6"
dependencies = [
"blake3",
"ed25519-dalek",

View file

@ -857,7 +857,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a"
[[package]]
name = "wai-quantum"
version = "0.6.5"
version = "0.6.6"
dependencies = [
"blake3",
"ed25519-dalek",

View file

@ -1079,7 +1079,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a"
[[package]]
name = "wai-quantum"
version = "0.6.5"
version = "0.6.6"
dependencies = [
"blake3",
"block2",

View file

@ -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"]

View file

@ -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

View file

@ -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}"),
}
}

View file

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

View file

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

View file

@ -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;

File diff suppressed because it is too large Load diff

View file

@ -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.

2
wai/wai-rs/Cargo.lock generated
View file

@ -4310,7 +4310,7 @@ version = "0.1.0"
[[package]]
name = "wai-quantum"
version = "0.6.5"
version = "0.6.6"
dependencies = [
"blake3",
"ed25519-dalek",

View file

@ -1241,7 +1241,7 @@ version = "0.1.0"
[[package]]
name = "wai-quantum"
version = "0.6.5"
version = "0.6.6"
dependencies = [
"blake3",
"ed25519-dalek",