From 7152e2a9cc8c609975884b4c6de7a2faee9f9c64 Mon Sep 17 00:00:00 2001 From: Transaction Science Date: Mon, 19 Oct 2026 06:31:58 -0400 Subject: [PATCH 1/3] wai-quantum 0.6.5 and the service: programs that measure in the middle MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The library: - Tableau::measure measures a qubit of the wide tableau and leaves the state the outcome selects; a random outcome is the caller's to choose, so a branch can be followed on purpose. The anticommuting stabilizer is multiplied into every other row at once, each row's phase summed mod 4 in two bit planes. Tableau::is_random answers in one column scan. - StateVector::probability_one (fixed-chunk sum) and StateVector::collapse. - Tests: against the state vector's projection over 60 random circuits of six measurements each; on wide tableaux, measuring every qubit lands in the support held before, one random outcome per dimension. Mutants caught: signs kept, Y phase swapped, destabilizer not moved. The service: - reads measurement and reset anywhere, if/else on measured bits (bit, !bit, ==, !=, OpenQASM 2's if(c==n)), and `bit b = measure q;`. A program that needs them is kept as steps (qasm::Instr); every other program is exactly what it was. - runs such a program on statevector or stabilizer (src/dynamic.rs): counts split their shots at each measurement and reset as independent shots would, each outcome running on in its own copy; expectation values are the exact mixture over the branches, of the state before the final measurements. Noise in such a program, amplitudes of one, and loops are refused with their reasons. - Tests: teleportation ( = cos θ, = sin θ exactly; the tableau and the state vector agree), reset and re-measurement, a bit-flip syndrome that corrects the flip, the step structure with its jumps, the planner. Mutants caught: terminal measurements made, a branch's outcome swapped, reset without its flip, condition bits reversed. Lockfiles: wai-quantum 0.6.5 in every workspace that pins it; cargo's rewrite of wai-rs's also re-picks two Windows-only edges (windows-sys 0.52, getrandom 0.3), both already in the lock. Co-Authored-By: Claude Opus 5.5 --- wai/CHANGELOG.md | 30 ++ 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-service/README.md | 22 +- wai/wai-quantum-service/src/catalog.rs | 76 +++++ wai/wai-quantum-service/src/discovery.rs | 2 +- wai/wai-quantum-service/src/dynamic.rs | 339 +++++++++++++++++++++++ wai/wai-quantum-service/src/main.rs | 1 + wai/wai-quantum-service/src/mcp.rs | 2 +- wai/wai-quantum-service/src/qasm.rs | 261 +++++++++++++++-- wai/wai-quantum-service/src/run.rs | 142 ++++++++++ wai/wai-quantum/Cargo.lock | 2 +- wai/wai-quantum/Cargo.toml | 2 +- wai/wai-quantum/src/quantum_sv.rs | 42 +++ wai/wai-quantum/src/quantum_tableau.rs | 183 ++++++++++++ wai/wai-rs/Cargo.lock | 14 +- wai/wai-web/Cargo.lock | 2 +- 21 files changed, 1089 insertions(+), 43 deletions(-) create mode 100644 wai/wai-quantum-service/src/dynamic.rs diff --git a/wai/CHANGELOG.md b/wai/CHANGELOG.md index 348697874..31e7830ee 100644 --- a/wai/CHANGELOG.md +++ b/wai/CHANGELOG.md @@ -6,6 +6,36 @@ wire format may change and such changes are recorded here. ## v0.4.0 — current +### 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 +runs exactly as before. + +- **`Tableau::measure`** (`quantum_tableau`) measures a qubit of the wide tableau and leaves + the state the outcome selects. A random outcome is the caller's to choose, so a branch + can be followed on purpose. The anticommuting stabilizer is multiplied into every other + row at once, each row's phase summed mod 4 in two bit planes. `Tableau::is_random` says + in one column scan whether an outcome is random. + - Tested against the state vector's projection, Pauli string by Pauli string, over 60 + random circuits of six measurements each. + - On wide tableaux, measuring every qubit lands inside the support held before, with one + random outcome per dimension. + - Three mutants are caught: signs kept, the Y phase swapped, the destabilizer not moved. +- **`StateVector::probability_one` and `StateVector::collapse`** (`quantum_sv`): a qubit's + probability of reading 1, summed in fixed chunks, and the state a reading leaves. +- **The service runs dynamic programs.** It reads measurement and `reset` anywhere, + `if (c[0]) …`, `if (!c[0]) …`, `if (c == 3) { … } else { … }`, OpenQASM 2's + `if(c==n) gate;`, and `bit b = measure q;`. + - Such a program runs on `statevector` or `stabilizer`. Its shots split at each + measurement and reset as independent shots would, and each outcome runs on in its own + copy of the state. + - Expectation values are the exact mixture over the branches, of the state before the + final measurements. + - Tests: teleportation gives ⟨Z⟩ = cos θ and ⟨X⟩ = sin θ exactly. It agrees between the + tableau and the state vector, reset and re-measurement behave, and a bit-flip syndrome + corrects the flip. Four mutants are caught. + - Loops, noise in such a program, and amplitudes of one are refused with their reasons. + ### wai-quantum-service: every job's energy, signed with its result - **What a result carries.** `energy` holds the CPU energy of the job's worker process, from diff --git a/wai/tools/quantum-meter/Cargo.lock b/wai/tools/quantum-meter/Cargo.lock index aed470f60..904f1b44e 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.4" +version = "0.6.5" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/tools/quantum-wgsl-check/Cargo.lock b/wai/tools/quantum-wgsl-check/Cargo.lock index 0398234a3..e9f467917 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.4" +version = "0.6.5" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-pack/Cargo.lock b/wai/wai-pack/Cargo.lock index 686cab658..a0356395f 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.4" +version = "0.6.5" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-quantum-qpu/Cargo.lock b/wai/wai-quantum-qpu/Cargo.lock index f7017f72d..143d92e94 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.4" +version = "0.6.5" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-quantum-serve/Cargo.lock b/wai/wai-quantum-serve/Cargo.lock index dbefe1eff..2ac649c7d 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.4" +version = "0.6.5" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-quantum-service/Cargo.lock b/wai/wai-quantum-service/Cargo.lock index 6a19a8d79..55d4fabf3 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.4" +version = "0.6.5" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-quantum-service/README.md b/wai/wai-quantum-service/README.md index 725935a10..cf6393608 100644 --- a/wai/wai-quantum-service/README.md +++ b/wai/wai-quantum-service/README.md @@ -28,7 +28,7 @@ A job is JSON: Counts are keyed by the program's classical bits, the highest bit first. With no measurement in the program, every qubit is measured into the bit of its own index. -The program may use the standard gates, user `gate` definitions, and the control modifiers `ctrl @`, `ctrl(n) @`, `negctrl @` and `negctrl(n) @` (on `x y z h s sdg t tdg p` and `cx cy cz ch cp ccx ccz`). Classical control, `inv @`, `pow @` and gates after a qubit's measurement are refused, each with its line and what to write instead (`src/qasm.rs` lists everything read). +The program may use the standard gates, user `gate` definitions, and the control modifiers `ctrl @`, `ctrl(n) @`, `negctrl @` and `negctrl(n) @` (on `x y z h s sdg t tdg p` and `cx cy cz ch cp ccx ccz`). Measurement and `reset` may come anywhere, and `if` acts on measured bits (`if (c[0]) x q[1];`, `if (c == 3) { … } else { … }`, OpenQASM 2's `if(c==1) x q;`). Loops, `inv @` and `pow @` are refused, each with its line and what to write instead (`src/qasm.rs` lists everything read). ## Methods @@ -51,6 +51,26 @@ How the methods run: - **stabilizer**, with noise, samples by Pauli frames against one noiseless outcome of the tableau, 64 shots to a machine word. Its counts have the law of independent noisy shots. Its expectation values are estimates with standard errors, from at least 10,000 shots unless `trajectories` sets the number. Amplitude damping is not a Pauli channel and is refused. - **trajectories** treats every shot as an independent noisy trajectory in law. Shots that took the same errors share one state, so the error-free part of a circuit, and the error-free majority of the shots, run once. `histories` in the result counts the distinct error histories that ran. +### Programs that measure in the middle + +A program that measures a qubit and then acts on it, resets a qubit after its first gate, or uses `if` on measured bits runs on `statevector` or `stabilizer`. `auto` picks the stabilizer tableau when every gate is Clifford, and the state vector otherwise. + +- **Counts** split their shots. The shots enter together, and at each measurement or reset they divide between the outcomes as independent shots would: binomially, in the outcome's probability. Each outcome then runs on in its own copy of the state. Shots that read the same outcomes share one state, and the counts have the law of independent runs. +- **Expectation values** are exact: the mixture over every branch of non-zero probability, at most 4096 branches. As for any program, they describe the state before the final measurements. A measurement whose qubit nothing touches afterwards, and whose bit no later condition reads, is not made for them. +- The result's `dynamic` member reports what the walks did: + - `measurements_and_resets`; + - `branches_run` and `states_at_once` for the counts; + - `branches_weighed` for the expectation values. + +Noise in such a program, and amplitudes of one, are refused, each with its reason. The state vector holds a copy for each branch still running, so it takes one qubit fewer than for a program that measures only at the end. + +```sh +# Teleportation: ⟨Z2⟩ = cos 0.7 and ⟨X2⟩ = sin 0.7, exactly. +curl -s https://wai.transaction.science/quantum-jobs/v1/run -H 'content-type: application/json' -d '{ + "qasm": "OPENQASM 3.0; qubit[3] q; bit[3] c; ry(0.7) q[0]; h q[1]; cx q[1], q[2]; cx q[0], q[1]; h q[0]; c[0] = measure q[0]; c[1] = measure q[1]; if (c[1]) x q[2]; if (c[0]) z q[2]; c[2] = measure q[2];", + "observables": ["Z2", "X2"], "shots": 1000 }' +``` + `auto` picks the first method that takes the job. With noise, it tries `density_matrix`, then `stabilizer`, then `trajectories`. For counts, it tries `stabilizer`, `exact`, `statevector`, `statevector_single`, then `mps`. For expectations alone, it tries `stabilizer`, `exact`, `statevector`, `statevector_single`, `pauli`, then `mps`. For amplitudes alone, it tries `exact`, `statevector`, `statevector_single`, `tensor_network`, then `mps`. The job's `why` names the method chosen and the ones that could not take it. The accuracy classes: diff --git a/wai/wai-quantum-service/src/catalog.rs b/wai/wai-quantum-service/src/catalog.rs index 90b816e75..eb3711fc8 100644 --- a/wai/wai-quantum-service/src/catalog.rs +++ b/wai/wai-quantum-service/src/catalog.rs @@ -261,6 +261,53 @@ pub fn mps_work(r: &Request, chi: u64) -> u64 { } /// Whether `method` can take the job, and the plan if so; otherwise why not. +/// A program that measures in the middle, resets after a qubit's first gate, +/// or acts on measured bits ([`crate::dynamic`]): the dense state vector and +/// the stabilizer tableau run it, splitting its shots at each measurement and +/// keeping a copy of the state for each branch that runs on. +fn try_dynamic(method: &str, r: &Request, clifford: &Result, String>, cap: u64, steps: &[crate::qasm::Instr]) -> Result { + let n = r.program.n_qubits; + let gates = r.program.ops.len() as u64; + let measured = crate::dynamic::measurements(steps); + if r.noise.is_some() { + return Err("noise in a program that measures in the middle, resets or acts on measured bits is not simulated yet".into()); + } + if !r.amplitudes.is_empty() { + return Err("a program that measures in the middle leaves no single state to give amplitudes of; ask for counts or expectation values".into()); + } + let plan = |memory_bytes: u64, work: u64, why: String| Plan { method: method.into(), why, memory_bytes, work, qubits: n }; + match method { + "statevector" => { + let one = statevector_bytes(n); + if n > wai_quantum::quantum_sv::MAX_QUBITS || 2 * one > cap { + return Err(format!("statevector runs a program that measures in the middle on at most {} qubits: a branch holds a second state", statevector_max_qubits(cap / 2))); + } + let states = (cap / one).clamp(2, 16); + Ok(plan(states * one + r.shots * 24, (gates << n).saturating_mul(r.shots.clamp(1, 64)), format!("dense double-precision state vector of {n} qubits, its shots split at each of {measured} mid-circuit measurements and resets"))) + } + "stabilizer" => { + if let Err(e) = clifford { + return Err(format!("stabilizer needs a Clifford circuit: {e}")); + } + if n > STABILIZER_MAX_QUBITS { + return Err(format!("stabilizer holds at most {STABILIZER_MAX_QUBITS} qubits")); + } + let nn = n as u64; + // A tableau: X and Z bits for 2n rows on each qubit. + let one = nn * nn / 2 + 4096; + let states = (cap.saturating_sub(64 * MIB) / one).clamp(2, 64); + // Each gate a pass over a column pair; a random measurement a pass + // over the tableau; a determined one an elimination. + let work = gates.saturating_mul(2 * nn.div_ceil(64)).saturating_mul(r.shots.clamp(1, 64)) + measured as u64 * nn * nn; + if work > STABILIZER_MAX_WORK { + return Err(format!("this program is more stabilizer work than one job may take ({work} word operations, at most {STABILIZER_MAX_WORK})")); + } + Ok(plan(states * one + 64 * MIB + r.shots.min(1 << 20) * (r.program.n_bits as u64 / 8 + 64), work, format!("Clifford circuit: the stabilizer tableau of {n} qubits, its shots split at each of {measured} mid-circuit measurements and resets"))) + } + _ => Err(format!("{method} does not run programs that measure in the middle, reset, or act on measured bits; statevector and stabilizer do")), + } +} + fn try_method(method: &str, r: &Request, dyadic: &Result, clifford: &Result, String>, cap: u64) -> Result { let n = r.program.n_qubits; let gates = r.program.ops.len() as u64; @@ -271,6 +318,9 @@ fn try_method(method: &str, r: &Request, dyadic: &Result { if let Err(e) = dyadic { @@ -513,6 +563,32 @@ mod tests { assert!(plan(&job(rot, json!({ "method": "nonsense" })), CAP).unwrap_err().contains("unknown method")); } + /// A program that measures in the middle goes to the tableau when it is + /// Clifford and to the state vector otherwise; every other method, noise + /// and amplitudes are refused with what would run it. + #[test] + fn programs_that_measure_in_the_middle_are_planned_for_the_methods_that_branch() { + let clifford = "qubit[2] q; bit[2] c; h q[0]; c[0] = measure q[0]; if (c[0]) x q[1]; c[1] = measure q[1];"; + let rotated = "qubit[2] q; bit[2] c; rx(0.3) q[0]; c[0] = measure q[0]; if (c[0]) x q[1]; c[1] = measure q[1];"; + let p = plan(&job(clifford, json!({})), CAP).unwrap(); + assert_eq!(p.method, "stabilizer"); + assert!(p.why.contains("split at each of 2 mid-circuit measurements"), "{}", p.why); + let p = plan(&job(rotated, json!({ "observables": ["Z1"] })), CAP).unwrap(); + assert_eq!(p.method, "statevector"); + assert!(p.memory_bytes >= 2 * statevector_bytes(2), "room for a branch's copy"); + for (extra, needle) in [ + (json!({ "method": "mps" }), "statevector and stabilizer do"), + (json!({ "method": "exact" }), "statevector and stabilizer do"), + (json!({ "noise": { "depolarizing_1q": 0.01 } }), "not simulated yet"), + (json!({ "amplitudes": ["00"] }), "no single state"), + ] { + let e = plan(&job(rotated, extra.clone()), CAP).unwrap_err(); + assert!(e.contains(needle), "{extra}: {e}"); + } + let e = plan(&job("qubit[31] q; bit c; h q; c = measure q[0]; x q[0];", json!({ "method": "statevector" })), 40 * GIB).unwrap_err(); + assert!(e.contains("a branch holds a second state"), "{e}"); + } + #[test] fn the_statevector_limit_follows_the_job_memory() { assert_eq!(statevector_max_qubits(72 * GIB), 32); diff --git a/wai/wai-quantum-service/src/discovery.rs b/wai/wai-quantum-service/src/discovery.rs index 3f0719fbf..80b2cdf9d 100644 --- a/wai/wai-quantum-service/src/discovery.rs +++ b/wai/wai-quantum-service/src/discovery.rs @@ -95,7 +95,7 @@ pub fn openapi(server: &str) -> Value { "type": "object", "required": ["qasm"], "properties": { - "qasm": { "type": "string", "description": "The program: OpenQASM 3, or OpenQASM 2. Registers, standard gates (h, x, y, z, s, sdg, t, tdg, sx, rx, ry, rz, p, u, cx, cy, cz, ch, swap, iswap, cp, crx, cry, crz, ccx, ccz, cswap, rxx, ryy, rzz), the control modifiers ctrl @, ctrl(n) @ and negctrl @ on x y z h s sdg t tdg p and the controlled names, user gate definitions and terminal measurement. Mid-circuit measurement, classical control, inv @ and pow @ are refused, with the line." }, + "qasm": { "type": "string", "description": "The program: OpenQASM 3, or OpenQASM 2. Registers, standard gates (h, x, y, z, s, sdg, t, tdg, sx, rx, ry, rz, p, u, cx, cy, cz, ch, swap, iswap, cp, crx, cry, crz, ccx, ccz, cswap, rxx, ryy, rzz), the control modifiers ctrl @, ctrl(n) @ and negctrl @ on x y z h s sdg t tdg p and the controlled names, user gate definitions, measurement anywhere, reset anywhere, and `if (c[0]) …` / `if (c == 3) { … } else { … }` on measured bits (OpenQASM 2's `if(c==n) gate;` too). A program that measures in the middle, resets after a gate or acts on measured bits runs on statevector or stabilizer: counts split their shots at each measurement, expectation values are the exact mixture over the branches, of the state before the final measurements. Loops, inv @ and pow @ are refused, with the line." }, "method": { "type": "string", "enum": ["auto", "exact", "statevector", "statevector_single", "stabilizer", "pauli", "mps", "tensor_network", "density_matrix", "trajectories"], "default": "auto" }, "shots": { "type": "integer", "minimum": 0, "maximum": 1000000, "description": "Measurement shots (default 1024 when nothing else is asked for)." }, "seed": { "type": "integer", "default": 1, "description": "The sampling seed: the same job and seed give the same counts." }, diff --git a/wai/wai-quantum-service/src/dynamic.rs b/wai/wai-quantum-service/src/dynamic.rs new file mode 100644 index 000000000..659e275a0 --- /dev/null +++ b/wai/wai-quantum-service/src/dynamic.rs @@ -0,0 +1,339 @@ +//! Dynamic programs: measurement in the middle of a circuit, reset after a +//! qubit's first gate, and gates that act on what was measured. +//! +//! The program runs as its steps ([`Instr`]). **Counts** come from shot +//! splitting: the shots enter together, and at each measurement or reset they +//! divide between the two outcomes as independent shots would, binomially in +//! the outcome's probability. Each outcome then runs on in its own copy of the +//! state, the one with more shots in the state itself. So shots that read the +//! same outcomes share one state, an outcome no shot reads costs nothing, and +//! the counts have the law of `shots` independent runs. **Expectation values** +//! are exact: the mixture over every branch of non-zero probability, each +//! branch's final state weighted by its probability. As for a program that +//! measures only at the end, they are of the state before the final +//! measurements: a measurement whose qubit nothing touches afterwards and +//! whose bit no later condition reads is not made for them ([`terminal`]). +//! +//! Two methods run them: the dense state vector (any gate) and the stabilizer +//! tableau (Clifford gates, where a random outcome has probability ½). + +use crate::qasm::{Instr, Op}; +use std::collections::HashMap; +use wai_quantum::gates::Gate as Matrix; +use wai_quantum::quantum_sv::StateVector; +use wai_quantum::quantum_tableau::{Clifford, Tableau}; + +/// The most branches the exact expectation walk follows. +pub const MAX_LEAVES: u64 = 4096; +/// Below this probability a state-vector outcome is rounding, not a branch, +/// for the exact expectation walk. (Shot splitting needs no threshold: a +/// negligible outcome draws no shots.) +const NEGLIGIBLE: f64 = 1e-12; + +/// A state a dynamic program runs on. +pub trait Engine: Clone { + /// Apply `ops[i]`. + fn gate(&mut self, i: usize) -> Result<(), String>; + /// The probability that measuring `q` reads 1. + fn prob_one(&mut self, q: u32) -> Result; + /// Leave the state a measurement of `q` reading `outcome` leaves. + fn collapse(&mut self, q: u32, outcome: bool) -> Result<(), String>; + /// `X` on `q` (a reset that read 1). + fn flip(&mut self, q: u32) -> Result<(), String>; + /// `⟨P⟩` in the current state. + fn expectation(&self, string: &[(u32, char)]) -> Result; +} + +/// The dense double-precision state vector. +#[derive(Clone)] +pub struct Dense<'a> { + pub sv: StateVector, + pub gates: &'a [Matrix], +} + +impl Engine for Dense<'_> { + fn gate(&mut self, i: usize) -> Result<(), String> { + self.sv.apply(&self.gates[i]).map_err(|e| format!("{e:?}")) + } + fn prob_one(&mut self, q: u32) -> Result { + Ok(self.sv.probability_one(q).clamp(0.0, 1.0)) + } + fn collapse(&mut self, q: u32, outcome: bool) -> Result<(), String> { + match self.sv.collapse(q, outcome).map_err(|e| format!("{e:?}"))? { + p if p > 0.0 => Ok(()), + _ => Err(format!("qubit {q} cannot read {}", u8::from(outcome))), + } + } + fn flip(&mut self, q: u32) -> Result<(), String> { + self.sv.apply(&Matrix::x(q)).map_err(|e| format!("{e:?}")) + } + fn expectation(&self, string: &[(u32, char)]) -> Result { + Ok(self.sv.expectation(&[(1.0, string.to_vec())])) + } +} + +/// The stabilizer tableau, each op as its Clifford gates. +#[derive(Clone)] +pub struct Stabilizer<'a> { + pub t: Tableau, + pub ops: &'a [Vec], +} + +impl Engine for Stabilizer<'_> { + fn gate(&mut self, i: usize) -> Result<(), String> { + for &c in &self.ops[i] { + self.t.apply(c).map_err(|e| e.to_string())?; + } + Ok(()) + } + fn prob_one(&mut self, q: u32) -> Result { + if self.t.is_random(q).map_err(|e| e.to_string())? { + return Ok(0.5); + } + Ok(if self.t.expectation(&[(q, 'Z')]).map_err(|e| e.to_string())? < 0 { 1.0 } else { 0.0 }) + } + fn collapse(&mut self, q: u32, outcome: bool) -> Result<(), String> { + let (got, _) = self.t.measure(q, || outcome).map_err(|e| e.to_string())?; + if got == outcome { Ok(()) } else { Err(format!("qubit {q} cannot read {}", u8::from(outcome))) } + } + fn flip(&mut self, q: u32) -> Result<(), String> { + self.t.apply(Clifford::X(q)).map_err(|e| e.to_string()) + } + fn expectation(&self, string: &[(u32, char)]) -> Result { + self.t.expectation(string).map(f64::from).map_err(|e| e.to_string()) + } +} + +/// What a sampling walk found. +pub struct Sampled { + /// Each outcome of the classical bits (bit `n_bits − 1` first) and its shots. + pub counts: HashMap, + /// Distinct branches run to the end. + pub branches: u64, + /// The most states alive at once. + pub peak_states: usize, +} + +struct Walk<'a> { + steps: &'a [Instr], + n_bits: usize, + rng: u64, + max_states: usize, + live: usize, + peak: usize, + /// Measurements the expectation walk does not make ([`terminal`]). + skip: Vec, +} + +impl Walk<'_> { + fn next(&mut self) -> u64 { + self.rng = self.rng.wrapping_add(0x9e37_79b9_7f4a_7c15); + let mut z = self.rng; + z = (z ^ (z >> 30)).wrapping_mul(0xbf58_476d_1ce4_e5b9); + z = (z ^ (z >> 27)).wrapping_mul(0x94d0_49bb_1331_11eb); + z ^ (z >> 31) + } + + /// How many of `shots` independent shots read 1, at probability `p`: + /// one draw a shot, so the law is exactly binomial. + fn binomial(&mut self, shots: u64, p: f64) -> u64 { + if p <= 0.0 { + return 0; + } + if p >= 1.0 { + return shots; + } + let mut k = 0; + for _ in 0..shots { + if ((self.next() >> 11) as f64 + 0.5) / 9_007_199_254_740_992.0 < p { + k += 1; + } + } + k + } + + fn copy(&mut self, state: &E) -> Result { + if self.live >= self.max_states { + return Err(format!( + "this program's measurements branch into more states at once than the job's memory holds ({}); measure fewer qubits in the middle of the circuit, or use fewer qubits", + self.max_states + )); + } + self.live += 1; + self.peak = self.peak.max(self.live); + Ok(state.clone()) + } + + /// Run `shots` from step `pc` to the end, in `state`, with the bits `bits`. + fn sample(&mut self, mut pc: usize, mut state: E, mut bits: Vec, mut shots: u64, out: &mut Sampled) -> Result<(), String> { + while pc < self.steps.len() { + match &self.steps[pc] { + Instr::Gate(i) => state.gate(*i)?, + Instr::Measure { qubit, .. } | Instr::Reset(qubit) => { + let q = *qubit; + let (bit, reset) = step_target(&self.steps[pc]); + let p1 = state.prob_one(q)?; + let ones = self.binomial(shots, p1); + // The outcome more shots read carries on in this state. + let (more, fewer) = if ones * 2 > shots { (true, shots - ones) } else { (false, ones) }; + if fewer > 0 { + // The smaller share runs first, in a copy. + let mut child = self.copy(&state)?; + settle(&mut child, q, !more, reset)?; + let mut child_bits = bits.clone(); + if let Some(b) = bit { + child_bits[b] = !more; + } + self.sample(pc + 1, child, child_bits, fewer, out)?; + self.live -= 1; + } + shots -= fewer; + settle(&mut state, q, more, reset)?; + if let Some(b) = bit { + bits[b] = more; + } + } + Instr::JumpUnless(cond, to) => { + if !cond.holds(&bits) { + pc = *to; + continue; + } + } + Instr::Jump(to) => { + pc = *to; + continue; + } + } + pc += 1; + } + let key: String = (0..self.n_bits).rev().map(|b| if bits[b] { '1' } else { '0' }).collect(); + *out.counts.entry(key).or_insert(0) += shots; + out.branches += 1; + Ok(()) + } + + /// The exact mixture's expectation values from step `pc`, the branch + /// weighted by `weight`, added into `sums`; `leaves` counts the branches. + fn expect(&mut self, mut pc: usize, mut state: E, mut bits: Vec, weight: f64, observables: &[Vec<(u32, char)>], sums: &mut [f64], leaves: &mut u64) -> Result<(), String> { + while pc < self.steps.len() { + match &self.steps[pc] { + Instr::Gate(i) => state.gate(*i)?, + Instr::Measure { .. } if self.skip.get(pc) == Some(&true) => {} + Instr::Measure { qubit, .. } | Instr::Reset(qubit) => { + let q = *qubit; + let (bit, reset) = step_target(&self.steps[pc]); + let p1 = state.prob_one(q)?; + let branches: Vec<(bool, f64)> = [(false, 1.0 - p1), (true, p1)].into_iter().filter(|&(_, p)| p > NEGLIGIBLE).collect(); + if let [(o, p), _] = branches[..] { + // Outcome 0 in a copy; outcome 1 in this state. + let mut child = self.copy(&state)?; + settle(&mut child, q, o, reset)?; + let mut child_bits = bits.clone(); + if let Some(b) = bit { + child_bits[b] = o; + } + self.expect(pc + 1, child, child_bits, weight * p, observables, sums, leaves)?; + self.live -= 1; + } + let (o, p) = *branches.last().ok_or("a measurement with no possible outcome")?; + settle(&mut state, q, o, reset)?; + if let Some(b) = bit { + bits[b] = o; + } + if branches.len() == 2 { + return self.expect(pc + 1, state, bits, weight * p, observables, sums, leaves); + } + } + Instr::JumpUnless(cond, to) => { + if !cond.holds(&bits) { + pc = *to; + continue; + } + } + Instr::Jump(to) => { + pc = *to; + continue; + } + } + pc += 1; + } + *leaves += 1; + if *leaves > MAX_LEAVES { + return Err(format!( + "this program's measurements branch {MAX_LEAVES}+ ways, more than the exact expectation values follow; ask for counts, whose shots split only where they fall" + )); + } + for (s, o) in sums.iter_mut().zip(observables) { + *s += weight * state.expectation(o)?; + } + Ok(()) + } +} + +/// `shots` shots of `steps` from `start`, split at each measurement and reset. +pub fn sample(steps: &[Instr], n_bits: u32, start: E, shots: u64, seed: u64, max_states: usize) -> Result { + let mut walk = Walk { steps, n_bits: n_bits as usize, rng: seed, max_states: max_states.max(1), live: 1, peak: 1, skip: Vec::new() }; + let mut out = Sampled { counts: HashMap::new(), branches: 0, peak_states: 1 }; + if shots > 0 { + walk.sample(0, start, vec![false; n_bits as usize], shots, &mut out)?; + } + out.peak_states = walk.peak; + Ok(out) +} + +/// The exact expectation values of `observables` over the program's +/// branches, before its final measurements, and how many branches there were. +pub fn expectations(steps: &[Instr], ops: &[Op], n_bits: u32, start: E, observables: &[Vec<(u32, char)>], max_states: usize) -> Result<(Vec, u64), String> { + let mut walk = Walk { steps, n_bits: n_bits as usize, rng: 0, max_states: max_states.max(1), live: 1, peak: 1, skip: terminal(steps, ops) }; + let mut sums = vec![0.0; observables.len()]; + let mut leaves = 0; + walk.expect(0, start, vec![false; n_bits as usize], 1.0, observables, &mut sums, &mut leaves)?; + Ok((sums, leaves)) +} + +/// For each step, whether it is a final measurement: no later step touches +/// its qubit (a gate, a reset, another measurement) or reads its bit (a +/// condition). Steps only jump forward, so "later" is "after it as written". +pub fn terminal(steps: &[Instr], ops: &[Op]) -> Vec { + let mut used = std::collections::HashSet::new(); + let mut read = std::collections::HashSet::new(); + let mut out = vec![false; steps.len()]; + for (i, step) in steps.iter().enumerate().rev() { + match step { + Instr::Gate(j) => used.extend(ops[*j].qubits.iter().copied()), + Instr::Reset(q) => { + used.insert(*q); + } + Instr::Measure { qubit, bit } => { + out[i] = !used.contains(qubit) && !read.contains(bit); + used.insert(*qubit); + } + Instr::JumpUnless(cond, _) => read.extend(cond.bits.iter().copied()), + Instr::Jump(_) => {} + } + } + out +} + +/// The bit a step writes (a measurement's), and whether it is a reset. +fn step_target(step: &Instr) -> (Option, bool) { + match step { + Instr::Measure { bit, .. } => (Some(*bit as usize), false), + _ => (None, true), + } +} + +/// The state a measurement (or reset) of `q` reading `outcome` leaves: a +/// reset that read 1 turns the qubit back to `|0⟩`. +fn settle(state: &mut E, q: u32, outcome: bool, reset: bool) -> Result<(), String> { + state.collapse(q, outcome)?; + if reset && outcome { + state.flip(q)?; + } + Ok(()) +} + +/// The measurements and resets a program makes, as written. +pub fn measurements(steps: &[Instr]) -> usize { + steps.iter().filter(|s| matches!(s, Instr::Measure { .. } | Instr::Reset(_))).count() +} diff --git a/wai/wai-quantum-service/src/main.rs b/wai/wai-quantum-service/src/main.rs index 05d98f35c..42608b028 100644 --- a/wai/wai-quantum-service/src/main.rs +++ b/wai/wai-quantum-service/src/main.rs @@ -15,6 +15,7 @@ mod catalog; mod discovery; +mod dynamic; mod energy; mod http; mod jobs; diff --git a/wai/wai-quantum-service/src/mcp.rs b/wai/wai-quantum-service/src/mcp.rs index 0cb4205e2..81b266551 100644 --- a/wai/wai-quantum-service/src/mcp.rs +++ b/wai/wai-quantum-service/src/mcp.rs @@ -24,7 +24,7 @@ const INSTRUCTIONS: &str = "Run quantum programs written in OpenQASM 3 (or 2). C fn tools() -> Value { let job_props = json!({ - "qasm": { "type": "string", "description": "The program, OpenQASM 3 or 2. Supported: qubit/bit registers and broadcast over them; gates h x y z s sdg t tdg sx sxdg rx ry rz p u cx cy cz ch swap iswap cp crx cry crz cu ccx ccz cswap rxx ryy rzz; the control modifiers `ctrl @`, `ctrl(n) @`, `negctrl @` on x y z h s sdg t tdg p and the controlled names; user `gate` definitions; angle expressions with pi; `reset` before a qubit's first gate; terminal measurement (`c = measure q;`). Not supported: mid-circuit measurement or reset, classical control, `inv @`, `pow @`." }, + "qasm": { "type": "string", "description": "The program, OpenQASM 3 or 2. Supported: qubit/bit registers and broadcast over them; gates h x y z s sdg t tdg sx sxdg rx ry rz p u cx cy cz ch swap iswap cp crx cry crz cu ccx ccz cswap rxx ryy rzz; the control modifiers `ctrl @`, `ctrl(n) @`, `negctrl @` on x y z h s sdg t tdg p and the controlled names; user `gate` definitions; angle expressions with pi; measurement and `reset` anywhere; classical control on measured bits (`if (c[0]) x q[1];`, `if (c == 3) { … } else { … }`, `if(c==1) x q;`). A program that measures in the middle or acts on what it measured runs on statevector or stabilizer: counts split their shots at each measurement; expectation values are exact, of the state before the final measurements. Not supported: loops, `inv @`, `pow @`." }, "shots": { "type": "integer", "minimum": 0, "maximum": 1000000, "description": "Measurement shots. Default 1024 when no observables or amplitudes are asked for." }, "observables": { "type": "array", "items": { "type": "string" }, "description": "Pauli strings to take expectation values of, e.g. [\"Z0\", \"Z0 Z1\", \"X2 Y3\"]. Qubits count from 0." }, "amplitudes": { "type": "array", "items": { "type": "string" }, "description": "Basis states to give amplitudes of: one character per qubit, qubit n-1 first, e.g. \"0011\"." }, diff --git a/wai/wai-quantum-service/src/qasm.rs b/wai/wai-quantum-service/src/qasm.rs index 456d6d574..1a4cbfb9d 100644 --- a/wai/wai-quantum-service/src/qasm.rs +++ b/wai/wai-quantum-service/src/qasm.rs @@ -6,19 +6,25 @@ //! gate calls with angle expressions (`pi`, `tau`, `euler`, `+ - * / ^`, and //! `sin cos tan exp ln sqrt arcsin arccos arctan`), register broadcast //! (`h q;`, `cx a, b;` over equal-length registers), user `gate` definitions -//! (inlined, with parameters), `barrier` (ignored), `reset` before a qubit's -//! first gate, and terminal measurement in every form (`c[0] = measure q[0];`, -//! `c = measure q;`, `measure q -> c;`, `measure q;`). With no measurement, -//! every qubit is measured into a bit of its own index. +//! (inlined, with parameters), `barrier` (ignored), `reset`, and measurement +//! in every form (`c[0] = measure q[0];`, `c = measure q;`, `measure q -> c;`, +//! `measure q;`, `bit b = measure q;`). With no measurement, every qubit is +//! measured into a bit of its own index. +//! +//! Classical control: `if (c[0]) …`, `if (!c[0]) …`, `if (c == 3) …`, +//! `if (c != 3) …`, with a `{ block }` or one statement, and `else`; and +//! OpenQASM 2's `if(c==n) gate;`. A program that measures a qubit and then +//! acts on it, resets one after its first gate, or uses `if` is dynamic: it +//! is kept as steps ([`Instr`]) as well as gates, and runs by +//! [`crate::dynamic`]. Every other program is exactly what it was before. //! //! Control modifiers, `ctrl @`, `ctrl(n) @`, `negctrl @` and `negctrl(n) @`, //! on `x y z h s sdg t tdg p` and on the controlled names `cx cy cz ch cp ccx //! ccz`: a gate the standard names is read as that gate, and any other as //! standard gates, exactly (see [`Parser::phase_n`]). //! -//! Refused, each with its line: classical control, the `inv` and `pow` -//! modifiers, a gate on a qubit after its measurement, and anything else not -//! listed. +//! Refused, each with its line: loops, the `inv` and `pow` modifiers, +//! declarations inside `if`, and anything else not listed. //! //! The standard gates: `id x y z h s sdg t tdg sx sxdg rx ry rz p phase u1 u2 //! u u3 U cx CX cnot cy cz ch cp cphase cu1 crx cry crz cu cu3 swap iswap rxx @@ -175,6 +181,45 @@ pub struct Program { /// `(qubit, bit)` for each measurement, in order. pub measured: Vec<(u32, u32)>, pub n_bits: u32, + /// The program as steps in order, when it measures in the middle, resets + /// a qubit after its first gate, or acts on what it measured: the gates + /// then run as [`Instr::Gate`] steps between the others. `None` for a + /// program whose measurements all come at the end, which runs as `ops` + /// followed by `measured`. + pub dynamic: Option>, +} + +/// A condition on classical bits: `bits` (the register's lowest first) read +/// as an unsigned integer, equal (or, with `equal` false, not equal) to +/// `value`. +#[derive(Clone, Debug, PartialEq, Eq)] +pub struct Cond { + pub bits: Vec, + pub value: u64, + pub equal: bool, +} + +impl Cond { + /// Whether the condition holds for the bits' current values. + pub fn holds(&self, bits: &[bool]) -> bool { + let v = self.bits.iter().enumerate().fold(0u64, |v, (i, &b)| v | (u64::from(bits[b as usize]) << i)); + (v == self.value) == self.equal + } +} + +/// One step of a dynamic program. +#[derive(Clone, Debug, PartialEq)] +pub enum Instr { + /// The gate `ops[i]`. + Gate(usize), + /// Measure `qubit` into `bit`. + Measure { qubit: u32, bit: u32 }, + /// Return `qubit` to `|0⟩`. + Reset(u32), + /// Go on to the next step if the condition holds, else to step `.1`. + JumpUnless(Cond, usize), + /// Go to step `.0`. + Jump(usize), } // --------------------------------------------------------------------------- @@ -255,6 +300,7 @@ fn lex(src: &str) -> Result, String> { "->" => Some("->"), "**" => Some("**"), "==" => Some("=="), + "!=" => Some("!="), _ => None, }; if let Some(s) = sym { @@ -279,6 +325,7 @@ fn lex(src: &str) -> Result, String> { '^' => "^", '@' => "@", ':' => ":", + '!' => "!", _ => return Err(format!("line {line}: unexpected character {c:?}")), }; out.push((Tok::Sym(s), line)); @@ -374,6 +421,11 @@ struct Parser { measured: Vec<(u32, u32)>, measured_qubits: Vec, touched: Vec, + /// The program as steps, kept for every program; used when it is dynamic. + steps: Vec, + dynamic: bool, + /// How many `if` blocks the parser is inside. + depth: usize, } impl Parser { @@ -562,6 +614,9 @@ impl Parser { "include" => { self.skip_to_semicolon(); } + "qubit" | "bit" | "qreg" | "creg" | "gate" if self.depth > 0 => { + return Err(format!("line {line}: declare `{word}` outside `if` blocks")); + } "qubit" | "bit" => { let size = if self.eat("[") { let n = self.integer()?; @@ -572,7 +627,16 @@ impl Parser { }; let name = self.ident()?; if self.eat("=") { - return self.err("initialised declarations are not supported"); + // `bit[n] c = measure q;` + if word != "bit" || !matches!(self.next(), Some(Tok::Ident(m)) if m == "measure") { + return Err(format!("line {line}: only a bit may be initialised, and only with `measure`")); + } + self.declare(name.clone(), size, true)?; + let qs = self.operand(false)?; + self.expect(";")?; + let (start, len) = self.cregs[&name]; + let bs: Vec = (start..start + len).collect(); + return self.measure(&qs, &bs, line); } self.expect(";")?; self.declare(name, size, word == "bit")?; @@ -590,9 +654,12 @@ impl Parser { "reset" => { let qs = self.operand(false)?; self.expect(";")?; - if qs.iter().any(|&q| self.touched[q as usize]) { - return Err(format!("line {line}: reset after a qubit's first gate is not supported (mid-circuit reset)")); + // Before a qubit's first gate a reset changes nothing; after + // it, or inside `if`, the program is dynamic. + if self.depth > 0 || qs.iter().any(|&q| self.touched[q as usize] || self.measured_qubits[q as usize]) { + self.dynamic = true; } + self.steps.extend(qs.iter().map(|&q| Instr::Reset(q))); } "measure" => { let qs = self.operand(false)?; @@ -615,10 +682,29 @@ impl Parser { } self.call_with(name, line, &polarity)?; } - "if" | "while" | "for" | "def" | "defcal" | "cal" | "box" | "let" | "const" | "input" | "output" | "inv" | "pow" | "delay" | "gphase" => { + "if" => { + // `if (c[0]) { … } else { … }`, `if (c == 3) …`, `if (!c[1]) …`. + self.dynamic = true; + self.expect("(")?; + let cond = self.condition(line)?; + self.expect(")")?; + let at = self.steps.len(); + self.steps.push(Instr::JumpUnless(cond.clone(), 0)); + self.block()?; + if matches!(self.peek(), Some(Tok::Ident(w)) if w == "else") { + self.pos += 1; + let jump = self.steps.len(); + self.steps.push(Instr::Jump(0)); + self.steps[at] = Instr::JumpUnless(cond, jump + 1); + self.block()?; + self.steps[jump] = Instr::Jump(self.steps.len()); + } else { + self.steps[at] = Instr::JumpUnless(cond, self.steps.len()); + } + } + "while" | "for" | "def" | "defcal" | "cal" | "box" | "let" | "const" | "input" | "output" | "inv" | "pow" | "delay" | "gphase" => { let hint = match word.as_str() { - "if" | "while" => "; classical control is not simulated here: measure every qubit at the end instead", - "for" => "; unroll the loop into gates", + "while" | "for" => "; unroll the loop into statements", "inv" | "pow" => "; write the inverse or powered gate out, or define it with `gate`", _ => "", }; @@ -645,16 +731,78 @@ impl Parser { Ok(()) } + /// The body of an `if` or `else`: `{ statements }` or one statement. + fn block(&mut self) -> Result<(), String> { + self.depth += 1; + if self.eat("{") { + while !self.eat("}") { + if self.pos >= self.toks.len() { + return self.err("expected `}`"); + } + self.statement()?; + } + } else { + self.statement()?; + } + self.depth -= 1; + Ok(()) + } + + /// An `if` condition: `c[0]`, `!c[0]`, `c[0] == 1`, `c == 5`, `c != 5`. + fn condition(&mut self, line: usize) -> Result { + let negate = self.eat("!"); + let name = self.ident()?; + let Some(&(start, len)) = self.cregs.get(&name) else { + return Err(format!("line {line}: `{name}` is not a bit register; a condition reads measured bits")); + }; + let bits: Vec = if self.eat("[") { + let i = self.integer()?; + self.expect("]")?; + if i >= len { + return Err(format!("line {line}: {name}[{i}] is outside {name}, which has {len}")); + } + vec![start + i] + } else { + (start..start + len).collect() + }; + if bits.len() > 64 { + return Err(format!("line {line}: a condition reads at most 64 bits")); + } + let compare = if self.eat("==") { + Some(true) + } else if self.eat("!=") { + Some(false) + } else { + None + }; + let (value, equal) = match compare { + Some(_) if negate => return Err(format!("line {line}: write `{name} != v` rather than negating a comparison")), + Some(eq) => { + let v = self.integer()? as u64; + if bits.len() < 64 && v >> bits.len() != 0 { + return Err(format!("line {line}: {v} does not fit in {} bits", bits.len())); + } + (v, eq) + } + None if bits.len() == 1 => (1, !negate), + None => return Err(format!("line {line}: compare the register `{name}` with a number, as `{name} == 3`")), + }; + Ok(Cond { bits, value, equal }) + } + fn measure(&mut self, qs: &[u32], bs: &[u32], line: usize) -> Result<(), String> { if qs.len() != bs.len() { return Err(format!("line {line}: {} qubits measured into {} bits", qs.len(), bs.len())); } for (&q, &b) in qs.iter().zip(bs) { if self.measured_qubits[q as usize] { - return Err(format!("line {line}: qubit {q} is measured twice")); + // Measured again: the program is dynamic. + self.dynamic = true; + } else { + self.measured_qubits[q as usize] = true; + self.measured.push((q, b)); } - self.measured_qubits[q as usize] = true; - self.measured.push((q, b)); + self.steps.push(Instr::Measure { qubit: q, bit: b }); } if self.cregs.is_empty() { self.n_bits = self.n_bits.max(bs.iter().map(|&b| b + 1).max().unwrap_or(0)); @@ -929,7 +1077,8 @@ impl Parser { return Err(format!("line {line}: qubit {q} named twice in one gate")); } if self.measured_qubits[*q as usize] { - return Err(format!("line {line}: a gate on qubit {q} after its measurement is not supported (mid-circuit measurement)")); + // A gate after a measurement: the program is dynamic. + self.dynamic = true; } } if let Some((kind, given)) = builtin(name) { @@ -949,6 +1098,7 @@ impl Parser { self.touched[q as usize] = true; } self.ops.push(Op { kind, qubits: qubits.to_vec(), params }); + self.steps.push(Instr::Gate(self.ops.len() - 1)); return Ok(()); } let Some(def) = self.defs.get(name).cloned() else { @@ -983,6 +1133,9 @@ pub fn parse(src: &str) -> Result { measured: Vec::new(), measured_qubits: Vec::new(), touched: Vec::new(), + steps: Vec::new(), + dynamic: false, + depth: 0, }; while p.pos < p.toks.len() { p.statement()?; @@ -991,11 +1144,14 @@ pub fn parse(src: &str) -> Result { return Err("the program declares no qubits".into()); } let (mut measured, mut n_bits) = (p.measured, p.n_bits); + let mut steps = p.steps; if measured.is_empty() { measured = (0..p.n_qubits).map(|q| (q, q)).collect(); n_bits = p.n_qubits; + steps.extend((0..p.n_qubits).map(|q| Instr::Measure { qubit: q, bit: q })); } - Ok(Program { n_qubits: p.n_qubits, ops: p.ops, measured, n_bits }) + let dynamic = p.dynamic.then_some(steps); + Ok(Program { n_qubits: p.n_qubits, ops: p.ops, measured, n_bits, dynamic }) } // --------------------------------------------------------------------------- @@ -1724,10 +1880,15 @@ mod tests { #[test] fn what_is_not_supported_is_refused_with_its_line() { for (src, needle) in [ - ("qubit[1] q;\nbit c;\nc = measure q;\nx q;", "line 4"), - ("qubit[1] q;\nif (c) x q;", "not supported"), - ("qubit[1] q;\nif (c) x q;", "measure every qubit at the end"), + ("qubit[1] q;\nif (c) x q;", "not a bit register"), ("qubit[1] q;\nfor uint i in [0:3] { x q; }", "unroll the loop"), + ("qubit[1] q;\nbit c;\nwhile (c) { x q; }", "unroll the loop"), + ("qubit[1] q;\nbit c;\nc = measure q;\nif (c) { qubit r; }", "outside `if` blocks"), + ("qubit[2] q;\nbit[2] c;\nc = measure q;\nif (c) x q[0];", "compare the register"), + ("qubit[2] q;\nbit[2] c;\nc = measure q;\nif (c == 4) x q[0];", "does not fit in 2 bits"), + ("qubit[1] q;\nbit c;\nif (!c == 1) x q;", "rather than negating"), + ("qubit[1] q;\nbit c;\nif (c) { x q;", "expected `}`"), + ("qubit[1] q;\nbit c = h q;", "only with `measure`"), ("qubit[1] q;\ninv @ x q[0];", "write the inverse"), ("qubit[2] q;\nctrl @ rz(0.1) q[0], q[1];", "a control modifier on `rz` is not supported"), ("qubit[3] q;\nctrl @ x q[0], q[1], q[2];", "takes 2 qubits"), @@ -1736,7 +1897,6 @@ mod tests { ("qubit[1] q;\nfoo q;", "unknown gate foo"), ("qubit[1] q;\nrx q;", "rx takes 1 parameters"), ("qubit[2] q;\ncx q[0], q[0];", "named twice"), - ("qubit[2] q;\nx q[0];\nreset q[0];", "mid-circuit reset"), ("qubit[2] q;\nh q[2];", "outside q"), ("bit c;", "no qubits"), ] { @@ -1745,6 +1905,59 @@ mod tests { } } + /// Measurement in the middle, reset after a gate and `if`/`else` make a + /// program dynamic, as steps with their jumps; measurements at the end, + /// even between gates on other qubits, keep it as it was. + #[test] + fn dynamic_programs_compile_to_steps() { + let fin = parse("qubit[2] q; bit[2] c; h q[0]; c[0] = measure q[0]; x q[1]; c[1] = measure q[1];").unwrap(); + assert_eq!(fin.dynamic, None); + assert_eq!(fin.measured, vec![(0, 0), (1, 1)]); + let b = parse("qubit q; bit b = measure q;").unwrap(); + assert_eq!((b.dynamic, b.measured, b.n_bits), (None, vec![(0, 0)], 1)); + + let p = parse("qubit[2] q; bit[2] c; h q[0]; c[0] = measure q[0]; if (c[0]) { x q[1]; } else { z q[1]; } c[1] = measure q[1];").unwrap(); + let one = Cond { bits: vec![0], value: 1, equal: true }; + assert_eq!( + p.dynamic.unwrap(), + vec![ + Instr::Gate(0), + Instr::Measure { qubit: 0, bit: 0 }, + Instr::JumpUnless(one, 5), + Instr::Gate(1), + Instr::Jump(6), + Instr::Gate(2), + Instr::Measure { qubit: 1, bit: 1 }, + ] + ); + assert_eq!(p.ops.iter().map(|o| o.kind).collect::>(), vec![Kind::H, Kind::X, Kind::Z]); + + // OpenQASM 2's `if (creg == n) gate;`, a register condition, and `!`. + let q2 = parse("OPENQASM 2.0; include \"qelib1.inc\"; qreg q[2]; creg c[2]; measure q -> c; if(c==2) x q[0]; measure q -> c;").unwrap(); + assert!(matches!(&q2.dynamic.unwrap()[2], Instr::JumpUnless(Cond { bits, value: 2, equal: true }, 4) if *bits == vec![0, 1])); + let not = parse("qubit[2] q; bit[2] c; c[1] = measure q[1]; if (!c[1]) x q[0]; c[0] = measure q[0];").unwrap(); + assert!(matches!(¬.dynamic.unwrap()[1], Instr::JumpUnless(Cond { bits, value: 1, equal: false }, 3) if *bits == vec![1])); + + // Reset after a gate, and a gate after a measurement. + assert!(parse("qubit[2] q; x q[0]; reset q[0];").unwrap().dynamic.is_some()); + assert!(parse("qubit[1] q; bit c; c = measure q; x q;").unwrap().dynamic.is_some()); + assert!(parse("qubit[1] q; bit[2] c; c[0] = measure q; c[1] = measure q;").unwrap().dynamic.is_some()); + // A reset before a qubit's first gate changes nothing. + assert_eq!(parse("qubit[1] q; reset q; h q;").unwrap().dynamic, None); + } + + #[test] + fn conditions_read_the_bits_lowest_first() { + let c = Cond { bits: vec![3, 1], value: 2, equal: true }; + let mut bits = vec![false; 4]; + assert!(!c.holds(&bits)); + bits[1] = true; + assert!(c.holds(&bits), "bit 1 is the condition's second bit, worth 2"); + bits[3] = true; + assert!(!c.holds(&bits)); + assert!(Cond { equal: false, ..c }.holds(&bits)); + } + #[test] fn dyadic_angles_are_recognised() { use std::f64::consts::PI; @@ -1783,7 +1996,7 @@ mod tests { ops.push(Op { kind, qubits, params }); } } - Program { n_qubits: n, ops, measured: (0..n).map(|q| (q, q)).collect(), n_bits: n } + Program { n_qubits: n, ops, measured: (0..n).map(|q| (q, q)).collect(), n_bits: n, dynamic: None } } fn dense(p: &Program) -> wai_quantum::quantum_sv::StateVector { @@ -1938,7 +2151,7 @@ mod tests { ops.push(Op { kind, qubits, params }); } } - Program { n_qubits: n, ops, measured: (0..n).map(|q| (q, q)).collect(), n_bits: n } + Program { n_qubits: n, ops, measured: (0..n).map(|q| (q, q)).collect(), n_bits: n, dynamic: None } } #[test] diff --git a/wai/wai-quantum-service/src/run.rs b/wai/wai-quantum-service/src/run.rs index 02150748f..f8b46883f 100644 --- a/wai/wai-quantum-service/src/run.rs +++ b/wai/wai-quantum-service/src/run.rs @@ -95,6 +95,11 @@ pub fn run(request: &Value, plan: &Plan) -> Result { "qubits": n, "gates": r.program.ops.len(), }); + if let Some(steps) = &r.program.dynamic { + run_dynamic(&r, plan, steps, &mut out)?; + out["elapsed_ms"] = json!(started.elapsed().as_millis() as u64); + return Ok(out); + } match plan.method.as_str() { "exact" => { let c = r.program.dyadic()?; @@ -433,6 +438,55 @@ pub fn run(request: &Value, plan: &Plan) -> Result { Ok(out) } +/// A program that measures in the middle, resets, or acts on measured bits +/// ([`crate::dynamic`]): counts by shot splitting, expectation values as the +/// exact mixture over its branches. +fn run_dynamic(r: &catalog::Request, plan: &Plan, steps: &[crate::qasm::Instr], out: &mut Value) -> Result<(), String> { + use crate::dynamic::{self, Dense, Stabilizer}; + let (n, nb) = (r.program.n_qubits, r.program.n_bits); + let strings: Vec> = r.observables.iter().map(|(_, _, s)| s.clone()).collect(); + let (sampled, exact) = match plan.method.as_str() { + "statevector" => { + let gates = r.program.matrices(); + let states = (plan.memory_bytes / catalog::statevector_bytes(n)).max(1) as usize; + let start = Dense { sv: wai_quantum::quantum_sv::StateVector::zero(n).map_err(|e| format!("{e:?}"))?, gates: &gates }; + let s = if r.shots > 0 { Some(dynamic::sample(steps, nb, start.clone(), r.shots, r.seed, states)?) } else { None }; + let e = if strings.is_empty() { None } else { Some(dynamic::expectations(steps, &r.program.ops, nb, start, &strings, states)?) }; + (s, e) + } + "stabilizer" => { + let ops: Vec> = r.program.clifford_ops()?.into_iter().map(|(g, _)| g).collect(); + let start = Stabilizer { t: wai_quantum::quantum_tableau::Tableau::new(n).map_err(|e| e.to_string())?, ops: &ops }; + let s = if r.shots > 0 { Some(dynamic::sample(steps, nb, start.clone(), r.shots, r.seed, 64)?) } else { None }; + let e = if strings.is_empty() { None } else { Some(dynamic::expectations(steps, &r.program.ops, nb, start, &strings, 64)?) }; + (s, e) + } + m => return Err(format!("{m} does not run programs that measure in the middle")), + }; + out["accuracy"] = json!({ + "class": "exact", + "note": "counts: the shots split at each mid-circuit measurement and reset as independent shots would, so they have the law of independent runs; expectation values: the exact mixture over every branch of non-zero probability, of the state before the final measurements" + }); + let mut d = json!({ "measurements_and_resets": dynamic::measurements(steps) }); + if let Some(s) = sampled { + d["branches_run"] = json!(s.branches); + d["states_at_once"] = json!(s.peak_states); + merge(out, counts(s.counts)); + } + if let Some((values, leaves)) = exact { + d["branches_weighed"] = json!(leaves); + out["expectations"] = Value::Array( + r.observables + .iter() + .zip(values) + .map(|((text, coeff, _), v)| json!({ "observable": text, "coeff": coeff, "value": coeff * v })) + .collect(), + ); + } + out["dynamic"] = d; + Ok(()) +} + fn merge(out: &mut Value, extra: Value) { if let (Some(o), Some(e)) = (out.as_object_mut(), extra.as_object()) { for (k, v) in e { @@ -733,4 +787,92 @@ mod tests { assert_eq!(a["counts"], b["counts"]); assert_eq!(a["state_hash"].as_str().unwrap().len(), 64); } + + + /// Teleportation: q0's state arrives on q2 through two measurements and + /// the corrections they select, so ⟨Z2⟩ is cos θ exactly and q2 reads 1 + /// with probability sin²(θ/2). + const TELEPORT: &str = "OPENQASM 3.0; qubit[3] q; bit[3] c; ry(THETA) q[0]; h q[1]; cx q[1], q[2]; cx q[0], q[1]; h q[0]; c[0] = measure q[0]; c[1] = measure q[1]; if (c[1]) x q[2]; if (c[0]) z q[2]; c[2] = measure q[2];"; + + #[test] + fn teleportation_delivers_the_state() { + for theta in [0.0, 0.7, 1.9, std::f64::consts::PI] { + let qasm = TELEPORT.replace("THETA", &format!("{theta}")); + let r = go(json!({ "qasm": qasm, "method": "statevector", "observables": ["Z2", "X2"], "shots": 20_000, "seed": 9 })); + assert_eq!(r["method"], "statevector", "{r}"); + let z = r["expectations"][0]["value"].as_f64().unwrap(); + let x = r["expectations"][1]["value"].as_f64().unwrap(); + assert!((z - theta.cos()).abs() < 1e-12 && (x - theta.sin()).abs() < 1e-12, "θ {theta}: ⟨Z2⟩ {z}, ⟨X2⟩ {x}"); + let ones: u64 = r["counts"].as_object().unwrap().iter().filter(|(k, _)| k.starts_with('1')).map(|(_, v)| v.as_u64().unwrap()).sum(); + let p = (theta / 2.0).sin().powi(2); + let sd = (20_000.0 * p * (1.0 - p)).sqrt().max(1.0); + assert!((ones as f64 - 20_000.0 * p).abs() < 5.0 * sd, "θ {theta}: {ones} of 20000 read 1, expected {:.0}", 20_000.0 * p); + assert_eq!(r["counts"].as_object().unwrap().values().map(|v| v.as_u64().unwrap()).sum::(), 20_000); + // c0 and c1 are uniform and independent of what arrives. + assert_eq!(r["dynamic"]["branches_weighed"], 4, "{r}"); + } + } + + #[test] + fn a_clifford_teleport_runs_on_the_tableau_as_on_the_state_vector() { + // |+i⟩ teleported: ⟨Y2⟩ = 1. + let qasm = "qubit[3] q; bit[3] c; h q[0]; s q[0]; h q[1]; cx q[1], q[2]; cx q[0], q[1]; h q[0]; c[0] = measure q[0]; c[1] = measure q[1]; if (c[1] == 1) { x q[2]; } if (c[0] == 1) { z q[2]; } c[2] = measure q[2];"; + let auto = go(json!({ "qasm": qasm, "observables": ["Y2", "Z2"] })); + assert_eq!(auto["method"], "stabilizer", "{auto}"); + let sv = go(json!({ "qasm": qasm, "method": "statevector", "observables": ["Y2", "Z2"] })); + for (a, b) in auto["expectations"].as_array().unwrap().iter().zip(sv["expectations"].as_array().unwrap()) { + assert!((a["value"].as_f64().unwrap() - b["value"].as_f64().unwrap()).abs() < 1e-12, "{a} vs {b}"); + } + assert_eq!(auto["expectations"][0]["value"], 1.0); + // The counts agree in law: c2 reads 0 or 1 evenly (Z on |+i⟩), c0 c1 uniform. + let st = go(json!({ "qasm": qasm, "method": "stabilizer", "shots": 40_000, "seed": 3 })); + assert_eq!(st["distinct_outcomes"], 8, "{st}"); + for (k, v) in st["counts"].as_object().unwrap() { + let v = v.as_u64().unwrap() as f64; + assert!((v - 5_000.0).abs() < 5.0 * (5_000.0f64 * 7.0 / 8.0).sqrt(), "{k}: {v}"); + } + } + + #[test] + fn reset_and_measurement_mid_circuit() { + // A reset returns the qubit to |0⟩ whatever it held. + let r = go(json!({ "qasm": "qubit[2] q; bit[2] c; h q[0]; cx q[0], q[1]; reset q[0]; c = measure q;", "method": "statevector", "shots": 4000, "observables": ["Z0", "Z1"] })); + let counts = r["counts"].as_object().unwrap(); + assert!(counts.keys().all(|k| k.ends_with('0')), "q0 always reads 0 after reset: {r}"); + assert_eq!(r["expectations"][0]["value"], 1.0); + assert!(r["expectations"][1]["value"].as_f64().unwrap().abs() < 1e-12, "{r}"); + // Measured, flipped, measured again: the second reading is the first's opposite. + for method in ["statevector", "stabilizer"] { + let r = go(json!({ "qasm": "qubit[1] q; bit[2] c; h q[0]; c[0] = measure q[0]; x q[0]; c[1] = measure q[0];", "method": method, "shots": 2000, "seed": 4 })); + let counts = r["counts"].as_object().unwrap(); + assert!(counts.keys().all(|k| k == "01" || k == "10") && counts.len() == 2, "{method}: {r}"); + } + } + + /// The three-qubit bit-flip code: an injected flip is found by two + /// ancilla measurements and undone by the correction they select. + #[test] + fn a_syndrome_measurement_corrects_the_flip() { + for flipped in 0..3 { + let qasm = format!( + "qubit[5] q; bit[2] s; bit r; ry(0.7) q[0]; cx q[0], q[1]; cx q[0], q[2]; x q[{flipped}]; cx q[0], q[3]; cx q[1], q[3]; cx q[1], q[4]; cx q[2], q[4]; s[0] = measure q[3]; s[1] = measure q[4]; if (s == 1) x q[0]; if (s == 3) x q[1]; if (s == 2) x q[2]; cx q[0], q[1]; cx q[0], q[2]; r = measure q[0];" + ); + let r = go(json!({ "qasm": qasm, "method": "statevector", "observables": ["Z0"], "shots": 1000 })); + assert!((r["expectations"][0]["value"].as_f64().unwrap() - 0.7f64.cos()).abs() < 1e-12, "flip on {flipped}: {r}"); + // The syndrome is determined: one branch, and its bits name the flip. + assert_eq!(r["dynamic"]["branches_weighed"], 1, "{r}"); + let want = ["01", "11", "10"][flipped]; + assert!(r["counts"].as_object().unwrap().keys().all(|k| k.ends_with(want)), "flip on {flipped}: {r}"); + } + } + + #[test] + fn dynamic_runs_repeat_with_their_seed() { + let qasm = TELEPORT.replace("THETA", "1.1"); + let a = go(json!({ "qasm": qasm, "method": "statevector", "shots": 5000, "seed": 21 })); + let b = go(json!({ "qasm": qasm, "method": "statevector", "shots": 5000, "seed": 21 })); + let c = go(json!({ "qasm": qasm, "method": "statevector", "shots": 5000, "seed": 22 })); + assert_eq!(a["counts"], b["counts"]); + assert_ne!(a["counts"], c["counts"]); + } } diff --git a/wai/wai-quantum/Cargo.lock b/wai/wai-quantum/Cargo.lock index 662dd7c4b..3294d0fa5 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.4" +version = "0.6.5" dependencies = [ "blake3", "block2", diff --git a/wai/wai-quantum/Cargo.toml b/wai/wai-quantum/Cargo.toml index 425891c35..66f7547ff 100644 --- a/wai/wai-quantum/Cargo.toml +++ b/wai/wai-quantum/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "wai-quantum" -version = "0.6.4" +version = "0.6.5" edition = "2024" rust-version = "1.97" license = "Apache-2.0" diff --git a/wai/wai-quantum/src/quantum_sv.rs b/wai/wai-quantum/src/quantum_sv.rs index 1e2b0741c..9693c801a 100644 --- a/wai/wai-quantum/src/quantum_sv.rs +++ b/wai/wai-quantum/src/quantum_sv.rs @@ -858,6 +858,48 @@ impl StateVector { self.chunked_sum(&|x| self.re[x] * self.re[x] + self.im[x] * self.im[x]) } + /// The probability that a measurement of qubit `q` reads 1: `Σ |amplitude|²` + /// over the indices whose bit for `q` is set, in fixed chunks, so the same + /// at any thread count. Assumes a normalised state; 0 for a qubit outside + /// the register. + pub fn probability_one(&self, q: u32) -> f64 { + if q >= self.n { + return 0.0; + } + let bit = 1usize << (self.n - 1 - q); + self.chunked_sum(&|x| if x & bit != 0 { self.re[x] * self.re[x] + self.im[x] * self.im[x] } else { 0.0 }) + } + + /// The state a measurement of qubit `q` leaves when it reads `outcome`: + /// the amplitudes that disagree are zeroed and the rest renormalised. + /// Returns the outcome's probability. An outcome of probability 0 cannot + /// be read, and leaves the state as it was, returning 0. + /// [`SvError::BadGate`] for a qubit outside the register. + pub fn collapse(&mut self, q: u32, outcome: bool) -> Result { + if q >= self.n { + return Err(SvError::BadGate(0)); + } + let bit = 1usize << (self.n - 1 - q); + let p = self.chunked_sum(&|x| if (x & bit != 0) == outcome { self.re[x] * self.re[x] + self.im[x] * self.im[x] } else { 0.0 }); + if p <= 0.0 { + return Ok(0.0); + } + let s = 1.0 / p.sqrt(); + let f = |re: &mut [f64], im: &mut [f64], off: usize| { + for (j, (a, b)) in re.iter_mut().zip(im.iter_mut()).enumerate() { + if ((off + j) & bit != 0) == outcome { + *a *= s; + *b *= s; + } else { + *a = 0.0; + *b = 0.0; + } + } + }; + par_slices(&mut self.re, &mut self.im, 1, &f); + Ok(p) + } + /// `Σ_x term(x)`, summed in [`CHUNKS`] fixed chunks (each in order, then /// the chunks in order), threaded where available. fn chunked_sum(&self, term: &(dyn Fn(usize) -> f64 + Sync)) -> f64 { diff --git a/wai/wai-quantum/src/quantum_tableau.rs b/wai/wai-quantum/src/quantum_tableau.rs index d6bac6f5b..66e9a7bab 100644 --- a/wai/wai-quantum/src/quantum_tableau.rs +++ b/wai/wai-quantum/src/quantum_tableau.rs @@ -118,6 +118,23 @@ fn set(v: &mut [u64], i: usize) { v[i / 64] |= 1 << (i % 64); } +/// The first set bit of `v` in `from..to`. +fn first_set(v: &[u64], from: usize, to: usize) -> Option { + let mut w = from / 64; + let mut word = v[w] & (u64::MAX << (from % 64)); + loop { + if word != 0 { + let i = w * 64 + word.trailing_zeros() as usize; + return (i < to).then_some(i); + } + w += 1; + if w * 64 >= to || w >= v.len() { + return None; + } + word = v[w]; + } +} + fn splitmix(s: &mut u64) -> u64 { *s = s.wrapping_add(0x9E37_79B9_7F4A_7C15); let mut z = *s; @@ -338,6 +355,99 @@ impl Tableau { Ok(()) } + /// 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. + pub fn is_random(&self, a: u32) -> Result { + let a = self.check(a)?; + Ok(first_set(&self.xs[a * self.rw..(a + 1) * self.rw], self.n, 2 * self.n).is_some()) + } + + /// Measures qubit `a` in the computational basis and leaves the state the + /// outcome selects: `(outcome, random)`, the outcome `true` for 1. + /// + /// A determined outcome (the state is an eigenstate of `Z_a`) is read + /// from the tableau and changes nothing. A random one has probability ½ + /// each way; it is `random()`'s, so a caller decides it (a seeded draw, or + /// a branch it wants), and the state becomes the matching eigenstate. This + /// is the Aaronson–Gottesman measurement: the first stabilizer that + /// anticommutes with `Z_a` is multiplied into every other row that does, + /// all rows at once in this layout (each row's phase summed mod 4 in two + /// bit planes), moved to its destabilizer, and replaced by `±Z_a`. + pub fn measure(&mut self, a: u32, random: impl FnOnce() -> bool) -> Result<(bool, bool), TableauError> { + let a = self.check(a)?; + let (n, rw) = (self.n, self.rw); + let col: Vec = self.xs[a * rw..(a + 1) * rw].to_vec(); + let Some(p) = first_set(&col, n, 2 * n) else { + // Z_a, up to sign, is a product of the stabilizers. + return Ok((self.expectation(&[(a as u32, 'Z')])? < 0, false)); + }; + let outcome = random(); + // Every other row whose X part holds qubit a: row ← row p · row. + let mut m = col; + m[p / 64] &= !(1u64 << (p % 64)); + let (mut lo, mut hi) = (vec![0u64; rw], vec![0u64; rw]); + for q in 0..n { + let (xp, zp) = (get(&self.xs[q * rw..(q + 1) * rw], p), get(&self.zs[q * rw..(q + 1) * rw], p)); + if !xp && !zp { + continue; + } + for w in 0..rw { + let mw = m[w]; + if mw == 0 { + continue; + } + let (x2, z2) = (self.xs[q * rw + w] & mw, self.zs[q * rw + w] & mw); + // The power of i row p's letter and this row's letter carry + // (`phase_word` with row p's letter fixed). + let (plus, minus) = match (xp, zp) { + (true, true) => (z2 & !x2, x2 & !z2), + (true, false) => (z2 & x2, z2 & !x2), + _ => (x2 & !z2, x2 & z2), + }; + // +1 where `plus`, +3 where `minus`, mod 4. + let c = lo[w] & plus; + lo[w] ^= plus; + hi[w] ^= c; + let c = lo[w] & minus; + lo[w] ^= minus; + hi[w] ^= c ^ minus; + if xp { + self.xs[q * rw + w] ^= mw; + } + if zp { + self.zs[q * rw + w] ^= mw; + } + } + } + // Each row's sign: 2·r + 2·r_p + Σ ≡ 2 (mod 4) is −. The rows that + // commute with row p sum evenly; destabilizer p − n, which does not, + // is overwritten next. + let d = p - n; + let rp = if get(&self.rs, p) { u64::MAX } else { 0 }; + for w in 0..rw { + let mw = m[w]; + let skip = if d / 64 == w { 1u64 << (d % 64) } else { 0 }; + debug_assert_eq!(lo[w] & mw & !skip, 0, "a product of commuting rows has a real phase"); + let negative = !lo[w] & (hi[w] ^ self.rs[w] ^ rp); + self.rs[w] = (self.rs[w] & !mw) | (negative & mw); + } + // Destabilizer p − n ← the old row p; row p ← ±Z_a. + for q in 0..n { + for bits in [&mut self.xs, &mut self.zs] { + let col = &mut bits[q * rw..(q + 1) * rw]; + let was = get(col, p); + col[d / 64] = (col[d / 64] & !(1u64 << (d % 64))) | (u64::from(was) << (d % 64)); + col[p / 64] &= !(1u64 << (p % 64)); + } + } + let was = get(&self.rs, p); + self.rs[d / 64] = (self.rs[d / 64] & !(1u64 << (d % 64))) | (u64::from(was) << (d % 64)); + set(&mut self.zs[a * rw..(a + 1) * rw], p); + self.rs[p / 64] = (self.rs[p / 64] & !(1u64 << (p % 64))) | (u64::from(outcome) << (p % 64)); + Ok((outcome, true)) + } + /// Rows `from..to` of the tableau, held row by row. fn rows(&self, from: usize, to: usize) -> Rows { let (n, rw) = (self.n, self.rw); @@ -1134,6 +1244,79 @@ mod tests { } } } + + /// Measurements between gates: each outcome has the state vector's + /// probability (½ when the tableau calls it random, 0 or 1 when it + /// reads it), and leaves the state the state vector's projection + /// leaves, Pauli string by Pauli string, round after round. + #[test] + fn mid_circuit_measurements_collapse_as_the_state_vector_projects() { + let (mut random, mut determined) = (0, 0); + for seed in 0..60u64 { + let n = 2 + (seed % 8) as u32; + let mut r = seed ^ 0x5eed; + let mut t = Tableau::new(n).unwrap(); + let mut sv = StateVector::zero(n).unwrap(); + for round in 0..6 { + for o in random_clifford(n, 10, seed * 31 + round) { + t.apply(o).unwrap(); + sv.apply(&matrix(o)).unwrap(); + } + let a = (next(&mut r) % n as u64) as u32; + let p1 = sv.probability_one(a); + let want = next(&mut r) & 1 == 1; + let said_random = t.is_random(a).unwrap(); + let (got, was_random) = t.measure(a, || want).unwrap(); + assert_eq!(said_random, was_random, "seed {seed}: is_random agrees with the measurement"); + if was_random { + random += 1; + assert!((p1 - 0.5).abs() < 1e-9, "seed {seed} round {round}: random, yet P(1) = {p1}"); + assert_eq!(got, want, "seed {seed}: a random outcome is the caller's"); + } else { + determined += 1; + assert!((p1 - f64::from(u8::from(got))).abs() < 1e-9, "seed {seed} round {round}: read {got}, P(1) = {p1}"); + } + assert!(sv.collapse(a, got).unwrap() > 0.4, "seed {seed}: the outcome was possible"); + assert_eq!(t.expectation(&[(a, 'Z')]).unwrap(), if got { -1 } else { 1 }); + for _ in 0..24 { + let string: Vec<(u32, char)> = (0..n).filter_map(|q| match next(&mut r) % 4 { 0 => None, 1 => Some((q, 'X')), 2 => Some((q, 'Y')), _ => Some((q, 'Z')) }).collect(); + let want = sv.expectation(&[(1.0, string.clone())]); + let got = t.expectation(&string).unwrap() as f64; + assert!((got - want).abs() < 1e-9, "seed {seed} round {round}: {string:?}: {got} vs {want}"); + } + } + } + assert!(random > 50 && determined > 50, "both kinds of outcome were exercised: {random} random, {determined} determined"); + } + } + + /// Wide tableaux (several words a column): measuring every qubit in turn, + /// with outcomes chosen where random, gives a bit string inside the + /// support the tableau held before, and measuring again repeats it. + #[test] + fn measuring_every_qubit_lands_inside_the_support() { + for seed in 0..12u64 { + let n = 40 + (seed * 11 % 110) as u32; + let mut t = Tableau::run(n, &random_clifford(n, 6 * n as usize, seed)).unwrap(); + let support = t.support(); + let mut r = seed ^ 0xfeed; + let mut bits = vec![0u64; (n as usize).div_ceil(64)]; + let mut random = 0; + for q in 0..n { + let pick = next(&mut r) & 1 == 1; + let (b, was_random) = t.measure(q, || pick).unwrap(); + random += usize::from(was_random); + if b { + bits[q as usize / 64] |= 1 << (q % 64); + } + } + assert_eq!(random, support.dimension(), "seed {seed}: one random outcome per dimension of the support"); + assert!(support.contains(&bits), "seed {seed}, {n} qubits: the outcome is outside the support"); + for q in 0..n { + let (b, was_random) = t.measure(q, || unreachable!("a measured qubit is determined")).unwrap(); + assert!(!was_random && b == (bits[q as usize / 64] >> (q % 64) & 1 == 1), "seed {seed}: qubit {q} read again differs"); + } + } } } diff --git a/wai/wai-rs/Cargo.lock b/wai/wai-rs/Cargo.lock index 47f996cd2..d27c01998 100644 --- a/wai/wai-rs/Cargo.lock +++ b/wai/wai-rs/Cargo.lock @@ -1295,7 +1295,7 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "39cab71617ae0d63f51a36d69f866391735b51691dbda63cf6f96d042b63efeb" dependencies = [ "libc", - "windows-sys 0.61.2", + "windows-sys 0.52.0", ] [[package]] @@ -3354,7 +3354,7 @@ dependencies = [ "errno", "libc", "linux-raw-sys", - "windows-sys 0.61.2", + "windows-sys 0.52.0", ] [[package]] @@ -3411,7 +3411,7 @@ dependencies = [ "security-framework", "security-framework-sys", "webpki-root-certs", - "windows-sys 0.61.2", + "windows-sys 0.52.0", ] [[package]] @@ -3934,10 +3934,10 @@ source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "32497e9a4c7b38532efcdebeef879707aa9f794296a4f0244f6f69e9bc8574bd" dependencies = [ "fastrand", - "getrandom 0.4.3", + "getrandom 0.3.4", "once_cell", "rustix", - "windows-sys 0.61.2", + "windows-sys 0.52.0", ] [[package]] @@ -4310,7 +4310,7 @@ version = "0.1.0" [[package]] name = "wai-quantum" -version = "0.6.4" +version = "0.6.5" dependencies = [ "blake3", "ed25519-dalek", @@ -4611,7 +4611,7 @@ version = "0.1.11" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "c2a7b1c03c876122aa43f3020e6c3c3ee5c05081c9a00739faf7503aeba10d22" dependencies = [ - "windows-sys 0.61.2", + "windows-sys 0.52.0", ] [[package]] diff --git a/wai/wai-web/Cargo.lock b/wai/wai-web/Cargo.lock index 2f02cf70e..a18a919fd 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.4" +version = "0.6.5" dependencies = [ "blake3", "ed25519-dalek", -- 2.54.0 (Apple Git-157) From d73788c5b0cdf25eacd153e939ec286ab4620317 Mon Sep 17 00:00:00 2001 From: Transaction Science Date: Mon, 19 Oct 2026 06:34:42 -0400 Subject: [PATCH 2/3] wai-quantum-service: the exact walk's eight parameters, as clippy allows the library's walker Co-Authored-By: Claude Opus 5.5 --- wai/wai-quantum-service/src/dynamic.rs | 1 + 1 file changed, 1 insertion(+) diff --git a/wai/wai-quantum-service/src/dynamic.rs b/wai/wai-quantum-service/src/dynamic.rs index 659e275a0..2de9314e5 100644 --- a/wai/wai-quantum-service/src/dynamic.rs +++ b/wai/wai-quantum-service/src/dynamic.rs @@ -214,6 +214,7 @@ impl Walk<'_> { /// The exact mixture's expectation values from step `pc`, the branch /// weighted by `weight`, added into `sums`; `leaves` counts the branches. + #[allow(clippy::too_many_arguments)] fn expect(&mut self, mut pc: usize, mut state: E, mut bits: Vec, weight: f64, observables: &[Vec<(u32, char)>], sums: &mut [f64], leaves: &mut u64) -> Result<(), String> { while pc < self.steps.len() { match &self.steps[pc] { -- 2.54.0 (Apple Git-157) From 33397f925ad2b1399360d7c79ea90f4026e0b34d Mon Sep 17 00:00:00 2001 From: Transaction Science Date: Mon, 19 Oct 2026 08:21:31 -0400 Subject: [PATCH 3/3] wai-quantum-service smoke: `if` is no longer the refusal example; teleportation checks the dynamic path The refusal check sent `if (c) x q;`, which now parses (and is refused for an undeclared register instead). It sends a `while` loop, still refused. A new check runs teleportation through the service: = cos 0.7 and = sin 0.7 exactly, four branches weighed, the shots all counted. All checks held against the branch's release build on the build host. Co-Authored-By: Claude Opus 5.5 --- wai/wai-quantum-service/smoke.sh | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/wai/wai-quantum-service/smoke.sh b/wai/wai-quantum-service/smoke.sh index 20f111d6e..39611dc5d 100644 --- a/wai/wai-quantum-service/smoke.sh +++ b/wai/wai-quantum-service/smoke.sh @@ -226,8 +226,18 @@ hits=r["counts"].get("10110", 0) assert d["status"]=="done" and hits > 980, d print("grover-5 (ctrl modifiers): ok,", r["method"], hits, "of 1000 on the marked item")' "$tmp/grover.out" +# A program that measures in the middle and acts on what it read: teleportation +# delivers ry(0.7)|0⟩ to q[2], so ⟨Z2⟩ = cos 0.7 and ⟨X2⟩ = sin 0.7 exactly. +py 'import json; print(json.dumps({"qasm": "OPENQASM 3.0; qubit[3] q; bit[3] c; ry(0.7) q[0]; h q[1]; cx q[1], q[2]; cx q[0], q[1]; h q[0]; c[0] = measure q[0]; c[1] = measure q[1]; if (c[1]) x q[2]; if (c[0]) z q[2]; c[2] = measure q[2];", "observables": ["Z2", "X2"], "shots": 2000}))' > "$tmp/tele.json" +post /v1/run "$tmp/tele.json" > "$tmp/tele.out" +py 'import json,math,sys +d=json.load(open(sys.argv[1])); r=d["result"]; e=[x["value"] for x in r["expectations"]] +assert d["status"]=="done" and abs(e[0]-math.cos(0.7))<1e-12 and abs(e[1]-math.sin(0.7))<1e-12, d +assert r["dynamic"]["branches_weighed"]==4 and sum(r["counts"].values())==2000, r +print("teleport (mid-circuit measurement, if): ok,", r["method"], e)' "$tmp/tele.out" + # A refusal says why. -py 'import json; print(json.dumps({"qasm": "qubit[2] q; if (c) x q;"}))' > "$tmp/bad.json" +py 'import json; print(json.dumps({"qasm": "qubit[1] q; bit c; c = measure q; while (c) { x q; c = measure q; }"}))' > "$tmp/bad.json" code=$(curl -sS -o "$tmp/bad.out" -w '%{http_code}' -X POST "$U/v1/jobs" -H 'content-type: application/json' --data-binary "@$tmp/bad.json") [ "$code" = 400 ] || { echo "a bad program got HTTP $code"; exit 1; } py 'import json,sys; d=json.load(open(sys.argv[1])); assert "not supported" in d["error"], d; print("refusal: ok,", d["error"])' "$tmp/bad.out" -- 2.54.0 (Apple Git-157)