diff --git a/wai/CHANGELOG.md b/wai/CHANGELOG.md index 5bed6ab66..2b36338ac 100644 --- a/wai/CHANGELOG.md +++ b/wai/CHANGELOG.md @@ -6,6 +6,36 @@ wire format may change and such changes are recorded here. ## v0.4.0 — current +### wai-quantum 0.6.9: cut-rank amplitudes with the circuit's own global phase + +Additive. Every pinned digest is unchanged. + +- **`PhasedTableau`** (`quantum_tableau`) is a stabilizer tableau with its global phase. + - It tracks one basis state `r` of the support and its exact amplitude + `⟨r|ψ⟩ = 2^{-k/2} e^{iπm/4}` through every gate. + - Paulis, phase gates, CX, CY, CZ and SWAP move `r` and the phase directly. + - For a Hadamard on qubit `a`, the amplitude of `r ⊕ e_a` comes from the stabilizer + whose X part is `e_a`, found by reducing the stabilizers (about `n³/64` word + operations). `√X` is `H S H`, exactly. + - Tested against the dense state vector after every gate of 40 random circuits + (2 to 8 qubits, over 2,000 gates). A mutant of the Hadamard update is caught. +- **`ProductMeasured::with_last`** keeps the graph form, order and steps with other closing + gates. **`ProductMeasured::scale`** multiplies every amplitude by a factor, folded into + one qubit's weights. Tested against fresh builds. +- **The service's `cutrank` method** returns amplitudes with the circuit's own global + phase. + - The factor is `e^{iΘ} a/b`. Here `a` is the tracked amplitude of the Clifford part, + `b` is the graph form's amplitude at the same basis state (closing gates set to the + identity), and `e^{iΘ}` is the phase the lowering to Clifford gates drops. That + phase is read gate by gate from each gate's matrix against its lowering's. + - Past about a minute of tracking (counted from the Hadamards), the result says + `global_phase` is not fixed, as before. + - `auto` now tries `cutrank` for amplitudes too, before `tensor_network`. + - Tests compare amplitudes with the state vector's directly (no phase aligned). The + circuits include phase-carrying lowerings: `rz`, `rx`, `ry`, `u`, `p`, `cp`, `rzz`, + `rxx`, `iswap`, `sx`, `sxdg`, `cy`, `swap`, `crz` and `cu`. Two mutants are caught: + the lowering phase dropped, and the tracked amplitude ignored. + ### wai-quantum 0.6.8: threads only where they pay, and a cut-rank search that recounts only what moved Additive. No result changes a bit: every kernel path gives the same bits at any thread count, diff --git a/wai/tools/quantum-meter/Cargo.lock b/wai/tools/quantum-meter/Cargo.lock index 6eec16195..985e84317 100644 --- a/wai/tools/quantum-meter/Cargo.lock +++ b/wai/tools/quantum-meter/Cargo.lock @@ -1017,7 +1017,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a" [[package]] name = "wai-quantum" -version = "0.6.8" +version = "0.6.9" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/tools/quantum-wgsl-check/Cargo.lock b/wai/tools/quantum-wgsl-check/Cargo.lock index db830d95f..f4b0bcb14 100644 --- a/wai/tools/quantum-wgsl-check/Cargo.lock +++ b/wai/tools/quantum-wgsl-check/Cargo.lock @@ -599,7 +599,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a" [[package]] name = "wai-quantum" -version = "0.6.8" +version = "0.6.9" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-pack/Cargo.lock b/wai/wai-pack/Cargo.lock index df27f67e6..0572fe23e 100644 --- a/wai/wai-pack/Cargo.lock +++ b/wai/wai-pack/Cargo.lock @@ -713,7 +713,7 @@ dependencies = [ [[package]] name = "wai-quantum" -version = "0.6.8" +version = "0.6.9" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-quantum-qpu/Cargo.lock b/wai/wai-quantum-qpu/Cargo.lock index 4e000db03..0b9cf7c2a 100644 --- a/wai/wai-quantum-qpu/Cargo.lock +++ b/wai/wai-quantum-qpu/Cargo.lock @@ -422,7 +422,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a" [[package]] name = "wai-quantum" -version = "0.6.8" +version = "0.6.9" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-quantum-serve/Cargo.lock b/wai/wai-quantum-serve/Cargo.lock index 6a9270c14..f9a3beec5 100644 --- a/wai/wai-quantum-serve/Cargo.lock +++ b/wai/wai-quantum-serve/Cargo.lock @@ -599,7 +599,7 @@ dependencies = [ [[package]] name = "wai-quantum" -version = "0.6.8" +version = "0.6.9" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-quantum-service/Cargo.lock b/wai/wai-quantum-service/Cargo.lock index 2b283d594..50058c247 100644 --- a/wai/wai-quantum-service/Cargo.lock +++ b/wai/wai-quantum-service/Cargo.lock @@ -857,7 +857,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a" [[package]] name = "wai-quantum" -version = "0.6.8" +version = "0.6.9" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-quantum-service/README.md b/wai/wai-quantum-service/README.md index 531274049..839c1a298 100644 --- a/wai/wai-quantum-service/README.md +++ b/wai/wai-quantum-service/README.md @@ -41,7 +41,7 @@ The program may use the standard gates, user `gate` definitions, and the control | `pauli` | any circuit, up to 1024 qubits | expectations | bounded | | `mps` | any circuit, up to 4096 qubits, at bond dimension up to 256; best when entanglement stays low | counts, expectations, amplitudes | exact when nothing is truncated, otherwise heuristic | | `tensor_network` | any circuit | amplitudes | exact | -| `cutrank` | a Clifford circuit, then one layer of one-qubit gates at any angles, up to 4096 qubits, as wide as the state's cut-rank allows | counts, amplitudes | exact (amplitudes up to one shared phase); single when only a single-precision table fits | +| `cutrank` | a Clifford circuit, then one layer of one-qubit gates at any angles, up to 4096 qubits, as wide as the state's cut-rank allows | counts, amplitudes | exact; single when only a single-precision table fits | | `density_matrix` | noisy circuits, up to 16 qubits within a job's memory | counts, expectations | exact | | `trajectories` | noisy circuits, up to 32 qubits | counts, expectations | estimate | @@ -52,7 +52,7 @@ How the methods run: - It writes the Clifford part's state as a graph state with local Cliffords, and sums the amplitude one qubit at a time along an order of the qubits, in a table of `2^cut` entries. Here `cut` is the state's cut-rank, its entanglement in bits, across the cut the sum has reached. - The planner makes one short search for an order with narrower cuts, enough to bound the job's memory and work. The run searches longer, from several seeds on threads of their own, and keeps the cheapest order that is no wider than the planner's. The cost is set by the widest cut, `max_cut` in the result, not by the number of qubits. - Counts are exact samples, drawn qubit by qubit from the table's conditional probabilities. - - Amplitudes are exact up to one phase shared by all of the job's amplitudes, because a stabilizer tableau does not fix the state's global phase; the result's `global_phase` says so. + - Amplitudes carry the circuit's own global phase. A stabilizer tableau does not fix it, so one amplitude of the Clifford part is tracked exactly, gate by gate, and the summed amplitudes are scaled to agree with it, together with the phase the lowering to Clifford gates drops. A circuit with so many Hadamards that the tracking would take more than about a minute reports `global_phase` as not fixed instead. - The table is held in double precision when it fits the job's memory, and in single precision otherwise (accuracy class `single`). - **pauli** reports each value with its truncation `bound`. It holds at most as many Pauli terms as the job's memory allows; a sum that would outgrow it raises its threshold until it fits (`threshold_reached`), and the bound covers everything that dropped. Unless `ladder` is `false`, it also reports a `ladder`: the value and bound at up to three thresholds a decade above the one asked for, so you can see the value settle as the threshold drops. - **stabilizer**, with noise, samples by Pauli frames against one noiseless outcome of the tableau, 64 shots to a machine word. Its counts have the law of independent noisy shots. Its expectation values are estimates with standard errors, from at least 10,000 shots unless `trajectories` sets the number. Amplitude damping is not a Pauli channel and is refused. @@ -78,7 +78,7 @@ curl -s https://wai.transaction.science/quantum-jobs/v1/run -H 'content-type: ap "observables": ["Z2", "X2"], "shots": 1000 }' ``` -`auto` picks the first method that takes the job. With noise, it tries `density_matrix`, then `stabilizer`, then `trajectories`. For counts, it tries `stabilizer`, `exact`, `statevector`, `statevector_single`, `cutrank`, then `mps`. For expectations alone, it tries `stabilizer`, `exact`, `statevector`, `statevector_single`, `pauli`, then `mps`. For amplitudes alone, it tries `exact`, `statevector`, `statevector_single`, `tensor_network`, then `mps`. The job's `why` names the method chosen and the ones that could not take it. +`auto` picks the first method that takes the job. With noise, it tries `density_matrix`, then `stabilizer`, then `trajectories`. For counts, it tries `stabilizer`, `exact`, `statevector`, `statevector_single`, `cutrank`, then `mps`. For expectations alone, it tries `stabilizer`, `exact`, `statevector`, `statevector_single`, `pauli`, then `mps`. For amplitudes alone, it tries `exact`, `statevector`, `statevector_single`, `cutrank`, `tensor_network`, then `mps`. The job's `why` names the method chosen and the ones that could not take it. The accuracy classes: diff --git a/wai/wai-quantum-service/src/catalog.rs b/wai/wai-quantum-service/src/catalog.rs index 4e03e9fc3..5a0a10e2e 100644 --- a/wai/wai-quantum-service/src/catalog.rs +++ b/wai/wai-quantum-service/src/catalog.rs @@ -14,7 +14,7 @@ //! | `pauli` | any, ≤ 1024 qubits | expectations | bounded | //! | `mps` | any, ≤ 4096 qubits, best on low-entanglement circuits | counts, expectations, amplitudes | exact, or heuristic with its fidelity estimate | //! | `tensor_network` | any | amplitudes | exact | -//! | `cutrank` | Clifford, then one layer of one-qubit gates at any angles, ≤ 4096 qubits, as wide as its cut-rank allows | counts, amplitudes | exact (amplitudes up to one shared phase), or single | +//! | `cutrank` | Clifford, then one layer of one-qubit gates at any angles, ≤ 4096 qubits, as wide as its cut-rank allows | counts, amplitudes | exact, or single | //! | `density_matrix` | any, noisy, ≤ 16 qubits in a job's memory | counts, expectations | exact | //! | `trajectories` | any, noisy, ≤ 32 qubits | counts, expectations | estimate | //! @@ -22,7 +22,7 @@ //! for counts, `stabilizer`, `exact`, `statevector`, `statevector_single`, //! `cutrank`, `mps`; for expectations only, `stabilizer`, `exact`, `statevector`, //! `statevector_single`, `pauli`, `mps`; for amplitudes only, `exact`, -//! `statevector`, `statevector_single`, `tensor_network`, `mps`. A method that cannot take the circuit is +//! `statevector`, `statevector_single`, `cutrank`, `tensor_network`, `mps`. A method that cannot take the circuit is //! skipped, and the plan says why the one chosen was chosen. use crate::qasm::{self, Program}; @@ -107,7 +107,7 @@ pub fn describe(job_bytes: u64) -> Value { { "method": "pauli", "accuracy": "bounded", "circuits": "any", "max_qubits": PAULI_MAX_QUBITS, "answers": ["expectations"], "options": { "threshold": "coefficients below it are dropped (default 1e-5); the bound covers what was dropped", "max_weight": "Pauli strings heavier than this are dropped", "ladder": "also report the value at up to three thresholds a decade apart above it, to show it settling (default true)", "max_terms": "the most Pauli terms to hold (default: what the job's memory holds); past it the threshold rises, and the bound covers what drops" } }, { "method": "mps", "accuracy": "exact when nothing is truncated, otherwise heuristic", "circuits": "any; best on low-entanglement circuits", "max_qubits": MPS_MAX_QUBITS, "answers": ["counts", "expectations", "amplitudes"], "options": { "max_bond": "the largest bond dimension kept (default 64, at most 256)", "cutoff": "singular values whose squared share is below it are dropped (default 1e-12)" } }, { "method": "tensor_network", "accuracy": "exact", "circuits": "any", "max_qubits": qasm::MAX_QUBITS, "answers": ["amplitudes"] }, - { "method": "cutrank", "accuracy": "exact (amplitudes up to one phase shared by the job's amplitudes); single when only a single-precision table fits", "circuits": "a Clifford circuit, then one layer of one-qubit gates at any angles: each qubit's gates from its first non-Clifford one on act on that qubit alone; the cost is set by the state's cut-rank along an order of the qubits, not by their number", "max_qubits": CUTRANK_MAX_QUBITS, "answers": ["counts", "amplitudes"] }, + { "method": "cutrank", "accuracy": "exact; single when only a single-precision table fits", "circuits": "a Clifford circuit, then one layer of one-qubit gates at any angles: each qubit's gates from its first non-Clifford one on act on that qubit alone; the cost is set by the state's cut-rank along an order of the qubits, not by their number", "max_qubits": CUTRANK_MAX_QUBITS, "answers": ["counts", "amplitudes"] }, { "method": "density_matrix", "accuracy": "exact", "circuits": "any, with noise", "max_qubits": density_max_qubits(job_bytes), "answers": ["counts", "expectations"] }, { "method": "trajectories", "accuracy": "estimate", "circuits": "any, with noise", "max_qubits": TRAJECTORY_MAX_QUBITS.min(statevector_max_qubits(job_bytes)), "answers": ["counts", "expectations"], "options": { "trajectories": "how many (default: the shots, or 1000 for expectations alone); shots that take the same errors share one state" } }, ]) @@ -291,7 +291,7 @@ pub fn product_measured(r: &Request, thorough: bool) -> Result CUTRANK_MAX_QUBITS { return Err(format!("cutrank holds at most {CUTRANK_MAX_QUBITS} qubits")); } - let (bulk, last) = r.program.product_measured()?; + let (bulk, last, _) = r.program.product_measured()?; let t = Tableau::run(n, &bulk).map_err(|e| e.to_string())?; let pm = ProductMeasured::new(&t, &last).map_err(|e| e.to_string())?; // Recounting every cut of an order takes about k · cut_k · words @@ -331,6 +331,49 @@ pub fn product_measured(r: &Request, thorough: bool) -> Result Result, String> { + use wai_quantum::gates::complex; + use wai_quantum::quantum_tableau::{Clifford, PhasedTableau}; + let n = r.program.n_qubits as u64; + let (bulk, _, theta) = r.program.product_measured()?; + let hadamards = bulk.iter().map(|g| match g { + Clifford::H(_) => 1, + Clifford::SX(_) | Clifford::SXdg(_) => 2, + _ => 0, + }).sum::(); + if hadamards.saturating_mul(n * n * n.div_ceil(64)) > CUTRANK_PHASE_OPS { + return Ok(None); + } + let tracked = PhasedTableau::run(n as u32, &bulk).map_err(|e| e.to_string())?; + let (bits, (are, aim)) = tracked.reference(); + let (o, one) = (complex(0.0, 0.0), complex(1.0, 0.0)); + let bulk_only = pm.with_last(&vec![[one, o, o, one]; n as usize]).map_err(|e| e.to_string())?; + let b = bulk_only.amplitude_in(&bits, cut, precision).map_err(|e| e.to_string())?; + let bb = b.re * b.re + b.im * b.im; + // a / b, then times e^{iΘ}. + let ratio = complex((are * b.re + aim * b.im) / bb, (aim * b.re - are * b.im) / bb); + let k = complex(theta.re * ratio.re - theta.im * ratio.im, theta.re * ratio.im + theta.im * ratio.re); + let off = (k.re.hypot(k.im) - 1.0).abs(); + let tolerance = if precision == wai_quantum::quantum_cutrank::Precision::Double { 1e-9 } else { 1e-5 }; + if off > tolerance { + return Err(format!("internal: the tracked and summed amplitudes differ in modulus by {off:e}")); + } + Ok(Some(k)) +} + /// The memory a `cutrank` job takes: its largest table with a margin, the /// graph form, and room to run. pub fn cutrank_bytes(pm: &wai_quantum::quantum_cutrank::ProductMeasured, precision: wai_quantum::quantum_cutrank::Precision) -> u64 { @@ -602,7 +645,7 @@ pub fn plan(v: &Value, job_bytes: u64) -> Result { } else if !r.observables.is_empty() && r.amplitudes.is_empty() { &["stabilizer", "exact", "statevector", "statevector_single", "pauli", "mps"] } else { - &["exact", "statevector", "statevector_single", "tensor_network", "mps"] + &["exact", "statevector", "statevector_single", "cutrank", "tensor_network", "mps"] }; let mut refusals = Vec::new(); for m in order { diff --git a/wai/wai-quantum-service/src/qasm.rs b/wai/wai-quantum-service/src/qasm.rs index b084e6ea6..b79f1730c 100644 --- a/wai/wai-quantum-service/src/qasm.rs +++ b/wai/wai-quantum-service/src/qasm.rs @@ -1244,6 +1244,73 @@ fn one_qubit(kind: Kind, p: &[f64]) -> [C; 4] { } } +/// A program split for the `cutrank` method: its Clifford gates, each qubit's +/// closing one-qubit matrix, and the global phase the lowering drops. +pub type ProductMeasuredParts = (Vec, Vec<[C; 4]>, C); + +fn cmul(a: C, b: C) -> C { + complex(a.re * b.re - a.im * b.im, a.re * b.im + a.im * b.re) +} + +/// A Clifford gate as a dense matrix on its qubits (the first the most +/// significant), as the stabilizer tableau reads it. +fn clifford_matrix(c: Clifford) -> Matrix { + let (z, o, i) = (complex(0.0, 0.0), complex(1.0, 0.0), complex(0.0, 1.0)); + match c { + Clifford::H(q) => Matrix::h(q), + Clifford::S(q) => Matrix::s(q), + Clifford::Sdg(q) => Matrix::sdg(q), + Clifford::X(q) => Matrix::x(q), + Clifford::Y(q) => Matrix::y(q), + Clifford::Z(q) => Matrix::z(q), + Clifford::SX(q) => Matrix::sqrt_x(q), + Clifford::SXdg(q) => { + let (a, b) = (complex(0.5, -0.5), complex(0.5, 0.5)); + Matrix { qubits: vec![q], matrix: vec![a, b, b, a] } + } + Clifford::CX(c, t) => Matrix::cx(c, t), + Clifford::CY(c, t) => Matrix { qubits: vec![c, t], matrix: vec![o, z, z, z, z, o, z, z, z, z, z, complex(0.0, -1.0), z, z, i, z] }, + Clifford::CZ(a, b) => Matrix::cz(a, b), + Clifford::Swap(a, b) => Matrix::swap(a, b), + } +} + +/// `e^{iθ}` with `op`'s matrix equal to `e^{iθ}` times the product of the +/// Clifford gates `gates` (in order) on the same qubits: the product is built +/// column by column on the op's own qubits, and the phase read off its +/// largest entry. +fn lowering_phase(op: &Op, gates: &[Clifford]) -> C { + let k = op.qubits.len() as u32; + let local = |q: u32| op.qubits.iter().position(|&p| p == q).expect("a lowered gate stays on its op's qubits") as u32; + let relabel = |c: Clifford| -> Matrix { + let mut m = clifford_matrix(c); + m.qubits = m.qubits.iter().map(|&q| local(q)).collect(); + m + }; + let local_gates: Vec = gates.iter().map(|&c| relabel(c)).collect(); + let d = 1usize << k; + let want = op.matrix().matrix; + let mut best = (0.0f64, complex(1.0, 0.0)); + for col in 0..d { + let mut sv = wai_quantum::quantum_sv::StateVector::zero(k).expect("at most three qubits"); + sv.re[0] = 0.0; + sv.re[col] = 1.0; + for g in &local_gates { + sv.apply(g).expect("a lowered gate on the op's qubits"); + } + for row in 0..d { + let (lr, li) = (sv.re[row], sv.im[row]); + let mag = lr * lr + li * li; + if mag > best.0 + 1e-12 { + // want / lowered + let w = want[row * d + col]; + best = (mag, complex((w.re * lr + w.im * li) / mag, (w.im * lr - w.re * li) / mag)); + } + } + } + best.1 +} + impl Op { /// The gate as a matrix on its qubits (the first the most significant). pub fn matrix(&self) -> Matrix { @@ -1737,7 +1804,11 @@ impl Program { /// into one matrix (the identity for a qubit with none). The input of the /// `cutrank` method. Refused when a gate on two or more qubits comes after /// a qubit's first non-Clifford gate, or is not Clifford itself. - pub fn product_measured(&self) -> Result<(Vec, Vec<[C; 4]>), String> { + /// + /// The third member is the global phase the Clifford lowering drops: each + /// gate's matrix is `e^{iθ}` times the product of the Clifford gates it is + /// written as, and this is the product of the `e^{iθ}`. + pub fn product_measured(&self) -> Result { let one = complex(1.0, 0.0); let zero = complex(0.0, 0.0); let mul = |a: &[C], b: &[C; 4]| -> [C; 4] { @@ -1746,6 +1817,7 @@ impl Program { [add(m(a[0], b[0]), m(a[1], b[2])), add(m(a[0], b[1]), m(a[1], b[3])), add(m(a[2], b[0]), m(a[3], b[2])), add(m(a[2], b[1]), m(a[3], b[3]))] }; let mut bulk = Vec::new(); + let mut phase = one; let mut tail: Vec> = vec![None; self.n_qubits as usize]; for op in &self.ops { if let [q] = op.qubits[..] { @@ -1754,7 +1826,10 @@ impl Program { tail[q] = Some(mul(&op.matrix().matrix, &t)); } else { match Self::lower_clifford(std::slice::from_ref(op)) { - Ok((gates, _)) => bulk.extend(gates), + Ok((gates, _)) => { + phase = cmul(phase, lowering_phase(op, &gates)); + bulk.extend(gates); + } Err(_) => tail[q] = Some(mul(&op.matrix().matrix, &[one, zero, zero, one])), } } @@ -1764,9 +1839,12 @@ impl Program { return Err(format!("{:?} on qubits {:?} comes after a gate on qubit {q} that is not Clifford; cutrank takes a Clifford circuit, then one layer of one-qubit gates", op.kind, op.qubits)); } let (gates, _) = Self::lower_clifford(std::slice::from_ref(op)).map_err(|e| format!("{e}; cutrank takes a Clifford circuit, then one layer of one-qubit gates"))?; + phase = cmul(phase, lowering_phase(op, &gates)); bulk.extend(gates); } - Ok((bulk, tail.into_iter().map(|t| t.unwrap_or([one, zero, zero, one])).collect())) + // Held to modulus 1 against the rounding of many products. + let m = phase.re.hypot(phase.im); + Ok((bulk, tail.into_iter().map(|t| t.unwrap_or([one, zero, zero, one])).collect(), complex(phase.re / m, phase.im / m))) } /// The Clifford lowering, and where each op's gates end in it. diff --git a/wai/wai-quantum-service/src/run.rs b/wai/wai-quantum-service/src/run.rs index 3bb8c00d4..d87ce1194 100644 --- a/wai/wai-quantum-service/src/run.rs +++ b/wai/wai-quantum-service/src/run.rs @@ -372,7 +372,7 @@ pub fn run(request: &Value, plan: &Plan) -> Result { "class": if precision == Precision::Double { "exact" } else { "single" }, "max_cut": cut, "order_work": pm.work(), - "note": "the state's sum taken qubit by qubit along an order of the qubits, in a table of 2^max_cut entries (single class: the table stored in single precision, its arithmetic in double); counts are exact samples, seeded; amplitudes are exact up to one phase shared by all of the job's amplitudes, since a stabilizer tableau does not fix the state's global phase" + "note": "the state's sum taken qubit by qubit along an order of the qubits, in a table of 2^max_cut entries (single class: the table stored in single precision, its arithmetic in double); counts are exact samples, seeded; amplitudes carry the circuit's own global phase, found by tracking one amplitude of the Clifford part exactly through the circuit" }); if r.shots > 0 { let mut map: HashMap = HashMap::new(); @@ -382,6 +382,14 @@ pub fn run(request: &Value, plan: &Plan) -> Result { merge(&mut out, counts(map)); } if !r.amplitudes.is_empty() { + // The circuit's global phase: the lowering's own, times the + // ratio of one exactly tracked amplitude of the Clifford part + // to the graph form's at the same basis state. + let mut pm = pm; + let fixed = catalog::cutrank_phase(&r, &pm, cut, precision)?; + if let Some(kappa) = fixed { + pm.scale(kappa); + } let mut entries = Vec::new(); for b in &r.amplitudes { let x: Vec = (0..n).map(|q| b.as_bytes()[(n - 1 - q) as usize] == b'1').collect(); @@ -389,7 +397,9 @@ pub fn run(request: &Value, plan: &Plan) -> Result { entries.push(json!({ "bits": b, "re": a.re, "im": a.im })); } out["amplitudes"] = Value::Array(entries); - out["global_phase"] = json!("not fixed: the amplitudes share one unknown phase; their magnitudes and relative phases are exact"); + if fixed.is_none() { + out["global_phase"] = json!("not fixed: tracking it exactly would take more than one job may; the amplitudes share one unknown phase, and their magnitudes and relative phases are exact"); + } } } "tensor_network" => { @@ -943,35 +953,34 @@ mod tests { q } - #[test] - fn cutrank_amplitudes_are_the_state_vectors_up_to_one_phase() { - for (n, seed) in [(6u32, 1u64), (9, 2), (11, 3)] { - let qasm = clifford_then_layer(n, 3 * n as usize, seed); - let bits: Vec = (0..24u64).map(|k| (0..n).map(|q| if (k.wrapping_mul(0x9e37_79b9) >> q) & 1 == 1 { '1' } else { '0' }).collect()).collect(); - let cr = go(json!({ "qasm": qasm, "method": "cutrank", "amplitudes": bits })); - let sv = go(json!({ "qasm": qasm, "method": "statevector", "amplitudes": bits })); - assert_eq!(cr["method"], "cutrank", "{cr}"); - assert_eq!(cr["accuracy"]["class"], "exact", "{cr}"); - assert!(cr["global_phase"].as_str().is_some_and(|g| g.starts_with("not fixed")), "{cr}"); - let amp = |v: &Value, i: usize| (v["amplitudes"][i]["re"].as_f64().unwrap(), v["amplitudes"][i]["im"].as_f64().unwrap()); - // The shared phase, from the largest amplitude. - let k = (0..bits.len()).max_by(|&a, &b| { - let (x, y) = (amp(&sv, a), amp(&sv, b)); - (x.0 * x.0 + x.1 * x.1).total_cmp(&(y.0 * y.0 + y.1 * y.1)) - }).unwrap(); - let ((cr_k_re, cr_k_im), (sv_k_re, sv_k_im)) = (amp(&cr, k), amp(&sv, k)); - let norm = sv_k_re * sv_k_re + sv_k_im * sv_k_im; - // phase = cr_k / sv_k - let (pr, pi) = ((cr_k_re * sv_k_re + cr_k_im * sv_k_im) / norm, (cr_k_im * sv_k_re - cr_k_re * sv_k_im) / norm); - assert!(((pr * pr + pi * pi) - 1.0).abs() < 1e-9, "n {n}: |phase| {}", (pr * pr + pi * pi).sqrt()); - for (i, b) in bits.iter().enumerate() { - let ((cr_re, cr_im), (sv_re, sv_im)) = (amp(&cr, i), amp(&sv, i)); - let (want_re, want_im) = (pr * sv_re - pi * sv_im, pr * sv_im + pi * sv_re); - assert!((cr_re - want_re).abs() < 1e-10 && (cr_im - want_im).abs() < 1e-10, "n {n}, {b}: cutrank {cr_re} {cr_im}, statevector × phase {want_re} {want_im}"); - } + /// cutrank's amplitudes, global phase and all, are the state vector's. + fn assert_amplitudes_match(qasm: &str, n: u32, label: &str) { + let bits: Vec = (0..24u64).map(|k| (0..n).map(|q| if (k.wrapping_mul(0x9e37_79b9) >> q) & 1 == 1 { '1' } else { '0' }).collect()).collect(); + let cr = go(json!({ "qasm": qasm, "method": "cutrank", "amplitudes": bits })); + let sv = go(json!({ "qasm": qasm, "method": "statevector", "amplitudes": bits })); + assert_eq!(cr["method"], "cutrank", "{label}: {cr}"); + assert_eq!(cr["accuracy"]["class"], "exact", "{label}: {cr}"); + assert!(cr.get("global_phase").is_none(), "{label}: the phase is fixed: {cr}"); + for (i, b) in bits.iter().enumerate() { + let get = |v: &Value| (v["amplitudes"][i]["re"].as_f64().unwrap(), v["amplitudes"][i]["im"].as_f64().unwrap()); + let ((cre, cim), (sre, sim)) = (get(&cr), get(&sv)); + assert!((cre - sre).abs() < 1e-10 && (cim - sim).abs() < 1e-10, "{label}, {b}: cutrank {cre} {cim}, statevector {sre} {sim}"); } } + #[test] + fn cutrank_amplitudes_are_the_state_vectors_phase_and_all() { + for (n, seed) in [(6u32, 1u64), (9, 2), (11, 3)] { + assert_amplitudes_match(&clifford_then_layer(n, 3 * n as usize, seed), n, &format!("n {n}")); + } + // Gates whose lowering to Clifford ones drops a phase, in the Clifford part. + let phased = "OPENQASM 3.0; qubit[5] q; h q; rz(pi/2) q[0]; rx(pi) q[1]; ry(pi/2) q[2]; u(pi/2, 0, pi) q[3]; p(pi/2) q[4]; \ + cp(pi) q[0], q[1]; rzz(pi/2) q[1], q[2]; iswap q[2], q[3]; sx q[4]; sxdg q[0]; cy q[3], q[4]; swap q[0], q[2]; \ + rxx(pi/2) q[3], q[1]; crz(pi) q[4], q[0]; cu(pi, 0, pi, pi/2) q[2], q[1]; \ + t q[0]; ry(0.4) q[1]; rz(1.3) q[2]; h q[3]; rx(0.7) q[3]; tdg q[4];"; + assert_amplitudes_match(phased, 5, "phase-carrying lowerings"); + } + #[test] fn cutrank_counts_follow_the_exact_distribution() { let n = 5u32; diff --git a/wai/wai-quantum/Cargo.lock b/wai/wai-quantum/Cargo.lock index cfe101006..47a7b5bcb 100644 --- a/wai/wai-quantum/Cargo.lock +++ b/wai/wai-quantum/Cargo.lock @@ -1079,7 +1079,7 @@ checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a" [[package]] name = "wai-quantum" -version = "0.6.8" +version = "0.6.9" dependencies = [ "blake3", "block2", diff --git a/wai/wai-quantum/Cargo.toml b/wai/wai-quantum/Cargo.toml index de6a36666..a86503a1b 100644 --- a/wai/wai-quantum/Cargo.toml +++ b/wai/wai-quantum/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "wai-quantum" -version = "0.6.8" +version = "0.6.9" edition = "2024" rust-version = "1.97" license = "Apache-2.0" diff --git a/wai/wai-quantum/README.md b/wai/wai-quantum/README.md index f60808a86..dfba3c2ed 100644 --- a/wai/wai-quantum/README.md +++ b/wai/wai-quantum/README.md @@ -467,6 +467,16 @@ the search takes one qubit to another place, and only the cuts between the two places are recounted, so 300 000 moves on the 70-qubit circuit below take a few seconds. `sample` prints each exact sample with `2^n · p`. +A stabilizer tableau fixes a state only up to a global phase, and the graph +form inherits that. `PhasedTableau` carries the phase: alongside the tableau it +tracks one basis state of the support and its exact amplitude, a power of +`2^{-1/2}` times an eighth root of unity, gate by gate. A Hadamard reads the +stabilizer whose X part is that one qubit, which is a reduction of the +stabilizers. `ProductMeasured::scale` then brings every amplitude to the +circuit's own phase: by the ratio of the tracked amplitude to the graph form's +at the same basis state (`with_last` with identities), times whatever phase +the caller's lowering to Clifford gates dropped. + A 70-qubit graph-state sampling circuit has a nearest-neighbour Clifford bulk, then a non-Clifford basis on every qubit. With an annealed order, every cut is 29 bits or fewer. On one workstation, one amplitude takes about half a second diff --git a/wai/wai-quantum/src/quantum_cutrank.rs b/wai/wai-quantum/src/quantum_cutrank.rs index 8627392f8..7a2c651fd 100644 --- a/wai/wai-quantum/src/quantum_cutrank.rs +++ b/wai/wai-quantum/src/quantum_cutrank.rs @@ -67,6 +67,12 @@ fn mul(a: &Unitary, b: &Unitary) -> Unitary { [add(m(a[0], b[0]), m(a[1], b[2])), add(m(a[0], b[1]), m(a[1], b[3])), add(m(a[2], b[0]), m(a[3], b[2])), add(m(a[2], b[1]), m(a[3], b[3]))] } +/// `W_q / √2 = V_q L_q / √2` for each qubit. +fn weights(local: &[Unitary], last: &[Unitary]) -> Vec { + let r = std::f64::consts::FRAC_1_SQRT_2; + local.iter().zip(last).map(|(l, v)| mul(v, l).map(|c| complex(c.re * r, c.im * r))).collect() +} + fn get(v: &[u64], i: usize) -> bool { v[i / 64] >> (i % 64) & 1 == 1 } @@ -78,6 +84,8 @@ pub struct ProductMeasured { adjacency: Vec>, /// `W_q / √2`. w: Vec, + /// The graph form's local Cliffords `L_q`, so the closing gates can change. + local: Vec, order: Vec, /// The cut-rank after each qubit of the order. cuts: Vec, @@ -105,7 +113,7 @@ impl ProductMeasured { let h = [complex(r, 0.0), complex(r, 0.0), complex(r, 0.0), complex(-r, 0.0)]; let s = [one, o, o, complex(0.0, 1.0)]; let z = [one, o, o, complex(-1.0, 0.0)]; - let w = (0..n) + let local: Vec = (0..n) .map(|q| { let l = form.local[q]; let mut m = [one, o, o, one]; @@ -118,16 +126,37 @@ impl ProductMeasured { if l.z { m = mul(&m, &z); } - let wq = mul(&last[q], &m); - wq.map(|c| complex(c.re * r, c.im * r)) + m }) .collect(); - let mut pm = ProductMeasured { n, adjacency: form.adjacency.clone(), w, order: (0..n).collect(), cuts: Vec::new(), steps: Vec::new() }; + let w = weights(&local, last); + let mut pm = ProductMeasured { n, adjacency: form.adjacency.clone(), w, local, order: (0..n).collect(), cuts: Vec::new(), steps: Vec::new() }; pm.cuts = pm.cut_ranks_of(&pm.order); pm.steps = plan(&pm.adjacency, n, &pm.order); pm } + /// The same graph form, order and steps with other closing gates: + /// `last[q]` on qubit q. + pub fn with_last(&self, last: &[Unitary]) -> Result { + if last.len() != self.n { + return Err(CutRankError::Gates { qubits: self.n, gates: last.len() }); + } + let mut pm = self.clone(); + pm.w = weights(&self.local, last); + Ok(pm) + } + + /// Every amplitude times `factor`. Each term of the sum holds exactly one + /// entry of each qubit's weights, so it is folded into one qubit's; a + /// factor of modulus 1 (a phase) leaves every probability as it was. + pub fn scale(&mut self, factor: C) { + let q = self.order.first().copied().unwrap_or(0); + if let Some(wq) = self.w.get_mut(q) { + *wq = wq.map(|c| complex(c.re * factor.re - c.im * factor.im, c.re * factor.im + c.im * factor.re)); + } + } + /// The same state, summed in `order`. pub fn with_order(mut self, order: &[usize]) -> Result { let mut seen = vec![false; self.n]; @@ -1164,6 +1193,34 @@ mod tests { assert!((a.re - b.re).hypot(a.im - b.im) < 1e-9 * a.re.hypot(a.im).max(1e-300), "{a:?} vs {b:?}"); } + /// Other closing gates on the same graph form and order give the same + /// amplitudes as a fresh build; a phase scales every amplitude by itself. + #[test] + fn closing_gates_change_and_a_phase_scales() { + let n = 10usize; + let mut s = 11u64; + let ops = random_clifford(n as u32, 6 * n, 23); + let t = Tableau::run(n as u32, &ops).unwrap(); + let first: Vec = (0..n).map(|_| random_unitary(&mut s)).collect(); + let second: Vec = (0..n).map(|_| random_unitary(&mut s)).collect(); + let pm = ProductMeasured::new(&t, &first).unwrap(); + let order = pm.anneal_order(200, 1); + let pm = pm.with_order(&order).unwrap(); + let swapped = pm.with_last(&second).unwrap(); + let fresh = ProductMeasured::new(&t, &second).unwrap().with_order(&order).unwrap(); + let phase = complex(0.6, 0.8); + let mut scaled = fresh.clone(); + scaled.scale(phase); + for k in 0..20 { + let x: Vec = (0..n).map(|q| (k * 7 + q * 3) % 5 < 2).collect(); + let (a, b, c) = (swapped.amplitude(&x, n).unwrap(), fresh.amplitude(&x, n).unwrap(), scaled.amplitude(&x, n).unwrap()); + assert!((a.re - b.re).abs() < 1e-14 && (a.im - b.im).abs() < 1e-14, "{a:?} vs {b:?}"); + let want = complex(b.re * phase.re - b.im * phase.im, b.re * phase.im + b.im * phase.re); + assert!((c.re - want.re).abs() < 1e-14 && (c.im - want.im).abs() < 1e-14, "{c:?} vs {want:?}"); + } + assert!(matches!(pm.with_last(&second[..3]), Err(CutRankError::Gates { .. }))); + } + #[test] fn cut_ranks_are_the_entanglement() { // A GHZ chain: one bit across every cut. A product state: none. diff --git a/wai/wai-quantum/src/quantum_tableau.rs b/wai/wai-quantum/src/quantum_tableau.rs index 3a8f8c4ac..13cfe8e7a 100644 --- a/wai/wai-quantum/src/quantum_tableau.rs +++ b/wai/wai-quantum/src/quantum_tableau.rs @@ -691,6 +691,211 @@ impl Tableau { } } +impl Tableau { + /// The stabilizer whose X part is exactly `e_a`, if the group holds one: + /// its Z part and its sign, as a row of letters (`Y` where both are set). + /// Found by reducing the stabilizers' X parts: in reduced form a vector + /// lies in their span only through the row whose pivot it sets. + fn stabilizer_with_x(&self, a: usize) -> Option<(Vec, bool)> { + let n = self.n; + let mut rows = self.rows(n, 2 * n); + let mut k = 0; + let mut pivot_row = None; + for q in 0..n { + let Some(p) = (k..n).find(|&i| get(rows.xr(i), q)) else { continue }; + rows.swap(k, p); + for i in 0..n { + if i != k && get(rows.xr(i), q) { + rows.mul_into(i, k); + } + } + if q == a { + pivot_row = Some(k); + } + k += 1; + } + let j = pivot_row?; + let x = rows.xr(j); + let only_a = x.iter().enumerate().all(|(w, &bits)| bits == if w == a / 64 { 1u64 << (a % 64) } else { 0 }); + only_a.then(|| (rows.zr(j).to_vec(), rows.r[j])) + } +} + +/// A stabilizer state with its global phase: the [`Tableau`], which fixes the +/// state only up to a phase, and one basis state `r` of its support with the +/// exact amplitude `⟨r|ψ⟩ = 2^{−k/2} e^{iπm/4}`. Every gate but `H` moves `r` +/// and the phase directly. For `H` on qubit `a` the amplitude of `r ⊕ e_a` is +/// needed: it is `c · ⟨r|ψ⟩` with `c` read from the stabilizer whose X part is +/// `e_a` (`c = 0` when there is none), which costs a reduction of the +/// stabilizers, about `n³/64` word operations. +#[derive(Clone, Debug, PartialEq, Eq)] +pub struct PhasedTableau { + t: Tableau, + r: Vec, + /// `|⟨r|ψ⟩| = 2^{−k/2}`. + k: u32, + /// The phase, in eighths of a turn. + m: u8, +} + +impl PhasedTableau { + /// `|0…0⟩`, with amplitude 1 on `|0…0⟩`. + pub fn new(n: u32) -> Result { + Ok(PhasedTableau { t: Tableau::new(n)?, r: vec![0; words(n as usize)], k: 0, m: 0 }) + } + + /// The circuit `ops` on `n` qubits from `|0…0⟩`. + pub fn run(n: u32, ops: &[Clifford]) -> Result { + let mut t = PhasedTableau::new(n)?; + for &op in ops { + t.apply(op)?; + } + Ok(t) + } + + /// The tableau, without the phase. + pub fn tableau(&self) -> &Tableau { + &self.t + } + + /// The reference basis state (qubit `q` at index `q`) and `⟨r|ψ⟩`. + pub fn reference(&self) -> (Vec, (f64, f64)) { + let bits = (0..self.t.n).map(|q| get(&self.r, q)).collect(); + let mag = (-(self.k as f64) / 2.0).exp2(); + // Exact for the eight angles, so the entries are the same bits everywhere. + let (c, s) = match self.m % 8 { + 0 => (1.0, 0.0), + 2 => (0.0, 1.0), + 4 => (-1.0, 0.0), + 6 => (0.0, -1.0), + 1 => (std::f64::consts::FRAC_1_SQRT_2, std::f64::consts::FRAC_1_SQRT_2), + 3 => (-std::f64::consts::FRAC_1_SQRT_2, std::f64::consts::FRAC_1_SQRT_2), + 5 => (-std::f64::consts::FRAC_1_SQRT_2, -std::f64::consts::FRAC_1_SQRT_2), + _ => (std::f64::consts::FRAC_1_SQRT_2, -std::f64::consts::FRAC_1_SQRT_2), + }; + (bits, (mag * c, mag * s)) + } + + fn bit(&self, q: usize) -> bool { + get(&self.r, q) + } + + fn flip(&mut self, q: usize) { + self.r[q / 64] ^= 1 << (q % 64); + } + + fn turn(&mut self, eighths: u8) { + self.m = (self.m + eighths) % 8; + } + + /// `H` on qubit `a`, the reference first (it needs the state before). + fn h(&mut self, a: usize) { + // c = ⟨r ⊕ e_a|ψ⟩ / ⟨r|ψ⟩ in quarter turns, or None when r ⊕ e_a is + // outside the support. With P = (−1)^s i^y X_a Z^z a stabilizer, + // P|r⟩ = (−1)^{s + z·r} i^y |r ⊕ e_a⟩. + let c = self.t.stabilizer_with_x(a).map(|(z, s)| { + let y = get(&z, a); + let parity = z.iter().zip(&self.r).map(|(a, b)| (a & b).count_ones()).sum::() % 2 == 1; + (u8::from(s ^ parity) * 2 + u8::from(y)) % 4 + }); + let b = self.bit(a); + // ⟨r|Hψ⟩ / ⟨r|ψ⟩ is (1 + c)/√2 for b = 0 and (c − 1)/√2 for b = 1; when + // that is zero, r ⊕ e_a takes over with (c + 1)/√2 or (1 − c)/√2. + match (c, b) { + (None, false) => self.k += 1, + (None, true) => { + self.k += 1; + self.turn(4); + } + (Some(0), false) => self.k -= 1, + (Some(0), true) => { + self.k -= 1; + self.flip(a); + } + (Some(2), false) => { + self.k -= 1; + self.flip(a); + } + (Some(2), true) => { + self.k -= 1; + self.turn(4); + } + (Some(1), false) => self.turn(1), + (Some(1), true) => self.turn(3), + (Some(_), false) => self.turn(7), + (Some(_), true) => self.turn(5), + } + self.t.h(a); + } + + /// Applies one gate. + pub fn apply(&mut self, op: Clifford) -> Result<(), TableauError> { + match op { + Clifford::H(a) => { + let a = self.t.check(a)?; + self.h(a); + return Ok(()); + } + Clifford::SX(a) | Clifford::SXdg(a) => { + // √X = H S H and √X† = H S† H, exactly. + let a = self.t.check(a)?; + self.h(a); + self.apply(if matches!(op, Clifford::SX(_)) { Clifford::S(a as u32) } else { Clifford::Sdg(a as u32) })?; + self.h(a); + return Ok(()); + } + Clifford::X(a) => { + let a = self.t.check(a)?; + self.flip(a); + } + Clifford::Y(a) => { + // Y|b⟩ = i(−1)^b |1 − b⟩. + let a = self.t.check(a)?; + self.turn(2 + 4 * u8::from(self.bit(a))); + self.flip(a); + } + Clifford::Z(a) => { + let a = self.t.check(a)?; + self.turn(4 * u8::from(self.bit(a))); + } + Clifford::S(a) => { + let a = self.t.check(a)?; + self.turn(2 * u8::from(self.bit(a))); + } + Clifford::Sdg(a) => { + let a = self.t.check(a)?; + self.turn(6 * u8::from(self.bit(a))); + } + Clifford::CX(a, b) => { + let (a, b) = self.t.pair(a, b)?; + if self.bit(a) { + self.flip(b); + } + } + Clifford::CY(a, b) => { + let (a, b) = self.t.pair(a, b)?; + if self.bit(a) { + self.turn(2 + 4 * u8::from(self.bit(b))); + self.flip(b); + } + } + Clifford::CZ(a, b) => { + let (a, b) = self.t.pair(a, b)?; + self.turn(4 * u8::from(self.bit(a) && self.bit(b))); + } + Clifford::Swap(a, b) => { + let (a, b) = self.t.pair(a, b)?; + let (x, y) = (self.bit(a), self.bit(b)); + if x != y { + self.flip(a); + self.flip(b); + } + } + } + self.t.apply(op) + } +} + /// The computational-basis outcomes of a stabilizer state: uniform over /// `x₀ + span(basis)`. Bit `q` of an outcome is qubit `q`. #[derive(Clone, Debug, PartialEq, Eq)] @@ -1593,3 +1798,81 @@ mod frame_tests { assert!(e[0].0 > 0.8 && e[0].0 < 1.0, "{e:?}"); } } + +#[cfg(all(test, feature = "quantum_sv"))] +mod phase_tests { + use super::*; + use crate::gates::{complex, Gate}; + use crate::quantum_sv::StateVector; + + fn gate(c: Clifford) -> Gate { + let (z, o, i) = (complex(0.0, 0.0), complex(1.0, 0.0), complex(0.0, 1.0)); + match c { + Clifford::H(q) => Gate::h(q), + Clifford::S(q) => Gate::s(q), + Clifford::Sdg(q) => Gate::sdg(q), + Clifford::X(q) => Gate::x(q), + Clifford::Y(q) => Gate::y(q), + Clifford::Z(q) => Gate::z(q), + Clifford::SX(q) => Gate::sqrt_x(q), + Clifford::SXdg(q) => { + let (a, b) = (complex(0.5, -0.5), complex(0.5, 0.5)); + Gate { qubits: vec![q], matrix: vec![a, b, b, a] } + } + Clifford::CX(c, t) => Gate::cx(c, t), + Clifford::CY(c, t) => Gate { qubits: vec![c, t], matrix: vec![o, z, z, z, z, o, z, z, z, z, z, complex(0.0, -1.0), z, z, i, z] }, + Clifford::CZ(a, b) => Gate::cz(a, b), + Clifford::Swap(a, b) => Gate::swap(a, b), + } + } + + fn circuit(n: u32, len: usize, seed: u64) -> Vec { + let mut s = seed; + (0..len) + .map(|_| { + let a = (splitmix(&mut s) % u64::from(n)) as u32; + let b = (a + 1 + (splitmix(&mut s) % u64::from(n - 1)) as u32) % n; + // Hadamards often, so every case of the reference's update is met. + match splitmix(&mut s) % 14 { + 0..=2 => Clifford::H(a), + 3 => Clifford::S(a), + 4 => Clifford::Sdg(a), + 5 => Clifford::SX(a), + 6 => Clifford::SXdg(a), + 7 => Clifford::X(a), + 8 => Clifford::Y(a), + 9 => Clifford::Z(a), + 10 => Clifford::CX(a, b), + 11 => Clifford::CY(a, b), + 12 => Clifford::CZ(a, b), + _ => Clifford::Swap(a, b), + } + }) + .collect() + } + + /// The reference amplitude is the state vector's, phase and all, after + /// every gate of random Clifford circuits. + #[test] + fn the_reference_amplitude_is_the_state_vectors() { + let mut checked = 0; + for seed in 0..40u64 { + let n = 2 + (seed % 7) as u32; + let ops = circuit(n, 12 * n as usize, seed); + let mut t = PhasedTableau::new(n).unwrap(); + let mut sv = StateVector::zero(n).unwrap(); + for (g, &op) in ops.iter().enumerate() { + t.apply(op).unwrap(); + sv.apply(&gate(op)).unwrap(); + let (bits, (re, im)) = t.reference(); + let index = bits.iter().enumerate().fold(0usize, |x, (q, &b)| x | (usize::from(b) << (n as usize - 1 - q))); + let a = sv.amplitude(index); + assert!((a.re - re).abs() < 1e-12 && (a.im - im).abs() < 1e-12, "seed {seed}, gate {g} ({op:?}): tableau {re} {im}, state vector {} {}", a.re, a.im); + assert!(re.hypot(im) > 1e-9, "the reference stays in the support"); + checked += 1; + } + assert_eq!(t.tableau(), &Tableau::run(n, &ops).unwrap(), "the tableau is the plain one"); + } + assert!(checked > 2000); + } +} diff --git a/wai/wai-rs/Cargo.lock b/wai/wai-rs/Cargo.lock index ae1093f44..890a68780 100644 --- a/wai/wai-rs/Cargo.lock +++ b/wai/wai-rs/Cargo.lock @@ -4310,7 +4310,7 @@ version = "0.1.0" [[package]] name = "wai-quantum" -version = "0.6.8" +version = "0.6.9" dependencies = [ "blake3", "ed25519-dalek", diff --git a/wai/wai-web/Cargo.lock b/wai/wai-web/Cargo.lock index 21070fae1..2bc7cea66 100644 --- a/wai/wai-web/Cargo.lock +++ b/wai/wai-web/Cargo.lock @@ -1241,7 +1241,7 @@ version = "0.1.0" [[package]] name = "wai-quantum" -version = "0.6.8" +version = "0.6.9" dependencies = [ "blake3", "ed25519-dalek",