wai-quantum 0.3.28: the dense state vector at any angle, the same bits threaded or not #147
17 changed files with 1304 additions and 44 deletions
2
wai/wai-quantum-serve/Cargo.lock
generated
2
wai/wai-quantum-serve/Cargo.lock
generated
|
|
@ -589,7 +589,7 @@ dependencies = [
|
|||
|
||||
[[package]]
|
||||
name = "wai-quantum"
|
||||
version = "0.3.27"
|
||||
version = "0.3.28"
|
||||
dependencies = [
|
||||
"blake3",
|
||||
"ed25519-dalek",
|
||||
|
|
|
|||
2
wai/wai-quantum/Cargo.lock
generated
2
wai/wai-quantum/Cargo.lock
generated
|
|
@ -415,7 +415,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a"
|
|||
|
||||
[[package]]
|
||||
name = "wai-quantum"
|
||||
version = "0.3.27"
|
||||
version = "0.3.28"
|
||||
dependencies = [
|
||||
"blake3",
|
||||
"ed25519-dalek",
|
||||
|
|
|
|||
|
|
@ -1,6 +1,6 @@
|
|||
[package]
|
||||
name = "wai-quantum"
|
||||
version = "0.3.27"
|
||||
version = "0.3.28"
|
||||
edition = "2024"
|
||||
rust-version = "1.97"
|
||||
license = "Apache-2.0"
|
||||
|
|
@ -62,6 +62,9 @@ quantum_tn = []
|
|||
# Gaussian boson sampling, exactly: hafnians, torontonians, lossy Gaussian states,
|
||||
# pattern probabilities, exact samples and the marginals any sampler must match.
|
||||
quantum_gbs = []
|
||||
# The dense state vector at any angle: gate fusion, diagonal fast paths, threads,
|
||||
# the same bits everywhere.
|
||||
quantum_sv = []
|
||||
# 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"]
|
||||
|
|
@ -101,13 +104,17 @@ quantum_phasor_meter = ["_signing"]
|
|||
|
||||
# Everything.
|
||||
full = [
|
||||
"quantum", "quantum_toolchain", "quantum_mps", "quantum_stabilizer", "quantum_pauli", "quantum_spd", "quantum_spd_receipt", "quantum_bptn", "quantum_resource", "quantum_resource_receipt", "quantum_frame", "quantum_tn", "quantum_gbs", "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_bptn", "quantum_resource", "quantum_resource_receipt", "quantum_frame", "quantum_tn", "quantum_gbs", "quantum_sv", "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 = "sv_run"
|
||||
required-features = ["quantum_sv"]
|
||||
|
||||
[[example]]
|
||||
name = "gbs_check"
|
||||
required-features = ["quantum_gbs"]
|
||||
|
|
|
|||
|
|
@ -72,6 +72,7 @@ never creates information.
|
|||
| `quantum_frame` | **noisy Clifford circuits at QEC scale**: reads the ecosystem's circuit text; samples detection events by bit-packed Pauli-frame simulation (64 shots per word); derives the detector error model by backward sensitivity, refusing non-deterministic detectors; generates the rotated surface-code memory experiment byte for byte as the reference layout; and decodes by **exact minimum-weight perfect matching** (Edmonds' blossom on integer weights). Measures logical error rates and fits the code model resource estimates assume |
|
||||
| `quantum_tn` | **tensor-network contraction**: circuits (any unitary on any qubits) become networks; a seeded hyper-optimiser finds contraction trees from randomised-greedy and multilevel-partition families, reconfigures them optimally in subtrees, and slices them to a memory cap with a second search around the slices; slices contract in a fixed order, threaded, to the same bits everywhere. Single amplitudes or open-output batches |
|
||||
| `quantum_gbs` | **Gaussian boson sampling, exactly**: hafnians (power traces, Hessenberg + characteristic polynomial, `O(n³ 2^{n/2})`), torontonians, lossy Gaussian states from squeezers and an interferometer, exact click and photon-number probabilities, exact samples drawn mode by mode, and a check of observed samples against the exact one- and two-mode click statistics |
|
||||
| `quantum_sv` | **the dense state vector at any angle**: one-, two- and `k`-qubit kernels with a diagonal fast path, gate fusion around multi-qubit gates, threads over independent blocks; every kernel path computes each amplitude in the same order, so threading never changes a bit. Agrees with an independent simulator to 2·10⁻¹⁶ |
|
||||
| `quantum_compile` | Clifford routing + stabilizer-tableau equivalence proof |
|
||||
| `quantum_atom` | neutral-atom register preparation (Hungarian / LSAP) |
|
||||
| `quantum_qir` | **QIR export + ingest** — emit and parse QIR (the QIR Alliance's LLVM-based IR); emitted modules validated with `llvm-as`, and 9 of 10 catalog algorithms re-ingest to a bit-identical statevector |
|
||||
|
|
@ -177,6 +178,20 @@ fitted code model comes out at `a ≈ 0.023`, `p* ≈ 0.0100`. That is a measure
|
|||
value for the constants a resource estimate (`quantum_resource`) otherwise takes
|
||||
on trust, and it can be fed straight into one.
|
||||
|
||||
## State vector
|
||||
|
||||
```sh
|
||||
wai-quantum sv --file circuit.jsonl --top 8 # likeliest outcomes
|
||||
wai-quantum sv --file circuit.jsonl --expect "Z0 Z1; X3"
|
||||
wai-quantum sv --file circuit.jsonl --sample 1000 --seed 7
|
||||
```
|
||||
|
||||
`quantum_sv` runs any gate at any angle on up to 32 qubits, in double
|
||||
precision. `gates` provides the standard gates by name (H, S, T, √X, √Y, √W,
|
||||
rotations, CX, CZ, iSWAP, fSim, …) for it and for the tensor networks. A
|
||||
26-qubit, depth-20 random circuit runs in about nine seconds on one machine,
|
||||
and the bits do not depend on how many threads ran it.
|
||||
|
||||
## Tensor networks
|
||||
|
||||
An amplitude of a circuit is a tensor network, and the order it is contracted
|
||||
|
|
|
|||
35
wai/wai-quantum/examples/sv_run.rs
Normal file
35
wai/wai-quantum/examples/sv_run.rs
Normal file
|
|
@ -0,0 +1,35 @@
|
|||
//! Run a circuit (JSON lines, one gate per line: {"q":[a,b],"m":[re,im,…]}) on
|
||||
//! the dense state vector and print the final state, or time it.
|
||||
//!
|
||||
//! ```text
|
||||
//! cargo run --release --example sv_run --features quantum_sv -- CIRCUIT N WIDTH [--state | --time]
|
||||
//! ```
|
||||
use wai_quantum::gates::{complex, Gate};
|
||||
use wai_quantum::quantum_sv::StateVector;
|
||||
|
||||
fn main() {
|
||||
let a: Vec<String> = std::env::args().collect();
|
||||
let text = std::fs::read_to_string(&a[1]).expect("read circuit");
|
||||
let n: u32 = a[2].parse().unwrap();
|
||||
let width: usize = a[3].parse().unwrap();
|
||||
let gates: Vec<Gate> = text
|
||||
.lines()
|
||||
.filter(|l| !l.trim().is_empty())
|
||||
.map(|l| {
|
||||
let v: serde_json::Value = serde_json::from_str(l).unwrap();
|
||||
let qubits = v["q"].as_array().unwrap().iter().map(|x| x.as_u64().unwrap() as u32).collect();
|
||||
let m: Vec<f64> = v["m"].as_array().unwrap().iter().map(|x| x.as_f64().unwrap()).collect();
|
||||
Gate { qubits, matrix: m.chunks(2).map(|c| complex(c[0], c[1])).collect() }
|
||||
})
|
||||
.collect();
|
||||
let mut sv = StateVector::zero(n).unwrap();
|
||||
let t0 = std::time::Instant::now();
|
||||
sv.run(&gates, width).unwrap();
|
||||
let dt = t0.elapsed().as_secs_f64();
|
||||
if a.iter().any(|x| x == "--state") {
|
||||
let amps: Vec<String> = (0..sv.re.len()).map(|x| format!("[{:e},{:e}]", sv.re[x], sv.im[x])).collect();
|
||||
println!("[{}]", amps.join(","));
|
||||
} else {
|
||||
println!("{n} qubits, {} gates, width {width}: {dt:.3} s, norm {:.15}", gates.len(), sv.norm2());
|
||||
}
|
||||
}
|
||||
139
wai/wai-quantum/src/gates.rs
Normal file
139
wai/wai-quantum/src/gates.rs
Normal file
|
|
@ -0,0 +1,139 @@
|
|||
//! Gates as matrices, shared by the dense state vector and the tensor
|
||||
//! networks: any unitary on any qubits, and the standard ones by name.
|
||||
//!
|
||||
//! Matrices are row-major, rows are outputs and columns inputs, and the first
|
||||
//! listed qubit is the most significant bit of the local index. Angles enter
|
||||
//! through the crate's own `sin`/`cos`, so a named gate is the same bits on
|
||||
//! every machine.
|
||||
|
||||
pub use crate::linalg::C;
|
||||
|
||||
/// A gate: a unitary on `qubits`, as a `2ᵏ × 2ᵏ` row-major matrix whose rows
|
||||
/// are outputs and columns inputs, with `qubits[0]` the most significant bit.
|
||||
#[derive(Clone, Debug)]
|
||||
pub struct Gate {
|
||||
pub qubits: Vec<u32>,
|
||||
pub matrix: Vec<C>,
|
||||
}
|
||||
|
||||
/// A complex number for building gates.
|
||||
pub fn complex(re: f64, im: f64) -> C {
|
||||
C { re, im }
|
||||
}
|
||||
|
||||
const Z: C = C { re: 0.0, im: 0.0 };
|
||||
const O: C = C { re: 1.0, im: 0.0 };
|
||||
const I: C = C { re: 0.0, im: 1.0 };
|
||||
|
||||
fn cis(theta: f64) -> C {
|
||||
let (s, c) = crate::repro::sin_cos(theta);
|
||||
C { re: c, im: s }
|
||||
}
|
||||
|
||||
fn one(q: u32, m: [C; 4]) -> Gate {
|
||||
Gate { qubits: vec![q], matrix: m.to_vec() }
|
||||
}
|
||||
|
||||
fn two(a: u32, b: u32, m: [C; 16]) -> Gate {
|
||||
Gate { qubits: vec![a, b], matrix: m.to_vec() }
|
||||
}
|
||||
|
||||
impl Gate {
|
||||
pub fn h(q: u32) -> Gate {
|
||||
let r = core::f64::consts::FRAC_1_SQRT_2;
|
||||
one(q, [complex(r, 0.0), complex(r, 0.0), complex(r, 0.0), complex(-r, 0.0)])
|
||||
}
|
||||
pub fn x(q: u32) -> Gate {
|
||||
one(q, [Z, O, O, Z])
|
||||
}
|
||||
pub fn y(q: u32) -> Gate {
|
||||
one(q, [Z, complex(0.0, -1.0), I, Z])
|
||||
}
|
||||
pub fn z(q: u32) -> Gate {
|
||||
one(q, [O, Z, Z, complex(-1.0, 0.0)])
|
||||
}
|
||||
pub fn s(q: u32) -> Gate {
|
||||
one(q, [O, Z, Z, I])
|
||||
}
|
||||
pub fn sdg(q: u32) -> Gate {
|
||||
one(q, [O, Z, Z, complex(0.0, -1.0)])
|
||||
}
|
||||
pub fn t(q: u32) -> Gate {
|
||||
let r = core::f64::consts::FRAC_1_SQRT_2;
|
||||
one(q, [O, Z, Z, complex(r, r)])
|
||||
}
|
||||
pub fn tdg(q: u32) -> Gate {
|
||||
let r = core::f64::consts::FRAC_1_SQRT_2;
|
||||
one(q, [O, Z, Z, complex(r, -r)])
|
||||
}
|
||||
/// √X = ½[[1+i, 1−i], [1−i, 1+i]].
|
||||
pub fn sqrt_x(q: u32) -> Gate {
|
||||
let (p, m) = (complex(0.5, 0.5), complex(0.5, -0.5));
|
||||
one(q, [p, m, m, p])
|
||||
}
|
||||
/// √Y = ½[[1+i, −1−i], [1+i, 1+i]].
|
||||
pub fn sqrt_y(q: u32) -> Gate {
|
||||
let p = complex(0.5, 0.5);
|
||||
one(q, [p, complex(-0.5, -0.5), p, p])
|
||||
}
|
||||
/// √W, with W = (X + Y)/√2: √X conjugated by diag(1, e^{iπ/4}).
|
||||
pub fn sqrt_w(q: u32) -> Gate {
|
||||
let (p, m) = (complex(0.5, 0.5), complex(0.5, -0.5));
|
||||
let (e, ec) = (cis(core::f64::consts::FRAC_PI_4), cis(-core::f64::consts::FRAC_PI_4));
|
||||
let mul = |a: C, b: C| complex(a.re * b.re - a.im * b.im, a.re * b.im + a.im * b.re);
|
||||
one(q, [p, mul(m, ec), mul(m, e), p])
|
||||
}
|
||||
/// `exp(−iθX/2)`.
|
||||
pub fn rx(q: u32, theta: f64) -> Gate {
|
||||
let (s, c) = crate::repro::sin_cos(theta / 2.0);
|
||||
one(q, [complex(c, 0.0), complex(0.0, -s), complex(0.0, -s), complex(c, 0.0)])
|
||||
}
|
||||
/// `exp(−iθY/2)`.
|
||||
pub fn ry(q: u32, theta: f64) -> Gate {
|
||||
let (s, c) = crate::repro::sin_cos(theta / 2.0);
|
||||
one(q, [complex(c, 0.0), complex(-s, 0.0), complex(s, 0.0), complex(c, 0.0)])
|
||||
}
|
||||
/// `exp(−iθZ/2)`.
|
||||
pub fn rz(q: u32, theta: f64) -> Gate {
|
||||
one(q, [cis(-theta / 2.0), Z, Z, cis(theta / 2.0)])
|
||||
}
|
||||
/// `diag(1, e^{iθ})`.
|
||||
pub fn phase(q: u32, theta: f64) -> Gate {
|
||||
one(q, [O, Z, Z, cis(theta)])
|
||||
}
|
||||
/// Controlled X, control first.
|
||||
pub fn cx(c: u32, t: u32) -> Gate {
|
||||
two(c, t, [O, Z, Z, Z, Z, O, Z, Z, Z, Z, Z, O, Z, Z, O, Z])
|
||||
}
|
||||
pub fn cz(a: u32, b: u32) -> Gate {
|
||||
two(a, b, [O, Z, Z, Z, Z, O, Z, Z, Z, Z, O, Z, Z, Z, Z, complex(-1.0, 0.0)])
|
||||
}
|
||||
pub fn swap(a: u32, b: u32) -> Gate {
|
||||
two(a, b, [O, Z, Z, Z, Z, Z, O, Z, Z, O, Z, Z, Z, Z, Z, O])
|
||||
}
|
||||
pub fn iswap(a: u32, b: u32) -> Gate {
|
||||
two(a, b, [O, Z, Z, Z, Z, Z, I, Z, Z, I, Z, Z, Z, Z, Z, O])
|
||||
}
|
||||
/// `diag(1, 1, 1, e^{iφ})`.
|
||||
pub fn cphase(a: u32, b: u32, phi: f64) -> Gate {
|
||||
two(a, b, [O, Z, Z, Z, Z, O, Z, Z, Z, Z, O, Z, Z, Z, Z, cis(phi)])
|
||||
}
|
||||
/// `exp(−iθ Z⊗Z / 2)`.
|
||||
pub fn rzz(a: u32, b: u32, theta: f64) -> Gate {
|
||||
let (m, p) = (cis(-theta / 2.0), cis(theta / 2.0));
|
||||
two(a, b, [m, Z, Z, Z, Z, p, Z, Z, Z, Z, p, Z, Z, Z, Z, m])
|
||||
}
|
||||
/// The fermionic simulation gate: an iSWAP-like rotation by `θ` and a
|
||||
/// conditional phase `e^{−iφ}` on `|11⟩`.
|
||||
pub fn fsim(a: u32, b: u32, theta: f64, phi: f64) -> Gate {
|
||||
let (s, c) = crate::repro::sin_cos(theta);
|
||||
let (cc, ms) = (complex(c, 0.0), complex(0.0, -s));
|
||||
two(a, b, [O, Z, Z, Z, Z, cc, ms, Z, Z, ms, cc, Z, Z, Z, Z, cis(-phi)])
|
||||
}
|
||||
|
||||
/// Whether the matrix is diagonal (exactly zero off the diagonal).
|
||||
pub fn is_diagonal(&self) -> bool {
|
||||
let d = 1usize << self.qubits.len();
|
||||
(0..d).all(|r| (0..d).all(|c| r == c || (self.matrix[r * d + c].re == 0.0 && self.matrix[r * d + c].im == 0.0)))
|
||||
}
|
||||
}
|
||||
|
|
@ -17,9 +17,13 @@
|
|||
#[cfg(feature = "_signing")]
|
||||
mod merkle;
|
||||
|
||||
#[cfg(any(feature = "quantum_mps", feature = "quantum_bptn", feature = "quantum_tn", feature = "quantum_gbs"))]
|
||||
#[cfg(any(feature = "quantum_mps", feature = "quantum_bptn", feature = "quantum_tn", feature = "quantum_gbs", feature = "quantum_sv"))]
|
||||
mod linalg;
|
||||
|
||||
// Gates as matrices, shared by the state vector and the tensor networks.
|
||||
#[cfg(any(feature = "quantum_tn", feature = "quantum_sv"))]
|
||||
pub mod gates;
|
||||
|
||||
#[cfg(any(
|
||||
feature = "quantum_spd",
|
||||
feature = "quantum_pauli",
|
||||
|
|
@ -30,7 +34,8 @@ mod linalg;
|
|||
feature = "quantum_resource",
|
||||
feature = "quantum_frame",
|
||||
feature = "quantum_tn",
|
||||
feature = "quantum_gbs"
|
||||
feature = "quantum_gbs",
|
||||
feature = "quantum_sv"
|
||||
))]
|
||||
mod repro;
|
||||
|
||||
|
|
@ -106,6 +111,9 @@ pub mod quantum_tn;
|
|||
#[cfg(feature = "quantum_gbs")]
|
||||
pub mod quantum_gbs;
|
||||
|
||||
#[cfg(feature = "quantum_sv")]
|
||||
pub mod quantum_sv;
|
||||
|
||||
#[cfg(feature = "quantum_info")]
|
||||
pub mod quantum_info;
|
||||
|
||||
|
|
|
|||
870
wai/wai-quantum/src/quantum_sv.rs
Normal file
870
wai/wai-quantum/src/quantum_sv.rs
Normal file
|
|
@ -0,0 +1,870 @@
|
|||
//! The dense state vector, any gate at any angle — `wai.quantum.sv`.
|
||||
//!
|
||||
//! The crate's byte-exact simulator ([`crate::quantum`]) is fixed-point and
|
||||
//! limited to the dyadic gate family. That limit is deliberate: it is what
|
||||
//! makes its amplitudes a portable hash. Most circuits are not dyadic.
|
||||
//! Rotations at arbitrary angles, √W, fSim and random two-qubit unitaries need
|
||||
//! a general state vector. It is also the ground truth every approximate
|
||||
//! method is measured against, up to about thirty qubits.
|
||||
//!
|
||||
//! This module is that state vector, reproducible bit for bit.
|
||||
//!
|
||||
//! - **Layout.** Real and imaginary parts are held in separate arrays, so the
|
||||
//! update loops vectorise. Qubit 0 is the most significant bit, as in
|
||||
//! [`crate::gates`] and the tensor networks.
|
||||
//! - **Kernels.** One-qubit, two-qubit and general `k`-qubit kernels, plus a
|
||||
//! diagonal fast path (phases only) for `Z`, `S`, `T`, `RZ`, `CZ`, phase
|
||||
//! gates and fused diagonal blocks.
|
||||
//! - **Fusion.** [`fuse`] merges consecutive gates into one matrix while
|
||||
//! their union spans at most `width` qubits. Fewer, denser passes over the
|
||||
//! state are what bound a state vector's speed. Fusion changes the rounding,
|
||||
//! so the width is part of a run's definition. A run's bits are fixed by its
|
||||
//! gates and its width.
|
||||
//! - **Threads.** Every amplitude group is updated independently, so a gate's
|
||||
//! result does not depend on how the work is split. Reductions (norms,
|
||||
//! expectations, sampling) are summed in a fixed number of chunks, each in
|
||||
//! order, and the chunks are added in order. Threaded natively, sequential
|
||||
//! under WASI: the same bits either way.
|
||||
//! - **One arithmetic per amplitude.** Gates on the top qubits leave too few
|
||||
//! independent blocks to share between threads, so they take other paths:
|
||||
//! split at the top target bit, or into quarters for adjacent pairs. Every
|
||||
//! path computes each amplitude through the same row sum, in the same
|
||||
//! order, so the path taken never changes a bit. A test forces one, two,
|
||||
//! three, eight and thirty-two threads and demands identical states.
|
||||
//!
|
||||
//! # Checked
|
||||
//!
|
||||
//! - **Final states** of random circuits (phased-XZ rotations, fSim, partial
|
||||
//! iSWAPs) on 8, 12 and 16 qubits agree with an independent simulator to
|
||||
//! 2·10⁻¹⁶ at every fusion width.
|
||||
//! - **Amplitudes** contracted as tensor networks agree with this state vector
|
||||
//! to 10⁻¹² on 18 qubits.
|
||||
//! - **Speed.** A 26-qubit, depth-20 brickwork circuit of 770 gates runs in
|
||||
//! about nine seconds at fusion width 2, in double precision, on one
|
||||
//! machine.
|
||||
|
||||
use crate::gates::Gate;
|
||||
pub use crate::gates::C;
|
||||
|
||||
/// The largest register this module allocates (2³² amplitudes, 64 GiB).
|
||||
pub const MAX_QUBITS: u32 = 32;
|
||||
|
||||
/// Reductions are summed in this many fixed chunks.
|
||||
const CHUNKS: usize = 64;
|
||||
|
||||
#[derive(Clone, Debug, PartialEq, Eq)]
|
||||
pub enum SvError {
|
||||
/// More qubits than [`MAX_QUBITS`].
|
||||
TooLarge,
|
||||
/// A gate whose matrix is not `4^k` entries, repeats a qubit, or acts
|
||||
/// outside the register.
|
||||
BadGate(usize),
|
||||
}
|
||||
|
||||
/// A pure state of `n` qubits.
|
||||
#[derive(Clone, Debug, PartialEq)]
|
||||
pub struct StateVector {
|
||||
pub n: u32,
|
||||
pub re: Vec<f64>,
|
||||
pub im: Vec<f64>,
|
||||
}
|
||||
|
||||
/// Run `f` over `re`/`im` split into contiguous slices of whole `unit`s,
|
||||
/// threaded where the target has threads; `f` gets each slice and its offset.
|
||||
/// Each thread owns its slices, so nothing is shared.
|
||||
/// A kernel over a slice of the real parts, the matching imaginary parts, and
|
||||
/// the slice's offset.
|
||||
type SliceFn<'a> = dyn Fn(&mut [f64], &mut [f64], usize) + Sync + 'a;
|
||||
|
||||
fn par_slices(re: &mut [f64], im: &mut [f64], unit: usize, f: &SliceFn) {
|
||||
#[cfg(not(target_arch = "wasm32"))]
|
||||
{
|
||||
let threads = available_threads();
|
||||
let units = re.len() / unit;
|
||||
if threads > 1 && units >= threads && re.len() >= 1 << 14 {
|
||||
let chunk = units.div_ceil(threads) * unit;
|
||||
std::thread::scope(|sc| {
|
||||
for (k, (r, i)) in re.chunks_mut(chunk).zip(im.chunks_mut(chunk)).enumerate() {
|
||||
sc.spawn(move || f(r, i, k * chunk));
|
||||
}
|
||||
});
|
||||
return;
|
||||
}
|
||||
}
|
||||
let _ = unit;
|
||||
f(re, im, 0);
|
||||
}
|
||||
|
||||
/// Insert a zero bit at each position in `pos` (ascending) into `t`.
|
||||
#[inline]
|
||||
fn spread(mut t: usize, pos: &[u32]) -> usize {
|
||||
for &p in pos {
|
||||
let low = t & ((1usize << p) - 1);
|
||||
t = ((t >> p) << (p + 1)) | low;
|
||||
}
|
||||
t
|
||||
}
|
||||
|
||||
/// One output amplitude of a gate on two or more qubits: the matrix row
|
||||
/// against the group's members in local-index order, accumulated from zero in
|
||||
/// that order. Every kernel for such gates goes through here, so each
|
||||
/// amplitude is the same sequence of operations whichever path computes it.
|
||||
#[inline(always)]
|
||||
fn row_dot(row: &[C], vr: &[f64], vi: &[f64]) -> (f64, f64) {
|
||||
let (mut sr, mut si) = (0.0, 0.0);
|
||||
for (c, e) in row.iter().enumerate() {
|
||||
sr += e.re * vr[c] - e.im * vi[c];
|
||||
si += e.re * vi[c] + e.im * vr[c];
|
||||
}
|
||||
(sr, si)
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
thread_local! {
|
||||
/// Tests force a thread count here to compare kernel paths.
|
||||
static FORCE_THREADS: core::cell::Cell<Option<usize>> = const { core::cell::Cell::new(None) };
|
||||
}
|
||||
|
||||
fn available_threads() -> usize {
|
||||
#[cfg(test)]
|
||||
if let Some(t) = FORCE_THREADS.with(|f| f.get()) {
|
||||
return t;
|
||||
}
|
||||
#[cfg(not(target_arch = "wasm32"))]
|
||||
{
|
||||
std::thread::available_parallelism().map_or(1, |t| t.get())
|
||||
}
|
||||
#[cfg(target_arch = "wasm32")]
|
||||
{
|
||||
1
|
||||
}
|
||||
}
|
||||
|
||||
/// Apply a gate to a slice holding whole groups (block-aligned). One- and
|
||||
/// two-qubit gates take unrolled loops whose inner stride is contiguous.
|
||||
fn kernel_block(re: &mut [f64], im: &mut [f64], m: &[C], offsets: &[usize], sorted: &[u32]) {
|
||||
let k = sorted.len();
|
||||
if k == 1 {
|
||||
let s = offsets[1];
|
||||
let (m00, m01, m10, m11) = (m[0], m[1], m[2], m[3]);
|
||||
let mut base = 0;
|
||||
while base < re.len() {
|
||||
for i in base..base + s {
|
||||
let (ar, ai, br, bi) = (re[i], im[i], re[i + s], im[i + s]);
|
||||
re[i] = m00.re * ar - m00.im * ai + m01.re * br - m01.im * bi;
|
||||
im[i] = m00.re * ai + m00.im * ar + m01.re * bi + m01.im * br;
|
||||
re[i + s] = m10.re * ar - m10.im * ai + m11.re * br - m11.im * bi;
|
||||
im[i + s] = m10.re * ai + m10.im * ar + m11.re * bi + m11.im * br;
|
||||
}
|
||||
base += 2 * s;
|
||||
}
|
||||
return;
|
||||
}
|
||||
if k == 2 {
|
||||
let (slo, shi) = (1usize << sorted[0], 1usize << sorted[1]);
|
||||
let mut b1 = 0;
|
||||
while b1 < re.len() {
|
||||
let mut b2 = b1;
|
||||
while b2 < b1 + shi {
|
||||
for i in b2..b2 + slo {
|
||||
let vr = [re[i + offsets[0]], re[i + offsets[1]], re[i + offsets[2]], re[i + offsets[3]]];
|
||||
let vi = [im[i + offsets[0]], im[i + offsets[1]], im[i + offsets[2]], im[i + offsets[3]]];
|
||||
for (r, &o) in offsets.iter().enumerate() {
|
||||
let (a, b) = row_dot(&m[r * 4..r * 4 + 4], &vr, &vi);
|
||||
re[i + o] = a;
|
||||
im[i + o] = b;
|
||||
}
|
||||
}
|
||||
b2 += 2 * slo;
|
||||
}
|
||||
b1 += 2 * shi;
|
||||
}
|
||||
return;
|
||||
}
|
||||
let d = offsets.len();
|
||||
let mut vr = vec![0.0; d];
|
||||
let mut vi = vec![0.0; d];
|
||||
for t in 0..re.len() >> k {
|
||||
let base = spread(t, sorted);
|
||||
for (j, &o) in offsets.iter().enumerate() {
|
||||
vr[j] = re[base | o];
|
||||
vi[j] = im[base | o];
|
||||
}
|
||||
for (r, &o) in offsets.iter().enumerate() {
|
||||
let (sr, si) = row_dot(&m[r * d..(r + 1) * d], &vr, &vi);
|
||||
re[base | o] = sr;
|
||||
im[base | o] = si;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Apply a gate whose top target bit splits the data into a lower and an
|
||||
/// upper half: `(lre, lim)` and `(hre, him)` are aligned chunks of the two
|
||||
/// halves, `members[j]` says which half local index `j` lives in and where.
|
||||
fn kernel_pair(lre: &mut [f64], lim: &mut [f64], hre: &mut [f64], him: &mut [f64], m: &[C], members: &[(bool, usize)], rest: &[u32]) {
|
||||
let d = members.len();
|
||||
if d == 2 {
|
||||
let (m00, m01, m10, m11) = (m[0], m[1], m[2], m[3]);
|
||||
// Which of the two is in the upper half decides the matrix orientation.
|
||||
let lower_first = !members[0].0;
|
||||
for i in 0..lre.len() {
|
||||
let (ar, ai, br, bi) = if lower_first { (lre[i], lim[i], hre[i], him[i]) } else { (hre[i], him[i], lre[i], lim[i]) };
|
||||
let nr0 = m00.re * ar - m00.im * ai + m01.re * br - m01.im * bi;
|
||||
let ni0 = m00.re * ai + m00.im * ar + m01.re * bi + m01.im * br;
|
||||
let nr1 = m10.re * ar - m10.im * ai + m11.re * br - m11.im * bi;
|
||||
let ni1 = m10.re * ai + m10.im * ar + m11.re * bi + m11.im * br;
|
||||
if lower_first {
|
||||
(lre[i], lim[i], hre[i], him[i]) = (nr0, ni0, nr1, ni1);
|
||||
} else {
|
||||
(hre[i], him[i], lre[i], lim[i]) = (nr0, ni0, nr1, ni1);
|
||||
}
|
||||
}
|
||||
return;
|
||||
}
|
||||
let k = rest.len() + 1;
|
||||
let mut vr = vec![0.0; d];
|
||||
let mut vi = vec![0.0; d];
|
||||
for t in 0..lre.len() >> (k - 1) {
|
||||
let base = spread(t, rest);
|
||||
for (j, &(up, o)) in members.iter().enumerate() {
|
||||
let (r, i) = if up { (&*hre, &*him) } else { (&*lre, &*lim) };
|
||||
vr[j] = r[base | o];
|
||||
vi[j] = i[base | o];
|
||||
}
|
||||
for (row, &(up, o)) in members.iter().enumerate() {
|
||||
let (sr, si) = row_dot(&m[row * d..(row + 1) * d], &vr, &vi);
|
||||
if up {
|
||||
hre[base | o] = sr;
|
||||
him[base | o] = si;
|
||||
} else {
|
||||
lre[base | o] = sr;
|
||||
lim[base | o] = si;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
impl StateVector {
|
||||
/// `|0…0⟩` on `n` qubits.
|
||||
pub fn zero(n: u32) -> Result<StateVector, SvError> {
|
||||
if n > MAX_QUBITS {
|
||||
return Err(SvError::TooLarge);
|
||||
}
|
||||
let dim = 1usize << n;
|
||||
let mut re = vec![0.0; dim];
|
||||
re[0] = 1.0;
|
||||
Ok(StateVector { n, re, im: vec![0.0; dim] })
|
||||
}
|
||||
|
||||
pub fn amplitude(&self, x: usize) -> C {
|
||||
C { re: self.re[x], im: self.im[x] }
|
||||
}
|
||||
|
||||
fn check(&self, g: &Gate, gi: usize) -> Result<(), SvError> {
|
||||
let k = g.qubits.len();
|
||||
let mut q = g.qubits.clone();
|
||||
q.sort_unstable();
|
||||
q.dedup();
|
||||
if k == 0 || q.len() != k || g.qubits.iter().any(|&x| x >= self.n) || g.matrix.len() != 1 << (2 * k) {
|
||||
return Err(SvError::BadGate(gi));
|
||||
}
|
||||
Ok(())
|
||||
}
|
||||
|
||||
/// Apply one gate.
|
||||
pub fn apply(&mut self, g: &Gate) -> Result<(), SvError> {
|
||||
self.check(g, 0)?;
|
||||
self.apply_unchecked(g);
|
||||
Ok(())
|
||||
}
|
||||
|
||||
/// Apply gates in order, fused to `width` qubits first (`width` 0 or 1
|
||||
/// applies them as given).
|
||||
pub fn run(&mut self, gates: &[Gate], width: usize) -> Result<(), SvError> {
|
||||
for (i, g) in gates.iter().enumerate() {
|
||||
self.check(g, i)?;
|
||||
}
|
||||
let fused;
|
||||
let list = if width >= 2 {
|
||||
fused = fuse(gates, width);
|
||||
&fused
|
||||
} else {
|
||||
gates
|
||||
};
|
||||
for g in list {
|
||||
self.apply_unchecked(g);
|
||||
}
|
||||
Ok(())
|
||||
}
|
||||
|
||||
fn apply_unchecked(&mut self, g: &Gate) {
|
||||
let n = self.n;
|
||||
let k = g.qubits.len();
|
||||
// Bit position of each gate qubit; local index bit (k−1−i) is qubits[i].
|
||||
let bitpos: Vec<u32> = g.qubits.iter().map(|&q| n - 1 - q).collect();
|
||||
let d = 1usize << k;
|
||||
let m = &g.matrix;
|
||||
if g.is_diagonal() {
|
||||
let diag: Vec<C> = (0..d).map(|i| m[i * d + i]).collect();
|
||||
let f = |re: &mut [f64], im: &mut [f64], off: usize| {
|
||||
for (j, (a, b)) in re.iter_mut().zip(im.iter_mut()).enumerate() {
|
||||
let x = off + j;
|
||||
let mut local = 0usize;
|
||||
for &p in &bitpos {
|
||||
local = (local << 1) | ((x >> p) & 1);
|
||||
}
|
||||
let c = diag[local];
|
||||
if c.re == 1.0 && c.im == 0.0 {
|
||||
continue;
|
||||
}
|
||||
let (ar, ai) = (*a, *b);
|
||||
*a = c.re * ar - c.im * ai;
|
||||
*b = c.re * ai + c.im * ar;
|
||||
}
|
||||
};
|
||||
par_slices(&mut self.re, &mut self.im, 1, &f);
|
||||
return;
|
||||
}
|
||||
let mut sorted = bitpos.clone();
|
||||
sorted.sort_unstable();
|
||||
let top = sorted[k - 1];
|
||||
// Every group lies inside a block of 2^(top + 1) amplitudes.
|
||||
let block = 1usize << (top + 1);
|
||||
let len = self.re.len();
|
||||
let offsets: Vec<usize> = (0..d)
|
||||
.map(|local| bitpos.iter().enumerate().fold(0usize, |o, (i, &p)| o | (((local >> (k - 1 - i)) & 1) << p)))
|
||||
.collect();
|
||||
let threads = available_threads();
|
||||
if len / block >= threads || threads == 1 || len < 1 << 14 {
|
||||
let f = |re: &mut [f64], im: &mut [f64], _off: usize| kernel_block(re, im, m, &offsets, &sorted);
|
||||
par_slices(&mut self.re, &mut self.im, block, &f);
|
||||
return;
|
||||
}
|
||||
// Few blocks: split each by the top target bit into its two halves, and
|
||||
// thread over aligned chunk pairs of the halves. A chunk of
|
||||
// 2^(next target + 1) amplitudes holds whole sub-groups.
|
||||
// Two adjacent target bits: each block is four contiguous quarters, one
|
||||
// per member of every group, so aligned chunks of the four zip.
|
||||
if k == 2 && sorted[0] + 1 == sorted[1] && block / 4 >= threads {
|
||||
let quarter = block / 4;
|
||||
let (slo, shi) = (1usize << sorted[0], 1usize << sorted[1]);
|
||||
// The quarter holding each local index's member.
|
||||
let q_of: Vec<usize> = offsets.iter().map(|&o| [0, slo, shi, slo + shi].iter().position(|&x| x == o).unwrap()).collect();
|
||||
let per = quarter.div_ceil(threads);
|
||||
for (bre, bim) in self.re.chunks_mut(block).zip(self.im.chunks_mut(block)) {
|
||||
let (r01, r23) = bre.split_at_mut(2 * quarter);
|
||||
let (r0, r1) = r01.split_at_mut(quarter);
|
||||
let (r2, r3) = r23.split_at_mut(quarter);
|
||||
let (i01, i23) = bim.split_at_mut(2 * quarter);
|
||||
let (i0, i1) = i01.split_at_mut(quarter);
|
||||
let (i2, i3) = i23.split_at_mut(quarter);
|
||||
std::thread::scope(|sc| {
|
||||
let parts = r0
|
||||
.chunks_mut(per)
|
||||
.zip(r1.chunks_mut(per))
|
||||
.zip(r2.chunks_mut(per))
|
||||
.zip(r3.chunks_mut(per))
|
||||
.zip(i0.chunks_mut(per).zip(i1.chunks_mut(per)).zip(i2.chunks_mut(per)).zip(i3.chunks_mut(per)));
|
||||
for ((((a0, a1), a2), a3), (((b0, b1), b2), b3)) in parts {
|
||||
let q_of = &q_of;
|
||||
sc.spawn(move || {
|
||||
let (re4, im4) = ([a0, a1, a2, a3], [b0, b1, b2, b3]);
|
||||
for i in 0..re4[0].len() {
|
||||
let vr = [re4[q_of[0]][i], re4[q_of[1]][i], re4[q_of[2]][i], re4[q_of[3]][i]];
|
||||
let vi = [im4[q_of[0]][i], im4[q_of[1]][i], im4[q_of[2]][i], im4[q_of[3]][i]];
|
||||
let mut out = [(0.0, 0.0); 4];
|
||||
for (r, o) in out.iter_mut().enumerate() {
|
||||
*o = row_dot(&m[r * 4..r * 4 + 4], &vr, &vi);
|
||||
}
|
||||
for (r, &(a, b)) in out.iter().enumerate() {
|
||||
re4[q_of[r]][i] = a;
|
||||
im4[q_of[r]][i] = b;
|
||||
}
|
||||
}
|
||||
});
|
||||
}
|
||||
});
|
||||
}
|
||||
return;
|
||||
}
|
||||
let half = block / 2;
|
||||
let chunk = if k == 1 { 1usize } else { 1usize << (sorted[k - 2] + 1) };
|
||||
if half / chunk < threads {
|
||||
kernel_block(&mut self.re, &mut self.im, m, &offsets, &sorted);
|
||||
return;
|
||||
}
|
||||
let top_bit = 1usize << top;
|
||||
// For each local index: in the upper half or not, and its offset there.
|
||||
let members: Vec<(bool, usize)> = offsets.iter().map(|&o| (o & top_bit != 0, o & !top_bit)).collect();
|
||||
let rest = &sorted[..k - 1];
|
||||
let per = (half / chunk).div_ceil(threads) * chunk;
|
||||
for (bre, bim) in self.re.chunks_mut(block).zip(self.im.chunks_mut(block)) {
|
||||
let (lre, hre) = bre.split_at_mut(half);
|
||||
let (lim, him) = bim.split_at_mut(half);
|
||||
std::thread::scope(|sc| {
|
||||
for (((a, b), c), e) in lre.chunks_mut(per).zip(lim.chunks_mut(per)).zip(hre.chunks_mut(per)).zip(him.chunks_mut(per)) {
|
||||
let members = &members;
|
||||
sc.spawn(move || kernel_pair(a, b, c, e, m, members, rest));
|
||||
}
|
||||
});
|
||||
}
|
||||
}
|
||||
|
||||
/// `Σ |amplitude|²`, in fixed chunks.
|
||||
pub fn norm2(&self) -> f64 {
|
||||
self.chunked_sum(&|x| self.re[x] * self.re[x] + self.im[x] * self.im[x])
|
||||
}
|
||||
|
||||
/// `Σ_x term(x)`, summed in [`CHUNKS`] fixed chunks (each in order, then
|
||||
/// the chunks in order), threaded where available.
|
||||
fn chunked_sum(&self, term: &(dyn Fn(usize) -> f64 + Sync)) -> f64 {
|
||||
let len = self.re.len();
|
||||
let chunks = len.min(CHUNKS);
|
||||
let per = len / chunks;
|
||||
let sum_chunk = |c: usize| -> f64 {
|
||||
let mut s = 0.0;
|
||||
for x in c * per..(c + 1) * per {
|
||||
s += term(x);
|
||||
}
|
||||
s
|
||||
};
|
||||
#[cfg(not(target_arch = "wasm32"))]
|
||||
let parts: Vec<f64> = {
|
||||
let threads = available_threads().min(chunks);
|
||||
if threads > 1 && len >= 1 << 14 {
|
||||
let size = chunks.div_ceil(threads);
|
||||
std::thread::scope(|sc| {
|
||||
let handles: Vec<_> = (0..chunks)
|
||||
.step_by(size)
|
||||
.map(|lo| {
|
||||
let sum_chunk = &sum_chunk;
|
||||
sc.spawn(move || (lo..(lo + size).min(chunks)).map(sum_chunk).collect::<Vec<f64>>())
|
||||
})
|
||||
.collect();
|
||||
handles.into_iter().flat_map(|h| h.join().expect("a chunk")).collect()
|
||||
})
|
||||
} else {
|
||||
(0..chunks).map(sum_chunk).collect()
|
||||
}
|
||||
};
|
||||
#[cfg(target_arch = "wasm32")]
|
||||
let parts: Vec<f64> = (0..chunks).map(sum_chunk).collect();
|
||||
parts.into_iter().fold(0.0, |a, b| a + b)
|
||||
}
|
||||
|
||||
/// `⟨ψ| Σ c·P |ψ⟩` for Pauli strings given as `(qubit, 'X'|'Y'|'Z')`.
|
||||
pub fn expectation(&self, observable: &[(f64, Vec<(u32, char)>)]) -> f64 {
|
||||
let n = self.n;
|
||||
let mut total = 0.0;
|
||||
for (coef, string) in observable {
|
||||
// P|x⟩ = phase(x) |x ⊕ flip⟩.
|
||||
let mut flip = 0usize;
|
||||
let mut zmask = 0usize;
|
||||
let mut ys = 0u32;
|
||||
for &(q, p) in string {
|
||||
let bit = 1usize << (n - 1 - q);
|
||||
match p {
|
||||
'X' => flip |= bit,
|
||||
'Y' => {
|
||||
flip |= bit;
|
||||
zmask |= bit;
|
||||
ys += 1;
|
||||
}
|
||||
'Z' => zmask |= bit,
|
||||
_ => {}
|
||||
}
|
||||
}
|
||||
// ⟨ψ|P|ψ⟩ = Σ_x conj(ψ[x⊕flip]) · i^{ys} (−1)^{popcount(x & zmask)} ψ[x].
|
||||
let v = self.chunked_sum(&|x| {
|
||||
let y = x ^ flip;
|
||||
let sign = if (x & zmask).count_ones().is_multiple_of(2) { 1.0 } else { -1.0 };
|
||||
// conj(ψ_y) ψ_x
|
||||
let (ar, ai) = (self.re[y], -self.im[y]);
|
||||
let (br, bi) = (self.re[x], self.im[x]);
|
||||
let (pr, pi) = (ar * br - ai * bi, ar * bi + ai * br);
|
||||
// times i^{ys}: only the real part survives in the sum for Hermitian P.
|
||||
let re = match ys % 4 {
|
||||
0 => pr,
|
||||
1 => -pi,
|
||||
2 => -pr,
|
||||
_ => pi,
|
||||
};
|
||||
sign * re
|
||||
});
|
||||
total += coef * v;
|
||||
}
|
||||
total
|
||||
}
|
||||
|
||||
/// `shots` measurement outcomes in the computational basis, from a seeded
|
||||
/// stream: sorted uniforms swept once through the cumulative distribution.
|
||||
pub fn sample(&self, shots: usize, seed: u64) -> Vec<u64> {
|
||||
let mut s = seed;
|
||||
let mut u: Vec<(f64, usize)> = (0..shots)
|
||||
.map(|i| {
|
||||
s = s.wrapping_add(0x9e37_79b9_7f4a_7c15);
|
||||
let mut z = s;
|
||||
z = (z ^ (z >> 30)).wrapping_mul(0xbf58_476d_1ce4_e5b9);
|
||||
z = (z ^ (z >> 27)).wrapping_mul(0x94d0_49bb_1331_11eb);
|
||||
((((z ^ (z >> 31)) >> 11) as f64 + 0.5) / 9_007_199_254_740_992.0, i)
|
||||
})
|
||||
.collect();
|
||||
u.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap().then(a.1.cmp(&b.1)));
|
||||
let total = self.norm2();
|
||||
let mut out = vec![0u64; shots];
|
||||
let mut cum = 0.0;
|
||||
let mut k = 0;
|
||||
let last = self.re.len() - 1;
|
||||
for x in 0..self.re.len() {
|
||||
cum += (self.re[x] * self.re[x] + self.im[x] * self.im[x]) / total;
|
||||
while k < shots && (u[k].0 < cum || x == last) {
|
||||
out[u[k].1] = x as u64;
|
||||
k += 1;
|
||||
}
|
||||
if k == shots {
|
||||
break;
|
||||
}
|
||||
}
|
||||
out
|
||||
}
|
||||
}
|
||||
|
||||
/// Fuse gates into fewer, denser ones:
|
||||
///
|
||||
/// 1. each qubit's pending one-qubit gates are absorbed into the next
|
||||
/// multi-qubit gate on that qubit;
|
||||
/// 2. one-qubit gates after a qubit's last multi-qubit gate are absorbed into
|
||||
/// that gate;
|
||||
/// 3. runs of one-qubit gates on a qubit that has no multi-qubit gate become one
|
||||
/// gate;
|
||||
/// 4. with `width > 2`, adjacent multi-qubit gates are merged while their union
|
||||
/// spans at most `width` qubits.
|
||||
///
|
||||
/// Every move only passes a gate across gates on other qubits, with which it
|
||||
/// commutes, so the fused circuit is the same unitary.
|
||||
pub fn fuse(gates: &[Gate], width: usize) -> Vec<Gate> {
|
||||
let nq = gates.iter().flat_map(|g| g.qubits.iter()).map(|&q| q as usize + 1).max().unwrap_or(0);
|
||||
let mut pending: Vec<Option<Gate>> = vec![None; nq];
|
||||
let mut out: Vec<Gate> = Vec::new();
|
||||
let mut last_multi: Vec<Option<usize>> = vec![None; nq];
|
||||
let compose = |later: &Gate, earlier: &Gate| -> Gate {
|
||||
let mut union = earlier.qubits.clone();
|
||||
for &q in &later.qubits {
|
||||
if !union.contains(&q) {
|
||||
union.push(q);
|
||||
}
|
||||
}
|
||||
let d = 1usize << union.len();
|
||||
Gate { qubits: union.clone(), matrix: matmul(&embed(later, &union), &embed(earlier, &union), d) }
|
||||
};
|
||||
for g in gates {
|
||||
if g.qubits.len() == 1 {
|
||||
let q = g.qubits[0] as usize;
|
||||
pending[q] = Some(match pending[q].take() {
|
||||
Some(p) => compose(g, &p),
|
||||
None => g.clone(),
|
||||
});
|
||||
continue;
|
||||
}
|
||||
// Absorb the pending one-qubit gates on this gate's qubits.
|
||||
let mut cur = g.clone();
|
||||
for &q in &g.qubits {
|
||||
if let Some(p) = pending[q as usize].take() {
|
||||
let embedded = Gate { qubits: cur.qubits.clone(), matrix: embed(&p, &cur.qubits) };
|
||||
let d = 1usize << cur.qubits.len();
|
||||
cur.matrix = matmul(&cur.matrix, &embedded.matrix, d);
|
||||
}
|
||||
}
|
||||
for &q in &cur.qubits {
|
||||
last_multi[q as usize] = Some(out.len());
|
||||
}
|
||||
out.push(cur);
|
||||
}
|
||||
// Trailing one-qubit gates join the last multi-qubit gate on their qubit,
|
||||
// or stand alone.
|
||||
let mut alone = Vec::new();
|
||||
for (q, p) in pending.into_iter().enumerate() {
|
||||
let Some(p) = p else { continue };
|
||||
match last_multi[q] {
|
||||
Some(i) => {
|
||||
let host = &out[i];
|
||||
let embedded = embed(&p, &host.qubits);
|
||||
let d = 1usize << host.qubits.len();
|
||||
let m = matmul(&embedded, &host.matrix, d);
|
||||
out[i].matrix = m;
|
||||
}
|
||||
None => alone.push(p),
|
||||
}
|
||||
}
|
||||
out.extend(alone);
|
||||
if width <= 2 {
|
||||
return out;
|
||||
}
|
||||
// Merge adjacent multi-qubit gates within the width.
|
||||
let mut merged: Vec<Gate> = Vec::new();
|
||||
for g in out {
|
||||
if let Some(prev) = merged.last() {
|
||||
let mut union = prev.qubits.clone();
|
||||
for &q in &g.qubits {
|
||||
if !union.contains(&q) {
|
||||
union.push(q);
|
||||
}
|
||||
}
|
||||
if union.len() <= width {
|
||||
let prev = merged.pop().unwrap();
|
||||
merged.push(compose(&g, &prev));
|
||||
continue;
|
||||
}
|
||||
}
|
||||
merged.push(g);
|
||||
}
|
||||
merged
|
||||
}
|
||||
|
||||
fn matmul(a: &[C], b: &[C], d: usize) -> Vec<C> {
|
||||
let mut out = vec![C { re: 0.0, im: 0.0 }; d * d];
|
||||
for i in 0..d {
|
||||
for k in 0..d {
|
||||
let x = a[i * d + k];
|
||||
if x.re == 0.0 && x.im == 0.0 {
|
||||
continue;
|
||||
}
|
||||
for j in 0..d {
|
||||
let y = b[k * d + j];
|
||||
let o = &mut out[i * d + j];
|
||||
o.re += x.re * y.re - x.im * y.im;
|
||||
o.im += x.re * y.im + x.im * y.re;
|
||||
}
|
||||
}
|
||||
}
|
||||
out
|
||||
}
|
||||
|
||||
/// A gate's matrix on the larger qubit list `onto` (identity elsewhere).
|
||||
fn embed(g: &Gate, onto: &[u32]) -> Vec<C> {
|
||||
let (k, m) = (g.qubits.len(), onto.len());
|
||||
let (dk, dm) = (1usize << k, 1usize << m);
|
||||
// Position of each gate qubit within `onto` (as a bit of the big index).
|
||||
let pos: Vec<usize> = g.qubits.iter().map(|q| m - 1 - onto.iter().position(|o| o == q).unwrap()).collect();
|
||||
let rest_mask: usize = (0..m).filter(|b| !pos.contains(b)).fold(0, |a, b| a | (1 << b));
|
||||
let mut out = vec![C { re: 0.0, im: 0.0 }; dm * dm];
|
||||
for r in 0..dm {
|
||||
for c in 0..dm {
|
||||
if (r & rest_mask) != (c & rest_mask) {
|
||||
continue;
|
||||
}
|
||||
let lr = pos.iter().enumerate().fold(0, |a, (i, &p)| a | (((r >> p) & 1) << (k - 1 - i)));
|
||||
let lc = pos.iter().enumerate().fold(0, |a, (i, &p)| a | (((c >> p) & 1) << (k - 1 - i)));
|
||||
out[r * dm + c] = g.matrix[lr * dk + lc];
|
||||
}
|
||||
}
|
||||
out
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
fn rng(seed: u64) -> impl FnMut() -> f64 {
|
||||
let mut s = seed;
|
||||
move || {
|
||||
s = s.wrapping_add(0x9e37_79b9_7f4a_7c15);
|
||||
let mut z = s;
|
||||
z = (z ^ (z >> 30)).wrapping_mul(0xbf58_476d_1ce4_e5b9);
|
||||
z = (z ^ (z >> 27)).wrapping_mul(0x94d0_49bb_1331_11eb);
|
||||
(((z ^ (z >> 31)) >> 11) as f64 + 0.5) / 9_007_199_254_740_992.0
|
||||
}
|
||||
}
|
||||
|
||||
/// A random circuit over the named gates and arbitrary angles.
|
||||
fn random_circuit(n: u32, depth: usize, seed: u64) -> Vec<Gate> {
|
||||
let mut u = rng(seed);
|
||||
let mut g = Vec::new();
|
||||
for layer in 0..depth {
|
||||
for q in 0..n {
|
||||
let a = u() * core::f64::consts::TAU;
|
||||
g.push(match (u() * 7.0) as u32 {
|
||||
0 => Gate::rx(q, a),
|
||||
1 => Gate::ry(q, a),
|
||||
2 => Gate::rz(q, a),
|
||||
3 => Gate::sqrt_w(q),
|
||||
4 => Gate::sqrt_x(q),
|
||||
5 => Gate::t(q),
|
||||
_ => Gate::h(q),
|
||||
});
|
||||
}
|
||||
let mut q = (layer % 2) as u32;
|
||||
while q + 1 < n {
|
||||
g.push(match (u() * 4.0) as u32 {
|
||||
0 => Gate::cx(q, q + 1),
|
||||
1 => Gate::fsim(q + 1, q, u() * 3.0, u() * 3.0),
|
||||
2 => Gate::cz(q, q + 1),
|
||||
_ => Gate::rzz(q, q + 1, u() * 3.0),
|
||||
});
|
||||
q += 2;
|
||||
}
|
||||
let (a, b) = ((u() * n as f64) as u32, (u() * n as f64) as u32);
|
||||
if a != b {
|
||||
g.push(Gate::iswap(a, b));
|
||||
}
|
||||
}
|
||||
g
|
||||
}
|
||||
|
||||
/// Reference: amplitude by amplitude, gate by gate, no kernels.
|
||||
fn reference(n: u32, gates: &[Gate]) -> Vec<C> {
|
||||
let dim = 1usize << n;
|
||||
let mut psi = vec![C { re: 0.0, im: 0.0 }; dim];
|
||||
psi[0] = C { re: 1.0, im: 0.0 };
|
||||
for g in gates {
|
||||
let k = g.qubits.len();
|
||||
let mut out = vec![C { re: 0.0, im: 0.0 }; dim];
|
||||
for (x, a) in psi.iter().enumerate() {
|
||||
let col = g.qubits.iter().fold(0usize, |c, &q| (c << 1) | ((x >> (n - 1 - q)) & 1));
|
||||
for row in 0..1usize << k {
|
||||
let mut y = x;
|
||||
for (i, &q) in g.qubits.iter().enumerate() {
|
||||
let s = n - 1 - q;
|
||||
y = (y & !(1 << s)) | (((row >> (k - 1 - i)) & 1) << s);
|
||||
}
|
||||
let e = g.matrix[row * (1 << k) + col];
|
||||
out[y].re += e.re * a.re - e.im * a.im;
|
||||
out[y].im += e.re * a.im + e.im * a.re;
|
||||
}
|
||||
}
|
||||
psi = out;
|
||||
}
|
||||
psi
|
||||
}
|
||||
|
||||
#[test]
|
||||
#[allow(clippy::needless_range_loop)]
|
||||
fn kernels_and_fusion_match_the_reference() {
|
||||
for seed in 0..5 {
|
||||
let n = 9;
|
||||
let gates = random_circuit(n, 8, seed);
|
||||
let want = reference(n, &gates);
|
||||
for width in [0usize, 2, 3, 4] {
|
||||
let mut sv = StateVector::zero(n).unwrap();
|
||||
sv.run(&gates, width).unwrap();
|
||||
for x in 0..1usize << n {
|
||||
let (a, b) = (sv.amplitude(x), want[x]);
|
||||
assert!((a.re - b.re).abs() < 1e-12 && (a.im - b.im).abs() < 1e-12, "seed {seed} width {width} x {x}");
|
||||
}
|
||||
assert!((sv.norm2() - 1.0).abs() < 1e-12);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Sixteen qubits: big enough that gates on the top qubits leave fewer
|
||||
/// blocks than threads and take the split-half path, for plain one- and
|
||||
/// two-qubit gates and for fused three- and four-qubit ones.
|
||||
#[test]
|
||||
fn the_split_half_path_matches_the_reference() {
|
||||
let n = 16;
|
||||
let gates = random_circuit(n, 5, 9);
|
||||
let want = reference(n, &gates);
|
||||
for width in [0usize, 2, 3, 4] {
|
||||
let mut sv = StateVector::zero(n).unwrap();
|
||||
sv.run(&gates, width).unwrap();
|
||||
let worst = (0..1usize << n).map(|x| (sv.re[x] - want[x].re).abs().max((sv.im[x] - want[x].im).abs())).fold(0.0, f64::max);
|
||||
assert!(worst < 1e-12, "width {width}: {worst}");
|
||||
// Threaded or not, the same bits: run again and compare exactly.
|
||||
let mut again = StateVector::zero(n).unwrap();
|
||||
again.run(&gates, width).unwrap();
|
||||
assert_eq!(sv, again);
|
||||
}
|
||||
}
|
||||
|
||||
/// One thread takes the block kernels everywhere (as on a target without
|
||||
/// threads); many threads take the split paths for the top qubits. The
|
||||
/// amplitudes and reductions must be the same bits either way.
|
||||
#[test]
|
||||
fn threading_never_changes_a_bit() {
|
||||
let n = 16;
|
||||
let gates = random_circuit(n, 6, 13);
|
||||
for width in [0usize, 2, 3, 4] {
|
||||
let run = |threads: Option<usize>| {
|
||||
FORCE_THREADS.with(|f| f.set(threads));
|
||||
let mut sv = StateVector::zero(n).unwrap();
|
||||
sv.run(&gates, width).unwrap();
|
||||
let e = sv.expectation(&[(1.0, vec![(0, 'X'), (9, 'Z')])]);
|
||||
let norm = sv.norm2();
|
||||
FORCE_THREADS.with(|f| f.set(None));
|
||||
(sv, e.to_bits(), norm.to_bits())
|
||||
};
|
||||
let one = run(Some(1));
|
||||
for t in [2usize, 3, 8, 32] {
|
||||
assert_eq!(run(Some(t)), one, "width {width}, {t} threads");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn named_gates_are_unitary_and_correct() {
|
||||
// √X² = X, √Y² = Y, √W² = W, T² = S, S² = Z, fSim(π/2, 0) = iSWAP up to sign.
|
||||
let sq = |g: Gate| {
|
||||
let mut sv = StateVector::zero(1).unwrap();
|
||||
sv.apply(&Gate::h(0)).unwrap();
|
||||
sv.apply(&Gate::rz(0, 0.3)).unwrap();
|
||||
let mut a = sv.clone();
|
||||
a.apply(&g).unwrap();
|
||||
a.apply(&g).unwrap();
|
||||
a
|
||||
};
|
||||
let once = |g: Gate| {
|
||||
let mut sv = StateVector::zero(1).unwrap();
|
||||
sv.apply(&Gate::h(0)).unwrap();
|
||||
sv.apply(&Gate::rz(0, 0.3)).unwrap();
|
||||
sv.apply(&g).unwrap();
|
||||
sv
|
||||
};
|
||||
let same = |a: &StateVector, b: &StateVector| (0..a.re.len()).all(|x| (a.re[x] - b.re[x]).abs() < 1e-15 && (a.im[x] - b.im[x]).abs() < 1e-15);
|
||||
assert!(same(&sq(Gate::sqrt_x(0)), &once(Gate::x(0))));
|
||||
assert!(same(&sq(Gate::sqrt_y(0)), &once(Gate::y(0))));
|
||||
assert!(same(&sq(Gate::t(0)), &once(Gate::s(0))));
|
||||
assert!(same(&sq(Gate::s(0)), &once(Gate::z(0))));
|
||||
// W = (X + Y)/√2.
|
||||
let r = core::f64::consts::FRAC_1_SQRT_2;
|
||||
let w = Gate { qubits: vec![0], matrix: vec![C { re: 0.0, im: 0.0 }, C { re: r, im: -r }, C { re: r, im: r }, C { re: 0.0, im: 0.0 }] };
|
||||
assert!(same(&sq(Gate::sqrt_w(0)), &once(w)));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn expectations_and_samples_are_right_and_reproducible() {
|
||||
// Bell state: ⟨ZZ⟩ = ⟨XX⟩ = 1, ⟨YY⟩ = −1, ⟨Z0⟩ = 0.
|
||||
let mut sv = StateVector::zero(2).unwrap();
|
||||
sv.run(&[Gate::h(0), Gate::cx(0, 1)], 0).unwrap();
|
||||
let e = |s: Vec<(u32, char)>| sv.expectation(&[(1.0, s)]);
|
||||
assert!((e(vec![(0, 'Z'), (1, 'Z')]) - 1.0).abs() < 1e-15);
|
||||
assert!((e(vec![(0, 'X'), (1, 'X')]) - 1.0).abs() < 1e-15);
|
||||
assert!((e(vec![(0, 'Y'), (1, 'Y')]) + 1.0).abs() < 1e-15);
|
||||
assert!(e(vec![(0, 'Z')]).abs() < 1e-15);
|
||||
let s = sv.sample(10_000, 3);
|
||||
assert!(s.iter().all(|&x| x == 0 || x == 3));
|
||||
let ones = s.iter().filter(|&&x| x == 3).count();
|
||||
assert!((ones as f64 / 10_000.0 - 0.5).abs() < 0.03);
|
||||
assert_eq!(s, sv.sample(10_000, 3));
|
||||
// A larger random state: expectation of Z_q equals Σ ±|ψ|².
|
||||
let n = 12;
|
||||
let mut big = StateVector::zero(n).unwrap();
|
||||
big.run(&random_circuit(n, 6, 4), 3).unwrap();
|
||||
for q in [0u32, 5, 11] {
|
||||
let direct: f64 = (0..1usize << n)
|
||||
.map(|x| {
|
||||
let p = big.re[x] * big.re[x] + big.im[x] * big.im[x];
|
||||
if (x >> (n - 1 - q)) & 1 == 0 { p } else { -p }
|
||||
})
|
||||
.sum();
|
||||
assert!((big.expectation(&[(1.0, vec![(q, 'Z')])]) - direct).abs() < 1e-12);
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn bad_input_is_refused() {
|
||||
assert_eq!(StateVector::zero(33).unwrap_err(), SvError::TooLarge);
|
||||
let mut sv = StateVector::zero(2).unwrap();
|
||||
assert_eq!(sv.apply(&Gate::cx(0, 2)).unwrap_err(), SvError::BadGate(0));
|
||||
assert_eq!(sv.run(&[Gate::h(0), Gate { qubits: vec![1, 1], matrix: vec![C { re: 1.0, im: 0.0 }; 16] }], 2).unwrap_err(), SvError::BadGate(1));
|
||||
}
|
||||
}
|
||||
|
|
@ -84,18 +84,7 @@ pub use crate::linalg::C;
|
|||
use crate::repro::ln;
|
||||
use std::collections::BTreeMap;
|
||||
|
||||
/// A complex number for building gates.
|
||||
pub fn complex(re: f64, im: f64) -> C {
|
||||
C { re, im }
|
||||
}
|
||||
|
||||
/// A gate: a unitary on `qubits`, as a `2ᵏ × 2ᵏ` row-major matrix whose rows
|
||||
/// are outputs and columns inputs, with `qubits[0]` the most significant bit.
|
||||
#[derive(Clone, Debug)]
|
||||
pub struct Gate {
|
||||
pub qubits: Vec<u32>,
|
||||
pub matrix: Vec<C>,
|
||||
}
|
||||
pub use crate::gates::{complex, Gate};
|
||||
|
||||
/// A dense tensor: row-major over `inds`, the first index most significant.
|
||||
#[derive(Clone, Debug)]
|
||||
|
|
|
|||
|
|
@ -13,7 +13,9 @@
|
|||
feature = "quantum_mps",
|
||||
feature = "quantum_vml",
|
||||
feature = "quantum_phasor",
|
||||
feature = "quantum_kernel"
|
||||
feature = "quantum_kernel",
|
||||
feature = "quantum_tn",
|
||||
feature = "quantum_sv"
|
||||
))]
|
||||
mod trig;
|
||||
#[cfg(any(
|
||||
|
|
@ -22,7 +24,9 @@ mod trig;
|
|||
feature = "quantum_mps",
|
||||
feature = "quantum_vml",
|
||||
feature = "quantum_phasor",
|
||||
feature = "quantum_kernel"
|
||||
feature = "quantum_kernel",
|
||||
feature = "quantum_tn",
|
||||
feature = "quantum_sv"
|
||||
))]
|
||||
pub use trig::sin_cos;
|
||||
|
||||
|
|
|
|||
|
|
@ -18,7 +18,8 @@
|
|||
feature = "quantum_resource",
|
||||
feature = "quantum_frame",
|
||||
feature = "quantum_tn",
|
||||
feature = "quantum_gbs"
|
||||
feature = "quantum_gbs",
|
||||
feature = "quantum_sv"
|
||||
))]
|
||||
|
||||
use wai_quantum::*;
|
||||
|
|
@ -251,3 +252,45 @@ fn gaussian_boson_sampling() {
|
|||
assert_eq!(digest(&v), 0x039d_e851_06c3_567d);
|
||||
}
|
||||
|
||||
/// A 16-qubit random circuit at fusion widths 0, 2 and 3: every kernel path
|
||||
/// (block-parallel, split by the top bit, split into quarters, general fused,
|
||||
/// diagonal), threaded natively and sequential under WASI, and the chunked
|
||||
/// reductions (norm, expectation, samples).
|
||||
#[test]
|
||||
fn dense_state_vector() {
|
||||
use gates::Gate;
|
||||
use quantum_sv::StateVector;
|
||||
let n = 16u32;
|
||||
let mut g = Vec::new();
|
||||
let mut a = 0.1f64;
|
||||
for layer in 0..6u32 {
|
||||
for q in 0..n {
|
||||
a += 0.731;
|
||||
g.push(match (q + layer) % 5 {
|
||||
0 => Gate::rx(q, a),
|
||||
1 => Gate::ry(q, a),
|
||||
2 => Gate::rz(q, a),
|
||||
3 => Gate::sqrt_w(q),
|
||||
_ => Gate::t(q),
|
||||
});
|
||||
}
|
||||
let mut q = layer % 2;
|
||||
while q + 1 < n {
|
||||
g.push(if q % 3 == 0 { Gate::fsim(q, q + 1, a, 0.5 * a) } else { Gate::cz(q, q + 1) });
|
||||
q += 2;
|
||||
}
|
||||
g.push(Gate::iswap(0, n - 1));
|
||||
}
|
||||
let mut v = Vec::new();
|
||||
for width in [0usize, 2, 3] {
|
||||
let mut sv = StateVector::zero(n).unwrap();
|
||||
sv.run(&g, width).unwrap();
|
||||
v.extend(sv.re.iter().step_by(97).copied());
|
||||
v.extend(sv.im.iter().step_by(89).copied());
|
||||
v.push(sv.norm2());
|
||||
v.push(sv.expectation(&[(0.5, vec![(0, 'Z'), (7, 'X')]), (1.5, vec![(3, 'Y'), (15, 'Y')])]));
|
||||
v.extend(sv.sample(64, width as u64).into_iter().map(|x| x as f64));
|
||||
}
|
||||
assert_eq!(digest(&v), 0xd0df_e437_5b28_903f);
|
||||
}
|
||||
|
||||
|
|
|
|||
62
wai/wai-quantum/tests/tn_vs_sv.rs
Normal file
62
wai/wai-quantum/tests/tn_vs_sv.rs
Normal file
|
|
@ -0,0 +1,62 @@
|
|||
//! Two independent methods, one answer: an amplitude contracted as a tensor
|
||||
//! network equals the state vector's, for every bitstring tried and for a
|
||||
//! batch of open outputs.
|
||||
#![cfg(all(feature = "quantum_tn", feature = "quantum_sv"))]
|
||||
|
||||
use wai_quantum::gates::Gate;
|
||||
use wai_quantum::quantum_sv::StateVector;
|
||||
use wai_quantum::quantum_tn::{contract_plan, find_path, Families, Network, PathOptions};
|
||||
|
||||
fn circuit(n: u32, depth: u32) -> Vec<Gate> {
|
||||
let mut g = Vec::new();
|
||||
let mut a = 0.2f64;
|
||||
for layer in 0..depth {
|
||||
for q in 0..n {
|
||||
a += 0.613;
|
||||
g.push(match (q * 7 + layer) % 4 {
|
||||
0 => Gate::rx(q, a),
|
||||
1 => Gate::sqrt_w(q),
|
||||
2 => Gate::ry(q, a),
|
||||
_ => Gate::sqrt_x(q),
|
||||
});
|
||||
}
|
||||
let mut q = layer % 2;
|
||||
while q + 1 < n {
|
||||
g.push(Gate::fsim(q, q + 1, 1.2 + 0.1 * q as f64, 0.5));
|
||||
q += 2;
|
||||
}
|
||||
g.push(Gate::cz(0, n - 1));
|
||||
}
|
||||
g
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn tensor_network_amplitudes_equal_the_state_vector() {
|
||||
let n = 18;
|
||||
let gates = circuit(n, 8);
|
||||
let mut sv = StateVector::zero(n).unwrap();
|
||||
sv.run(&gates, 2).unwrap();
|
||||
let opts = PathOptions { trials: 16, seed: 1, max_size: Some(4096.0), reconfigure: 6, families: Families::Both };
|
||||
for x in [0usize, 1, 77_777, 131_071, 262_143] {
|
||||
let bits: Vec<u8> = (0..n).map(|q| ((x >> (n - 1 - q)) & 1) as u8).collect();
|
||||
let mut net = Network::amplitude(n, &gates, &bits, &[]).unwrap();
|
||||
net.simplify();
|
||||
let plan = find_path(&net.shape(), &net.dims, &net.open, &opts);
|
||||
let t = contract_plan(&net, &plan).unwrap().data[0];
|
||||
let s = sv.amplitude(x);
|
||||
assert!((t.re - s.re).abs() < 1e-12 && (t.im - s.im).abs() < 1e-12, "x {x}: {t:?} vs {s:?}");
|
||||
}
|
||||
// A batch: qubits 4 and 11 open, the rest fixed to x = 77 777.
|
||||
let x = 77_777usize;
|
||||
let bits: Vec<u8> = (0..n).map(|q| ((x >> (n - 1 - q)) & 1) as u8).collect();
|
||||
let mut net = Network::amplitude(n, &gates, &bits, &[4, 11]).unwrap();
|
||||
net.simplify();
|
||||
let plan = find_path(&net.shape(), &net.dims, &net.open, &opts);
|
||||
let t = contract_plan(&net, &plan).unwrap();
|
||||
for j in 0..4usize {
|
||||
let (b4, b11) = ((j >> 1) & 1, j & 1);
|
||||
let y = (x & !(1 << (n - 1 - 4)) & !(1 << (n - 1 - 11))) | (b4 << (n - 1 - 4)) | (b11 << (n - 1 - 11));
|
||||
let s = sv.amplitude(y);
|
||||
assert!((t.data[j].re - s.re).abs() < 1e-12 && (t.data[j].im - s.im).abs() < 1e-12, "batch {j}");
|
||||
}
|
||||
}
|
||||
2
wai/wai-rs/Cargo.lock
generated
2
wai/wai-rs/Cargo.lock
generated
|
|
@ -3900,7 +3900,7 @@ dependencies = [
|
|||
|
||||
[[package]]
|
||||
name = "wai-quantum"
|
||||
version = "0.3.27"
|
||||
version = "0.3.28"
|
||||
dependencies = [
|
||||
"blake3",
|
||||
"ed25519-dalek",
|
||||
|
|
|
|||
|
|
@ -355,6 +355,8 @@ quantum_frame = ["dep:wai-quantum", "wai-quantum/quantum_frame"]
|
|||
quantum_tn = ["dep:wai-quantum", "wai-quantum/quantum_tn"]
|
||||
# Gaussian boson sampling: hafnians, torontonians, exact probabilities, samples, checks.
|
||||
quantum_gbs = ["dep:wai-quantum", "wai-quantum/quantum_gbs"]
|
||||
# The dense state vector at any angle: fusion, threads, the same bits everywhere.
|
||||
quantum_sv = ["dep:wai-quantum", "wai-quantum/quantum_sv"]
|
||||
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
|
||||
|
|
@ -486,6 +488,7 @@ required-features = [
|
|||
"quantum_frame",
|
||||
"quantum_tn",
|
||||
"quantum_gbs",
|
||||
"quantum_sv",
|
||||
]
|
||||
|
||||
[[bin]]
|
||||
|
|
|
|||
|
|
@ -165,6 +165,8 @@ fn usage() -> ! {
|
|||
wai-quantum gbs --squeezing r1,r2,… [--eta E] [--seed S | --unitary U.json] [--json]\n \
|
||||
wai-quantum gbs … --sample K (exact click patterns)\n \
|
||||
wai-quantum gbs … --check SAMPLES.txt (observed patterns vs exact marginals)\n\n\
|
||||
STATE VECTOR (any gate, any angle):\n \
|
||||
wai-quantum sv --file CIRCUIT.jsonl [-n N] [--fuse W] [--top K | --sample S [--seed X] | --expect \"Z0 Z1; X3\"] [--json]\n\n\
|
||||
TENSOR NETWORKS:\n \
|
||||
wai-quantum tn --file CIRCUIT.jsonl [-n N] [--bits B] [--open q,q] [--trials T] [--seed S] [--cap-log2 L] [--plan-only] [--json]\n\n\
|
||||
Channels: identity, bitflip, dephasing, depolarizing, damping.\n\
|
||||
|
|
@ -258,6 +260,34 @@ fn human_time(ns: f64) -> String {
|
|||
}
|
||||
}
|
||||
|
||||
/// A circuit file: JSON lines, one gate per line, `{"q":[a,b],"m":[re,im, …]}`
|
||||
/// (rows are outputs, the first qubit most significant). Returns the gates
|
||||
/// and the register size they need.
|
||||
fn read_gate_file(path: &str) -> (Vec<wai::gates::Gate>, u32) {
|
||||
use wai::gates::{complex, Gate};
|
||||
let text = std::fs::read_to_string(path).unwrap_or_else(|e| die(&format!("{path}: {e}")));
|
||||
let mut gates = Vec::new();
|
||||
let mut n: u32 = 0;
|
||||
for (i, line) in text.lines().enumerate().filter(|(_, l)| !l.trim().is_empty()) {
|
||||
let v: serde_json::Value = serde_json::from_str(line).unwrap_or_else(|_| die(&format!("line {}: not JSON", i + 1)));
|
||||
let qubits: Vec<u32> = v["q"]
|
||||
.as_array()
|
||||
.unwrap_or_else(|| die(&format!("line {}: missing \"q\"", i + 1)))
|
||||
.iter()
|
||||
.map(|x| x.as_u64().unwrap_or_else(|| die("qubits are integers")) as u32)
|
||||
.collect();
|
||||
let m: Vec<f64> = v["m"]
|
||||
.as_array()
|
||||
.unwrap_or_else(|| die(&format!("line {}: missing \"m\"", i + 1)))
|
||||
.iter()
|
||||
.map(|x| x.as_f64().unwrap_or_else(|| die("matrix entries are numbers")))
|
||||
.collect();
|
||||
n = n.max(qubits.iter().max().map_or(0, |q| q + 1));
|
||||
gates.push(Gate { qubits, matrix: m.chunks(2).map(|c| complex(c[0], *c.get(1).unwrap_or(&0.0))).collect() });
|
||||
}
|
||||
(gates, n)
|
||||
}
|
||||
|
||||
/// A program for `estimate`: a named example, or counts from flags.
|
||||
/// Returns the counts and the error budget the example is stated with.
|
||||
fn estimate_program(args: &[String], name: Option<&String>) -> (LogicalCounts, f64) {
|
||||
|
|
@ -771,28 +801,9 @@ fn main() {
|
|||
// {"q":[0,1],"m":[re,im, …]} (rows are outputs, first qubit most
|
||||
// significant).
|
||||
"tn" => {
|
||||
use wai::quantum_tn::{complex, contract_plan, find_path, Families, Gate, Network, PathOptions};
|
||||
use wai::quantum_tn::{contract_plan, find_path, Families, Network, PathOptions};
|
||||
let path = opt(args, "--file", "").unwrap_or_else(|| die("--file CIRCUIT.jsonl is required"));
|
||||
let text = std::fs::read_to_string(&path).unwrap_or_else(|e| die(&format!("{path}: {e}")));
|
||||
let mut gates = Vec::new();
|
||||
let mut n: u32 = 0;
|
||||
for (i, line) in text.lines().enumerate().filter(|(_, l)| !l.trim().is_empty()) {
|
||||
let v: serde_json::Value = serde_json::from_str(line).unwrap_or_else(|_| die(&format!("line {}: not JSON", i + 1)));
|
||||
let qubits: Vec<u32> = v["q"]
|
||||
.as_array()
|
||||
.unwrap_or_else(|| die(&format!("line {}: missing \"q\"", i + 1)))
|
||||
.iter()
|
||||
.map(|x| x.as_u64().unwrap_or_else(|| die("qubits are integers")) as u32)
|
||||
.collect();
|
||||
let m: Vec<f64> = v["m"]
|
||||
.as_array()
|
||||
.unwrap_or_else(|| die(&format!("line {}: missing \"m\"", i + 1)))
|
||||
.iter()
|
||||
.map(|x| x.as_f64().unwrap_or_else(|| die("matrix entries are numbers")))
|
||||
.collect();
|
||||
n = n.max(qubits.iter().max().map_or(0, |q| q + 1));
|
||||
gates.push(Gate { qubits, matrix: m.chunks(2).map(|c| complex(c[0], *c.get(1).unwrap_or(&0.0))).collect() });
|
||||
}
|
||||
let (gates, n) = read_gate_file(&path);
|
||||
let n: u32 = num(args, "--qubits", "-n", n);
|
||||
let bits: Vec<u8> = match opt(args, "--bits", "") {
|
||||
Some(b) => b.chars().map(|c| if c == '1' { 1 } else if c == '0' { 0 } else { die("--bits is a 0/1 string") }).collect(),
|
||||
|
|
@ -946,6 +957,74 @@ fn main() {
|
|||
}
|
||||
}
|
||||
|
||||
// The dense state vector, any gate at any angle: run a circuit file
|
||||
// (fused to --fuse qubits), then show the likeliest outcomes, draw
|
||||
// samples, or evaluate Pauli observables.
|
||||
"sv" => {
|
||||
use wai::quantum_sv::StateVector;
|
||||
let path = opt(args, "--file", "").unwrap_or_else(|| die("--file CIRCUIT.jsonl is required"));
|
||||
let (gates, n) = read_gate_file(&path);
|
||||
let n: u32 = num(args, "--qubits", "-n", n);
|
||||
let width: usize = num(args, "--fuse", "", 2);
|
||||
let mut sv = StateVector::zero(n).unwrap_or_else(|e| die(&format!("{e:?} (at most 32 qubits)")));
|
||||
sv.run(&gates, width).unwrap_or_else(|e| die(&format!("{e:?}")));
|
||||
let json = flag(args, "--json");
|
||||
if let Some(k) = opt(args, "--sample", "") {
|
||||
let k: usize = k.parse().unwrap_or_else(|_| die("bad --sample count"));
|
||||
for x in sv.sample(k, num(args, "--seed", "", 0)) {
|
||||
println!("{:0width$b}", x, width = n as usize);
|
||||
}
|
||||
return;
|
||||
}
|
||||
if let Some(spec) = opt(args, "--expect", "") {
|
||||
// "Z0 X3; Y1 Y2": terms separated by ';', each a product of Pauli factors.
|
||||
let mut values = Vec::new();
|
||||
for term in spec.split(';').map(str::trim).filter(|t| !t.is_empty()) {
|
||||
let string: Vec<(u32, char)> = term
|
||||
.split_whitespace()
|
||||
.map(|f| {
|
||||
let (p, q) = f.split_at(1);
|
||||
let p = p.chars().next().unwrap().to_ascii_uppercase();
|
||||
if !"XYZ".contains(p) {
|
||||
die("Pauli factors are X, Y or Z followed by a qubit, as in Z0");
|
||||
}
|
||||
(q.parse().unwrap_or_else(|_| die("bad qubit in --expect")), p)
|
||||
})
|
||||
.collect();
|
||||
if string.iter().any(|&(q, _)| q >= n) {
|
||||
die("--expect names a qubit outside the register");
|
||||
}
|
||||
values.push((term.to_string(), sv.expectation(&[(1.0, string)])));
|
||||
}
|
||||
if json {
|
||||
let items: Vec<String> = values.iter().map(|(t, v)| format!("{}:{v}", serde_json::to_string(t).unwrap())).collect();
|
||||
println!("{{{}}}", items.join(","));
|
||||
} else {
|
||||
for (t, v) in values {
|
||||
println!("<{t}> = {v:+.12}");
|
||||
}
|
||||
}
|
||||
return;
|
||||
}
|
||||
let top: usize = num(args, "--top", "", 8);
|
||||
let mut idx: Vec<usize> = (0..sv.re.len()).collect();
|
||||
let p = |x: usize| sv.re[x] * sv.re[x] + sv.im[x] * sv.im[x];
|
||||
idx.sort_by(|&a, &b| p(b).partial_cmp(&p(a)).unwrap().then(a.cmp(&b)));
|
||||
if json {
|
||||
let rows: Vec<String> = idx
|
||||
.iter()
|
||||
.take(top)
|
||||
.map(|&x| format!("{{\"amplitude\":[{},{}],\"bits\":\"{:0w$b}\",\"p\":{}}}", sv.re[x], sv.im[x], x, p(x), w = n as usize))
|
||||
.collect();
|
||||
println!("{{\"gates\":{},\"norm\":{},\"qubits\":{n},\"top\":[{}]}}", gates.len(), sv.norm2(), rows.join(","));
|
||||
} else {
|
||||
println!("{n} qubits, {} gates (fused to {width}), norm {:.15}", gates.len(), sv.norm2());
|
||||
for &x in idx.iter().take(top) {
|
||||
println!("|{:0w$b}> p = {:.10} amp = {:+.10} {:+.10}i", x, p(x), sv.re[x], sv.im[x], w = n as usize);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
"entropy" => {
|
||||
let c = circuit_from(args, algo);
|
||||
let sv = c.simulate().unwrap_or_else(|e| die(&format!("simulate: {e:?}")));
|
||||
|
|
|
|||
|
|
@ -630,6 +630,12 @@ pub use wai_quantum::quantum_tn;
|
|||
|
||||
#[cfg(feature = "quantum_gbs")]
|
||||
pub use wai_quantum::quantum_gbs;
|
||||
|
||||
#[cfg(feature = "quantum_sv")]
|
||||
pub use wai_quantum::quantum_sv;
|
||||
|
||||
#[cfg(any(feature = "quantum_tn", feature = "quantum_sv"))]
|
||||
pub use wai_quantum::gates;
|
||||
#[cfg(feature = "quantum_qir")]
|
||||
pub use wai_quantum::quantum_qir;
|
||||
|
||||
|
|
|
|||
2
wai/wai-web/Cargo.lock
generated
2
wai/wai-web/Cargo.lock
generated
|
|
@ -734,7 +734,7 @@ dependencies = [
|
|||
|
||||
[[package]]
|
||||
name = "wai-quantum"
|
||||
version = "0.3.27"
|
||||
version = "0.3.28"
|
||||
dependencies = [
|
||||
"blake3",
|
||||
"ed25519-dalek",
|
||||
|
|
|
|||
Loading…
Reference in a new issue