From 617075d3877e3abdd2831c383619494e0894a3f4 Mon Sep 17 00:00:00 2001 From: Transaction Science Date: Thu, 1 Oct 2026 15:09:35 -0400 Subject: [PATCH 1/2] wai-quantum 0.3.22: sparse Pauli dynamics at utility scale, and identical bits on every platform MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit quantum_spd computes observables of circuits up to 1024 qubits with arbitrary rotation angles, in the Heisenberg picture. Until now every backend took the dyadic circuit type, which allows only Clifford, T and 2π/2^k phases, and sparse-Pauli propagation stopped at 64 qubits. The experiments that classical methods are known for reproducing were out of reach. - Clifford gates and quarter-turn rotations are exact. - Continuous rotations merge in place, over an insertion-ordered term table, so every coefficient is added in the same order on every platform. - Every dropped coefficient goes into a truncation bound, a rigorous limit on how far truncation can have moved the result. - On the 127-qubit heavy-hex kicked-Ising circuits, it reproduces the exactly computed 5-step magnetization (arXiv:2308.05077 exact.csv) on all 127 qubits: 0.995131 against 0.995131 at θ_h = π/32, exact propagation. - 20-step ⟨Z62⟩ at 2π/32 falls inside the published converged classical values. quantum_spd_receipt seals a run under a meter. The result is computed inside the seal, floats are hashed and serialized as exact bits, and verify_computation re-runs the job anywhere and demands the same bits. The energy is attested, and can carry an energy-class label. Defects fixed: - quantum_pauli was not reproducible from run to run. It summed coefficients in std HashMap order, which is seeded afresh for every map, so six runs gave six bit patterns. It now visits terms in key order. - Four f64 layers (vml, phasor, qfhrr, kernel) gave different bits natively and as a WASI component, through the platform's cos, sin, ln and atan2. That falsified the crate's identical-native-browser-edge claim. src/repro.rs computes them in-crate (fdlibm's algorithms, + − × ÷ only, within an ulp). The MPS backend, whose contract excludes transcendentals, used them too. tests/bit_identity.rs pins a digest per layer, confirmed natively and under wasm32-wasip2. wai-rs: - passthrough features quantum_spd and quantum_spd_receipt; - `wai-quantum utility [--theta --steps --threshold -q --json --receipt]`; - `verify` reads both receipt kinds. The utility receipt is byte-identical natively and under wasm32-wasip2. Dependent lockfiles were moved. Every feature now builds on its own with no warnings. The merkle module is compiled only for the receipt features that use it, and the QEC encoder's row reduction only with the simulator it serves. Co-Authored-By: Claude Opus 5.5 --- wai/wai-quantum-serve/Cargo.lock | 2 +- wai/wai-quantum/Cargo.lock | 2 +- wai/wai-quantum/Cargo.toml | 20 +- wai/wai-quantum/README.md | 30 +- wai/wai-quantum/examples/utility.rs | 31 + wai/wai-quantum/examples/utility_mz.rs | 37 + wai/wai-quantum/src/lib.rs | 18 + wai/wai-quantum/src/quantum_kernel.rs | 9 +- wai/wai-quantum/src/quantum_mps.rs | 9 +- wai/wai-quantum/src/quantum_pauli.rs | 60 +- wai/wai-quantum/src/quantum_phasor.rs | 5 +- wai/wai-quantum/src/quantum_qec.rs | 1 + wai/wai-quantum/src/quantum_qfhrr.rs | 9 +- wai/wai-quantum/src/quantum_spd.rs | 1268 ++++++++++++++++++++ wai/wai-quantum/src/quantum_spd_receipt.rs | 427 +++++++ wai/wai-quantum/src/quantum_vml.rs | 6 +- wai/wai-quantum/src/repro.rs | 387 ++++++ wai/wai-quantum/tests/bit_identity.rs | 109 ++ wai/wai-quantum/tests/utility_reference.rs | 37 + wai/wai-rs/Cargo.lock | 2 +- wai/wai-rs/Cargo.toml | 5 + wai/wai-rs/src/bin/wai_quantum.rs | 72 ++ wai/wai-rs/src/lib.rs | 6 + wai/wai-web/Cargo.lock | 2 +- 24 files changed, 2516 insertions(+), 38 deletions(-) create mode 100644 wai/wai-quantum/examples/utility.rs create mode 100644 wai/wai-quantum/examples/utility_mz.rs create mode 100644 wai/wai-quantum/src/quantum_spd.rs create mode 100644 wai/wai-quantum/src/quantum_spd_receipt.rs create mode 100644 wai/wai-quantum/src/repro.rs create mode 100644 wai/wai-quantum/tests/bit_identity.rs create mode 100644 wai/wai-quantum/tests/utility_reference.rs diff --git a/wai/wai-quantum-serve/Cargo.lock b/wai/wai-quantum-serve/Cargo.lock index 8bc1dd68d..5d62fd3c3 100644 --- a/wai/wai-quantum-serve/Cargo.lock +++ b/wai/wai-quantum-serve/Cargo.lock @@ -589,7 +589,7 @@ dependencies = [ [[package]] name = "wai-quantum" -version = "0.3.21" +version = "0.3.22" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-quantum/Cargo.lock b/wai/wai-quantum/Cargo.lock index a4e98601e..dfcb512b1 100644 --- a/wai/wai-quantum/Cargo.lock +++ b/wai/wai-quantum/Cargo.lock @@ -415,7 +415,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a" [[package]] name = "wai-quantum" -version = "0.3.21" +version = "0.3.22" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-quantum/Cargo.toml b/wai/wai-quantum/Cargo.toml index 548109693..783681c9a 100644 --- a/wai/wai-quantum/Cargo.toml +++ b/wai/wai-quantum/Cargo.toml @@ -1,10 +1,10 @@ [package] name = "wai-quantum" -version = "0.3.21" +version = "0.3.22" edition = "2024" rust-version = "1.97" license = "Apache-2.0" -description = "A deterministic quantum stack in pure Rust: byte-exact circuit simulation (statevector / stabilizer / tensor-network MPS / sparse-Pauli backends), error mitigation, qLDPC decoding, noise learning, circuit-equivalence proofs, a phasor interference-ML layer, information-theoretic limits, noisy channels and state tomography, and signed energy-accounted receipts. No QPU, no cloud, no system libraries — identical results native, in the browser, and as a WASI component at the edge." +description = "A deterministic quantum stack in pure Rust: byte-exact circuit simulation (statevector / stabilizer / tensor-network MPS / sparse-Pauli backends), sparse Pauli dynamics at utility scale (arbitrary angles, 1024 qubits), error mitigation, qLDPC decoding, noise learning, circuit-equivalence proofs, a phasor interference-ML layer, information-theoretic limits, noisy channels and state tomography, and signed energy-accounted receipts. No QPU, no cloud, no system libraries — identical results native, in the browser, and as a WASI component at the edge." repository = "https://git.transaction.science/Transaction-Science/open-standards" homepage = "https://wai.transaction.science/quantum-ide" documentation = "https://docs.rs/wai-quantum" @@ -39,6 +39,12 @@ quantum_toolchain = ["quantum"] quantum_mps = ["quantum"] quantum_stabilizer = ["quantum"] quantum_pauli = ["quantum"] +# Sparse Pauli dynamics at utility scale: arbitrary-angle Pauli rotations on up +# to 1024 qubits, reproducible bit for bit, with a rigorous truncation bound. +# Independent of the dyadic core; `from_dyadic` appears when `quantum` is on. +quantum_spd = [] +# Signed, energy-accounted receipts for it, re-checkable bit for bit anywhere. +quantum_spd_receipt = ["quantum_spd", "_signing"] # Mixed states and the observables that distinguish quantum theory from a classical # account: density matrices, the partial trace, exact Pauli expectations, CHSH. quantum_info = ["quantum"] @@ -78,13 +84,21 @@ quantum_phasor_meter = ["_signing"] # Everything. full = [ - "quantum", "quantum_toolchain", "quantum_mps", "quantum_stabilizer", "quantum_pauli", "quantum_qir", "quantum_info", "quantum_qudit", "quantum_comm", "quantum_source", "quantum_channel", + "quantum", "quantum_toolchain", "quantum_mps", "quantum_stabilizer", "quantum_pauli", "quantum_spd", "quantum_spd_receipt", "quantum_qir", "quantum_info", "quantum_qudit", "quantum_comm", "quantum_source", "quantum_channel", "quantum_receipt", "quantum_ops", "quantum_cal", "quantum_mitigate", "quantum_qec", "quantum_control", "quantum_noise", "quantum_compile", "quantum_atom", "quantum_qbom", "quantum_vml", "quantum_phasor", "quantum_qfhrr", "quantum_kernel", "quantum_qdata", "quantum_phasor_meter", ] +[[example]] +name = "utility" +required-features = ["quantum_spd"] + +[[example]] +name = "utility_mz" +required-features = ["quantum_spd"] + [[example]] name = "bell" required-features = ["quantum", "quantum_toolchain"] diff --git a/wai/wai-quantum/README.md b/wai/wai-quantum/README.md index 32aed7eac..282470b08 100644 --- a/wai/wai-quantum/README.md +++ b/wai/wai-quantum/README.md @@ -58,7 +58,9 @@ never creates information. | `quantum_toolchain` | algorithm library, backend recommendation, **OpenQASM 2/3 ingest of what the mainstream toolchains actually export** — both register and both measurement spellings, `rx`/`ry`/`rz`/`sx`/`u`/`u2`/`u3`/`crz`, and exact rational-multiple-of-pi angles (non-dyadic angles are refused, never rounded) — Bloch / entanglement / purity analysis | | `quantum_stabilizer` | stabilizer (CHP) tableau — `O(n²)`, scales far past statevector | | `quantum_mps` | tensor-network (matrix-product-state) backend | -| `quantum_pauli` | sparse-Pauli / Heisenberg observable propagation | +| `quantum_pauli` | sparse-Pauli / Heisenberg observable propagation over the dyadic gate set (up to 64 qubits) | +| `quantum_spd` | **sparse Pauli dynamics at utility scale**: arbitrary-angle Pauli rotations on up to 1024 qubits, exact Clifford and quarter-turn rotations, a rigorous bound on what truncation can move a result, and the 127-qubit heavy-hex kicked-Ising circuits — whose exactly computed 5-step magnetization it reproduces on all 127 qubits | +| `quantum_spd_receipt` | a signed, energy-accounted receipt over a sparse-Pauli-dynamics run, re-checkable bit for bit on any machine | | `quantum_receipt` | signed, energy-accounted receipt over a reconstruction | | `quantum_ops` | the operations/attestation layer | | `quantum_cal`, `quantum_control`, `quantum_noise` | calibration, filter-function robust control, DD noise spectroscopy, Cycle-Benchmarking noise learning | @@ -93,13 +95,33 @@ the imaginary cross terms are multiplications by exactly zero. At 20 qubits that is 8.9x on `Z`, 7.1x on `X`, 4.6x on `CZ`, 4.3x on `CX`, 3.4x on `T`, 1.9x on `Y` and 1.8x on `H`. +## Utility scale + +`quantum_spd` evolves an observable instead of a state, so a circuit costs what +its non-Clifford rotations make it cost rather than `2ⁿ`: + +```sh +cargo run --release --example utility --features quantum_spd -- 0.2 20 1e-4 +``` + +computes ⟨Z₆₂⟩ after 20 Trotter steps of the 127-qubit heavy-hex kicked-Ising +circuit in a fraction of a second, with the bound on how far truncation can have moved +it. `--example utility_mz` checks the 5-step magnetization against the exactly +computed values for every angle, and `tests/utility_reference.rs` holds the +published references. + ## Determinism Simulation is **byte-exact**: amplitudes are dyadic fixed-point, so a circuit reconstructs to the same `statevector_hash` on every machine and every target. -The approximate methods (MPS truncation, the phasor/ML layer) are -*reproducible f64* — IEEE-754-strict, identical given the same inputs and seed — -and are documented as such rather than claimed byte-exact. +The approximate methods (MPS truncation, sparse Pauli dynamics, the phasor/ML +layer) are *reproducible f64*: IEEE-754-strict, and identical given the same +inputs and seed on every platform. They are documented as such rather than +claimed byte-exact. Their `sin`, `cos`, `ln` and `atan2` are computed in-crate +from `+ − × ÷`, because platform math libraries differ in the last place — enough +that four of these layers once gave different bits natively and as a WASI +component. `tests/bit_identity.rs` pins one result from each, so any platform +whose arithmetic differs fails there. A receipt separates the two halves of a cost honestly: **work** is portable-exact (`n_ops · 2^n` amplitude updates, recomputed by the verifier, so a diff --git a/wai/wai-quantum/examples/utility.rs b/wai/wai-quantum/examples/utility.rs new file mode 100644 index 000000000..d2a0db053 --- /dev/null +++ b/wai/wai-quantum/examples/utility.rs @@ -0,0 +1,31 @@ +//! The 127-qubit heavy-hex kicked-Ising experiment, on a laptop. +//! +//! `cargo run --release --example utility --features quantum_spd -- [θ_h] [steps] [threshold]` +//! +//! Computes ⟨Z₆₂⟩ after `steps` Trotter steps of `RX(θ_h)` on every qubit and +//! `RZZ(−π/2)` on every edge of the heavy-hex map. Prints the value, the bound +//! on what truncation can have moved it, and the work done. + +use std::time::Instant; +use wai_quantum::quantum_spd::{SpdConfig, expect_z, heavy_hex_127, kicked_ising}; + +fn main() { + let args: Vec = std::env::args().skip(1).map(|a| a.parse().expect("a number")).collect(); + let theta_h = args.first().copied().unwrap_or(std::f64::consts::FRAC_PI_4); + let steps = args.get(1).copied().unwrap_or(20.0) as u32; + let threshold = args.get(2).copied().unwrap_or(1e-3); + + let circuit = kicked_ising(127, &heavy_hex_127(), theta_h, steps); + let started = Instant::now(); + let r = expect_z(&circuit, 62, &SpdConfig { threshold, max_weight: None }).expect("a valid circuit"); + let secs = started.elapsed().as_secs_f64(); + + println!("theta_h {theta_h}"); + println!("steps {steps}"); + println!("threshold {threshold:e}"); + println!(" {:.6}", r.value); + println!("bound {:.3e} (|exact - value| is at most this)", r.truncation_bound); + println!("peak terms {}", r.peak_terms); + println!("updates {}", r.term_updates); + println!("seconds {secs:.2}"); +} diff --git a/wai/wai-quantum/examples/utility_mz.rs b/wai/wai-quantum/examples/utility_mz.rs new file mode 100644 index 000000000..cd1768671 --- /dev/null +++ b/wai/wai-quantum/examples/utility_mz.rs @@ -0,0 +1,37 @@ +//! The 5-step magnetization of the heavy-hex kicked-Ising experiment, against +//! the exactly computed values published for it (θ_h = kπ/32). +//! +//! `cargo run --release --example utility_mz --features quantum_spd -- [threshold]` + +use std::time::Instant; +use wai_quantum::quantum_spd::{Pauli, SpdConfig, expectation, heavy_hex_127, kicked_ising}; + +/// Exact Mz after 5 steps at θ_h = kπ/32, k = 1…16. +const EXACT: [f64; 16] = [ + 0.995131, 0.979836, 0.951490, 0.904843, 0.833134, 0.731464, 0.601104, 0.452804, 0.305680, 0.180317, 0.090301, + 0.037263, 0.012260, 0.002936, 0.000324, 0.0, +]; + +fn main() { + let threshold: f64 = std::env::args().nth(1).map_or(0.0, |a| a.parse().expect("a number")); + let kmax: u32 = std::env::args().nth(2).map_or(16, |a| a.parse().expect("a number")); + let edges = heavy_hex_127(); + let mz: Vec<(f64, Vec<(u32, Pauli)>)> = (0..127).map(|q| (1.0 / 127.0, vec![(q, Pauli::Z)])).collect(); + let mut worst: f64 = 0.0; + for (k, exact) in (1..=kmax).zip(EXACT) { + let theta_h = f64::from(k) * std::f64::consts::PI / 32.0; + let started = Instant::now(); + let r = expectation(&kicked_ising(127, &edges, theta_h, 5), &mz, &[], &SpdConfig { threshold, max_weight: None }) + .expect("a valid circuit"); + let err = (r.value - exact).abs(); + worst = worst.max(err); + println!( + "k={k:2} Mz={:.6} exact={exact:.6} |err|={err:.1e} bound={:.1e} peak={} {:.2}s", + r.value, + r.truncation_bound, + r.peak_terms, + started.elapsed().as_secs_f64() + ); + } + println!("worst |err| = {worst:.2e} (the published values carry 6 decimals)"); +} diff --git a/wai/wai-quantum/src/lib.rs b/wai/wai-quantum/src/lib.rs index 482d44c64..3b77cb02c 100644 --- a/wai/wai-quantum/src/lib.rs +++ b/wai/wai-quantum/src/lib.rs @@ -13,8 +13,20 @@ //! //! Everything is behind a feature; enable only what you use. +// Every receipt seals over a Merkle root; nothing else needs one. +#[cfg(feature = "_signing")] mod merkle; +#[cfg(any( + feature = "quantum_spd", + feature = "quantum_pauli", + feature = "quantum_mps", + feature = "quantum_vml", + feature = "quantum_phasor", + feature = "quantum_kernel" +))] +mod repro; + /// Re-exported so consumers get a **version-matched** signing type: the receipt /// sealing APIs take an `ed25519_dalek::SigningKey`, and a caller who depended on /// their own copy could easily pick an incompatible version. @@ -60,6 +72,12 @@ pub mod quantum_ops; #[cfg(feature = "quantum_pauli")] pub mod quantum_pauli; +#[cfg(feature = "quantum_spd")] +pub mod quantum_spd; + +#[cfg(feature = "quantum_spd_receipt")] +pub mod quantum_spd_receipt; + #[cfg(feature = "quantum_info")] pub mod quantum_info; diff --git a/wai/wai-quantum/src/quantum_kernel.rs b/wai/wai-quantum/src/quantum_kernel.rs index d14dbaebf..fd76ccd05 100644 --- a/wai/wai-quantum/src/quantum_kernel.rs +++ b/wai/wai-quantum/src/quantum_kernel.rs @@ -34,7 +34,7 @@ fn u01(s: &mut u64) -> f64 { fn gauss(s: &mut u64) -> f64 { let u1 = u01(s).max(1e-300); let u2 = u01(s); - (-2.0 * u1.ln()).sqrt() * (2.0 * PI * u2).cos() + (-2.0 * crate::repro::ln(u1)).sqrt() * crate::repro::sin_cos(2.0 * PI * u2).1 } /// A random-Fourier-feature (phasor) map approximating a Gaussian kernel @@ -71,7 +71,7 @@ pub fn features(rff: &Rff, x: &[f64]) -> Vec { .zip(&rff.phases) .map(|(w, &b)| { let dot: f64 = w.iter().zip(x).map(|(wi, xi)| wi * xi).sum(); - rff.scale * (dot + b).cos() + rff.scale * crate::repro::sin_cos(dot + b).1 }) .collect() } @@ -86,7 +86,7 @@ pub fn features_quantized(rff: &Rff, x: &[f64], bits: u32) -> Vec { let dot: f64 = w.iter().zip(x).map(|(wi, xi)| wi * xi).sum(); let t = dot + b; let tq = (t / (2.0 * PI) * levels).round() / levels * 2.0 * PI; - rff.scale * tq.cos() + rff.scale * crate::repro::sin_cos(tq).1 }) .collect() } @@ -202,7 +202,8 @@ pub fn circles(n: usize, seed: u64) -> (Vec>, Vec) { let inner = i % 2 == 0; let ang = u01(&mut st) * 2.0 * PI; let r = if inner { 0.35 * u01(&mut st) } else { 0.9 + 0.35 * u01(&mut st) }; - x.push(vec![r * ang.cos() + 0.03 * gauss(&mut st), r * ang.sin() + 0.03 * gauss(&mut st)]); + let (sin, cos) = crate::repro::sin_cos(ang); + x.push(vec![r * cos + 0.03 * gauss(&mut st), r * sin + 0.03 * gauss(&mut st)]); y.push(if inner { -1.0 } else { 1.0 }); } (x, y) diff --git a/wai/wai-quantum/src/quantum_mps.rs b/wai/wai-quantum/src/quantum_mps.rs index d2368f297..fcf815120 100644 --- a/wai/wai-quantum/src/quantum_mps.rs +++ b/wai/wai-quantum/src/quantum_mps.rs @@ -10,9 +10,9 @@ //! bounded `χ` is **approximate** — the whole point is to discard small Schmidt //! weight — so byte-exact-vs-analytic is not the promise here. What IS preserved is //! **determinism**: the contraction and truncation use only IEEE-754 `+ − × ÷ √` -//! (no FMA). The phase gate's angle is computed with the platform `cos` and `sin`, -//! so the result is bit-reproducible only where the math library returns the same -//! bits for those calls. The engine reports the **max bond dimension** it used and the **retained weight** +//! (no FMA), and the phase gate's `cos` and `sin` are computed in-crate rather than +//! by the platform's math library, so the result is bit-identical on every machine +//! and target. The engine reports the **max bond dimension** it used and the **retained weight** //! (the standard truncation-fidelity heuristic). Verified against the byte-exact dense //! simulator: exact agreement when `χ` is unbounded, graceful, measured degradation //! when it is squeezed. @@ -367,7 +367,8 @@ fn m1(base: BaseGate, param: u16) -> [[C; 2]; 2] { let no = C::new(-1.0, 0.0); let phase = |k: u16| { let th = 2.0 * std::f64::consts::PI / (1u64 << k) as f64; - C::new(th.cos(), th.sin()) + let (s, c) = crate::repro::sin_cos(th); + C::new(c, s) }; match base { BaseGate::I => [[o, z], [z, o]], diff --git a/wai/wai-quantum/src/quantum_pauli.rs b/wai/wai-quantum/src/quantum_pauli.rs index a63bc5278..e97c24ce4 100644 --- a/wai/wai-quantum/src/quantum_pauli.rs +++ b/wai/wai-quantum/src/quantum_pauli.rs @@ -15,10 +15,15 @@ //! Real arithmetic throughout in the Hermitian-Pauli basis (`X, Y, Z` real; the //! rotation coefficients `cos φ, sin φ` real), and the final expectation is //! `Σ_{x-support = 0} c_P` (only all-`I`/`Z` strings survive `⟨0|·|0⟩`). The -//! coefficients are `f64` and the truncation is deterministic, so the result is -//! reproducible; its honest contract, like the MPS backend, is determinism + a -//! stated truncation threshold, not the dense sim's byte-exact guarantee. Verified -//! against the dense simulator's exact expectation values. Supports `n ≤ 64`. +//! coefficients are `f64`. Terms are visited in sorted order wherever +//! coefficients are added, and `cos φ, sin φ` are computed in-crate rather than by +//! the platform's math library, so the result is bit-identical from run to run +//! and across platforms. Its contract, like the MPS backend's, is that +//! reproducibility plus a stated truncation threshold, not the dense sim's +//! byte-exact guarantee. Verified against the dense simulator's exact +//! expectation values. Supports `n ≤ 64` and the dyadic gate set; arbitrary +//! angles and registers up to 1024 qubits are +//! [`quantum_spd`](https://docs.rs/wai-quantum/latest/wai_quantum/quantum_spd/)'s. use crate::quantum::{BaseGate, Circuit, Gate}; use std::collections::HashMap; @@ -31,6 +36,15 @@ pub struct Unsupported(pub String); /// A qubit's factor is `I` (00), `X` (10), `Z` (01), `Y` (11). pub type PauliSum = HashMap<(u64, u64), f64>; +/// The terms of `sum` in key order. A `HashMap` iterates in an order seeded +/// afresh for every map; adding floating-point coefficients in that order +/// would change the last bits of the result from one run to the next. +fn sorted(sum: &PauliSum) -> Vec<((u64, u64), f64)> { + let mut v: Vec<((u64, u64), f64)> = sum.iter().map(|(&k, &c)| (k, c)).collect(); + v.sort_unstable_by_key(|&(k, _)| k); + v +} + #[inline] fn bit(v: u64, q: u8) -> u64 { (v >> q) & 1 @@ -118,8 +132,8 @@ fn conjugate(sum: &PauliSum, g: &Gate, out: &mut PauliSum) -> Result<(), Unsuppo // diagonal rotation: split if let Some(phi) = is_diag_rotation(g) { let q = g.target; - let (cs, sn) = (phi.cos(), phi.sin()); - for (&(x, z), &c) in sum { + let (sn, cs) = crate::repro::sin_cos(phi); + for ((x, z), c) in sorted(sum) { if bit(x, q) == 0 { add(x, z, c, out); // commutes } else { @@ -135,14 +149,14 @@ fn conjugate(sum: &PauliSum, g: &Gate, out: &mut PauliSum) -> Result<(), Unsuppo // gates by control count match (g.base, g.controls.as_slice()) { (BaseGate::I, []) => { - for (&k, &c) in sum { + for (k, c) in sorted(sum) { add(k.0, k.1, c, out); } } (BaseGate::H | BaseGate::S | BaseGate::Sdg | BaseGate::X | BaseGate::Y | BaseGate::Z, []) | (BaseGate::P, []) => { let q = g.target; - for (&(x, z), &c) in sum { + for ((x, z), c) in sorted(sum) { let (nx, nz, s) = conj1(g.base, g.param, bit(x, q), bit(z, q)); let (x2, z2) = set_factor(x, z, q, nx, nz); add(x2, z2, c * s, out); @@ -151,7 +165,7 @@ fn conjugate(sum: &PauliSum, g: &Gate, out: &mut PauliSum) -> Result<(), Unsuppo (BaseGate::X, [ctrl]) => { // CX(c,t): x_t ^= x_c; z_c ^= z_t; sign (-1)^{x_c z_t (x_t ⊕ z_c ⊕ 1)} let (c, t) = (*ctrl, g.target); - for (&(x, z), &co) in sum { + for ((x, z), co) in sorted(sum) { let (xc, zc, xt, zt) = (bit(x, c), bit(z, c), bit(x, t), bit(z, t)); let s = if (xc & zt & (xt ^ zc ^ 1)) == 1 { -1.0 } else { 1.0 }; let (x2, _) = set_factor(x, z, t, xt ^ xc, zt); @@ -199,7 +213,7 @@ pub fn propagate( /// `⟨0…0| (sum) |0…0⟩` — the sum of coefficients of all-`I`/`Z` strings. pub fn expectation_zero(sum: &PauliSum) -> f64 { - sum.iter().filter(|((x, _), _)| *x == 0).map(|(_, c)| c).sum() + sorted(sum).into_iter().filter(|((x, _), _)| *x == 0).map(|(_, c)| c).sum() } /// `⟨0…0| U† Z_q U |0…0⟩` — the single-qubit magnetization, via propagation. @@ -341,6 +355,32 @@ mod tests { assert!(peak < 1 << 20, "Pauli sum stayed bounded: {peak}"); } + #[test] + fn results_are_bit_identical_from_run_to_run() { + // Each HashMap is seeded afresh, so a sum accumulated in a map's + // iteration order would already differ between calls in one process. + let n = 20u8; + let mut c = Circuit::new(n); + for layer in 0..4u8 { + for q in 0..n { + c.h(q); + c.t(q); + } + for q in (layer % 2..n - 1).step_by(2) { + c.cx(q, q + 1); + } + for q in 0..n { + c.t(q); + } + } + let (a, pa) = expect_z(&c, 10, 1e-4).unwrap(); + assert!(pa > 500, "enough terms to merge in many orders: {pa}"); + for _ in 0..6 { + let (b, pb) = expect_z(&c, 10, 1e-4).unwrap(); + assert_eq!((a.to_bits(), pa), (b.to_bits(), pb)); + } + } + #[test] fn rejects_unsupported() { // a non-Clifford controlled-phase (CP(2)) isn't handled. diff --git a/wai/wai-quantum/src/quantum_phasor.rs b/wai/wai-quantum/src/quantum_phasor.rs index 2fa46a309..8d4544bc7 100644 --- a/wai/wai-quantum/src/quantum_phasor.rs +++ b/wai/wai-quantum/src/quantum_phasor.rs @@ -46,7 +46,8 @@ impl C { } #[inline] fn expi(t: f64) -> C { - C::new(t.cos(), t.sin()) + let (s, c) = crate::repro::sin_cos(t); + C::new(c, s) } // --------------------------------------------------------------------------- @@ -56,7 +57,7 @@ fn expi(t: f64) -> C { /// Trainable block `W = RZ(φ)·RY(θ)` applied to a single-qubit amplitude pair. fn apply_w(a: &mut [C; 2], theta: f64, phi: f64) { // RY(θ) - let (c, s) = ((theta * 0.5).cos(), (theta * 0.5).sin()); + let (s, c) = crate::repro::sin_cos(theta * 0.5); let (a0, a1) = (a[0], a[1]); let r0 = C::new(a0.re * c - a1.re * s, a0.im * c - a1.im * s); let r1 = C::new(a0.re * s + a1.re * c, a0.im * s + a1.im * c); diff --git a/wai/wai-quantum/src/quantum_qec.rs b/wai/wai-quantum/src/quantum_qec.rs index 7a2170b94..2652f2c5b 100644 --- a/wai/wai-quantum/src/quantum_qec.rs +++ b/wai/wai-quantum/src/quantum_qec.rs @@ -222,6 +222,7 @@ fn is_zero(row: &[u64]) -> bool { /// [`Gf2Basis`] stops at row-echelon form — a pivot may still appear in earlier /// rows — which is enough to test membership but not to read a row off as an /// independent generator. Encoding needs the stronger form. +#[cfg(feature = "quantum")] fn rref(rows_cols: &[Vec], n: usize) -> Vec<(usize, Vec)> { let mut m: Vec> = rows_cols .iter() diff --git a/wai/wai-quantum/src/quantum_qfhrr.rs b/wai/wai-quantum/src/quantum_qfhrr.rs index 8b3466c08..c132bebae 100644 --- a/wai/wai-quantum/src/quantum_qfhrr.rs +++ b/wai/wai-quantum/src/quantum_qfhrr.rs @@ -42,14 +42,15 @@ fn qfhrr(a: C, bits: u32) -> C { if m <= 0.0 { return C { re: 0.0, im: 0.0 }; } - let theta = a.im.atan2(a.re); + let theta = crate::repro::atan2(a.im, a.re); let levels = (1u64 << bits) as f64; let tq = (theta / (2.0 * PI) * levels).round() / levels * 2.0 * PI; - C { re: m * tq.cos(), im: m * tq.sin() } + let (s, c) = crate::repro::sin_cos(tq); + C { re: m * c, im: m * s } } fn apply_ry(amps: &mut [C], q: u8, theta: f64, bits: u32) { - let (c, s) = ((theta * 0.5).cos(), (theta * 0.5).sin()); + let (s, c) = crate::repro::sin_cos(theta * 0.5); let bit = 1usize << q; let mut i = 0; while i < amps.len() { @@ -63,7 +64,7 @@ fn apply_ry(amps: &mut [C], q: u8, theta: f64, bits: u32) { } } fn apply_rz(amps: &mut [C], q: u8, theta: f64, bits: u32) { - let (cm, sm) = ((theta * 0.5).cos(), (theta * 0.5).sin()); + let (sm, cm) = crate::repro::sin_cos(theta * 0.5); let bit = 1usize << q; for (i, a) in amps.iter_mut().enumerate() { let (er, ei) = if i & bit == 0 { (cm, -sm) } else { (cm, sm) }; diff --git a/wai/wai-quantum/src/quantum_spd.rs b/wai/wai-quantum/src/quantum_spd.rs new file mode 100644 index 000000000..49b844e79 --- /dev/null +++ b/wai/wai-quantum/src/quantum_spd.rs @@ -0,0 +1,1268 @@ +//! Sparse Pauli dynamics at utility scale (`wai.quantum.spd`). +//! +//! Computes `⟨ψ₀| U† O U |ψ₀⟩` for circuits of up to 1024 qubits by evolving the +//! *observable* in the Heisenberg picture as a sum of Pauli strings, +//! `O = Σ c_P P`. A Clifford gate maps each string to one signed string. A Pauli +//! rotation `exp(−iθP/2)` leaves strings that commute with `P` alone and splits +//! the rest in two: `Q ↦ cos θ·Q + sin θ·(−iQP)`. Strings whose coefficient falls +//! below a threshold are dropped. A circuit therefore costs what its +//! non-Clifford rotations and the observable's light cone make it cost, not +//! `2ⁿ`. This is the method that computed the 127-qubit kicked-Ising +//! observables of the heavy-hex "utility" experiment on a laptop. +//! +//! **Arbitrary angles.** The dyadic circuit core is byte-exact because every +//! phase it allows is a power-of-two root of unity. This module takes any real +//! angle, so its contract is the one the crate's other `f64` layers keep: +//! +//! - **Reproducible.** Every result is bit-identical across runs, machines and +//! targets. Terms live in an insertion-ordered table, so the order in which +//! coefficients are added never depends on hashing or on the platform's word +//! size. `sin` and `cos` are computed in-crate from `+ − × ÷`, because the +//! platform's libm may differ in the last place. +//! - **Bounded.** Every dropped coefficient is added to +//! [`SpdResult::truncation_bound`]. A dropped term `δ·P` would have evolved +//! under the remaining gates into `δ` times a unitary conjugate of `P`, whose +//! expectation is at most `|δ|` in magnitude. So `|exact − value|` is at most +//! the bound, up to floating-point rounding. Clifford gates and exact +//! quarter-turn rotations drop nothing. + +use std::fmt; + +pub use crate::repro::sin_cos; + +// --------------------------------------------------------------------------- +// Circuits +// --------------------------------------------------------------------------- + +/// The largest register this module takes: 16 words of 64 qubits. +pub const MAX_QUBITS: u32 = 1024; +/// The largest angle magnitude a rotation may carry. +pub const MAX_ANGLE: f64 = 1_048_576.0; + +/// A single-qubit Pauli factor. +#[derive(Clone, Copy, Debug, PartialEq, Eq, PartialOrd, Ord, Hash)] +pub enum Pauli { + X, + Y, + Z, +} + +/// A Pauli string as `(qubit, factor)` pairs; qubits not named carry `I`. +pub type PauliString = Vec<(u32, Pauli)>; + +/// One instruction. +#[derive(Clone, Debug, PartialEq)] +pub enum Op { + H(u32), + S(u32), + Sdg(u32), + /// √X. + SX(u32), + /// √X†. + SXdg(u32), + X(u32), + Y(u32), + Z(u32), + /// `CX(control, target)`. + CX(u32, u32), + CZ(u32, u32), + Swap(u32, u32), + /// `exp(−iθP/2)` for any Pauli string `P`: `RX`, `RY`, `RZ` for one factor, + /// `RZZ` and the like for two. + Rot { axis: PauliString, theta: f64 }, + /// `exp(−i(kπ/2)P/2)` for `k ∈ {1, 2, 3}`: a Pauli rotation through a + /// multiple of π/2. These are Clifford, and applied exactly. + QuarterRot { axis: PauliString, k: u8 }, +} + +/// A circuit on `n` qubits. +#[derive(Clone, Debug, PartialEq)] +pub struct RotCircuit { + pub n: u32, + pub ops: Vec, +} + +/// Why a circuit, observable or initial state was refused. +#[derive(Clone, Debug, PartialEq)] +pub enum SpdError { + TooManyQubits(u32), + QubitOutOfRange { qubit: u32, n: u32 }, + /// A two-qubit gate or Pauli string naming one qubit twice. + RepeatedQubit(u32), + EmptyAxis, + BadAngle(f64), + BadQuarter(u8), + /// A dyadic gate with no Pauli-rotation form here (a controlled `H`). + Unsupported(String), + /// A receipt's meter returned without running the job. + NotRun, +} + +impl fmt::Display for SpdError { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + match self { + SpdError::TooManyQubits(n) => write!(f, "{n} qubits; at most {MAX_QUBITS}"), + SpdError::QubitOutOfRange { qubit, n } => write!(f, "qubit {qubit} of a {n}-qubit register"), + SpdError::RepeatedQubit(q) => write!(f, "qubit {q} named twice"), + SpdError::EmptyAxis => write!(f, "a rotation about the identity"), + SpdError::BadAngle(a) => write!(f, "angle {a} is not finite or exceeds {MAX_ANGLE}"), + SpdError::BadQuarter(k) => write!(f, "quarter turns {k}; expected 1, 2 or 3"), + SpdError::Unsupported(g) => write!(f, "unsupported gate: {g}"), + SpdError::NotRun => write!(f, "the meter did not run the job"), + } + } +} + +impl std::error::Error for SpdError {} + +impl RotCircuit { + pub fn new(n: u32) -> Self { + RotCircuit { n, ops: Vec::new() } + } + + pub fn push(&mut self, op: Op) -> &mut Self { + self.ops.push(op); + self + } + + pub fn rx(&mut self, q: u32, theta: f64) -> &mut Self { + self.push(Op::Rot { axis: vec![(q, Pauli::X)], theta }) + } + + pub fn ry(&mut self, q: u32, theta: f64) -> &mut Self { + self.push(Op::Rot { axis: vec![(q, Pauli::Y)], theta }) + } + + pub fn rz(&mut self, q: u32, theta: f64) -> &mut Self { + self.push(Op::Rot { axis: vec![(q, Pauli::Z)], theta }) + } + + /// `exp(−iθ Z⊗Z/2)`. + pub fn rzz(&mut self, a: u32, b: u32, theta: f64) -> &mut Self { + self.push(Op::Rot { axis: vec![(a, Pauli::Z), (b, Pauli::Z)], theta }) + } + + /// `exp(−iθ X⊗X/2)`. + pub fn rxx(&mut self, a: u32, b: u32, theta: f64) -> &mut Self { + self.push(Op::Rot { axis: vec![(a, Pauli::X), (b, Pauli::X)], theta }) + } + + /// Check every qubit index, axis and angle. + pub fn validate(&self) -> Result<(), SpdError> { + if self.n == 0 || self.n > MAX_QUBITS { + return Err(SpdError::TooManyQubits(self.n)); + } + let q_ok = |q: u32| { + if q < self.n { + Ok(()) + } else { + Err(SpdError::QubitOutOfRange { qubit: q, n: self.n }) + } + }; + let pair = |a: u32, b: u32| { + q_ok(a)?; + q_ok(b)?; + if a == b { Err(SpdError::RepeatedQubit(a)) } else { Ok(()) } + }; + for op in &self.ops { + match op { + Op::H(q) | Op::S(q) | Op::Sdg(q) | Op::SX(q) | Op::SXdg(q) | Op::X(q) | Op::Y(q) | Op::Z(q) => { + q_ok(*q)? + } + Op::CX(a, b) | Op::CZ(a, b) | Op::Swap(a, b) => pair(*a, *b)?, + Op::Rot { axis, theta } => { + check_string(axis, self.n)?; + if axis.is_empty() { + return Err(SpdError::EmptyAxis); + } + if !theta.is_finite() || theta.abs() > MAX_ANGLE { + return Err(SpdError::BadAngle(*theta)); + } + } + Op::QuarterRot { axis, k } => { + check_string(axis, self.n)?; + if axis.is_empty() { + return Err(SpdError::EmptyAxis); + } + if !(1..=3).contains(k) { + return Err(SpdError::BadQuarter(*k)); + } + } + } + } + Ok(()) + } +} + +fn check_string(s: &[(u32, Pauli)], n: u32) -> Result<(), SpdError> { + for (i, (q, _)) in s.iter().enumerate() { + if *q >= n { + return Err(SpdError::QubitOutOfRange { qubit: *q, n }); + } + if s[..i].iter().any(|(r, _)| r == q) { + return Err(SpdError::RepeatedQubit(*q)); + } + } + Ok(()) +} + +// --------------------------------------------------------------------------- +// Pauli strings as bit words +// --------------------------------------------------------------------------- + +/// A Hermitian Pauli string on up to `64·W` qubits: factor `X^x Z^z` at each +/// qubit, with `Y` (both bits set) meaning `i·XZ`, so every string is Hermitian. +#[derive(Clone, Copy, PartialEq, Eq)] +struct Key { + x: [u64; W], + z: [u64; W], +} + +impl Key { + const ZERO: Self = Key { x: [0; W], z: [0; W] }; + + fn from_string(s: &[(u32, Pauli)]) -> Self { + let mut k = Self::ZERO; + for &(q, p) in s { + let (w, b) = ((q / 64) as usize, 1u64 << (q % 64)); + match p { + Pauli::X => k.x[w] |= b, + Pauli::Z => k.z[w] |= b, + Pauli::Y => { + k.x[w] |= b; + k.z[w] |= b; + } + } + } + k + } + + #[inline] + fn get(&self, q: u32) -> (bool, bool) { + let (w, b) = ((q / 64) as usize, q % 64); + ((self.x[w] >> b) & 1 == 1, (self.z[w] >> b) & 1 == 1) + } + + #[inline] + fn set(&mut self, q: u32, x: bool, z: bool) { + let (w, b) = ((q / 64) as usize, 1u64 << (q % 64)); + self.x[w] = (self.x[w] & !b) | if x { b } else { 0 }; + self.z[w] = (self.z[w] & !b) | if z { b } else { 0 }; + } + + /// Whether two strings commute: the symplectic product is even. + #[inline] + fn commutes(&self, o: &Self) -> bool { + let mut p = 0u32; + for w in 0..W { + p += (self.x[w] & o.z[w]).count_ones() + (self.z[w] & o.x[w]).count_ones(); + } + p & 1 == 0 + } + + /// `self · o = i^e · (self ⊕ o)`: the product string and the exponent `e` + /// (mod 4) of its phase. + #[inline] + fn mul(&self, o: &Self) -> (Self, u32) { + let mut out = Self::ZERO; + let mut e = 0u32; + for w in 0..W { + let (x, z) = (self.x[w] ^ o.x[w], self.z[w] ^ o.z[w]); + out.x[w] = x; + out.z[w] = z; + e = e.wrapping_add((self.x[w] & self.z[w]).count_ones()); + e = e.wrapping_add((o.x[w] & o.z[w]).count_ones()); + e = e.wrapping_add(2 * (self.z[w] & o.x[w]).count_ones()); + e = e.wrapping_sub((x & z).count_ones()); + } + (out, e & 3) + } + + fn weight(&self) -> u32 { + (0..W).map(|w| (self.x[w] | self.z[w]).count_ones()).sum() + } + + fn hash(&self) -> u64 { + let mut h = 0x9e37_79b9_7f4a_7c15u64; + for w in 0..W { + h = mix(h ^ self.x[w]); + h = mix(h ^ self.z[w]); + } + h + } +} + +#[inline] +fn mix(mut z: u64) -> u64 { + z = (z ^ (z >> 30)).wrapping_mul(0xbf58_476d_1ce4_e5b9); + z = (z ^ (z >> 27)).wrapping_mul(0x94d0_49bb_1331_11eb); + z ^ (z >> 31) +} + +/// Pauli terms in insertion order, with a hash index on the side. Iteration +/// follows the order terms were first inserted, so the order coefficients are +/// added in — and with it every bit of the result — is the same on every +/// platform. +/// +/// A dropped term keeps its slot with coefficient `0.0`, so dropping costs +/// nothing to the index, and a term that reappears is revived in place. The +/// table is compacted, in order, when dead slots outnumber live ones, and +/// before a Clifford gate rewrites every key anyway. When it happens depends +/// only on the computation, so the result stays the same everywhere. +struct Terms { + keys: Vec>, + coef: Vec, + index: Vec, + stale: bool, + /// Slots whose coefficient is `0.0`. + dead: usize, +} + +const EMPTY: u32 = u32::MAX; + +impl Terms { + fn new() -> Self { + Terms { keys: Vec::new(), coef: Vec::new(), index: vec![EMPTY; 16], stale: false, dead: 0 } + } + + /// Terms with a nonzero coefficient. + fn live(&self) -> usize { + self.keys.len() - self.dead + } + + /// Remove dead slots, keeping the order of the rest. + fn compact(&mut self) { + let mut keep = 0; + for i in 0..self.keys.len() { + if self.coef[i] != 0.0 { + self.keys[keep] = self.keys[i]; + self.coef[keep] = self.coef[i]; + keep += 1; + } + } + self.keys.truncate(keep); + self.coef.truncate(keep); + self.dead = 0; + self.stale = true; + } + + fn reindex(&mut self) { + let mut cap = 16usize; + while cap < 2 * self.keys.len() + 2 { + cap *= 2; + } + self.index.clear(); + self.index.resize(cap, EMPTY); + let mask = (cap - 1) as u64; + for (i, k) in self.keys.iter().enumerate() { + let mut slot = (k.hash() & mask) as usize; + while self.index[slot] != EMPTY { + slot = (slot + 1) & (cap - 1); + } + self.index[slot] = i as u32; + } + self.stale = false; + } + + /// Add `c` to the coefficient of `k`, inserting it at the end if absent. + fn add(&mut self, k: Key, c: f64) { + debug_assert!(!self.stale); + if 2 * (self.keys.len() + 1) > self.index.len() { + self.reindex(); + } + let cap = self.index.len(); + let mut slot = (k.hash() & (cap as u64 - 1)) as usize; + loop { + let i = self.index[slot]; + if i == EMPTY { + self.index[slot] = self.keys.len() as u32; + self.keys.push(k); + self.coef.push(c); + if c == 0.0 { + self.dead += 1; + } + return; + } + if self.keys[i as usize] == k { + let before = self.coef[i as usize]; + let after = before + c; + self.coef[i as usize] = after; + match (before == 0.0, after == 0.0) { + (true, false) => self.dead -= 1, + (false, true) => self.dead += 1, + _ => {} + } + return; + } + slot = (slot + 1) & (cap - 1); + } + } + + /// Drop terms with `|c| < threshold`, or heavier than `max_weight`. + /// Returns the sum of dropped `|c|`. + fn truncate(&mut self, threshold: f64, max_weight: Option) -> f64 { + let mut dropped = 0.0; + for i in 0..self.keys.len() { + let c = self.coef[i]; + if c == 0.0 { + continue; + } + if c.abs() < threshold || max_weight.is_some_and(|m| self.keys[i].weight() > m) { + dropped += c.abs(); + self.coef[i] = 0.0; + self.dead += 1; + } + } + if 2 * self.dead > self.keys.len() { + self.compact(); + } + dropped + } +} + +// --------------------------------------------------------------------------- +// Conjugation by one gate +// --------------------------------------------------------------------------- + +/// `G† (X^x Z^z) G` for a single-qubit Clifford `G` on one factor: the new +/// factor and whether the sign flips. +fn conj1(op: &Op, x: bool, z: bool) -> (bool, bool, bool) { + // factors: I (0,0), X (1,0), Z (0,1), Y (1,1) + match (op, x, z) { + (_, false, false) => (false, false, false), + (Op::H(_), true, false) => (false, true, false), // X → Z + (Op::H(_), false, true) => (true, false, false), // Z → X + (Op::H(_), true, true) => (true, true, true), // Y → −Y + (Op::S(_), true, false) => (true, true, true), // X → −Y + (Op::S(_), true, true) => (true, false, false), // Y → X + (Op::Sdg(_), true, false) => (true, true, false), // X → Y + (Op::Sdg(_), true, true) => (true, false, true), // Y → −X + (Op::SX(_), false, true) => (true, true, false), // Z → Y + (Op::SX(_), true, true) => (false, true, true), // Y → −Z + (Op::SXdg(_), false, true) => (true, true, true), // Z → −Y + (Op::SXdg(_), true, true) => (false, true, false), // Y → Z + (Op::X(_), _, true) => (x, z, true), // Z, Y → − + (Op::Y(_), true, false) | (Op::Y(_), false, true) => (x, z, true), // X, Z → − + (Op::Z(_), true, _) => (x, z, true), // X, Y → − + _ => (x, z, false), + } +} + +/// Apply a Clifford gate to every term in place. Conjugation by a unitary is a +/// bijection on Pauli strings, so no two terms collide. +fn clifford(t: &mut Terms, op: &Op) { + if t.dead > 0 { + t.compact(); + } + for i in 0..t.keys.len() { + let k = &mut t.keys[i]; + let neg = match *op { + Op::H(q) | Op::S(q) | Op::Sdg(q) | Op::SX(q) | Op::SXdg(q) | Op::X(q) | Op::Y(q) | Op::Z(q) => { + let (x, z) = k.get(q); + let (nx, nz, neg) = conj1(op, x, z); + k.set(q, nx, nz); + neg + } + Op::CX(c, tq) => { + let ((xc, zc), (xt, zt)) = (k.get(c), k.get(tq)); + k.set(tq, xt ^ xc, zt); + k.set(c, xc, zc ^ zt); + xc && zt && !(xt ^ zc) + } + Op::CZ(a, b) => { + let ((xa, za), (xb, zb)) = (k.get(a), k.get(b)); + k.set(a, xa, za ^ xb); + k.set(b, xb, zb ^ xa); + xa && xb && (za ^ zb) + } + Op::Swap(a, b) => { + let (fa, fb) = (k.get(a), k.get(b)); + k.set(a, fb.0, fb.1); + k.set(b, fa.0, fa.1); + false + } + _ => unreachable!("not a fixed Clifford"), + }; + if neg { + t.coef[i] = -t.coef[i]; + } + } + t.stale = true; +} + +/// Conjugate by `exp(−i(kπ/2)P/2)`, exactly. A string anticommuting with `P` +/// becomes `±(−iQP)` for `k = 1` or `3`, and `−Q` for `k = 2`. +fn quarter(t: &mut Terms, p: &Key, k: u8) { + if t.dead > 0 { + t.compact(); + } + for i in 0..t.keys.len() { + let q = t.keys[i]; + if q.commutes(p) { + continue; + } + if k == 2 { + t.coef[i] = -t.coef[i]; + continue; + } + let (r, e) = q.mul(p); + // −i·QP = i^{e+3}·R; e is odd here, so the factor is ±1. + let mut neg = (e + 3) & 3 == 2; + if k == 3 { + neg = !neg; + } + t.keys[i] = r; + if neg { + t.coef[i] = -t.coef[i]; + } + } + t.stale = true; +} + +/// Conjugate by `exp(−iθP/2)`: `Q ↦ cos θ·Q + sin θ·(−iQP)` for every `Q` that +/// anticommutes with `P`. The new terms are merged after every old one has been +/// scaled, so each old coefficient is read before it changes. +fn rotate(t: &mut Terms, p: &Key, sin: f64, cos: f64, pending: &mut Vec<(Key, f64)>) { + if t.stale { + t.reindex(); + } + pending.clear(); + for i in 0..t.keys.len() { + let c = t.coef[i]; + if c == 0.0 { + continue; + } + let q = t.keys[i]; + if q.commutes(p) { + continue; + } + let (r, e) = q.mul(p); + let v = c * sin; + pending.push((r, if (e + 3) & 3 == 2 { -v } else { v })); + t.coef[i] = c * cos; + } + for &(r, v) in pending.iter() { + t.add(r, v); + } +} + +// --------------------------------------------------------------------------- +// Propagation +// --------------------------------------------------------------------------- + +/// How hard to truncate. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct SpdConfig { + /// Drop terms whose coefficient magnitude falls below this after a rotation. + /// `0.0` drops only exact zeros: the computation is then exact up to + /// rounding, and its cost may be exponential. + pub threshold: f64, + /// Also drop strings acting on more than this many qubits. + pub max_weight: Option, +} + +impl Default for SpdConfig { + fn default() -> Self { + SpdConfig { threshold: 1e-4, max_weight: None } + } +} + +/// An expectation value, what it cost, and how far truncation can have moved it. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct SpdResult { + /// `⟨ψ₀| U† O U |ψ₀⟩` as computed. + pub value: f64, + /// The sum of every dropped coefficient's magnitude: a bound on + /// `|exact − value|`, up to floating-point rounding. + pub truncation_bound: f64, + /// The most terms the sum held after any gate. + pub peak_terms: usize, + /// Terms left when propagation finished. + pub final_terms: usize, + /// Term-gate updates performed: a portable measure of the work, equal on + /// every machine. + pub term_updates: u64, +} + +/// `⟨b| U† O U |b⟩` for the computational basis state `|b⟩` whose set qubits +/// are `ones` (empty for `|0…0⟩`), with `O` a weighted sum of Pauli strings. +pub fn expectation( + circuit: &RotCircuit, + observable: &[(f64, PauliString)], + ones: &[u32], + cfg: &SpdConfig, +) -> Result { + circuit.validate()?; + for (_, s) in observable { + check_string(s, circuit.n)?; + } + for &q in ones { + if q >= circuit.n { + return Err(SpdError::QubitOutOfRange { qubit: q, n: circuit.n }); + } + } + Ok(match circuit.n.div_ceil(64) { + 1 => run::<1>(circuit, observable, ones, cfg), + 2 => run::<2>(circuit, observable, ones, cfg), + 3 | 4 => run::<4>(circuit, observable, ones, cfg), + 5..=8 => run::<8>(circuit, observable, ones, cfg), + _ => run::<16>(circuit, observable, ones, cfg), + }) +} + +/// `⟨0…0| U† Z_q U |0…0⟩`. +pub fn expect_z(circuit: &RotCircuit, q: u32, cfg: &SpdConfig) -> Result { + expectation(circuit, &[(1.0, vec![(q, Pauli::Z)])], &[], cfg) +} + +fn run( + circuit: &RotCircuit, + observable: &[(f64, PauliString)], + ones: &[u32], + cfg: &SpdConfig, +) -> SpdResult { + let mut t = Terms::::new(); + for (c, s) in observable { + t.add(Key::from_string(s), *c); + } + let mut bound = t.truncate(0.0, None); + let mut peak = t.live(); + let mut work = 0u64; + let mut pending = Vec::new(); + for op in circuit.ops.iter().rev() { + work += t.live() as u64; + match op { + Op::Rot { axis, theta } => { + let (s, c) = sin_cos(*theta); + rotate(&mut t, &Key::from_string(axis), s, c, &mut pending); + bound += t.truncate(cfg.threshold, cfg.max_weight); + } + Op::QuarterRot { axis, k } => quarter(&mut t, &Key::from_string(axis), *k), + other => clifford(&mut t, other), + } + peak = peak.max(t.live()); + } + let b = Key::::from_string(&ones.iter().map(|&q| (q, Pauli::Z)).collect::>()); + let mut value = 0.0; + for i in 0..t.keys.len() { + let k = &t.keys[i]; + if t.coef[i] != 0.0 && k.x.iter().all(|&w| w == 0) { + let flips: u32 = (0..W).map(|w| (k.z[w] & b.z[w]).count_ones()).sum(); + value += if flips & 1 == 1 { -t.coef[i] } else { t.coef[i] }; + } + } + SpdResult { value, truncation_bound: bound, peak_terms: peak, final_terms: t.live(), term_updates: work } +} + +// --------------------------------------------------------------------------- +// The heavy-hex kicked-Ising experiment +// --------------------------------------------------------------------------- + +/// The 127-qubit heavy-hex coupling map: seven rows (14, 15, 15, 15, 15, 15 +/// and 14 qubits) joined by six rows of four bridge qubits. Returns the 144 +/// edges, each with the smaller qubit first. +pub fn heavy_hex_127() -> Vec<(u32, u32)> { + let rows: [(u32, u32); 7] = [(0, 13), (18, 32), (37, 51), (56, 70), (75, 89), (94, 108), (113, 126)]; + // Each bridge row: (first bridge qubit, column offset into the row above, + // column offset into the row below). + let bridges: [(u32, u32, u32); 6] = [(14, 0, 0), (33, 2, 2), (52, 0, 0), (71, 2, 2), (90, 0, 0), (109, 2, 2)]; + let mut edges = Vec::with_capacity(144); + for &(a, b) in &rows { + for q in a..b { + edges.push((q, q + 1)); + } + } + for (r, &(first, up, down)) in bridges.iter().enumerate() { + let (above, below) = (rows[r].0, rows[r + 1].0); + for j in 0..4 { + let top = above + up + 4 * j; + // Row 6 starts one column right of the rows above it, so its + // bridge columns are one index earlier. + let bottom = below + down + 4 * j - if r == 5 { 1 } else { 0 }; + edges.push((top.min(first + j), top.max(first + j))); + edges.push(((first + j).min(bottom), (first + j).max(bottom))); + } + } + edges.sort_unstable(); + edges +} + +/// The kicked-Ising circuit of the heavy-hex utility experiment: `steps` +/// Trotter steps, each `RX(θ_h)` on every qubit and then `RZZ(−π/2)` on every +/// edge (applied exactly, as quarter turns). +pub fn kicked_ising(n: u32, edges: &[(u32, u32)], theta_h: f64, steps: u32) -> RotCircuit { + let mut c = RotCircuit::new(n); + for _ in 0..steps { + for q in 0..n { + c.rx(q, theta_h); + } + for &(a, b) in edges { + // RZZ(−π/2) = exp(+iπ/4·ZZ) = exp(−i(3π/2)ZZ/2) up to a global phase. + c.push(Op::QuarterRot { axis: vec![(a, Pauli::Z), (b, Pauli::Z)], k: 3 }); + } + } + c +} + +// --------------------------------------------------------------------------- +// From the dyadic core +// --------------------------------------------------------------------------- + +/// The dyadic circuit as Pauli rotations. Fixed Cliffords map one to one; `T` +/// and `P(k)` become `Z` rotations; a controlled phase on any number of +/// controls becomes rotations about `Z` strings, from +/// `|1…1⟩⟨1…1| = Π (I − Z)/2`; controlled `X` and `Y` beyond one control are +/// that, conjugated by `H` (and `S`) on the target. Global phases drop out of a +/// Heisenberg-picture expectation and are not tracked. +#[cfg(feature = "quantum")] +pub fn from_dyadic(c: &crate::quantum::Circuit) -> Result { + use crate::quantum::BaseGate as B; + let mut out = RotCircuit::new(u32::from(c.n_qubits)); + let tau = 2.0 * std::f64::consts::PI; + for g in &c.ops { + let t = u32::from(g.target); + let ctrls: Vec = g.controls.iter().map(|&q| u32::from(q)).collect(); + // The phase a diagonal base applies to |1⟩, in turns of 2π. + let phase = |base: B| -> Option { + match base { + B::Z => Some(0.5), + B::S => Some(0.25), + B::Sdg => Some(-0.25), + B::T => Some(0.125), + B::Tdg => Some(-0.125), + B::P => Some(1.0 / (1u64 << g.param) as f64), + _ => None, + } + }; + match (g.base, ctrls.as_slice()) { + (B::I, _) => {} + (B::H, []) => _ = out.push(Op::H(t)), + (B::X, []) => _ = out.push(Op::X(t)), + (B::Y, []) => _ = out.push(Op::Y(t)), + (B::Z, []) => _ = out.push(Op::Z(t)), + (B::S, []) => _ = out.push(Op::S(t)), + (B::Sdg, []) => _ = out.push(Op::Sdg(t)), + (B::X, [a]) => _ = out.push(Op::CX(*a, t)), + (B::Z, [a]) => _ = out.push(Op::CZ(*a, t)), + (B::Y, [a]) => { + out.push(Op::Sdg(t)).push(Op::CX(*a, t)).push(Op::S(t)); + } + (B::X, _) => { + out.push(Op::H(t)); + controlled_phase(&mut out, &ctrls, t, 0.5 * tau); + out.push(Op::H(t)); + } + (B::Y, _) => { + out.push(Op::Sdg(t)).push(Op::H(t)); + controlled_phase(&mut out, &ctrls, t, 0.5 * tau); + out.push(Op::H(t)).push(Op::S(t)); + } + (base, _) => match phase(base) { + Some(turns) => controlled_phase(&mut out, &ctrls, t, turns * tau), + None => return Err(SpdError::Unsupported(format!("{base:?} on {} controls", ctrls.len()))), + }, + } + } + Ok(out) +} + +/// `diag(…, e^{iφ})` on `controls ∪ {target}`: the phase `φ` on `|1…1⟩` only. +/// With `m` qubits, `|1…1⟩⟨1…1| = 2⁻ᵐ Σ_S (−1)^{|S|} Z_S`, so the gate is +/// `Π_S exp(iφ(−1)^{|S|} Z_S / 2ᵐ)` over nonempty subsets `S` (the empty one is +/// a global phase): a rotation `exp(−iθZ_S/2)` with `θ = −2φ(−1)^{|S|}/2ᵐ`. +#[cfg(feature = "quantum")] +fn controlled_phase(out: &mut RotCircuit, controls: &[u32], target: u32, phi: f64) { + let mut qs = controls.to_vec(); + qs.push(target); + let m = qs.len() as u32; + for mask in 1u32..(1 << m) { + let axis: PauliString = (0..m).filter(|i| mask >> i & 1 == 1).map(|i| (qs[i as usize], Pauli::Z)).collect(); + let sign = if axis.len() % 2 == 0 { 1.0 } else { -1.0 }; + out.push(Op::Rot { axis, theta: -2.0 * phi * sign / f64::from(1u32 << m) }); + } +} + +#[cfg(test)] +mod tests { + use super::*; + + // ---- an independent Schrödinger-picture reference ---------------------- + + type C = (f64, f64); + + fn cmul(a: C, b: C) -> C { + (a.0 * b.0 - a.1 * b.1, a.0 * b.1 + a.1 * b.0) + } + + /// `P|ψ⟩` for a Pauli string, on a dense state. + fn apply_pauli(psi: &[C], s: &[(u32, Pauli)]) -> Vec { + let mut out = vec![(0.0, 0.0); psi.len()]; + for (j, &a) in psi.iter().enumerate() { + let mut k = j; + let mut ph: C = (1.0, 0.0); + for &(q, p) in s { + let bit = (j >> q) & 1; + match p { + Pauli::X => k ^= 1 << q, + Pauli::Y => { + k ^= 1 << q; + // Y|0⟩ = i|1⟩, Y|1⟩ = −i|0⟩ + ph = cmul(ph, if bit == 0 { (0.0, 1.0) } else { (0.0, -1.0) }); + } + Pauli::Z => { + if bit == 1 { + ph = (-ph.0, -ph.1); + } + } + } + } + let v = cmul(ph, a); + out[k].0 += v.0; + out[k].1 += v.1; + } + out + } + + fn one_qubit(psi: &mut [C], q: u32, m: [[C; 2]; 2]) { + for j in 0..psi.len() { + if (j >> q) & 1 == 0 { + let k = j | (1 << q); + let (a, b) = (psi[j], psi[k]); + let add = |x: C, y: C| (x.0 + y.0, x.1 + y.1); + psi[j] = add(cmul(m[0][0], a), cmul(m[0][1], b)); + psi[k] = add(cmul(m[1][0], a), cmul(m[1][1], b)); + } + } + } + + /// Run `circuit` on `|ones⟩` as a dense state, with the platform's `sin` + /// and `cos`: nothing shared with the code under test but the `Op` type. + fn dense(circuit: &RotCircuit, ones: &[u32]) -> Vec { + let n = circuit.n; + let mut psi = vec![(0.0, 0.0); 1 << n]; + psi[ones.iter().fold(0usize, |a, &q| a | 1 << q)] = (1.0, 0.0); + let r = std::f64::consts::FRAC_1_SQRT_2; + let (o, z, i, ni) = ((1.0, 0.0), (0.0, 0.0), (0.0, 1.0), (0.0, -1.0)); + for op in &circuit.ops { + match op { + Op::H(q) => one_qubit(&mut psi, *q, [[(r, 0.0), (r, 0.0)], [(r, 0.0), (-r, 0.0)]]), + Op::S(q) => one_qubit(&mut psi, *q, [[o, z], [z, i]]), + Op::Sdg(q) => one_qubit(&mut psi, *q, [[o, z], [z, ni]]), + Op::SX(q) => one_qubit(&mut psi, *q, [[(0.5, 0.5), (0.5, -0.5)], [(0.5, -0.5), (0.5, 0.5)]]), + Op::SXdg(q) => one_qubit(&mut psi, *q, [[(0.5, -0.5), (0.5, 0.5)], [(0.5, 0.5), (0.5, -0.5)]]), + Op::X(q) => psi = apply_pauli(&psi, &[(*q, Pauli::X)]), + Op::Y(q) => psi = apply_pauli(&psi, &[(*q, Pauli::Y)]), + Op::Z(q) => psi = apply_pauli(&psi, &[(*q, Pauli::Z)]), + Op::CX(c, t) => { + for j in 0..psi.len() { + if (j >> c) & 1 == 1 && (j >> t) & 1 == 0 { + psi.swap(j, j | 1 << t); + } + } + } + Op::CZ(a, b) => { + for (j, v) in psi.iter_mut().enumerate() { + if (j >> a) & 1 == 1 && (j >> b) & 1 == 1 { + *v = (-v.0, -v.1); + } + } + } + Op::Swap(a, b) => { + for j in 0..psi.len() { + if (j >> a) & 1 == 1 && (j >> b) & 1 == 0 { + psi.swap(j, (j & !(1 << a)) | 1 << b); + } + } + } + Op::Rot { axis, theta } => { + let p = apply_pauli(&psi, axis); + let (c, s) = ((theta / 2.0).cos(), (theta / 2.0).sin()); + for (v, w) in psi.iter_mut().zip(&p) { + // cos(θ/2)ψ − i sin(θ/2) Pψ + *v = (c * v.0 + s * w.1, c * v.1 - s * w.0); + } + } + Op::QuarterRot { axis, k } => { + let theta = f64::from(*k) * std::f64::consts::FRAC_PI_2; + let p = apply_pauli(&psi, axis); + let (c, s) = ((theta / 2.0).cos(), (theta / 2.0).sin()); + for (v, w) in psi.iter_mut().zip(&p) { + *v = (c * v.0 + s * w.1, c * v.1 - s * w.0); + } + } + } + } + psi + } + + fn dense_expect(circuit: &RotCircuit, obs: &[(u32, Pauli)], ones: &[u32]) -> f64 { + let psi = dense(circuit, ones); + let p = apply_pauli(&psi, obs); + psi.iter().zip(&p).map(|(a, b)| a.0 * b.0 + a.1 * b.1).sum() + } + + /// A small deterministic generator for test circuits. + struct Rng(u64); + impl Rng { + fn next(&mut self) -> u64 { + self.0 = self.0.wrapping_add(0x9e37_79b9_7f4a_7c15); + mix(self.0) + } + fn below(&mut self, n: u64) -> u64 { + self.next() % n + } + fn angle(&mut self) -> f64 { + (self.next() >> 11) as f64 / (1u64 << 53) as f64 * 8.0 - 4.0 + } + fn pauli(&mut self) -> Pauli { + [Pauli::X, Pauli::Y, Pauli::Z][self.below(3) as usize] + } + fn string(&mut self, n: u32, max: usize) -> PauliString { + let mut qs: Vec = (0..n).collect(); + let len = 1 + self.below(max as u64) as usize; + let mut s = Vec::new(); + for _ in 0..len.min(n as usize) { + let i = self.below(qs.len() as u64) as usize; + s.push((qs.swap_remove(i), self.pauli())); + } + s + } + } + + fn random_circuit(rng: &mut Rng, n: u32, len: usize) -> RotCircuit { + let mut c = RotCircuit::new(n); + for _ in 0..len { + let q = rng.below(u64::from(n)) as u32; + let mut r = rng.below(u64::from(n - 1)) as u32; + if r >= q { + r += 1; + } + let op = match rng.below(14) { + 0 => Op::H(q), + 1 => Op::S(q), + 2 => Op::Sdg(q), + 3 => Op::SX(q), + 4 => Op::SXdg(q), + 5 => Op::X(q), + 6 => Op::Y(q), + 7 => Op::Z(q), + 8 => Op::CX(q, r), + 9 => Op::CZ(q, r), + 10 => Op::Swap(q, r), + 11 => Op::QuarterRot { axis: rng.string(n, 3), k: 1 + rng.below(3) as u8 }, + _ => Op::Rot { axis: rng.string(n, 3), theta: rng.angle() }, + }; + c.push(op); + } + c + } + + // ---- trigonometry -------------------------------------------------------- + + #[test] + fn sin_cos_is_within_an_ulp_of_the_platform() { + let mut rng = Rng(7); + let ulp = |a: f64, b: f64| (a.to_bits() as i64 - b.to_bits() as i64).unsigned_abs(); + for i in 0..200_000 { + let x = match i % 4 { + 0 => rng.angle(), + 1 => rng.angle() * 1e-6, + 2 => rng.angle() * 1e4, + _ => rng.angle() * 0.2, + }; + let (s, c) = sin_cos(x); + // Both are faithful roundings; they may differ by one ulp at most, + // and only where the true value is close to a rounding boundary. + assert!(ulp(s, x.sin()) <= 1 || (s - x.sin()).abs() < 1e-300, "sin {x}: {s} vs {}", x.sin()); + assert!(ulp(c, x.cos()) <= 1, "cos {x}: {c} vs {}", x.cos()); + } + } + + #[test] + fn sin_cos_bits_are_pinned() { + // Bit patterns this implementation must produce on every target: they + // are what the cross-architecture run (wasm32-wasip2) is checked against. + // Each is within 0.51 ulp of the true value, checked against an + // 80-digit Taylor series when it was pinned. + let pins: [(f64, u64, u64); 4] = [ + (0.3, 0x3fd2_e9cd_95ba_ba33, 0x3fee_921d_d42f_09ba), + (1.0, 0x3fea_ed54_8f09_0cee, 0x3fe1_4a28_0fb5_068c), + (-2.5, 0xbfe3_26af_0dcf_cab0, 0xbfe9_a2f7_ef85_8b7d), + (100.0, 0xbfe0_3425_b78c_4db8, 0x3feb_981d_bf66_5fdf), + ]; + for (x, sb, cb) in pins { + let (s, c) = sin_cos(x); + assert_eq!((s.to_bits(), c.to_bits()), (sb, cb), "x = {x}: sin {s:e} cos {c:e}"); + } + } + + #[test] + fn quarter_turns_are_exact() { + for (k, x) in [(1, std::f64::consts::FRAC_PI_2), (2, std::f64::consts::PI)] { + let (s, c) = sin_cos(x); + // Not exactly 0 and ±1 — π/2 is not representable — which is why the + // Clifford rotations have their own exact op. + assert!(c.abs() < 1e-15 || k == 2, "k = {k}"); + assert!((s.abs() - if k == 1 { 1.0 } else { 0.0 }).abs() < 1e-15); + } + } + + // ---- Pauli algebra against matrices -------------------------------------- + + #[test] + fn every_clifford_conjugates_every_pauli_as_its_matrix_does() { + // For each gate G and each Pauli P on two qubits, the Heisenberg rule's + // G†PG must give the same expectation on random states as conjugating + // the state instead: ⟨ψ|G†PG|ψ⟩ = ⟨Gψ|P|Gψ⟩. + let gates = [ + Op::H(0), Op::S(0), Op::Sdg(0), Op::SX(0), Op::SXdg(0), Op::X(0), Op::Y(0), Op::Z(0), + Op::H(1), Op::SX(1), Op::CX(0, 1), Op::CX(1, 0), Op::CZ(0, 1), Op::Swap(0, 1), + ]; + let mut rng = Rng(11); + let paulis = [None, Some(Pauli::X), Some(Pauli::Y), Some(Pauli::Z)]; + for g in &gates { + for p0 in paulis { + for p1 in paulis { + let obs: PauliString = [(0, p0), (1, p1)].iter().filter_map(|&(q, p)| p.map(|p| (q, p))).collect(); + if obs.is_empty() { + continue; + } + // a random state prepared by rotations, then the gate + let mut c = RotCircuit::new(2); + for _ in 0..6 { + c.push(Op::Rot { axis: rng.string(2, 2), theta: rng.angle() }); + } + c.push(g.clone()); + let want = dense_expect(&c, &obs, &[]); + let got = expectation(&c, &[(1.0, obs.clone())], &[], &SpdConfig { threshold: 0.0, max_weight: None }).unwrap(); + assert!((got.value - want).abs() < 1e-12, "{g:?} {obs:?}: {} vs {want}", got.value); + } + } + } + } + + #[test] + fn random_circuits_match_the_dense_reference() { + let mut rng = Rng(2026); + for trial in 0..300 { + let n = 2 + rng.below(6) as u32; + let len = 6 + rng.below(30) as usize; + let c = random_circuit(&mut rng, n, len); + let obs = rng.string(n, 3); + let ones: Vec = (0..n).filter(|_| rng.below(2) == 1).collect(); + let want = dense_expect(&c, &obs, &ones); + let got = expectation(&c, &[(1.0, obs.clone())], &ones, &SpdConfig { threshold: 0.0, max_weight: None }).unwrap(); + assert!((got.value - want).abs() < 1e-10, "trial {trial}: {} vs {want}", got.value); + assert_eq!(got.truncation_bound, 0.0, "nothing is dropped at threshold 0 but exact zeros"); + } + } + + #[test] + fn the_truncation_bound_holds() { + // At every threshold the computed value is within the reported bound of + // the exact one — the guarantee, checked on circuits where the bound is + // not vacuous. + let mut rng = Rng(99); + let mut nontrivial = 0; + for _ in 0..120 { + let n = 6; + let c = random_circuit(&mut rng, n, 60); + let obs = vec![(rng.below(6) as u32, Pauli::Z)]; + let exact = dense_expect(&c, &obs, &[]); + for threshold in [1e-3, 1e-2, 5e-2, 0.2] { + let r = expectation(&c, &[(1.0, obs.clone())], &[], &SpdConfig { threshold, max_weight: None }).unwrap(); + assert!((r.value - exact).abs() <= r.truncation_bound + 1e-12, "|{} − {exact}| > {}", r.value, r.truncation_bound); + if r.truncation_bound > 0.0 && r.truncation_bound < 1.0 { + nontrivial += 1; + } + } + } + assert!(nontrivial > 50, "the bound was exercised {nontrivial} times"); + } + + #[test] + fn results_are_bit_identical_from_run_to_run() { + // A fresh std HashMap is seeded differently every time, so a term table + // that summed in its order would already differ here, in one process. + let c = kicked_ising(127, &heavy_hex_127(), 0.6, 5); + let cfg = SpdConfig { threshold: 1e-4, max_weight: None }; + let a = expect_z(&c, 62, &cfg).unwrap(); + assert!(a.peak_terms > 1000, "the run merges many terms: {}", a.peak_terms); + for _ in 0..4 { + let b = expect_z(&c, 62, &cfg).unwrap(); + assert_eq!(a.value.to_bits(), b.value.to_bits()); + assert_eq!(a.truncation_bound.to_bits(), b.truncation_bound.to_bits()); + assert_eq!((a.peak_terms, a.term_updates), (b.peak_terms, b.term_updates)); + } + } + + #[test] + fn word_width_does_not_change_the_answer() { + // The same circuit, placed in a 30-qubit and a 900-qubit register, is + // computed with one and with sixteen words per string: the bits agree. + let chain: Vec<(u32, u32)> = (0..29).map(|q| (q, q + 1)).collect(); + let small = kicked_ising(30, &chain, 0.5, 6); + let mut big = small.clone(); + big.n = 900; + let cfg = SpdConfig { threshold: 1e-6, max_weight: None }; + let a = expect_z(&small, 7, &cfg).unwrap(); + let b = expect_z(&big, 7, &cfg).unwrap(); + assert_eq!(a.value.to_bits(), b.value.to_bits()); + assert_eq!(a.term_updates, b.term_updates); + } + + #[test] + fn invalid_circuits_are_refused() { + let mut c = RotCircuit::new(3); + c.push(Op::CX(1, 1)); + assert_eq!(c.validate(), Err(SpdError::RepeatedQubit(1))); + let mut c = RotCircuit::new(3); + c.rx(3, 0.1); + assert_eq!(c.validate(), Err(SpdError::QubitOutOfRange { qubit: 3, n: 3 })); + let mut c = RotCircuit::new(3); + c.rx(0, f64::NAN); + assert!(matches!(c.validate(), Err(SpdError::BadAngle(_)))); + let mut c = RotCircuit::new(3); + c.push(Op::QuarterRot { axis: vec![(0, Pauli::Z)], k: 4 }); + assert_eq!(c.validate(), Err(SpdError::BadQuarter(4))); + assert_eq!(RotCircuit::new(1025).validate(), Err(SpdError::TooManyQubits(1025))); + } + + // ---- the heavy-hex experiment --------------------------------------------- + + #[test] + fn the_heavy_hex_map_has_the_published_shape() { + let e = heavy_hex_127(); + assert_eq!(e.len(), 144); + let mut deg = [0u32; 127]; + for &(a, b) in &e { + assert!(a < b && b < 127); + deg[a as usize] += 1; + deg[b as usize] += 1; + } + assert!(deg.iter().all(|&d| (1..=3).contains(&d))); + // Heavy-hex: every bridge qubit has degree 2, and no two degree-3 + // qubits are adjacent. + for &(a, b) in &e { + assert!(!(deg[a as usize] == 3 && deg[b as usize] == 3), "({a},{b})"); + } + for b in [14u32, 15, 16, 17, 33, 34, 35, 36, 52, 53, 54, 55, 71, 72, 73, 74, 90, 91, 92, 93, 109, 110, 111, 112] { + assert_eq!(deg[b as usize], 2, "bridge {b}"); + } + // connected + let mut seen = [false; 127]; + let mut stack = vec![0u32]; + while let Some(q) = stack.pop() { + if !std::mem::replace(&mut seen[q as usize], true) { + stack.extend(e.iter().filter(|&&(a, b)| a == q || b == q).map(|&(a, b)| a + b - q)); + } + } + assert!(seen.iter().all(|&s| s)); + // It is bipartite, so three matchings cover it, as the experiment's + // three RZZ layers do. + let mut colour = [u8::MAX; 127]; + colour[0] = 0; + let mut stack = vec![0u32]; + while let Some(q) = stack.pop() { + for &(a, b) in e.iter().filter(|&&(a, b)| a == q || b == q) { + let o = (a + b - q) as usize; + if colour[o] == u8::MAX { + colour[o] = 1 - colour[q as usize]; + stack.push(o as u32); + } + assert_ne!(colour[o], colour[q as usize]); + } + } + } + + #[test] + fn the_kicked_ising_clifford_points_are_exact() { + // θ_h = 0: no rotation, every ⟨Z_q⟩ = 1. θ_h = π/2 written as quarter + // turns: a Clifford circuit, so the Heisenberg sum stays one term and + // ⟨Z_q⟩ is exactly 0 or ±1 with nothing dropped. + let e = heavy_hex_127(); + let still = kicked_ising(127, &e, 0.0, 20); + let r = expect_z(&still, 62, &SpdConfig::default()).unwrap(); + assert_eq!((r.value, r.truncation_bound, r.peak_terms), (1.0, 0.0, 1)); + let mut cliff = RotCircuit::new(127); + for _ in 0..20 { + for q in 0..127 { + cliff.push(Op::QuarterRot { axis: vec![(q, Pauli::X)], k: 1 }); + } + for &(a, b) in &e { + cliff.push(Op::QuarterRot { axis: vec![(a, Pauli::Z), (b, Pauli::Z)], k: 3 }); + } + } + for q in [0, 13, 62, 126] { + let r = expect_z(&cliff, q, &SpdConfig::default()).unwrap(); + assert!(r.value == 0.0 || r.value.abs() == 1.0, "q {q}: {}", r.value); + assert_eq!((r.truncation_bound, r.peak_terms), (0.0, 1)); + } + } + + #[test] + fn kicked_ising_on_a_fragment_matches_the_dense_reference() { + // The same circuit family on a 15-qubit heavy-hex fragment, where a + // dense state fits: exact propagation agrees with it at every angle. + let e: Vec<(u32, u32)> = heavy_hex_127().into_iter().filter(|&(a, b)| a < 15 && b < 15).collect(); + assert!(e.len() >= 13); + for theta_h in [0.1, 0.4, 0.7854, 1.2] { + let c = kicked_ising(15, &e, theta_h, 4); + for q in [0, 4, 14] { + let want = dense_expect(&c, &[(q, Pauli::Z)], &[]); + let got = expect_z(&c, q, &SpdConfig { threshold: 0.0, max_weight: None }).unwrap(); + assert!((got.value - want).abs() < 1e-10, "θ_h {theta_h} q {q}: {} vs {want}", got.value); + } + } + } + + // ---- the dyadic adapter ------------------------------------------------------ + + #[cfg(feature = "quantum")] + #[test] + fn every_dyadic_algorithm_agrees_with_the_byte_exact_core() { + use crate::quantum::Circuit; + // Random dyadic circuits with every base gate and up to three controls: + // ⟨Z_q⟩ from the byte-exact simulator's probabilities and from the + // rotation form agree. + let mut rng = Rng(31); + for trial in 0..150 { + let n = 3 + rng.below(4) as u8; + let mut c = Circuit::new(n); + for _ in 0..25 { + let t = rng.below(u64::from(n)) as u8; + let mut ctrls = Vec::new(); + for _ in 0..rng.below(4) { + let q = rng.below(u64::from(n)) as u8; + if q != t && !ctrls.contains(&q) { + ctrls.push(q); + } + } + let (base, param) = match rng.below(9) { + 0 => (crate::quantum::BaseGate::H, 0), + 1 => (crate::quantum::BaseGate::X, 0), + 2 => (crate::quantum::BaseGate::Y, 0), + 3 => (crate::quantum::BaseGate::Z, 0), + 4 => (crate::quantum::BaseGate::S, 0), + 5 => (crate::quantum::BaseGate::Sdg, 0), + 6 => (crate::quantum::BaseGate::T, 0), + 7 => (crate::quantum::BaseGate::Tdg, 0), + _ => (crate::quantum::BaseGate::P, 1 + rng.below(6) as u16), + }; + if base == crate::quantum::BaseGate::H { + ctrls.clear(); + } + c.ops.push(crate::quantum::Gate { base, controls: ctrls, target: t, param }); + } + let sv = c.simulate().unwrap(); + let w = sv.prob_weights(); + let total: f64 = w.iter().map(|&x| x as f64).sum(); + let rc = from_dyadic(&c).unwrap(); + for q in 0..n { + let want: f64 = w.iter().enumerate().map(|(j, &x)| if (j >> q) & 1 == 0 { x as f64 } else { -(x as f64) }).sum::() / total; + let got = expect_z(&rc, u32::from(q), &SpdConfig { threshold: 0.0, max_weight: None }).unwrap(); + assert!((got.value - want).abs() < 1e-6, "trial {trial} q {q}: {} vs {want}", got.value); + } + } + } +} diff --git a/wai/wai-quantum/src/quantum_spd_receipt.rs b/wai/wai-quantum/src/quantum_spd_receipt.rs new file mode 100644 index 000000000..a8c1064c3 --- /dev/null +++ b/wai/wai-quantum/src/quantum_spd_receipt.rs @@ -0,0 +1,427 @@ +//! Signed, energy-accounted receipts for sparse Pauli dynamics +//! (`wai.quantum.spd`). +//! +//! A receipt binds three things: +//! +//! - **the job**: the circuit, the observable, the initial basis state and the +//! truncation settings. Angles and coefficients are hashed as their exact +//! IEEE-754 bits. +//! - **the result**: the expectation value, the truncation bound, the term +//! counts and the work done. +//! - **the energy** the sink measured while computing it (`joules_micro`, `0` +//! if unmetered). +//! +//! The job is run inside [`SpdReceipt::seal`], under a meter the caller +//! supplies, so the sealed result is never a caller's claim. Because the +//! computation is reproducible bit for bit on every platform, +//! [`SpdReceipt::verify_computation`] re-runs the job anywhere and demands the +//! identical result: value, bound, term counts and work. The energy stays what +//! it always is, attested by the signer's key. The receipt is a JWP profile like +//! the crate's others (Merkle root over leaf hashes, Ed25519), and can carry an +//! energy-class label (`quantum_energy::Labelled`). + +use crate::merkle::merkle_root; +use crate::quantum_spd::{Op, Pauli, PauliString, RotCircuit, SpdConfig, SpdError, SpdResult, expectation}; +use ed25519_dalek::{Signature, Signer, SigningKey, Verifier, VerifyingKey}; + +const DOMAIN_RECEIPT: &[u8] = b"wai:spd-receipt\x01"; +const DOMAIN_RECEIPT_ID: &[u8] = b"wai:spd-receipt-id\x01"; +const DOMAIN_CIRCUIT: &[u8] = b"wai:spd-circuit\x01"; +const DOMAIN_QUERY: &[u8] = b"wai:spd-query\x01"; +const DOMAIN_RESULT: &[u8] = b"wai:spd-result\x01"; + +/// Everything that determines a result: the circuit, the observable `Σ c·P`, +/// the initial basis state `|ones⟩`, and the truncation settings. +#[derive(Clone, Debug, PartialEq)] +pub struct SpdJob { + pub circuit: RotCircuit, + pub observable: Vec<(f64, PauliString)>, + pub ones: Vec, + pub config: SpdConfig, +} + +impl SpdJob { + pub fn run(&self) -> Result { + expectation(&self.circuit, &self.observable, &self.ones, &self.config) + } + + /// BLAKE3 over the circuit's canonical encoding. + pub fn circuit_hash(&self) -> [u8; 32] { + let mut o = DOMAIN_CIRCUIT.to_vec(); + o.extend_from_slice(&self.circuit.n.to_be_bytes()); + o.extend_from_slice(&(self.circuit.ops.len() as u64).to_be_bytes()); + for op in &self.circuit.ops { + encode_op(op, &mut o); + } + *blake3::hash(&o).as_bytes() + } + + /// BLAKE3 over the observable, the initial state and the settings. + pub fn query_hash(&self) -> [u8; 32] { + let mut o = DOMAIN_QUERY.to_vec(); + o.extend_from_slice(&(self.observable.len() as u64).to_be_bytes()); + for (c, s) in &self.observable { + o.extend_from_slice(&c.to_bits().to_be_bytes()); + encode_string(s, &mut o); + } + o.extend_from_slice(&(self.ones.len() as u64).to_be_bytes()); + for q in &self.ones { + o.extend_from_slice(&q.to_be_bytes()); + } + o.extend_from_slice(&self.config.threshold.to_bits().to_be_bytes()); + match self.config.max_weight { + Some(w) => { + o.push(1); + o.extend_from_slice(&w.to_be_bytes()); + } + None => o.push(0), + } + *blake3::hash(&o).as_bytes() + } +} + +fn encode_string(s: &[(u32, Pauli)], o: &mut Vec) { + o.extend_from_slice(&(s.len() as u64).to_be_bytes()); + for (q, p) in s { + o.extend_from_slice(&q.to_be_bytes()); + o.push(match p { + Pauli::X => b'X', + Pauli::Y => b'Y', + Pauli::Z => b'Z', + }); + } +} + +fn encode_op(op: &Op, o: &mut Vec) { + let one = |tag: u8, q: u32, o: &mut Vec| { + o.push(tag); + o.extend_from_slice(&q.to_be_bytes()); + }; + let two = |tag: u8, a: u32, b: u32, o: &mut Vec| { + o.push(tag); + o.extend_from_slice(&a.to_be_bytes()); + o.extend_from_slice(&b.to_be_bytes()); + }; + match op { + Op::H(q) => one(1, *q, o), + Op::S(q) => one(2, *q, o), + Op::Sdg(q) => one(3, *q, o), + Op::SX(q) => one(4, *q, o), + Op::SXdg(q) => one(5, *q, o), + Op::X(q) => one(6, *q, o), + Op::Y(q) => one(7, *q, o), + Op::Z(q) => one(8, *q, o), + Op::CX(a, b) => two(9, *a, *b, o), + Op::CZ(a, b) => two(10, *a, *b, o), + Op::Swap(a, b) => two(11, *a, *b, o), + Op::Rot { axis, theta } => { + o.push(12); + encode_string(axis, o); + o.extend_from_slice(&theta.to_bits().to_be_bytes()); + } + Op::QuarterRot { axis, k } => { + o.push(13); + encode_string(axis, o); + o.push(*k); + } + } +} + +fn result_hash(r: &SpdResult) -> [u8; 32] { + let mut o = DOMAIN_RESULT.to_vec(); + o.extend_from_slice(&r.value.to_bits().to_be_bytes()); + o.extend_from_slice(&r.truncation_bound.to_bits().to_be_bytes()); + o.extend_from_slice(&(r.peak_terms as u64).to_be_bytes()); + o.extend_from_slice(&(r.final_terms as u64).to_be_bytes()); + o.extend_from_slice(&r.term_updates.to_be_bytes()); + *blake3::hash(&o).as_bytes() +} + +fn hx(b: &[u8]) -> String { + b.iter().map(|x| format!("{x:02x}")).collect() +} + +fn unhex(s: &str) -> Option<[u8; N]> { + if s.len() != 2 * N || !s.bytes().all(|b| b.is_ascii_hexdigit()) { + return None; + } + let mut out = [0u8; N]; + for (i, c) in s.as_bytes().chunks(2).enumerate() { + out[i] = u8::from_str_radix(std::str::from_utf8(c).ok()?, 16).ok()?; + } + Some(out) +} + +/// A signed record of one sparse-Pauli-dynamics computation and its energy. +#[derive(Clone, Debug, PartialEq)] +pub struct SpdReceipt { + pub circuit_hash: [u8; 32], + pub query_hash: [u8; 32], + pub result: SpdResult, + pub result_hash: [u8; 32], + /// Merkle root over (circuit, query, result) hashes. + pub merkle_root: [u8; 32], + /// Measured microjoules the sink spent on the run (`0` = unmetered). + /// Attested by the signature, not verifiable by anyone else. + pub joules_micro: u64, + pub parent_receipt_hash: Option<[u8; 32]>, + pub signer_pubkey: [u8; 32], + pub signer_id: String, + pub sig: [u8; 64], +} + +impl SpdReceipt { + fn leaves(&self) -> Vec<[u8; 32]> { + vec![self.circuit_hash, self.query_hash, self.result_hash] + } + + fn payload_with(&self, label: Option<&crate::quantum_energy::EnergyClass>) -> Vec { + let mut o = crate::quantum_energy::labelled_domain(DOMAIN_RECEIPT, label); + o.extend_from_slice(&self.circuit_hash); + o.extend_from_slice(&self.query_hash); + o.extend_from_slice(&self.result_hash); + o.extend_from_slice(&self.merkle_root); + o.extend_from_slice(&self.result.term_updates.to_be_bytes()); + o.extend_from_slice(&self.joules_micro.to_be_bytes()); + crate::quantum_energy::write_label(label, &mut o); + match self.parent_receipt_hash { + Some(h) => { + o.push(1); + o.extend_from_slice(&h); + } + None => o.push(0), + } + o.extend_from_slice(&self.signer_pubkey); + o.extend_from_slice(self.signer_id.as_bytes()); + o + } + + /// Run `job` under `meter` and seal what it computed. `meter` is handed a + /// closure that performs the run; it must call it once, and returns the + /// microjoules it measured (`|run| { run(); 0 }` for an unmetered seal). + /// The result sealed is the one that closure produced. + pub fn seal( + signer: &SigningKey, + signer_id: impl Into, + job: &SpdJob, + parent_receipt_hash: Option<[u8; 32]>, + meter: impl FnOnce(&mut dyn FnMut()) -> u64, + ) -> Result { + let mut out = None; + let joules_micro = meter(&mut || out = Some(job.run())); + let result = out.ok_or(SpdError::NotRun)??; + let mut r = SpdReceipt { + circuit_hash: job.circuit_hash(), + query_hash: job.query_hash(), + result, + result_hash: result_hash(&result), + merkle_root: [0; 32], + joules_micro, + parent_receipt_hash, + signer_pubkey: signer.verifying_key().to_bytes(), + signer_id: signer_id.into(), + sig: [0; 64], + }; + r.merkle_root = merkle_root(&r.leaves()); + r.sig = signer.sign(&r.payload_with(None)).to_bytes(); + Ok(r) + } + + fn body_verifies(&self) -> bool { + result_hash(&self.result) == self.result_hash && merkle_root(&self.leaves()) == self.merkle_root + } + + /// The receipt is internally consistent and its signature is valid. Does + /// not re-run the job: see [`Self::verify_computation`]. + pub fn verify(&self) -> bool { + if !self.body_verifies() { + return false; + } + let Ok(k) = VerifyingKey::from_bytes(&self.signer_pubkey) else { + return false; + }; + k.verify(&self.payload_with(None), &Signature::from_bytes(&self.sig)).is_ok() + } + + /// The re-run half of [`Self::verify_computation`]. + pub(crate) fn computation_matches(&self, job: &SpdJob) -> bool { + job.circuit_hash() == self.circuit_hash + && job.query_hash() == self.query_hash + && job.run().is_ok_and(|r| result_hash(&r) == self.result_hash) + } + + /// The strong check: `job` is the job this receipt names, re-running it + /// here reproduces the sealed result bit for bit, and the receipt verifies. + pub fn verify_computation(&self, job: &SpdJob) -> bool { + self.computation_matches(job) && self.verify() + } + + /// Stable id, for chaining a derived receipt via `parent_receipt_hash`. + pub fn receipt_hash(&self) -> [u8; 32] { + let mut h = blake3::Hasher::new(); + h.update(DOMAIN_RECEIPT_ID); + h.update(&self.payload_with(None)); + h.update(&self.sig); + *h.finalize().as_bytes() + } + + /// Canonical JSON, keys sorted. Floating-point figures appear twice: as a + /// readable number, and as their exact bits (`*_bits`), which are what a + /// verifier reads back. + pub fn to_json(&self) -> String { + let parent = match self.parent_receipt_hash { + Some(h) => format!("\"{}\"", hx(&h)), + None => "null".into(), + }; + let r = &self.result; + format!( + "{{\"kind\":\"quantum-spd\",\"circuit_hash\":\"{}\",\"final_terms\":{},\"joules_micro\":{},\ + \"parent_receipt_hash\":{},\"peak_terms\":{},\"query_hash\":\"{}\",\"receipt_hash\":\"{}\",\ + \"result_hash\":\"{}\",\"root_hash\":\"{}\",\"sig\":\"{}\",\"signer_id\":{},\ + \"signer_pubkey\":\"{}\",\"term_updates\":{},\"truncation_bound\":{},\ + \"truncation_bound_bits\":\"{:016x}\",\"value\":{},\"value_bits\":\"{:016x}\"}}", + hx(&self.circuit_hash), + r.final_terms, + self.joules_micro, + parent, + r.peak_terms, + hx(&self.query_hash), + hx(&self.receipt_hash()), + hx(&self.result_hash), + hx(&self.merkle_root), + hx(&self.sig), + serde_json::to_string(&self.signer_id).unwrap(), + hx(&self.signer_pubkey), + r.term_updates, + serde_json::to_string(&r.truncation_bound).unwrap(), + r.truncation_bound.to_bits(), + serde_json::to_string(&r.value).unwrap(), + r.value.to_bits(), + ) + } + + pub fn from_json(s: &str) -> Option { + let v: serde_json::Value = serde_json::from_str(s).ok()?; + let o = v.as_object()?; + if o.get("kind")?.as_str()? != "quantum-spd" { + return None; + } + let u = |k: &str| o.get(k).and_then(|x| x.as_u64()); + let s = |k: &str| o.get(k).and_then(|x| x.as_str()); + let bits = |k: &str| u64::from_str_radix(s(k)?, 16).ok().map(f64::from_bits); + let parent = match o.get("parent_receipt_hash") { + Some(serde_json::Value::String(h)) => Some(unhex::<32>(h)?), + _ => None, + }; + Some(SpdReceipt { + circuit_hash: unhex(s("circuit_hash")?)?, + query_hash: unhex(s("query_hash")?)?, + result: SpdResult { + value: bits("value_bits")?, + truncation_bound: bits("truncation_bound_bits")?, + peak_terms: usize::try_from(u("peak_terms")?).ok()?, + final_terms: usize::try_from(u("final_terms")?).ok()?, + term_updates: u("term_updates")?, + }, + result_hash: unhex(s("result_hash")?)?, + merkle_root: unhex(s("root_hash")?)?, + joules_micro: u("joules_micro")?, + parent_receipt_hash: parent, + signer_pubkey: unhex(s("signer_pubkey")?)?, + signer_id: s("signer_id")?.to_string(), + sig: unhex(s("sig")?)?, + }) + } +} + +crate::quantum_energy::labellable!(SpdReceipt, DOMAIN_RECEIPT_ID); + +#[cfg(test)] +mod tests { + use super::*; + use crate::quantum_spd::{heavy_hex_127, kicked_ising}; + + fn key(b: u8) -> SigningKey { + SigningKey::from_bytes(&[b; 32]) + } + + fn job() -> SpdJob { + SpdJob { + circuit: kicked_ising(127, &heavy_hex_127(), 0.4, 4), + observable: vec![(1.0, vec![(62, Pauli::Z)])], + ones: vec![], + config: SpdConfig { threshold: 1e-4, max_weight: None }, + } + } + + fn unmetered(run: &mut dyn FnMut()) -> u64 { + run(); + 0 + } + + #[test] + fn a_sealed_receipt_verifies_and_reproduces() { + let j = job(); + let r = SpdReceipt::seal(&key(1), "sink-a", &j, None, |run| { + run(); + 4_200 + }) + .unwrap(); + assert!(r.verify()); + assert!(r.verify_computation(&j)); + assert_eq!(r.joules_micro, 4_200); + assert_eq!(r.result, j.run().unwrap()); + } + + #[test] + fn the_json_round_trips_to_the_same_bytes() { + let r = SpdReceipt::seal(&key(2), "sink \"b\"", &job(), Some([7; 32]), unmetered).unwrap(); + let back = SpdReceipt::from_json(&r.to_json()).unwrap(); + assert_eq!(back, r); + assert_eq!(back.to_json(), r.to_json()); + assert!(back.verify()); + } + + #[test] + fn tampering_is_caught() { + let j = job(); + let r = SpdReceipt::seal(&key(3), "sink", &j, None, unmetered).unwrap(); + let mut a = r.clone(); + a.joules_micro = 1; + assert!(!a.verify(), "the energy is signed"); + let mut b = r.clone(); + b.result.value += 1e-12; + assert!(!b.verify(), "the result is hashed"); + let mut c = r.clone(); + c.result.term_updates -= 1; + assert!(!c.verify(), "the work is hashed"); + // A different job — one angle moved by one ulp — is not the job sealed. + let mut other = j.clone(); + if let Op::Rot { theta, .. } = &mut other.circuit.ops[5] { + *theta = f64::from_bits(theta.to_bits() + 1); + } + assert!(!r.verify_computation(&other)); + let mut looser = j.clone(); + looser.config.threshold = 2e-4; + assert!(!r.verify_computation(&looser), "the settings are part of the job"); + } + + #[test] + fn a_meter_that_never_runs_the_job_cannot_seal() { + assert_eq!(SpdReceipt::seal(&key(4), "sink", &job(), None, |_run| 99), Err(SpdError::NotRun)); + } + + #[test] + fn an_energy_class_label_seals_over_it() { + use crate::quantum_energy::{EnergyClass, Labelled}; + let r = SpdReceipt::seal(&key(5), "sink", &job(), None, |run| { + run(); + 1_000 + }) + .unwrap(); + let class = EnergyClass::ModelBased; + let labelled = Labelled::seal(r, class, &key(5)).unwrap(); + assert!(labelled.verify()); + assert!(!labelled.receipt.verify(), "the inner signature now covers the label"); + } +} diff --git a/wai/wai-quantum/src/quantum_vml.rs b/wai/wai-quantum/src/quantum_vml.rs index dc77659d0..a6e8ab9ae 100644 --- a/wai/wai-quantum/src/quantum_vml.rs +++ b/wai/wai-quantum/src/quantum_vml.rs @@ -65,7 +65,7 @@ fn rand_angle(s: &mut u64) -> f64 { #[inline] fn apply_ry(amps: &mut [C], n: u8, q: u8, theta: f64) { - let (c, s) = ((theta * 0.5).cos(), (theta * 0.5).sin()); + let (s, c) = crate::repro::sin_cos(theta * 0.5); let bit = 1usize << q; let dim = 1usize << n; let mut i = 0; @@ -81,7 +81,7 @@ fn apply_ry(amps: &mut [C], n: u8, q: u8, theta: f64) { } #[inline] fn apply_rz(amps: &mut [C], q: u8, theta: f64) { - let (cm, sm) = ((theta * 0.5).cos(), (theta * 0.5).sin()); + let (sm, cm) = crate::repro::sin_cos(theta * 0.5); let e0 = C { re: cm, im: -sm }; // e^{-iθ/2} let e1 = C { re: cm, im: sm }; // e^{+iθ/2} let bit = 1usize << q; @@ -148,7 +148,7 @@ fn kl(target: &[f64], model: &[f64]) -> f64 { let mut s = 0.0; for (&q, &p) in target.iter().zip(model) { if q > 0.0 { - s += q * (q / p.max(1e-12)).ln(); + s += q * crate::repro::ln(q / p.max(1e-12)); } } s diff --git a/wai/wai-quantum/src/repro.rs b/wai/wai-quantum/src/repro.rs new file mode 100644 index 000000000..5b50dc8c1 --- /dev/null +++ b/wai/wai-quantum/src/repro.rs @@ -0,0 +1,387 @@ +//! Transcendental functions that give the same bits on every platform. +//! +//! The platform's libm may differ from another's in the last place. That is +//! enough to break a result promised bit-identical natively, in the browser and +//! at the edge, and it did: four layers of this crate computed different bits +//! natively and as a WASI component until they moved here. Everything below is +//! fdlibm's algorithm, as musl carries it, using only `+ − × ÷`, which Rust +//! never fuses. + +/// `sin` and `cos` of `x`, bit-identical on every IEEE-754 platform. The +/// algorithm is fdlibm's (as musl carries it): Cody–Waite reduction by π/2 in +/// three parts, then minimax polynomials on [−π/4, π/4]. Only `+ − × ÷` are +/// used, and Rust never fuses them, so the bits do not depend on the target. +/// Accurate to within an ulp of the true value. +/// +/// For `|x| > 2²⁰` (and for non-finite `x`) both are NaN: the reduction used +/// here is exact only below that, and a NaN is better than a value that +/// depends on the platform. +pub fn sin_cos(x: f64) -> (f64, f64) { + if !(x.abs() <= 1_048_576.0) { + return (f64::NAN, f64::NAN); + } + let ix = (x.to_bits() >> 32) as u32 & 0x7fff_ffff; + if ix <= 0x3fe9_21fb { + // |x| ≤ π/4 + if ix < 0x3e40_0000 { + // |x| < 2⁻²⁷: sin x = x and cos x = 1 to the last bit. + return (x, 1.0); + } + return (k_sin(x, 0.0, false), k_cos(x, 0.0)); + } + let (n, y0, y1) = rem_pio2(x, ix); + let (s, c) = (k_sin(y0, y1, true), k_cos(y0, y1)); + match n & 3 { + 0 => (s, c), + 1 => (c, -s), + 2 => (-s, -c), + _ => (-c, s), + } +} + +const S1: f64 = -1.666_666_666_666_663_243_48e-1; +const S2: f64 = 8.333_333_333_322_489_461_24e-3; +const S3: f64 = -1.984_126_982_985_794_931_34e-4; +const S4: f64 = 2.755_731_370_707_006_767_89e-6; +const S5: f64 = -2.505_076_025_340_686_341_95e-8; +const S6: f64 = 1.589_690_995_211_550_102_21e-10; + +fn k_sin(x: f64, y: f64, tail: bool) -> f64 { + let z = x * x; + let w = z * z; + let r = S2 + z * (S3 + z * S4) + z * w * (S5 + z * S6); + let v = z * x; + if !tail { + x + v * (S1 + z * r) + } else { + x - ((z * (0.5 * y - v * r) - y) - v * S1) + } +} + +const C1: f64 = 4.166_666_666_666_660_190_37e-2; +const C2: f64 = -1.388_888_888_887_410_957_49e-3; +const C3: f64 = 2.480_158_728_947_672_941_78e-5; +const C4: f64 = -2.755_731_435_139_066_330_35e-7; +const C5: f64 = 2.087_572_321_298_174_827_90e-9; +const C6: f64 = -1.135_964_755_778_819_482_65e-11; + +fn k_cos(x: f64, y: f64) -> f64 { + let z = x * x; + let w = z * z; + let r = z * (C1 + z * (C2 + z * C3)) + w * w * (C4 + z * (C5 + z * C6)); + let hz = 0.5 * z; + let w = 1.0 - hz; + w + (((1.0 - w) - hz) + (z * r - x * y)) +} + +const TOINT: f64 = 1.5 / f64::EPSILON; +const PIO4: f64 = 7.853_981_633_974_482_789_99e-1; +const INVPIO2: f64 = 6.366_197_723_675_813_824_33e-1; +const PIO2_1: f64 = 1.570_796_326_734_125_614_17; +const PIO2_1T: f64 = 6.077_100_506_506_192_249_32e-11; +const PIO2_2: f64 = 6.077_100_506_303_965_976_60e-11; +const PIO2_2T: f64 = 2.022_266_248_795_950_631_54e-21; +const PIO2_3: f64 = 2.022_266_248_711_166_455_80e-21; +const PIO2_3T: f64 = 8.478_427_660_368_899_569_97e-32; + +/// `x = n·π/2 + (y0 + y1)` with `|y0 + y1| ≤ π/4`, for `|x| ≤ 2²⁰`. +fn rem_pio2(x: f64, ix: u32) -> (i64, f64, f64) { + let mut f = (x * INVPIO2 + TOINT) - TOINT; + let mut r = x - f * PIO2_1; + let mut w = f * PIO2_1T; + // Matters only under directed rounding; kept so the port is exact. + if r - w < -PIO4 { + f -= 1.0; + r = x - f * PIO2_1; + w = f * PIO2_1T; + } else if r - w > PIO4 { + f += 1.0; + r = x - f * PIO2_1; + w = f * PIO2_1T; + } + let mut y0 = r - w; + let ex = (ix >> 20) as i32; + let ey = ((y0.to_bits() >> 52) & 0x7ff) as i32; + if ex - ey > 16 { + let t = r; + w = f * PIO2_2; + r = t - w; + w = f * PIO2_2T - ((t - r) - w); + y0 = r - w; + let ey = ((y0.to_bits() >> 52) & 0x7ff) as i32; + if ex - ey > 49 { + let t = r; + w = f * PIO2_3; + r = t - w; + w = f * PIO2_3T - ((t - r) - w); + y0 = r - w; + } + } + let y1 = (r - y0) - w; + (f as i64, y0, y1) +} + +// --------------------------------------------------------------------------- +// Natural logarithm +// --------------------------------------------------------------------------- + +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))] +const LN2_HI: f64 = 6.931_471_803_691_238_164_90e-1; +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))] +const LN2_LO: f64 = 1.908_214_929_270_587_700_02e-10; +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))] +const LG1: f64 = 6.666_666_666_666_735_130e-1; +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))] +const LG2: f64 = 3.999_999_999_940_941_908e-1; +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))] +const LG3: f64 = 2.857_142_874_366_239_149e-1; +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))] +const LG4: f64 = 2.222_219_843_214_978_396e-1; +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))] +const LG5: f64 = 1.818_357_216_161_805_012e-1; +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))] +const LG6: f64 = 1.531_383_769_920_937_332e-1; +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))] +const LG7: f64 = 1.479_819_860_511_658_591e-1; + +/// The natural logarithm, within an ulp of the true value: `−∞` at zero, NaN +/// below it. +#[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))] +pub fn ln(x: f64) -> f64 { + let mut bits = x.to_bits(); + let mut hx = (bits >> 32) as u32; + let mut k: i32 = 0; + let mut x = x; + if hx < 0x0010_0000 || hx >> 31 != 0 { + if bits << 1 == 0 { + return f64::NEG_INFINITY; + } + if hx >> 31 != 0 { + return f64::NAN; + } + // subnormal: scale up + k -= 54; + x *= 18_014_398_509_481_984.0; // 2⁵⁴ + bits = x.to_bits(); + hx = (bits >> 32) as u32; + } else if hx >= 0x7ff0_0000 { + return x; + } else if hx == 0x3ff0_0000 && bits << 32 == 0 { + return 0.0; + } + // reduce x into [√2/2, √2] + hx = hx.wrapping_add(0x3ff0_0000 - 0x3fe6_a09e); + k += (hx >> 20) as i32 - 0x3ff; + hx = (hx & 0x000f_ffff) + 0x3fe6_a09e; + let x = f64::from_bits((u64::from(hx) << 32) | (bits & 0xffff_ffff)); + let f = x - 1.0; + let hfsq = 0.5 * f * f; + let s = f / (2.0 + f); + let z = s * s; + let w = z * z; + let t1 = w * (LG2 + w * (LG4 + w * LG6)); + let t2 = z * (LG1 + w * (LG3 + w * (LG5 + w * LG7))); + let r = t2 + t1; + let dk = f64::from(k); + s * (hfsq + r) + dk * LN2_LO - hfsq + f + dk * LN2_HI +} + +// --------------------------------------------------------------------------- +// Arctangent +// --------------------------------------------------------------------------- + +#[cfg(feature = "quantum_qfhrr")] +const ATANHI: [f64; 4] = [ + 4.636_476_090_008_060_935_15e-1, + 7.853_981_633_974_482_789_99e-1, + 9.827_937_232_473_290_540_82e-1, + 1.570_796_326_794_896_558_00, +]; +#[cfg(feature = "quantum_qfhrr")] +const ATANLO: [f64; 4] = [ + 2.269_877_745_296_168_709_24e-17, + 3.061_616_997_868_383_017_93e-17, + 1.390_331_103_123_099_845_16e-17, + 6.123_233_995_736_766_035_87e-17, +]; +#[cfg(feature = "quantum_qfhrr")] +const AT: [f64; 11] = [ + 3.333_333_333_333_293_180_27e-1, + -1.999_999_999_987_648_324_76e-1, + 1.428_571_427_250_346_637_11e-1, + -1.111_111_040_546_235_578_80e-1, + 9.090_887_133_436_506_561_96e-2, + -7.691_876_205_044_829_994_95e-2, + 6.661_073_137_387_531_206_69e-2, + -5.833_570_133_790_573_486_45e-2, + 4.976_877_994_615_932_360_17e-2, + -3.653_157_274_421_691_552_70e-2, + 1.628_582_011_536_578_236_23e-2, +]; + +/// The arctangent, within an ulp of the true value. +#[cfg(feature = "quantum_qfhrr")] +pub fn atan(x: f64) -> f64 { + let ix = (x.to_bits() >> 32) as u32; + let negative = ix >> 31 != 0; + let ix = ix & 0x7fff_ffff; + if ix >= 0x4410_0000 { + // |x| ≥ 2⁶⁶ + if x.is_nan() { + return x; + } + let z = ATANHI[3]; + return if negative { -z } else { z }; + } + let mut x = x; + let id: i32; + if ix < 0x3fdc_0000 { + // |x| < 0.4375 + if ix < 0x3e40_0000 { + return x; + } + id = -1; + } else { + x = x.abs(); + if ix < 0x3ff3_0000 { + if ix < 0x3fe6_0000 { + id = 0; + x = (2.0 * x - 1.0) / (2.0 + x); + } else { + id = 1; + x = (x - 1.0) / (x + 1.0); + } + } else if ix < 0x4003_8000 { + id = 2; + x = (x - 1.5) / (1.0 + 1.5 * x); + } else { + id = 3; + x = -1.0 / x; + } + } + let z = x * x; + let w = z * z; + let s1 = z * (AT[0] + w * (AT[2] + w * (AT[4] + w * (AT[6] + w * (AT[8] + w * AT[10]))))); + let s2 = w * (AT[1] + w * (AT[3] + w * (AT[5] + w * (AT[7] + w * AT[9])))); + if id < 0 { + return x - x * (s1 + s2); + } + let i = id as usize; + let z = ATANHI[i] - (x * (s1 + s2) - ATANLO[i] - x); + if negative { -z } else { z } +} + +#[cfg(feature = "quantum_qfhrr")] +const PI: f64 = 3.141_592_653_589_793_116; +#[cfg(feature = "quantum_qfhrr")] +const PI_LO: f64 = 1.224_646_799_147_353_177_2e-16; + +/// The angle of the point `(x, y)`, in `(−π, π]`, within an ulp of the true +/// value. +#[cfg(feature = "quantum_qfhrr")] +pub fn atan2(y: f64, x: f64) -> f64 { + if x.is_nan() || y.is_nan() { + return x + y; + } + let (xb, yb) = (x.to_bits(), y.to_bits()); + let (ix, lx) = ((xb >> 32) as u32, xb as u32); + let (iy, ly) = ((yb >> 32) as u32, yb as u32); + if (ix.wrapping_sub(0x3ff0_0000) | lx) == 0 { + return atan(y); // x = 1 + } + let m = ((iy >> 31) & 1) | ((ix >> 30) & 2); + let (ix, iy) = (ix & 0x7fff_ffff, iy & 0x7fff_ffff); + if (iy | ly) == 0 { + return match m { + 0 | 1 => y, + 2 => PI, + _ => -PI, + }; + } + if (ix | lx) == 0 { + return if m & 1 != 0 { -PI / 2.0 } else { PI / 2.0 }; + } + if ix == 0x7ff0_0000 { + return if iy == 0x7ff0_0000 { + match m { + 0 => PI / 4.0, + 1 => -PI / 4.0, + 2 => 3.0 * PI / 4.0, + _ => -3.0 * PI / 4.0, + } + } else { + match m { + 0 => 0.0, + 1 => -0.0, + 2 => PI, + _ => -PI, + } + }; + } + if ix.wrapping_add(64 << 20) < iy || iy == 0x7ff0_0000 { + return if m & 1 != 0 { -PI / 2.0 } else { PI / 2.0 }; + } + let z = if m & 2 != 0 && iy.wrapping_add(64 << 20) < ix { 0.0 } else { atan((y / x).abs()) }; + match m { + 0 => z, + 1 => -z, + 2 => PI - (z - PI_LO), + _ => (z - PI_LO) - PI, + } +} + +#[cfg(test)] +mod tests { + #[allow(unused_imports)] + use super::*; + + #[allow(dead_code)] + fn ulps(a: f64, b: f64) -> u64 { + if a == b { + return 0; + } + (a.to_bits() as i64).abs_diff(b.to_bits() as i64) + } + + /// A deterministic spread of inputs over many magnitudes. + #[allow(dead_code)] + fn inputs() -> impl Iterator { + let mut s = 0x2545_f491_4f6c_dd1du64; + (0..200_000).map(move |i| { + s ^= s << 13; + s ^= s >> 7; + s ^= s << 17; + let u = (s >> 11) as f64 / (1u64 << 53) as f64; + let scale = [1e-300, 1e-8, 1e-3, 0.5, 1.0, 3.0, 1e3, 1e12, 1e300][i % 9]; + u * scale + }) + } + + #[cfg(any(feature = "quantum_vml", feature = "quantum_kernel"))] + #[test] + fn ln_is_within_an_ulp_of_the_platform() { + for x in inputs().filter(|&x| x > 0.0) { + assert!(ulps(ln(x), x.ln()) <= 1, "ln({x:e}) = {:e} vs {:e}", ln(x), x.ln()); + } + assert_eq!(ln(1.0), 0.0); + assert_eq!(ln(0.0), f64::NEG_INFINITY); + assert!(ln(-1.0).is_nan()); + assert_eq!(ln(f64::INFINITY), f64::INFINITY); + assert!(ulps(ln(f64::MIN_POSITIVE / 8.0), (f64::MIN_POSITIVE / 8.0).ln()) <= 1, "subnormal"); + } + + #[cfg(feature = "quantum_qfhrr")] + #[test] + fn atan2_is_within_an_ulp_of_the_platform() { + let xs: Vec = inputs().take(4000).collect(); + for (i, &a) in xs.iter().enumerate() { + let b = xs[(i * 7 + 3) % xs.len()]; + for (y, x) in [(a, b), (-a, b), (a, -b), (-a, -b)] { + assert!(ulps(atan2(y, x), y.atan2(x)) <= 1, "atan2({y:e}, {x:e}) = {:e} vs {:e}", atan2(y, x), y.atan2(x)); + } + } + assert_eq!(atan2(0.0, -1.0), PI); + assert_eq!(atan2(1.0, 0.0), PI / 2.0); + assert_eq!(atan(f64::INFINITY), ATANHI[3]); + } +} diff --git a/wai/wai-quantum/tests/bit_identity.rs b/wai/wai-quantum/tests/bit_identity.rs new file mode 100644 index 000000000..b93dc549d --- /dev/null +++ b/wai/wai-quantum/tests/bit_identity.rs @@ -0,0 +1,109 @@ +//! The crate promises identical results natively, in the browser and at the +//! edge. Before every transcendental moved in-crate, four of these layers gave +//! different bits natively and under `wasm32-wasip2`. These digests were taken +//! natively and confirmed under `wasm32-wasip2`; any platform whose math +//! differs fails here, not in a user's diff. +//! +//! `CARGO_TARGET_WASM32_WASIP2_RUNNER=wasmtime cargo test --release --target wasm32-wasip2 --features full` +#![cfg(all( + feature = "quantum", + feature = "quantum_vml", + feature = "quantum_phasor", + feature = "quantum_qfhrr", + feature = "quantum_kernel", + feature = "quantum_mps", + feature = "quantum_pauli", + feature = "quantum_spd" +))] + +use wai_quantum::*; + +/// FNV-1a over every value's IEEE-754 bits. +fn digest(xs: &[f64]) -> u64 { + let mut h: u64 = 0xcbf2_9ce4_8422_2325; + for x in xs { + for b in x.to_bits().to_le_bytes() { + h = (h ^ u64::from(b)).wrapping_mul(0x100_0000_01b3); + } + } + h +} + +fn trained() -> quantum_vml::BornFit { + quantum_vml::train_born(&quantum_vml::parity_target(3), 3, 2, 40, 0.1, 7) +} + +#[test] +fn born_machine_training() { + let fit = trained(); + let mut v = fit.params.clone(); + v.extend(&fit.loss_curve); + v.extend(&fit.final_probs); + v.push(fit.final_kl); + assert_eq!(digest(&v), 0xd8e5_bd02_2201_2e5d); +} + +#[test] +fn phasor_bridge() { + let p = quantum_phasor::random_params(3, 5); + let (a, b, c, d) = quantum_phasor::compare(3, &p, 64); + let mut v = a; + v.extend(b); + v.extend(c); + v.push(d); + assert_eq!(digest(&v), 0x08be_30dd_7ce1_9561); +} + +#[test] +fn qfhrr_born_machine() { + let fit = trained(); + let mut v = quantum_qfhrr::born_probs(3, 2, &fit.params, 4); + v.push(quantum_qfhrr::tv_vs_exact(3, 2, &fit.params, 4)); + assert_eq!(digest(&v), 0x5407_7c4d_1be0_d985); +} + +#[test] +fn phasor_kernel() { + let (x, y) = quantum_kernel::circles(200, 3); + let rff = quantum_kernel::sample_rff(2, 64, 1.0, 9); + let w = quantum_kernel::fit_ridge(&rff, &x, &y, 1e-3); + let mut v: Vec = x.iter().flatten().copied().collect(); + v.extend(&w); + assert_eq!(digest(&v), 0x58ba_1666_5ef6_fd9e); +} + +fn brickwork() -> quantum::Circuit { + let mut c = quantum::Circuit::new(8); + for l in 0..4u8 { + for q in 0..8 { + c.h(q); + c.t(q); + } + for q in (l % 2..7).step_by(2) { + c.cx(q, q + 1); + } + } + c +} + +#[test] +fn mps_truncation() { + let c = brickwork(); + let (mps, _) = quantum_mps::run(&c, 4).unwrap(); + assert_eq!(digest(&[quantum_mps::fidelity(&mps, &c.simulate().unwrap())]), 0xa807_44ca_a97f_bff8); +} + +#[test] +fn sparse_pauli_propagation() { + let (e, _) = quantum_pauli::expect_z(&brickwork(), 3, 1e-6).unwrap(); + assert_eq!(digest(&[e]), 0x8cc6_5829_1fb1_8909); +} + +#[test] +fn sparse_pauli_dynamics() { + let chain: Vec<(u32, u32)> = (0..29).map(|q| (q, q + 1)).collect(); + let circuit = quantum_spd::kicked_ising(30, &chain, 0.5, 6); + let cfg = quantum_spd::SpdConfig { threshold: 1e-6, max_weight: None }; + let r = quantum_spd::expect_z(&circuit, 7, &cfg).unwrap(); + assert_eq!(digest(&[r.value, r.truncation_bound]), 0xfc44_f4b6_0e70_68b2); +} diff --git a/wai/wai-quantum/tests/utility_reference.rs b/wai/wai-quantum/tests/utility_reference.rs new file mode 100644 index 000000000..d5c5d87a2 --- /dev/null +++ b/wai/wai-quantum/tests/utility_reference.rs @@ -0,0 +1,37 @@ +//! The 127-qubit heavy-hex kicked-Ising experiment against published values. +//! +//! - The 5-step magnetization `Mz = (1/127) Σ ⟨Z_q⟩` at `θ_h = kπ/32` was +//! computed exactly and published as `exact.csv` with arXiv:2308.05077. +//! With the truncation threshold at 0 this crate reproduces `k = 1` to +//! 3.5e-7 (0.995131 against 0.995131, 357M terms); here a threshold of +//! 1e-7 keeps the test fast, and the error is checked against the bound. +//! - ⟨Z₆₂⟩ after 20 steps at `θ_h = 2π/32`: the converged classical values +//! published with the same paper (sparse Pauli dynamics, a mixed method, +//! PEPS, PEPO) span 0.9776–0.9796. +#![cfg(feature = "quantum_spd")] + +use wai_quantum::quantum_spd::{Pauli, SpdConfig, expect_z, expectation, heavy_hex_127, kicked_ising}; + +fn theta(k: u32) -> f64 { + f64::from(k) * std::f64::consts::PI / 32.0 +} + +#[test] +fn five_step_magnetization_matches_the_exact_values() { + let mz: Vec<(f64, Vec<(u32, Pauli)>)> = (0..127).map(|q| (1.0 / 127.0, vec![(q, Pauli::Z)])).collect(); + for (k, exact) in [(1, 0.995131), (2, 0.979836)] { + let c = kicked_ising(127, &heavy_hex_127(), theta(k), 5); + let r = expectation(&c, &mz, &[], &SpdConfig { threshold: 1e-7, max_weight: None }).unwrap(); + let err = (r.value - exact).abs(); + // The published values carry six decimals: allow their rounding. + assert!(err <= r.truncation_bound + 5e-7, "k = {k}: error {err:e} exceeds the bound {:e}", r.truncation_bound); + assert!(err < 3e-5, "k = {k}: Mz = {} against {exact}", r.value); + } +} + +#[test] +fn twenty_step_z62_lies_among_the_converged_classical_values() { + let c = kicked_ising(127, &heavy_hex_127(), theta(2), 20); + let r = expect_z(&c, 62, &SpdConfig { threshold: 1e-4, max_weight: None }).unwrap(); + assert!((0.9776..=0.9796).contains(&r.value), "⟨Z₆₂⟩ = {}", r.value); +} diff --git a/wai/wai-rs/Cargo.lock b/wai/wai-rs/Cargo.lock index ddbd51d31..8f783e1d6 100644 --- a/wai/wai-rs/Cargo.lock +++ b/wai/wai-rs/Cargo.lock @@ -3900,7 +3900,7 @@ dependencies = [ [[package]] name = "wai-quantum" -version = "0.3.21" +version = "0.3.22" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-rs/Cargo.toml b/wai/wai-rs/Cargo.toml index 3cd82e8cf..552ccadff 100644 --- a/wai/wai-rs/Cargo.toml +++ b/wai/wai-rs/Cargo.toml @@ -335,6 +335,10 @@ quantum_stabilizer = ["quantum", "wai-quantum/quantum_stabilizer"] # Pauli dynamics). Computes ⟨O⟩ on large # low-magic circuits where no state vector exists; grows with magic, not qubits. quantum_pauli = ["quantum", "wai-quantum/quantum_pauli"] +# Sparse Pauli dynamics at utility scale (arbitrary angles, up to 1024 qubits), +# and its signed receipts. +quantum_spd = ["quantum", "wai-quantum/quantum_spd"] +quantum_spd_receipt = ["quantum_spd", "wai-quantum/quantum_spd_receipt"] quantum_qir = ["quantum_toolchain", "wai-quantum/quantum_qir"] # Variational quantum machine learning (extensions/quantum-sim § VML, # `wai.quantum.vml`): the interference-based ML the thermodynamic-computing pivot @@ -460,6 +464,7 @@ required-features = [ "quantum_channel", "quantum_comm", "quantum_qudit", + "quantum_spd_receipt", ] [[bin]] diff --git a/wai/wai-rs/src/bin/wai_quantum.rs b/wai/wai-rs/src/bin/wai_quantum.rs index d10953db5..a48f7cb95 100644 --- a/wai/wai-rs/src/bin/wai_quantum.rs +++ b/wai/wai-rs/src/bin/wai_quantum.rs @@ -27,6 +27,14 @@ //! wai-quantum qutrit # d = 3, not just qubits //! ``` //! +//! At utility scale, sparse Pauli dynamics computes observables of circuits far +//! past any state vector, with arbitrary rotation angles: +//! +//! ```text +//! wai-quantum utility --theta 0.2 --steps 20 # the 127-qubit heavy-hex experiment +//! wai-quantum utility --receipt | wai-quantum verify +//! ``` +//! //! The `channel` command prints the Holevo bound before and after, which can only //! fall: noise never creates information. Two results there are easy to miss — //! fully damping `|1⟩` yields a *pure* output that still carries nothing (both @@ -42,6 +50,7 @@ use wai::quantum_info::{chsh, DensityMatrix}; use wai::quantum_qudit::QuditState; use wai::quantum_energy::{EnergyClass, MaybeLabelled}; use wai::quantum_receipt::QuantumReceipt; +use wai::quantum_spd_receipt::{SpdJob, SpdReceipt}; use wai::quantum_source::{holevo_bound, schumacher_limit, von_neumann_entropy}; use wai::quantum_toolchain as qt; @@ -115,6 +124,8 @@ fn usage() -> ! { wai-quantum bb84 [--rounds N] [--seed S] [--eavesdrop] [--json]\n \ wai-quantum chsh [--json]\n \ wai-quantum qutrit [-n QUTRITS] [--json]\n\n\ + UTILITY SCALE:\n \ + wai-quantum utility [--theta T] [--steps N] [--threshold D] [-q QUBIT] [--json] [--receipt]\n\n\ Channels: identity, bitflip, dephasing, depolarizing, damping.\n\ Run `wai-quantum list` for the algorithm catalog." ); @@ -257,6 +268,24 @@ fn main() { // carries its acquisition class inside the signature // (`quantum_energy::Labelled`, as `wai_quantum_meter` emits). let src = read_stdin(); + if let Some(any) = MaybeLabelled::::from_json(&src) { + let ok = any.verify(); + let r = any.receipt(); + println!("signature {}", if ok { "VALID" } else { "INVALID" }); + println!("kind sparse Pauli dynamics"); + println!("circuit-hash {}", hx(&r.circuit_hash)); + println!("query-hash {}", hx(&r.query_hash)); + println!("value {}", r.result.value); + println!("bound {}", r.result.truncation_bound); + println!("work {} term-updates (peak {} terms)", r.result.term_updates, r.result.peak_terms); + println!("joules {} uJ", r.joules_micro); + println!("energy class {}", energy_class_line(any.energy_class(), r.joules_micro)); + println!("signer {}", r.signer_id); + if !ok { + std::process::exit(1); + } + return; + } match MaybeLabelled::::from_json(&src) { None => { println!("INVALID (could not parse receipt JSON, or its energy class is malformed)"); @@ -285,6 +314,49 @@ fn main() { // cannot answer this for a subsystem — half a Bell pair has no state // vector — so `--keep` traces the rest out and reports the entropy that // only a density matrix can see. + // The 127-qubit heavy-hex kicked-Ising experiment, by sparse Pauli + // dynamics: ⟨Z_q⟩ after `--steps` Trotter steps of RX(θ_h) on every + // qubit and RZZ(−π/2) on every edge, and the bound on how far + // truncation can have moved it. `--receipt` seals the run instead. + "utility" => { + use wai::quantum_spd::{heavy_hex_127, kicked_ising, Pauli, SpdConfig}; + let theta: f64 = num(args, "--theta", "", std::f64::consts::FRAC_PI_4); + let steps: u32 = num(args, "--steps", "", 20); + let threshold: f64 = num(args, "--threshold", "", 1e-3); + let qubit: u32 = num(args, "--qubit", "-q", 62); + let job = SpdJob { + circuit: kicked_ising(127, &heavy_hex_127(), theta, steps), + observable: vec![(1.0, vec![(qubit, Pauli::Z)])], + ones: vec![], + config: SpdConfig { threshold, max_weight: None }, + }; + if flag(args, "--receipt") { + // The same deterministic demo key as `receipt`: the run is + // bit-identical everywhere, so is its receipt. + let key = SigningKey::from_bytes(&[42u8; 32]); + let r = SpdReceipt::seal(&key, "did:key:wai-quantum-cli", &job, None, |run| { + run(); + 0 + }) + .unwrap_or_else(|e| die(&format!("seal: {e}"))); + println!("{}", r.to_json()); + return; + } + let r = job.run().unwrap_or_else(|e| die(&e.to_string())); + if flag(args, "--json") { + println!( + "{{\"bound\":{},\"final_terms\":{},\"peak_terms\":{},\"qubit\":{qubit},\"steps\":{steps},\ + \"term_updates\":{},\"theta_h\":{theta},\"threshold\":{threshold},\"value\":{}}}", + r.truncation_bound, r.final_terms, r.peak_terms, r.term_updates, r.value + ); + } else { + println!(" {:.6}", r.value); + println!("bound {:.3e} (|exact - value| is at most this)", r.truncation_bound); + println!("peak terms {}", r.peak_terms); + println!("work {} term-updates", r.term_updates); + } + } + "entropy" => { let c = circuit_from(args, algo); let sv = c.simulate().unwrap_or_else(|e| die(&format!("simulate: {e:?}"))); diff --git a/wai/wai-rs/src/lib.rs b/wai/wai-rs/src/lib.rs index 3922b7d44..ff3e83f15 100644 --- a/wai/wai-rs/src/lib.rs +++ b/wai/wai-rs/src/lib.rs @@ -581,6 +581,12 @@ pub use wai_quantum::quantum_stabilizer; // large low-magic circuits where no 2^n state vector exists. Verified vs dense. #[cfg(feature = "quantum_pauli")] pub use wai_quantum::quantum_pauli; + +#[cfg(feature = "quantum_spd")] +pub use wai_quantum::quantum_spd; + +#[cfg(feature = "quantum_spd_receipt")] +pub use wai_quantum::quantum_spd_receipt; #[cfg(feature = "quantum_qir")] pub use wai_quantum::quantum_qir; diff --git a/wai/wai-web/Cargo.lock b/wai/wai-web/Cargo.lock index 8e1dfad7c..c75264d66 100644 --- a/wai/wai-web/Cargo.lock +++ b/wai/wai-web/Cargo.lock @@ -734,7 +734,7 @@ dependencies = [ [[package]] name = "wai-quantum" -version = "0.3.21" +version = "0.3.22" dependencies = [ "blake3", "ed25519-dalek", -- 2.54.0 (Apple Git-157) From ba653dac59a6818769fe20db8126d6b6e1e744ea Mon Sep 17 00:00:00 2001 From: Transaction Science Date: Thu, 1 Oct 2026 15:15:24 -0400 Subject: [PATCH 2/2] ci: run wai-quantum's integration tests on every architecture in the matrix The energy-receipts workflow ran wai-quantum with --lib only, so the integration tests never ran in CI: tests/bit_identity.rs, which pins one result of each f64 layer and exists to make every architecture produce the same bits, and tests/utility_reference.rs, which holds the 127-qubit results to published values. --tests adds them on linux, linux-arm64 and macos. Co-Authored-By: Claude Opus 5.5 --- .github/workflows/wai-energy-receipts.yml | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/.github/workflows/wai-energy-receipts.yml b/.github/workflows/wai-energy-receipts.yml index c994a08a6..ff466c560 100644 --- a/.github/workflows/wai-energy-receipts.yml +++ b/.github/workflows/wai-energy-receipts.yml @@ -19,9 +19,13 @@ name: WAI energy receipts # the set can build (one whose required features are all on), which also # proves those bins compile; and the one integration test file # tests/video_receipt.rs. -# 2. wai-quantum: `cargo test --locked --features full --lib`: the crate's -# unit tests, among them the pins and label tests of every quantum -# receipt. The tests already marked #[ignore] stay ignored. +# 2. wai-quantum: `cargo test --locked --features full --lib --tests`: the +# crate's unit tests, among them the pins and label tests of every quantum +# receipt, and its integration tests — tests/bit_identity.rs, which pins +# one result of each f64 layer so every architecture here must produce +# the same bits, and tests/utility_reference.rs, which holds the 127-qubit +# results to published values. The tests already marked #[ignore] stay +# ignored. # 3. wai-quantum-serve: `cargo test --locked --lib`, native: the verdict # logic behind POST /verify, for a receipt in either form, and that an # unmetered receipt's 0 is rendered as the unmetered marker (null), never @@ -129,11 +133,13 @@ jobs: --lib --bins --test video_receipt # wai-quantum is published separately; its receipts' legacy pins and the - # labelled forms are tested in its own workspace, against its own lock. - - name: wai-quantum — receipt pins and labels + # labelled forms are tested in its own workspace, against its own lock, + # with the integration tests that pin its f64 layers' bits on every + # architecture this matrix runs. + - name: wai-quantum — receipt pins, labels and bit identity working-directory: wai/wai-quantum shell: bash - run: cargo test --locked --features full --lib + run: cargo test --locked --features full --lib --tests # The HTTP handler's POST /verify must accept a receipt in either form — # as sealed, or labelled as wai_quantum_meter emits it — and report the -- 2.54.0 (Apple Git-157)