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