e546a87a2f
QuantumReceipt gains gateset_of, which names the gateset under which a circuit's ops hash to the receipt's circuit_hash, and gateset_note, which turns that into the text to show when a reconstruction fails for that reason. verify_reconstruction stays strict. from_json refuses an n_qubits that does not fit a u8 and an n_ops that does not fit a u32 instead of truncating them. The fabric names its gate kernel and random corpus once per gateset: wai.quantum.circuit2/apply and wai.quantum.circuit2/random-v1 join the circuit/1 names. certify_lanes and lanes::random_circuit take a gateset and draw the same ops from a seed in both. detect_local offers the reference kernel of each gateset. The lanes32 corpus digests for seed 7 and 200 cases are pinned: d18454c6 for circuit/1, as 0.3.41 certified it, and 3dbde276 for circuit2. noisy_expect_z and route_line build their circuits with empty_like, so a derived circuit keeps its source's gateset. from_qasm refuses a gate that names one qubit twice as an invalid circuit when it reads it, in every gateset. quantum_spd's controlled_phase refuses a phase on 32 or more qubits instead of shifting out of range. The QIR and QASM export docs say that a circuit/1 circuit exports ideal angles, not its own numerics. The bell and receipt examples print the capability and handle the fallible calls. Co-Authored-By: Claude Opus 5.5 <[email protected]>
1408 lines
55 KiB
Rust
1408 lines
55 KiB
Rust
//! Sparse Pauli dynamics at utility scale (`wai.quantum.spd`).
|
||
//!
|
||
//! Computes `⟨ψ₀| U† O U |ψ₀⟩` for circuits of up to 1024 qubits by evolving the
|
||
//! *observable* in the Heisenberg picture as a sum of Pauli strings,
|
||
//! `O = Σ c_P P`. A Clifford gate maps each string to one signed string. A Pauli
|
||
//! rotation `exp(−iθP/2)` leaves strings that commute with `P` alone and splits
|
||
//! the rest in two: `Q ↦ cos θ·Q + sin θ·(−iQP)`. Strings whose coefficient falls
|
||
//! below a threshold are dropped. A circuit therefore costs what its
|
||
//! non-Clifford rotations and the observable's light cone make it cost, not
|
||
//! `2ⁿ`. This is the method that computed the 127-qubit kicked-Ising
|
||
//! observables of the heavy-hex "utility" experiment on a laptop.
|
||
//!
|
||
//! **Arbitrary angles.** The dyadic circuit core is byte-exact because every
|
||
//! phase it allows is a power-of-two root of unity. This module takes any real
|
||
//! angle, so its contract is the one the crate's other `f64` layers keep:
|
||
//!
|
||
//! - **Reproducible.** Every result is bit-identical across runs, machines and
|
||
//! targets. Terms live in an insertion-ordered table, so the order in which
|
||
//! coefficients are added never depends on hashing or on the platform's word
|
||
//! size. `sin` and `cos` are computed in-crate from `+ − × ÷`, because the
|
||
//! platform's libm may differ in the last place.
|
||
//! - **Bounded.** Every dropped coefficient is added to
|
||
//! [`SpdResult::truncation_bound`]. A dropped term `δ·P` would have evolved
|
||
//! under the remaining gates into `δ` times a unitary conjugate of `P`, whose
|
||
//! expectation is at most `|δ|` in magnitude. So `|exact − value|` is at most
|
||
//! the bound, up to floating-point rounding. Clifford gates and exact
|
||
//! quarter-turn rotations drop nothing.
|
||
|
||
use std::fmt;
|
||
|
||
pub use crate::repro::sin_cos;
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Circuits
|
||
// ---------------------------------------------------------------------------
|
||
|
||
/// The largest register this module takes: 16 words of 64 qubits.
|
||
pub const MAX_QUBITS: u32 = 1024;
|
||
/// The largest angle magnitude a rotation may carry.
|
||
pub const MAX_ANGLE: f64 = 1_048_576.0;
|
||
|
||
/// A single-qubit Pauli factor.
|
||
#[derive(Clone, Copy, Debug, PartialEq, Eq, PartialOrd, Ord, Hash)]
|
||
pub enum Pauli {
|
||
X,
|
||
Y,
|
||
Z,
|
||
}
|
||
|
||
/// A Pauli string as `(qubit, factor)` pairs; qubits not named carry `I`.
|
||
pub type PauliString = Vec<(u32, Pauli)>;
|
||
|
||
/// One instruction.
|
||
#[derive(Clone, Debug, PartialEq)]
|
||
pub enum Op {
|
||
H(u32),
|
||
S(u32),
|
||
Sdg(u32),
|
||
/// √X.
|
||
SX(u32),
|
||
/// √X†.
|
||
SXdg(u32),
|
||
X(u32),
|
||
Y(u32),
|
||
Z(u32),
|
||
/// `CX(control, target)`.
|
||
CX(u32, u32),
|
||
CZ(u32, u32),
|
||
Swap(u32, u32),
|
||
/// `exp(−iθP/2)` for any Pauli string `P`: `RX`, `RY`, `RZ` for one factor,
|
||
/// `RZZ` and the like for two.
|
||
Rot { axis: PauliString, theta: f64 },
|
||
/// `exp(−i(kπ/2)P/2)` for `k ∈ {1, 2, 3}`: a Pauli rotation through a
|
||
/// multiple of π/2. These are Clifford, and applied exactly.
|
||
QuarterRot { axis: PauliString, k: u8 },
|
||
}
|
||
|
||
/// A circuit on `n` qubits.
|
||
#[derive(Clone, Debug, PartialEq)]
|
||
pub struct RotCircuit {
|
||
pub n: u32,
|
||
pub ops: Vec<Op>,
|
||
}
|
||
|
||
/// Why a circuit, observable or initial state was refused.
|
||
#[derive(Clone, Debug, PartialEq)]
|
||
pub enum SpdError {
|
||
TooManyQubits(u32),
|
||
QubitOutOfRange { qubit: u32, n: u32 },
|
||
/// A two-qubit gate or Pauli string naming one qubit twice.
|
||
RepeatedQubit(u32),
|
||
EmptyAxis,
|
||
BadAngle(f64),
|
||
BadQuarter(u8),
|
||
/// A dyadic gate with no Pauli-rotation form here (a controlled `H`).
|
||
Unsupported(String),
|
||
/// A receipt's meter returned without running the job.
|
||
NotRun,
|
||
}
|
||
|
||
impl fmt::Display for SpdError {
|
||
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
|
||
match self {
|
||
SpdError::TooManyQubits(n) => write!(f, "{n} qubits; at most {MAX_QUBITS}"),
|
||
SpdError::QubitOutOfRange { qubit, n } => write!(f, "qubit {qubit} of a {n}-qubit register"),
|
||
SpdError::RepeatedQubit(q) => write!(f, "qubit {q} named twice"),
|
||
SpdError::EmptyAxis => write!(f, "a rotation about the identity"),
|
||
SpdError::BadAngle(a) => write!(f, "angle {a} is not finite or exceeds {MAX_ANGLE}"),
|
||
SpdError::BadQuarter(k) => write!(f, "quarter turns {k}; expected 1, 2 or 3"),
|
||
SpdError::Unsupported(g) => write!(f, "unsupported gate: {g}"),
|
||
SpdError::NotRun => write!(f, "the meter did not run the job"),
|
||
}
|
||
}
|
||
}
|
||
|
||
impl std::error::Error for SpdError {}
|
||
|
||
impl RotCircuit {
|
||
pub fn new(n: u32) -> Self {
|
||
RotCircuit { n, ops: Vec::new() }
|
||
}
|
||
|
||
pub fn push(&mut self, op: Op) -> &mut Self {
|
||
self.ops.push(op);
|
||
self
|
||
}
|
||
|
||
pub fn rx(&mut self, q: u32, theta: f64) -> &mut Self {
|
||
self.push(Op::Rot { axis: vec![(q, Pauli::X)], theta })
|
||
}
|
||
|
||
pub fn ry(&mut self, q: u32, theta: f64) -> &mut Self {
|
||
self.push(Op::Rot { axis: vec![(q, Pauli::Y)], theta })
|
||
}
|
||
|
||
pub fn rz(&mut self, q: u32, theta: f64) -> &mut Self {
|
||
self.push(Op::Rot { axis: vec![(q, Pauli::Z)], theta })
|
||
}
|
||
|
||
/// `exp(−iθ Z⊗Z/2)`.
|
||
pub fn rzz(&mut self, a: u32, b: u32, theta: f64) -> &mut Self {
|
||
self.push(Op::Rot { axis: vec![(a, Pauli::Z), (b, Pauli::Z)], theta })
|
||
}
|
||
|
||
/// `exp(−iθ X⊗X/2)`.
|
||
pub fn rxx(&mut self, a: u32, b: u32, theta: f64) -> &mut Self {
|
||
self.push(Op::Rot { axis: vec![(a, Pauli::X), (b, Pauli::X)], theta })
|
||
}
|
||
|
||
/// Check every qubit index, axis and angle.
|
||
pub fn validate(&self) -> Result<(), SpdError> {
|
||
if self.n == 0 || self.n > MAX_QUBITS {
|
||
return Err(SpdError::TooManyQubits(self.n));
|
||
}
|
||
let q_ok = |q: u32| {
|
||
if q < self.n {
|
||
Ok(())
|
||
} else {
|
||
Err(SpdError::QubitOutOfRange { qubit: q, n: self.n })
|
||
}
|
||
};
|
||
let pair = |a: u32, b: u32| {
|
||
q_ok(a)?;
|
||
q_ok(b)?;
|
||
if a == b { Err(SpdError::RepeatedQubit(a)) } else { Ok(()) }
|
||
};
|
||
for op in &self.ops {
|
||
match op {
|
||
Op::H(q) | Op::S(q) | Op::Sdg(q) | Op::SX(q) | Op::SXdg(q) | Op::X(q) | Op::Y(q) | Op::Z(q) => {
|
||
q_ok(*q)?
|
||
}
|
||
Op::CX(a, b) | Op::CZ(a, b) | Op::Swap(a, b) => pair(*a, *b)?,
|
||
Op::Rot { axis, theta } => {
|
||
check_string(axis, self.n)?;
|
||
if axis.is_empty() {
|
||
return Err(SpdError::EmptyAxis);
|
||
}
|
||
if !theta.is_finite() || theta.abs() > MAX_ANGLE {
|
||
return Err(SpdError::BadAngle(*theta));
|
||
}
|
||
}
|
||
Op::QuarterRot { axis, k } => {
|
||
check_string(axis, self.n)?;
|
||
if axis.is_empty() {
|
||
return Err(SpdError::EmptyAxis);
|
||
}
|
||
if !(1..=3).contains(k) {
|
||
return Err(SpdError::BadQuarter(*k));
|
||
}
|
||
}
|
||
}
|
||
}
|
||
Ok(())
|
||
}
|
||
}
|
||
|
||
fn check_string(s: &[(u32, Pauli)], n: u32) -> Result<(), SpdError> {
|
||
for (i, (q, _)) in s.iter().enumerate() {
|
||
if *q >= n {
|
||
return Err(SpdError::QubitOutOfRange { qubit: *q, n });
|
||
}
|
||
if s[..i].iter().any(|(r, _)| r == q) {
|
||
return Err(SpdError::RepeatedQubit(*q));
|
||
}
|
||
}
|
||
Ok(())
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Pauli strings as bit words
|
||
// ---------------------------------------------------------------------------
|
||
|
||
/// A Hermitian Pauli string on up to `64·W` qubits: factor `X^x Z^z` at each
|
||
/// qubit, with `Y` (both bits set) meaning `i·XZ`, so every string is Hermitian.
|
||
#[derive(Clone, Copy, PartialEq, Eq)]
|
||
struct Key<const W: usize> {
|
||
x: [u64; W],
|
||
z: [u64; W],
|
||
}
|
||
|
||
impl<const W: usize> Key<W> {
|
||
const ZERO: Self = Key { x: [0; W], z: [0; W] };
|
||
|
||
fn from_string(s: &[(u32, Pauli)]) -> Self {
|
||
let mut k = Self::ZERO;
|
||
for &(q, p) in s {
|
||
let (w, b) = ((q / 64) as usize, 1u64 << (q % 64));
|
||
match p {
|
||
Pauli::X => k.x[w] |= b,
|
||
Pauli::Z => k.z[w] |= b,
|
||
Pauli::Y => {
|
||
k.x[w] |= b;
|
||
k.z[w] |= b;
|
||
}
|
||
}
|
||
}
|
||
k
|
||
}
|
||
|
||
#[inline]
|
||
fn get(&self, q: u32) -> (bool, bool) {
|
||
let (w, b) = ((q / 64) as usize, q % 64);
|
||
((self.x[w] >> b) & 1 == 1, (self.z[w] >> b) & 1 == 1)
|
||
}
|
||
|
||
#[inline]
|
||
fn set(&mut self, q: u32, x: bool, z: bool) {
|
||
let (w, b) = ((q / 64) as usize, 1u64 << (q % 64));
|
||
self.x[w] = (self.x[w] & !b) | if x { b } else { 0 };
|
||
self.z[w] = (self.z[w] & !b) | if z { b } else { 0 };
|
||
}
|
||
|
||
/// Whether two strings commute: the symplectic product is even.
|
||
#[inline]
|
||
fn commutes(&self, o: &Self) -> bool {
|
||
let mut p = 0u32;
|
||
for w in 0..W {
|
||
p += (self.x[w] & o.z[w]).count_ones() + (self.z[w] & o.x[w]).count_ones();
|
||
}
|
||
p & 1 == 0
|
||
}
|
||
|
||
/// `self · o = i^e · (self ⊕ o)`: the product string and the exponent `e`
|
||
/// (mod 4) of its phase.
|
||
#[inline]
|
||
fn mul(&self, o: &Self) -> (Self, u32) {
|
||
let mut out = Self::ZERO;
|
||
let mut e = 0u32;
|
||
for w in 0..W {
|
||
let (x, z) = (self.x[w] ^ o.x[w], self.z[w] ^ o.z[w]);
|
||
out.x[w] = x;
|
||
out.z[w] = z;
|
||
e = e.wrapping_add((self.x[w] & self.z[w]).count_ones());
|
||
e = e.wrapping_add((o.x[w] & o.z[w]).count_ones());
|
||
e = e.wrapping_add(2 * (self.z[w] & o.x[w]).count_ones());
|
||
e = e.wrapping_sub((x & z).count_ones());
|
||
}
|
||
(out, e & 3)
|
||
}
|
||
|
||
fn weight(&self) -> u32 {
|
||
(0..W).map(|w| (self.x[w] | self.z[w]).count_ones()).sum()
|
||
}
|
||
|
||
fn hash(&self) -> u64 {
|
||
let mut h = 0x9e37_79b9_7f4a_7c15u64;
|
||
for w in 0..W {
|
||
h = mix(h ^ self.x[w]);
|
||
h = mix(h ^ self.z[w]);
|
||
}
|
||
h
|
||
}
|
||
}
|
||
|
||
#[inline]
|
||
fn mix(mut z: u64) -> u64 {
|
||
z = (z ^ (z >> 30)).wrapping_mul(0xbf58_476d_1ce4_e5b9);
|
||
z = (z ^ (z >> 27)).wrapping_mul(0x94d0_49bb_1331_11eb);
|
||
z ^ (z >> 31)
|
||
}
|
||
|
||
/// Pauli terms in insertion order, with a hash index on the side. Iteration
|
||
/// follows the order terms were first inserted, so the order coefficients are
|
||
/// added in — and with it every bit of the result — is the same on every
|
||
/// platform.
|
||
///
|
||
/// A dropped term keeps its slot with coefficient `0.0`, so dropping costs
|
||
/// nothing to the index, and a term that reappears is revived in place. The
|
||
/// table is compacted, in order, when dead slots outnumber live ones, and
|
||
/// before a Clifford gate rewrites every key anyway. When it happens depends
|
||
/// only on the computation, so the result stays the same everywhere.
|
||
struct Terms<const W: usize> {
|
||
keys: Vec<Key<W>>,
|
||
coef: Vec<f64>,
|
||
index: Vec<u32>,
|
||
stale: bool,
|
||
/// Slots whose coefficient is `0.0`.
|
||
dead: usize,
|
||
}
|
||
|
||
const EMPTY: u32 = u32::MAX;
|
||
|
||
impl<const W: usize> Terms<W> {
|
||
fn new() -> Self {
|
||
Terms { keys: Vec::new(), coef: Vec::new(), index: vec![EMPTY; 16], stale: false, dead: 0 }
|
||
}
|
||
|
||
/// Terms with a nonzero coefficient.
|
||
fn live(&self) -> usize {
|
||
self.keys.len() - self.dead
|
||
}
|
||
|
||
/// Remove dead slots, keeping the order of the rest.
|
||
fn compact(&mut self) {
|
||
let mut keep = 0;
|
||
for i in 0..self.keys.len() {
|
||
if self.coef[i] != 0.0 {
|
||
self.keys[keep] = self.keys[i];
|
||
self.coef[keep] = self.coef[i];
|
||
keep += 1;
|
||
}
|
||
}
|
||
self.keys.truncate(keep);
|
||
self.coef.truncate(keep);
|
||
self.dead = 0;
|
||
self.stale = true;
|
||
}
|
||
|
||
fn reindex(&mut self) {
|
||
let mut cap = 16usize;
|
||
while cap < 2 * self.keys.len() + 2 {
|
||
cap *= 2;
|
||
}
|
||
self.index.clear();
|
||
self.index.resize(cap, EMPTY);
|
||
let mask = (cap - 1) as u64;
|
||
for (i, k) in self.keys.iter().enumerate() {
|
||
let mut slot = (k.hash() & mask) as usize;
|
||
while self.index[slot] != EMPTY {
|
||
slot = (slot + 1) & (cap - 1);
|
||
}
|
||
self.index[slot] = i as u32;
|
||
}
|
||
self.stale = false;
|
||
}
|
||
|
||
/// Add `c` to the coefficient of `k`, inserting it at the end if absent.
|
||
fn add(&mut self, k: Key<W>, c: f64) {
|
||
debug_assert!(!self.stale);
|
||
if 2 * (self.keys.len() + 1) > self.index.len() {
|
||
self.reindex();
|
||
}
|
||
let cap = self.index.len();
|
||
let mut slot = (k.hash() & (cap as u64 - 1)) as usize;
|
||
loop {
|
||
let i = self.index[slot];
|
||
if i == EMPTY {
|
||
self.index[slot] = self.keys.len() as u32;
|
||
self.keys.push(k);
|
||
self.coef.push(c);
|
||
if c == 0.0 {
|
||
self.dead += 1;
|
||
}
|
||
return;
|
||
}
|
||
if self.keys[i as usize] == k {
|
||
let before = self.coef[i as usize];
|
||
let after = before + c;
|
||
self.coef[i as usize] = after;
|
||
match (before == 0.0, after == 0.0) {
|
||
(true, false) => self.dead -= 1,
|
||
(false, true) => self.dead += 1,
|
||
_ => {}
|
||
}
|
||
return;
|
||
}
|
||
slot = (slot + 1) & (cap - 1);
|
||
}
|
||
}
|
||
|
||
/// Drop terms with `|c| < threshold`, or heavier than `max_weight`.
|
||
/// Returns the sum of dropped `|c|`.
|
||
fn truncate(&mut self, threshold: f64, max_weight: Option<u32>) -> f64 {
|
||
let mut dropped = 0.0;
|
||
for i in 0..self.keys.len() {
|
||
let c = self.coef[i];
|
||
if c == 0.0 {
|
||
continue;
|
||
}
|
||
if c.abs() < threshold || max_weight.is_some_and(|m| self.keys[i].weight() > m) {
|
||
dropped += c.abs();
|
||
self.coef[i] = 0.0;
|
||
self.dead += 1;
|
||
}
|
||
}
|
||
if 2 * self.dead > self.keys.len() {
|
||
self.compact();
|
||
}
|
||
dropped
|
||
}
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Conjugation by one gate
|
||
// ---------------------------------------------------------------------------
|
||
|
||
/// `G† (X^x Z^z) G` for a single-qubit Clifford `G` on one factor: the new
|
||
/// factor and whether the sign flips.
|
||
fn conj1(op: &Op, x: bool, z: bool) -> (bool, bool, bool) {
|
||
// factors: I (0,0), X (1,0), Z (0,1), Y (1,1)
|
||
match (op, x, z) {
|
||
(_, false, false) => (false, false, false),
|
||
(Op::H(_), true, false) => (false, true, false), // X → Z
|
||
(Op::H(_), false, true) => (true, false, false), // Z → X
|
||
(Op::H(_), true, true) => (true, true, true), // Y → −Y
|
||
(Op::S(_), true, false) => (true, true, true), // X → −Y
|
||
(Op::S(_), true, true) => (true, false, false), // Y → X
|
||
(Op::Sdg(_), true, false) => (true, true, false), // X → Y
|
||
(Op::Sdg(_), true, true) => (true, false, true), // Y → −X
|
||
(Op::SX(_), false, true) => (true, true, false), // Z → Y
|
||
(Op::SX(_), true, true) => (false, true, true), // Y → −Z
|
||
(Op::SXdg(_), false, true) => (true, true, true), // Z → −Y
|
||
(Op::SXdg(_), true, true) => (false, true, false), // Y → Z
|
||
(Op::X(_), _, true) => (x, z, true), // Z, Y → −
|
||
(Op::Y(_), true, false) | (Op::Y(_), false, true) => (x, z, true), // X, Z → −
|
||
(Op::Z(_), true, _) => (x, z, true), // X, Y → −
|
||
_ => (x, z, false),
|
||
}
|
||
}
|
||
|
||
/// Apply a Clifford gate to every term in place. Conjugation by a unitary is a
|
||
/// bijection on Pauli strings, so no two terms collide.
|
||
fn clifford<const W: usize>(t: &mut Terms<W>, op: &Op) {
|
||
if t.dead > 0 {
|
||
t.compact();
|
||
}
|
||
for i in 0..t.keys.len() {
|
||
let k = &mut t.keys[i];
|
||
let neg = match *op {
|
||
Op::H(q) | Op::S(q) | Op::Sdg(q) | Op::SX(q) | Op::SXdg(q) | Op::X(q) | Op::Y(q) | Op::Z(q) => {
|
||
let (x, z) = k.get(q);
|
||
let (nx, nz, neg) = conj1(op, x, z);
|
||
k.set(q, nx, nz);
|
||
neg
|
||
}
|
||
Op::CX(c, tq) => {
|
||
let ((xc, zc), (xt, zt)) = (k.get(c), k.get(tq));
|
||
k.set(tq, xt ^ xc, zt);
|
||
k.set(c, xc, zc ^ zt);
|
||
xc && zt && !(xt ^ zc)
|
||
}
|
||
Op::CZ(a, b) => {
|
||
let ((xa, za), (xb, zb)) = (k.get(a), k.get(b));
|
||
k.set(a, xa, za ^ xb);
|
||
k.set(b, xb, zb ^ xa);
|
||
xa && xb && (za ^ zb)
|
||
}
|
||
Op::Swap(a, b) => {
|
||
let (fa, fb) = (k.get(a), k.get(b));
|
||
k.set(a, fb.0, fb.1);
|
||
k.set(b, fa.0, fa.1);
|
||
false
|
||
}
|
||
_ => unreachable!("not a fixed Clifford"),
|
||
};
|
||
if neg {
|
||
t.coef[i] = -t.coef[i];
|
||
}
|
||
}
|
||
t.stale = true;
|
||
}
|
||
|
||
/// Conjugate by `exp(−i(kπ/2)P/2)`, exactly. A string anticommuting with `P`
|
||
/// becomes `±(−iQP)` for `k = 1` or `3`, and `−Q` for `k = 2`.
|
||
fn quarter<const W: usize>(t: &mut Terms<W>, p: &Key<W>, k: u8) {
|
||
if t.dead > 0 {
|
||
t.compact();
|
||
}
|
||
for i in 0..t.keys.len() {
|
||
let q = t.keys[i];
|
||
if q.commutes(p) {
|
||
continue;
|
||
}
|
||
if k == 2 {
|
||
t.coef[i] = -t.coef[i];
|
||
continue;
|
||
}
|
||
let (r, e) = q.mul(p);
|
||
// −i·QP = i^{e+3}·R; e is odd here, so the factor is ±1.
|
||
let mut neg = (e + 3) & 3 == 2;
|
||
if k == 3 {
|
||
neg = !neg;
|
||
}
|
||
t.keys[i] = r;
|
||
if neg {
|
||
t.coef[i] = -t.coef[i];
|
||
}
|
||
}
|
||
t.stale = true;
|
||
}
|
||
|
||
/// Conjugate by `exp(−iθP/2)`: `Q ↦ cos θ·Q + sin θ·(−iQP)` for every `Q` that
|
||
/// anticommutes with `P`. The new terms are merged after every old one has been
|
||
/// scaled, so each old coefficient is read before it changes.
|
||
fn rotate<const W: usize>(t: &mut Terms<W>, p: &Key<W>, sin: f64, cos: f64, pending: &mut Vec<(Key<W>, f64)>) {
|
||
if t.stale {
|
||
t.reindex();
|
||
}
|
||
pending.clear();
|
||
for i in 0..t.keys.len() {
|
||
let c = t.coef[i];
|
||
if c == 0.0 {
|
||
continue;
|
||
}
|
||
let q = t.keys[i];
|
||
if q.commutes(p) {
|
||
continue;
|
||
}
|
||
let (r, e) = q.mul(p);
|
||
let v = c * sin;
|
||
pending.push((r, if (e + 3) & 3 == 2 { -v } else { v }));
|
||
t.coef[i] = c * cos;
|
||
}
|
||
for &(r, v) in pending.iter() {
|
||
t.add(r, v);
|
||
}
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// Propagation
|
||
// ---------------------------------------------------------------------------
|
||
|
||
/// How hard to truncate.
|
||
#[derive(Clone, Copy, Debug, PartialEq)]
|
||
pub struct SpdConfig {
|
||
/// Drop terms whose coefficient magnitude falls below this after a rotation.
|
||
/// `0.0` drops only exact zeros: the computation is then exact up to
|
||
/// rounding, and its cost may be exponential.
|
||
pub threshold: f64,
|
||
/// Also drop strings acting on more than this many qubits.
|
||
pub max_weight: Option<u32>,
|
||
}
|
||
|
||
impl Default for SpdConfig {
|
||
fn default() -> Self {
|
||
SpdConfig { threshold: 1e-4, max_weight: None }
|
||
}
|
||
}
|
||
|
||
/// An expectation value, what it cost, and how far truncation can have moved it.
|
||
#[derive(Clone, Copy, Debug, PartialEq)]
|
||
pub struct SpdResult {
|
||
/// `⟨ψ₀| U† O U |ψ₀⟩` as computed.
|
||
pub value: f64,
|
||
/// The sum of every dropped coefficient's magnitude: a bound on
|
||
/// `|exact − value|`, up to floating-point rounding.
|
||
pub truncation_bound: f64,
|
||
/// The most terms the sum held after any gate.
|
||
pub peak_terms: usize,
|
||
/// Terms left when propagation finished.
|
||
pub final_terms: usize,
|
||
/// Term-gate updates performed: a portable measure of the work, equal on
|
||
/// every machine.
|
||
pub term_updates: u64,
|
||
}
|
||
|
||
/// `⟨b| U† O U |b⟩` for the computational basis state `|b⟩` whose set qubits
|
||
/// are `ones` (empty for `|0…0⟩`), with `O` a weighted sum of Pauli strings.
|
||
pub fn expectation(
|
||
circuit: &RotCircuit,
|
||
observable: &[(f64, PauliString)],
|
||
ones: &[u32],
|
||
cfg: &SpdConfig,
|
||
) -> Result<SpdResult, SpdError> {
|
||
circuit.validate()?;
|
||
for (_, s) in observable {
|
||
check_string(s, circuit.n)?;
|
||
}
|
||
for &q in ones {
|
||
if q >= circuit.n {
|
||
return Err(SpdError::QubitOutOfRange { qubit: q, n: circuit.n });
|
||
}
|
||
}
|
||
Ok(match circuit.n.div_ceil(64) {
|
||
1 => run::<1>(circuit, observable, ones, cfg),
|
||
2 => run::<2>(circuit, observable, ones, cfg),
|
||
3 | 4 => run::<4>(circuit, observable, ones, cfg),
|
||
5..=8 => run::<8>(circuit, observable, ones, cfg),
|
||
_ => run::<16>(circuit, observable, ones, cfg),
|
||
})
|
||
}
|
||
|
||
/// `⟨0…0| U† Z_q U |0…0⟩`.
|
||
pub fn expect_z(circuit: &RotCircuit, q: u32, cfg: &SpdConfig) -> Result<SpdResult, SpdError> {
|
||
expectation(circuit, &[(1.0, vec![(q, Pauli::Z)])], &[], cfg)
|
||
}
|
||
|
||
fn run<const W: usize>(
|
||
circuit: &RotCircuit,
|
||
observable: &[(f64, PauliString)],
|
||
ones: &[u32],
|
||
cfg: &SpdConfig,
|
||
) -> SpdResult {
|
||
let mut t = Terms::<W>::new();
|
||
for (c, s) in observable {
|
||
t.add(Key::from_string(s), *c);
|
||
}
|
||
let mut bound = t.truncate(0.0, None);
|
||
let mut peak = t.live();
|
||
let mut work = 0u64;
|
||
let mut pending = Vec::new();
|
||
for op in circuit.ops.iter().rev() {
|
||
work += t.live() as u64;
|
||
match op {
|
||
Op::Rot { axis, theta } => {
|
||
let (s, c) = sin_cos(*theta);
|
||
rotate(&mut t, &Key::from_string(axis), s, c, &mut pending);
|
||
bound += t.truncate(cfg.threshold, cfg.max_weight);
|
||
}
|
||
Op::QuarterRot { axis, k } => quarter(&mut t, &Key::from_string(axis), *k),
|
||
other => clifford(&mut t, other),
|
||
}
|
||
peak = peak.max(t.live());
|
||
}
|
||
let b = Key::<W>::from_string(&ones.iter().map(|&q| (q, Pauli::Z)).collect::<Vec<_>>());
|
||
let mut value = 0.0;
|
||
for i in 0..t.keys.len() {
|
||
let k = &t.keys[i];
|
||
if t.coef[i] != 0.0 && k.x.iter().all(|&w| w == 0) {
|
||
let flips: u32 = (0..W).map(|w| (k.z[w] & b.z[w]).count_ones()).sum();
|
||
value += if flips & 1 == 1 { -t.coef[i] } else { t.coef[i] };
|
||
}
|
||
}
|
||
SpdResult { value, truncation_bound: bound, peak_terms: peak, final_terms: t.live(), term_updates: work }
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// The heavy-hex kicked-Ising experiment
|
||
// ---------------------------------------------------------------------------
|
||
|
||
/// The 127-qubit heavy-hex coupling map: seven rows (14, 15, 15, 15, 15, 15
|
||
/// and 14 qubits) joined by six rows of four bridge qubits. Returns the 144
|
||
/// edges, each with the smaller qubit first.
|
||
pub fn heavy_hex_127() -> Vec<(u32, u32)> {
|
||
let rows: [(u32, u32); 7] = [(0, 13), (18, 32), (37, 51), (56, 70), (75, 89), (94, 108), (113, 126)];
|
||
// Each bridge row: (first bridge qubit, column offset into the row above,
|
||
// column offset into the row below).
|
||
let bridges: [(u32, u32, u32); 6] = [(14, 0, 0), (33, 2, 2), (52, 0, 0), (71, 2, 2), (90, 0, 0), (109, 2, 2)];
|
||
let mut edges = Vec::with_capacity(144);
|
||
for &(a, b) in &rows {
|
||
for q in a..b {
|
||
edges.push((q, q + 1));
|
||
}
|
||
}
|
||
for (r, &(first, up, down)) in bridges.iter().enumerate() {
|
||
let (above, below) = (rows[r].0, rows[r + 1].0);
|
||
for j in 0..4 {
|
||
let top = above + up + 4 * j;
|
||
// Row 6 starts one column right of the rows above it, so its
|
||
// bridge columns are one index earlier.
|
||
let bottom = below + down + 4 * j - if r == 5 { 1 } else { 0 };
|
||
edges.push((top.min(first + j), top.max(first + j)));
|
||
edges.push(((first + j).min(bottom), (first + j).max(bottom)));
|
||
}
|
||
}
|
||
edges.sort_unstable();
|
||
edges
|
||
}
|
||
|
||
/// The edges as layers of vertex-disjoint edges: an edge colouring. A graph
|
||
/// whose maximum degree is `Δ` and which is bipartite, such as heavy-hex, gets
|
||
/// exactly `Δ` layers (König's theorem): each edge takes a colour free at both
|
||
/// ends, and when none is, an alternating path is recoloured to free one. A
|
||
/// graph with an odd cycle may need a layer more, which it is given. Within a
|
||
/// layer the edges are in their input order, and the layers in colour order.
|
||
pub fn edge_layers(n: u32, edges: &[(u32, u32)]) -> Vec<Vec<(u32, u32)>> {
|
||
let n = n as usize;
|
||
let mut deg = vec![0usize; n];
|
||
for &(a, b) in edges {
|
||
deg[a as usize] += 1;
|
||
deg[b as usize] += 1;
|
||
}
|
||
let mut ncol = deg.iter().copied().max().unwrap_or(0).max(1);
|
||
// at[v][c]: the edge of colour c at vertex v, if any
|
||
let mut at: Vec<Vec<Option<usize>>> = vec![vec![None; ncol]; n];
|
||
let mut colour = vec![usize::MAX; edges.len()];
|
||
let other = |f: usize, x: usize| {
|
||
let (a, b) = edges[f];
|
||
if a as usize == x { b as usize } else { a as usize }
|
||
};
|
||
for (e, &(u, v)) in edges.iter().enumerate() {
|
||
let (u, v) = (u as usize, v as usize);
|
||
let free_u = (0..ncol).find(|&c| at[u][c].is_none());
|
||
let free_v = (0..ncol).find(|&c| at[v][c].is_none());
|
||
let chosen = match (free_u, free_v) {
|
||
(Some(a), Some(_)) if at[v][a].is_none() => Some(a),
|
||
(Some(a), Some(b)) => {
|
||
// Walk the a/b alternating path from v; flipping it frees a at v.
|
||
let mut path = Vec::new();
|
||
let (mut x, mut c) = (v, a);
|
||
let mut reaches_u = false;
|
||
while let Some(f) = at[x][c] {
|
||
path.push(f);
|
||
x = other(f, x);
|
||
if x == u {
|
||
reaches_u = true;
|
||
break;
|
||
}
|
||
c = if c == a { b } else { a };
|
||
}
|
||
if reaches_u {
|
||
None
|
||
} else {
|
||
for &f in &path {
|
||
let (p, q) = (edges[f].0 as usize, edges[f].1 as usize);
|
||
at[p][colour[f]] = None;
|
||
at[q][colour[f]] = None;
|
||
}
|
||
for &f in &path {
|
||
let (p, q) = (edges[f].0 as usize, edges[f].1 as usize);
|
||
colour[f] = if colour[f] == a { b } else { a };
|
||
at[p][colour[f]] = Some(f);
|
||
at[q][colour[f]] = Some(f);
|
||
}
|
||
Some(a)
|
||
}
|
||
}
|
||
_ => None,
|
||
};
|
||
let c = chosen.unwrap_or_else(|| {
|
||
ncol += 1;
|
||
for row in at.iter_mut() {
|
||
row.push(None);
|
||
}
|
||
ncol - 1
|
||
});
|
||
colour[e] = c;
|
||
at[u][c] = Some(e);
|
||
at[v][c] = Some(e);
|
||
}
|
||
(0..ncol)
|
||
.map(|c| edges.iter().zip(&colour).filter(|&(_, &k)| k == c).map(|(&e, _)| e).collect::<Vec<_>>())
|
||
.filter(|layer| !layer.is_empty())
|
||
.collect()
|
||
}
|
||
|
||
/// The kicked-Ising circuit of the heavy-hex utility experiment: `steps`
|
||
/// Trotter steps, each `RX(θ_h)` on every qubit and then `RZZ(−π/2)` on every
|
||
/// edge, layer by layer as [`edge_layers`] partitions them: three layers of
|
||
/// vertex-disjoint edges on heavy-hex, as on the device. The `RZZ` gates are applied
|
||
/// exactly, as quarter turns.
|
||
pub fn kicked_ising(n: u32, edges: &[(u32, u32)], theta_h: f64, steps: u32) -> RotCircuit {
|
||
let layers = edge_layers(n, edges);
|
||
let mut c = RotCircuit::new(n);
|
||
for _ in 0..steps {
|
||
for q in 0..n {
|
||
c.rx(q, theta_h);
|
||
}
|
||
for layer in &layers {
|
||
for &(a, b) in layer {
|
||
// RZZ(−π/2) = exp(+iπ/4·ZZ) = exp(−i(3π/2)ZZ/2) up to a global phase.
|
||
c.push(Op::QuarterRot { axis: vec![(a, Pauli::Z), (b, Pauli::Z)], k: 3 });
|
||
}
|
||
}
|
||
}
|
||
c
|
||
}
|
||
|
||
// ---------------------------------------------------------------------------
|
||
// From the dyadic core
|
||
// ---------------------------------------------------------------------------
|
||
|
||
/// The dyadic circuit as Pauli rotations. Fixed Cliffords map one to one; `T`
|
||
/// and `P(k)` become `Z` rotations; a controlled phase on any number of
|
||
/// controls becomes rotations about `Z` strings, from
|
||
/// `|1…1⟩⟨1…1| = Π (I − Z)/2`; controlled `X` and `Y` beyond one control are
|
||
/// that, conjugated by `H` (and `S`) on the target. Global phases drop out of a
|
||
/// Heisenberg-picture expectation and are not tracked.
|
||
#[cfg(feature = "quantum")]
|
||
pub fn from_dyadic(c: &crate::quantum::Circuit) -> Result<RotCircuit, SpdError> {
|
||
use crate::quantum::BaseGate as B;
|
||
let mut out = RotCircuit::new(u32::from(c.n_qubits));
|
||
let tau = 2.0 * std::f64::consts::PI;
|
||
for g in &c.ops {
|
||
let t = u32::from(g.target);
|
||
let ctrls: Vec<u32> = g.controls.iter().map(|&q| u32::from(q)).collect();
|
||
// The phase a diagonal base applies to |1⟩, in turns of 2π.
|
||
let phase = |base: B| -> Option<f64> {
|
||
match base {
|
||
B::Z => Some(0.5),
|
||
B::S => Some(0.25),
|
||
B::Sdg => Some(-0.25),
|
||
B::T => Some(0.125),
|
||
B::Tdg => Some(-0.125),
|
||
B::P => Some(1.0 / (1u64 << g.param) as f64),
|
||
_ => None,
|
||
}
|
||
};
|
||
match (g.base, ctrls.as_slice()) {
|
||
(B::I, _) => {}
|
||
(B::H, []) => _ = out.push(Op::H(t)),
|
||
(B::X, []) => _ = out.push(Op::X(t)),
|
||
(B::Y, []) => _ = out.push(Op::Y(t)),
|
||
(B::Z, []) => _ = out.push(Op::Z(t)),
|
||
(B::S, []) => _ = out.push(Op::S(t)),
|
||
(B::Sdg, []) => _ = out.push(Op::Sdg(t)),
|
||
(B::X, [a]) => _ = out.push(Op::CX(*a, t)),
|
||
(B::Z, [a]) => _ = out.push(Op::CZ(*a, t)),
|
||
(B::Y, [a]) => {
|
||
out.push(Op::Sdg(t)).push(Op::CX(*a, t)).push(Op::S(t));
|
||
}
|
||
(B::X, _) => {
|
||
out.push(Op::H(t));
|
||
controlled_phase(&mut out, &ctrls, t, 0.5 * tau)?;
|
||
out.push(Op::H(t));
|
||
}
|
||
(B::Y, _) => {
|
||
out.push(Op::Sdg(t)).push(Op::H(t));
|
||
controlled_phase(&mut out, &ctrls, t, 0.5 * tau)?;
|
||
out.push(Op::H(t)).push(Op::S(t));
|
||
}
|
||
(base, _) => match phase(base) {
|
||
Some(turns) => controlled_phase(&mut out, &ctrls, t, turns * tau)?,
|
||
None => return Err(SpdError::Unsupported(format!("{base:?} on {} controls", ctrls.len()))),
|
||
},
|
||
}
|
||
}
|
||
Ok(out)
|
||
}
|
||
|
||
/// `diag(…, e^{iφ})` on `controls ∪ {target}`: the phase `φ` on `|1…1⟩` only.
|
||
/// With `m` qubits, `|1…1⟩⟨1…1| = 2⁻ᵐ Σ_S (−1)^{|S|} Z_S`, so the gate is
|
||
/// `Π_S exp(iφ(−1)^{|S|} Z_S / 2ᵐ)` over nonempty subsets `S` (the empty one is
|
||
/// a global phase): a rotation `exp(−iθZ_S/2)` with `θ = −2φ(−1)^{|S|}/2ᵐ`.
|
||
///
|
||
/// The subsets are counted in a `u32`, so `m` must be below 32; a gate on more
|
||
/// qubits (a circuit2 gate may carry up to 254 controls) is refused as
|
||
/// unsupported rather than shifted out of range.
|
||
#[cfg(feature = "quantum")]
|
||
fn controlled_phase(out: &mut RotCircuit, controls: &[u32], target: u32, phi: f64) -> Result<(), SpdError> {
|
||
let mut qs = controls.to_vec();
|
||
qs.push(target);
|
||
let subsets = u32::try_from(qs.len()).ok().and_then(|m| 1u32.checked_shl(m).map(|s| (m, s)));
|
||
let Some((m, subsets)) = subsets else {
|
||
return Err(SpdError::Unsupported(format!("a phase on {} qubits; at most 31", qs.len())));
|
||
};
|
||
for mask in 1u32..subsets {
|
||
let axis: PauliString = (0..m).filter(|i| mask >> i & 1 == 1).map(|i| (qs[i as usize], Pauli::Z)).collect();
|
||
let sign = if axis.len() % 2 == 0 { 1.0 } else { -1.0 };
|
||
out.push(Op::Rot { axis, theta: -2.0 * phi * sign / f64::from(subsets) });
|
||
}
|
||
Ok(())
|
||
}
|
||
|
||
#[cfg(test)]
|
||
mod tests {
|
||
use super::*;
|
||
|
||
// ---- an independent Schrödinger-picture reference ----------------------
|
||
|
||
type C = (f64, f64);
|
||
|
||
fn cmul(a: C, b: C) -> C {
|
||
(a.0 * b.0 - a.1 * b.1, a.0 * b.1 + a.1 * b.0)
|
||
}
|
||
|
||
/// `P|ψ⟩` for a Pauli string, on a dense state.
|
||
fn apply_pauli(psi: &[C], s: &[(u32, Pauli)]) -> Vec<C> {
|
||
let mut out = vec![(0.0, 0.0); psi.len()];
|
||
for (j, &a) in psi.iter().enumerate() {
|
||
let mut k = j;
|
||
let mut ph: C = (1.0, 0.0);
|
||
for &(q, p) in s {
|
||
let bit = (j >> q) & 1;
|
||
match p {
|
||
Pauli::X => k ^= 1 << q,
|
||
Pauli::Y => {
|
||
k ^= 1 << q;
|
||
// Y|0⟩ = i|1⟩, Y|1⟩ = −i|0⟩
|
||
ph = cmul(ph, if bit == 0 { (0.0, 1.0) } else { (0.0, -1.0) });
|
||
}
|
||
Pauli::Z => {
|
||
if bit == 1 {
|
||
ph = (-ph.0, -ph.1);
|
||
}
|
||
}
|
||
}
|
||
}
|
||
let v = cmul(ph, a);
|
||
out[k].0 += v.0;
|
||
out[k].1 += v.1;
|
||
}
|
||
out
|
||
}
|
||
|
||
fn one_qubit(psi: &mut [C], q: u32, m: [[C; 2]; 2]) {
|
||
for j in 0..psi.len() {
|
||
if (j >> q) & 1 == 0 {
|
||
let k = j | (1 << q);
|
||
let (a, b) = (psi[j], psi[k]);
|
||
let add = |x: C, y: C| (x.0 + y.0, x.1 + y.1);
|
||
psi[j] = add(cmul(m[0][0], a), cmul(m[0][1], b));
|
||
psi[k] = add(cmul(m[1][0], a), cmul(m[1][1], b));
|
||
}
|
||
}
|
||
}
|
||
|
||
/// Run `circuit` on `|ones⟩` as a dense state, with the platform's `sin`
|
||
/// and `cos`: nothing shared with the code under test but the `Op` type.
|
||
fn dense(circuit: &RotCircuit, ones: &[u32]) -> Vec<C> {
|
||
let n = circuit.n;
|
||
let mut psi = vec![(0.0, 0.0); 1 << n];
|
||
psi[ones.iter().fold(0usize, |a, &q| a | 1 << q)] = (1.0, 0.0);
|
||
let r = std::f64::consts::FRAC_1_SQRT_2;
|
||
let (o, z, i, ni) = ((1.0, 0.0), (0.0, 0.0), (0.0, 1.0), (0.0, -1.0));
|
||
for op in &circuit.ops {
|
||
match op {
|
||
Op::H(q) => one_qubit(&mut psi, *q, [[(r, 0.0), (r, 0.0)], [(r, 0.0), (-r, 0.0)]]),
|
||
Op::S(q) => one_qubit(&mut psi, *q, [[o, z], [z, i]]),
|
||
Op::Sdg(q) => one_qubit(&mut psi, *q, [[o, z], [z, ni]]),
|
||
Op::SX(q) => one_qubit(&mut psi, *q, [[(0.5, 0.5), (0.5, -0.5)], [(0.5, -0.5), (0.5, 0.5)]]),
|
||
Op::SXdg(q) => one_qubit(&mut psi, *q, [[(0.5, -0.5), (0.5, 0.5)], [(0.5, 0.5), (0.5, -0.5)]]),
|
||
Op::X(q) => psi = apply_pauli(&psi, &[(*q, Pauli::X)]),
|
||
Op::Y(q) => psi = apply_pauli(&psi, &[(*q, Pauli::Y)]),
|
||
Op::Z(q) => psi = apply_pauli(&psi, &[(*q, Pauli::Z)]),
|
||
Op::CX(c, t) => {
|
||
for j in 0..psi.len() {
|
||
if (j >> c) & 1 == 1 && (j >> t) & 1 == 0 {
|
||
psi.swap(j, j | 1 << t);
|
||
}
|
||
}
|
||
}
|
||
Op::CZ(a, b) => {
|
||
for (j, v) in psi.iter_mut().enumerate() {
|
||
if (j >> a) & 1 == 1 && (j >> b) & 1 == 1 {
|
||
*v = (-v.0, -v.1);
|
||
}
|
||
}
|
||
}
|
||
Op::Swap(a, b) => {
|
||
for j in 0..psi.len() {
|
||
if (j >> a) & 1 == 1 && (j >> b) & 1 == 0 {
|
||
psi.swap(j, (j & !(1 << a)) | 1 << b);
|
||
}
|
||
}
|
||
}
|
||
Op::Rot { axis, theta } => {
|
||
let p = apply_pauli(&psi, axis);
|
||
let (c, s) = ((theta / 2.0).cos(), (theta / 2.0).sin());
|
||
for (v, w) in psi.iter_mut().zip(&p) {
|
||
// cos(θ/2)ψ − i sin(θ/2) Pψ
|
||
*v = (c * v.0 + s * w.1, c * v.1 - s * w.0);
|
||
}
|
||
}
|
||
Op::QuarterRot { axis, k } => {
|
||
let theta = f64::from(*k) * std::f64::consts::FRAC_PI_2;
|
||
let p = apply_pauli(&psi, axis);
|
||
let (c, s) = ((theta / 2.0).cos(), (theta / 2.0).sin());
|
||
for (v, w) in psi.iter_mut().zip(&p) {
|
||
*v = (c * v.0 + s * w.1, c * v.1 - s * w.0);
|
||
}
|
||
}
|
||
}
|
||
}
|
||
psi
|
||
}
|
||
|
||
fn dense_expect(circuit: &RotCircuit, obs: &[(u32, Pauli)], ones: &[u32]) -> f64 {
|
||
let psi = dense(circuit, ones);
|
||
let p = apply_pauli(&psi, obs);
|
||
psi.iter().zip(&p).map(|(a, b)| a.0 * b.0 + a.1 * b.1).sum()
|
||
}
|
||
|
||
/// A small deterministic generator for test circuits.
|
||
struct Rng(u64);
|
||
impl Rng {
|
||
fn next(&mut self) -> u64 {
|
||
self.0 = self.0.wrapping_add(0x9e37_79b9_7f4a_7c15);
|
||
mix(self.0)
|
||
}
|
||
fn below(&mut self, n: u64) -> u64 {
|
||
self.next() % n
|
||
}
|
||
fn angle(&mut self) -> f64 {
|
||
(self.next() >> 11) as f64 / (1u64 << 53) as f64 * 8.0 - 4.0
|
||
}
|
||
fn pauli(&mut self) -> Pauli {
|
||
[Pauli::X, Pauli::Y, Pauli::Z][self.below(3) as usize]
|
||
}
|
||
fn string(&mut self, n: u32, max: usize) -> PauliString {
|
||
let mut qs: Vec<u32> = (0..n).collect();
|
||
let len = 1 + self.below(max as u64) as usize;
|
||
let mut s = Vec::new();
|
||
for _ in 0..len.min(n as usize) {
|
||
let i = self.below(qs.len() as u64) as usize;
|
||
s.push((qs.swap_remove(i), self.pauli()));
|
||
}
|
||
s
|
||
}
|
||
}
|
||
|
||
fn random_circuit(rng: &mut Rng, n: u32, len: usize) -> RotCircuit {
|
||
let mut c = RotCircuit::new(n);
|
||
for _ in 0..len {
|
||
let q = rng.below(u64::from(n)) as u32;
|
||
let mut r = rng.below(u64::from(n - 1)) as u32;
|
||
if r >= q {
|
||
r += 1;
|
||
}
|
||
let op = match rng.below(14) {
|
||
0 => Op::H(q),
|
||
1 => Op::S(q),
|
||
2 => Op::Sdg(q),
|
||
3 => Op::SX(q),
|
||
4 => Op::SXdg(q),
|
||
5 => Op::X(q),
|
||
6 => Op::Y(q),
|
||
7 => Op::Z(q),
|
||
8 => Op::CX(q, r),
|
||
9 => Op::CZ(q, r),
|
||
10 => Op::Swap(q, r),
|
||
11 => Op::QuarterRot { axis: rng.string(n, 3), k: 1 + rng.below(3) as u8 },
|
||
_ => Op::Rot { axis: rng.string(n, 3), theta: rng.angle() },
|
||
};
|
||
c.push(op);
|
||
}
|
||
c
|
||
}
|
||
|
||
// ---- trigonometry --------------------------------------------------------
|
||
|
||
#[test]
|
||
fn sin_cos_is_within_an_ulp_of_the_platform() {
|
||
let mut rng = Rng(7);
|
||
let ulp = |a: f64, b: f64| (a.to_bits() as i64 - b.to_bits() as i64).unsigned_abs();
|
||
for i in 0..200_000 {
|
||
let x = match i % 4 {
|
||
0 => rng.angle(),
|
||
1 => rng.angle() * 1e-6,
|
||
2 => rng.angle() * 1e4,
|
||
_ => rng.angle() * 0.2,
|
||
};
|
||
let (s, c) = sin_cos(x);
|
||
// Both are faithful roundings; they may differ by one ulp at most,
|
||
// and only where the true value is close to a rounding boundary.
|
||
assert!(ulp(s, x.sin()) <= 1 || (s - x.sin()).abs() < 1e-300, "sin {x}: {s} vs {}", x.sin());
|
||
assert!(ulp(c, x.cos()) <= 1, "cos {x}: {c} vs {}", x.cos());
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn sin_cos_bits_are_pinned() {
|
||
// Bit patterns this implementation must produce on every target: they
|
||
// are what the cross-architecture run (wasm32-wasip2) is checked against.
|
||
// Each is within 0.51 ulp of the true value, checked against an
|
||
// 80-digit Taylor series when it was pinned.
|
||
let pins: [(f64, u64, u64); 4] = [
|
||
(0.3, 0x3fd2_e9cd_95ba_ba33, 0x3fee_921d_d42f_09ba),
|
||
(1.0, 0x3fea_ed54_8f09_0cee, 0x3fe1_4a28_0fb5_068c),
|
||
(-2.5, 0xbfe3_26af_0dcf_cab0, 0xbfe9_a2f7_ef85_8b7d),
|
||
(100.0, 0xbfe0_3425_b78c_4db8, 0x3feb_981d_bf66_5fdf),
|
||
];
|
||
for (x, sb, cb) in pins {
|
||
let (s, c) = sin_cos(x);
|
||
assert_eq!((s.to_bits(), c.to_bits()), (sb, cb), "x = {x}: sin {s:e} cos {c:e}");
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn quarter_turns_are_exact() {
|
||
for (k, x) in [(1, std::f64::consts::FRAC_PI_2), (2, std::f64::consts::PI)] {
|
||
let (s, c) = sin_cos(x);
|
||
// Not exactly 0 and ±1 — π/2 is not representable — which is why the
|
||
// Clifford rotations have their own exact op.
|
||
assert!(c.abs() < 1e-15 || k == 2, "k = {k}");
|
||
assert!((s.abs() - if k == 1 { 1.0 } else { 0.0 }).abs() < 1e-15);
|
||
}
|
||
}
|
||
|
||
// ---- Pauli algebra against matrices --------------------------------------
|
||
|
||
#[test]
|
||
fn every_clifford_conjugates_every_pauli_as_its_matrix_does() {
|
||
// For each gate G and each Pauli P on two qubits, the Heisenberg rule's
|
||
// G†PG must give the same expectation on random states as conjugating
|
||
// the state instead: ⟨ψ|G†PG|ψ⟩ = ⟨Gψ|P|Gψ⟩.
|
||
let gates = [
|
||
Op::H(0), Op::S(0), Op::Sdg(0), Op::SX(0), Op::SXdg(0), Op::X(0), Op::Y(0), Op::Z(0),
|
||
Op::H(1), Op::SX(1), Op::CX(0, 1), Op::CX(1, 0), Op::CZ(0, 1), Op::Swap(0, 1),
|
||
];
|
||
let mut rng = Rng(11);
|
||
let paulis = [None, Some(Pauli::X), Some(Pauli::Y), Some(Pauli::Z)];
|
||
for g in &gates {
|
||
for p0 in paulis {
|
||
for p1 in paulis {
|
||
let obs: PauliString = [(0, p0), (1, p1)].iter().filter_map(|&(q, p)| p.map(|p| (q, p))).collect();
|
||
if obs.is_empty() {
|
||
continue;
|
||
}
|
||
// a random state prepared by rotations, then the gate
|
||
let mut c = RotCircuit::new(2);
|
||
for _ in 0..6 {
|
||
c.push(Op::Rot { axis: rng.string(2, 2), theta: rng.angle() });
|
||
}
|
||
c.push(g.clone());
|
||
let want = dense_expect(&c, &obs, &[]);
|
||
let got = expectation(&c, &[(1.0, obs.clone())], &[], &SpdConfig { threshold: 0.0, max_weight: None }).unwrap();
|
||
assert!((got.value - want).abs() < 1e-12, "{g:?} {obs:?}: {} vs {want}", got.value);
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn random_circuits_match_the_dense_reference() {
|
||
let mut rng = Rng(2026);
|
||
for trial in 0..300 {
|
||
let n = 2 + rng.below(6) as u32;
|
||
let len = 6 + rng.below(30) as usize;
|
||
let c = random_circuit(&mut rng, n, len);
|
||
let obs = rng.string(n, 3);
|
||
let ones: Vec<u32> = (0..n).filter(|_| rng.below(2) == 1).collect();
|
||
let want = dense_expect(&c, &obs, &ones);
|
||
let got = expectation(&c, &[(1.0, obs.clone())], &ones, &SpdConfig { threshold: 0.0, max_weight: None }).unwrap();
|
||
assert!((got.value - want).abs() < 1e-10, "trial {trial}: {} vs {want}", got.value);
|
||
assert_eq!(got.truncation_bound, 0.0, "nothing is dropped at threshold 0 but exact zeros");
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn the_truncation_bound_holds() {
|
||
// At every threshold the computed value is within the reported bound of
|
||
// the exact one — the guarantee, checked on circuits where the bound is
|
||
// not vacuous.
|
||
let mut rng = Rng(99);
|
||
let mut nontrivial = 0;
|
||
for _ in 0..120 {
|
||
let n = 6;
|
||
let c = random_circuit(&mut rng, n, 60);
|
||
let obs = vec![(rng.below(6) as u32, Pauli::Z)];
|
||
let exact = dense_expect(&c, &obs, &[]);
|
||
for threshold in [1e-3, 1e-2, 5e-2, 0.2] {
|
||
let r = expectation(&c, &[(1.0, obs.clone())], &[], &SpdConfig { threshold, max_weight: None }).unwrap();
|
||
assert!((r.value - exact).abs() <= r.truncation_bound + 1e-12, "|{} − {exact}| > {}", r.value, r.truncation_bound);
|
||
if r.truncation_bound > 0.0 && r.truncation_bound < 1.0 {
|
||
nontrivial += 1;
|
||
}
|
||
}
|
||
}
|
||
assert!(nontrivial > 50, "the bound was exercised {nontrivial} times");
|
||
}
|
||
|
||
#[test]
|
||
fn results_are_bit_identical_from_run_to_run() {
|
||
// A fresh std HashMap is seeded differently every time, so a term table
|
||
// that summed in its order would already differ here, in one process.
|
||
let c = kicked_ising(127, &heavy_hex_127(), 0.6, 5);
|
||
let cfg = SpdConfig { threshold: 1e-4, max_weight: None };
|
||
let a = expect_z(&c, 62, &cfg).unwrap();
|
||
assert!(a.peak_terms > 1000, "the run merges many terms: {}", a.peak_terms);
|
||
for _ in 0..4 {
|
||
let b = expect_z(&c, 62, &cfg).unwrap();
|
||
assert_eq!(a.value.to_bits(), b.value.to_bits());
|
||
assert_eq!(a.truncation_bound.to_bits(), b.truncation_bound.to_bits());
|
||
assert_eq!((a.peak_terms, a.term_updates), (b.peak_terms, b.term_updates));
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn word_width_does_not_change_the_answer() {
|
||
// The same circuit, placed in a 30-qubit and a 900-qubit register, is
|
||
// computed with one and with sixteen words per string: the bits agree.
|
||
let chain: Vec<(u32, u32)> = (0..29).map(|q| (q, q + 1)).collect();
|
||
let small = kicked_ising(30, &chain, 0.5, 6);
|
||
let mut big = small.clone();
|
||
big.n = 900;
|
||
let cfg = SpdConfig { threshold: 1e-6, max_weight: None };
|
||
let a = expect_z(&small, 7, &cfg).unwrap();
|
||
let b = expect_z(&big, 7, &cfg).unwrap();
|
||
assert_eq!(a.value.to_bits(), b.value.to_bits());
|
||
assert_eq!(a.term_updates, b.term_updates);
|
||
}
|
||
|
||
#[test]
|
||
fn invalid_circuits_are_refused() {
|
||
let mut c = RotCircuit::new(3);
|
||
c.push(Op::CX(1, 1));
|
||
assert_eq!(c.validate(), Err(SpdError::RepeatedQubit(1)));
|
||
let mut c = RotCircuit::new(3);
|
||
c.rx(3, 0.1);
|
||
assert_eq!(c.validate(), Err(SpdError::QubitOutOfRange { qubit: 3, n: 3 }));
|
||
let mut c = RotCircuit::new(3);
|
||
c.rx(0, f64::NAN);
|
||
assert!(matches!(c.validate(), Err(SpdError::BadAngle(_))));
|
||
let mut c = RotCircuit::new(3);
|
||
c.push(Op::QuarterRot { axis: vec![(0, Pauli::Z)], k: 4 });
|
||
assert_eq!(c.validate(), Err(SpdError::BadQuarter(4)));
|
||
assert_eq!(RotCircuit::new(1025).validate(), Err(SpdError::TooManyQubits(1025)));
|
||
}
|
||
|
||
// ---- the heavy-hex experiment ---------------------------------------------
|
||
|
||
#[test]
|
||
fn the_heavy_hex_map_has_the_published_shape() {
|
||
let e = heavy_hex_127();
|
||
assert_eq!(e.len(), 144);
|
||
let mut deg = [0u32; 127];
|
||
for &(a, b) in &e {
|
||
assert!(a < b && b < 127);
|
||
deg[a as usize] += 1;
|
||
deg[b as usize] += 1;
|
||
}
|
||
assert!(deg.iter().all(|&d| (1..=3).contains(&d)));
|
||
// Heavy-hex: every bridge qubit has degree 2, and no two degree-3
|
||
// qubits are adjacent.
|
||
for &(a, b) in &e {
|
||
assert!(!(deg[a as usize] == 3 && deg[b as usize] == 3), "({a},{b})");
|
||
}
|
||
for b in [14u32, 15, 16, 17, 33, 34, 35, 36, 52, 53, 54, 55, 71, 72, 73, 74, 90, 91, 92, 93, 109, 110, 111, 112] {
|
||
assert_eq!(deg[b as usize], 2, "bridge {b}");
|
||
}
|
||
// connected
|
||
let mut seen = [false; 127];
|
||
let mut stack = vec![0u32];
|
||
while let Some(q) = stack.pop() {
|
||
if !std::mem::replace(&mut seen[q as usize], true) {
|
||
stack.extend(e.iter().filter(|&&(a, b)| a == q || b == q).map(|&(a, b)| a + b - q));
|
||
}
|
||
}
|
||
assert!(seen.iter().all(|&s| s));
|
||
// It is bipartite, so three matchings cover it, as the experiment's
|
||
// three RZZ layers do.
|
||
let mut colour = [u8::MAX; 127];
|
||
colour[0] = 0;
|
||
let mut stack = vec![0u32];
|
||
while let Some(q) = stack.pop() {
|
||
for &(a, b) in e.iter().filter(|&&(a, b)| a == q || b == q) {
|
||
let o = (a + b - q) as usize;
|
||
if colour[o] == u8::MAX {
|
||
colour[o] = 1 - colour[q as usize];
|
||
stack.push(o as u32);
|
||
}
|
||
assert_ne!(colour[o], colour[q as usize]);
|
||
}
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn edge_layers_are_proper_colourings_with_as_few_layers_as_the_degree() {
|
||
let check = |n: u32, edges: &[(u32, u32)], want: usize| {
|
||
let layers = edge_layers(n, edges);
|
||
assert_eq!(layers.len(), want);
|
||
let mut all: Vec<(u32, u32)> = layers.iter().flatten().copied().collect();
|
||
all.sort_unstable();
|
||
let mut sorted = edges.to_vec();
|
||
sorted.sort_unstable();
|
||
assert_eq!(all, sorted, "every edge exactly once");
|
||
for layer in &layers {
|
||
let mut seen = std::collections::BTreeSet::new();
|
||
for &(a, b) in layer {
|
||
assert!(seen.insert(a) && seen.insert(b), "vertex-disjoint");
|
||
}
|
||
}
|
||
};
|
||
check(127, &heavy_hex_127(), 3);
|
||
let chain: Vec<(u32, u32)> = (0..9).map(|q| (q, q + 1)).collect();
|
||
check(10, &chain, 2);
|
||
// A triangle is not bipartite: three edges, degree two, three layers.
|
||
check(3, &[(0, 1), (1, 2), (0, 2)], 3);
|
||
}
|
||
|
||
#[test]
|
||
fn the_kicked_ising_clifford_points_are_exact() {
|
||
// θ_h = 0: no rotation, every ⟨Z_q⟩ = 1. θ_h = π/2 written as quarter
|
||
// turns: a Clifford circuit, so the Heisenberg sum stays one term and
|
||
// ⟨Z_q⟩ is exactly 0 or ±1 with nothing dropped.
|
||
let e = heavy_hex_127();
|
||
let still = kicked_ising(127, &e, 0.0, 20);
|
||
let r = expect_z(&still, 62, &SpdConfig::default()).unwrap();
|
||
assert_eq!((r.value, r.truncation_bound, r.peak_terms), (1.0, 0.0, 1));
|
||
let mut cliff = RotCircuit::new(127);
|
||
for _ in 0..20 {
|
||
for q in 0..127 {
|
||
cliff.push(Op::QuarterRot { axis: vec![(q, Pauli::X)], k: 1 });
|
||
}
|
||
for &(a, b) in &e {
|
||
cliff.push(Op::QuarterRot { axis: vec![(a, Pauli::Z), (b, Pauli::Z)], k: 3 });
|
||
}
|
||
}
|
||
for q in [0, 13, 62, 126] {
|
||
let r = expect_z(&cliff, q, &SpdConfig::default()).unwrap();
|
||
assert!(r.value == 0.0 || r.value.abs() == 1.0, "q {q}: {}", r.value);
|
||
assert_eq!((r.truncation_bound, r.peak_terms), (0.0, 1));
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn kicked_ising_on_a_fragment_matches_the_dense_reference() {
|
||
// The same circuit family on a 15-qubit heavy-hex fragment, where a
|
||
// dense state fits: exact propagation agrees with it at every angle.
|
||
let e: Vec<(u32, u32)> = heavy_hex_127().into_iter().filter(|&(a, b)| a < 15 && b < 15).collect();
|
||
assert!(e.len() >= 13);
|
||
for theta_h in [0.1, 0.4, 0.7854, 1.2] {
|
||
let c = kicked_ising(15, &e, theta_h, 4);
|
||
for q in [0, 4, 14] {
|
||
let want = dense_expect(&c, &[(q, Pauli::Z)], &[]);
|
||
let got = expect_z(&c, q, &SpdConfig { threshold: 0.0, max_weight: None }).unwrap();
|
||
assert!((got.value - want).abs() < 1e-10, "θ_h {theta_h} q {q}: {} vs {want}", got.value);
|
||
}
|
||
}
|
||
}
|
||
|
||
// ---- the dyadic adapter ------------------------------------------------------
|
||
|
||
/// A controlled phase counts its `2^m` subsets in a `u32`. Forty controls
|
||
/// would shift out of range; the gate is refused instead.
|
||
#[cfg(feature = "quantum")]
|
||
#[test]
|
||
fn a_phase_on_forty_controls_does_not_overflow() {
|
||
use crate::quantum::{BaseGate, Circuit, Gate};
|
||
let mut out = RotCircuit::new(41);
|
||
let controls: Vec<u32> = (0..40).collect();
|
||
assert!(matches!(controlled_phase(&mut out, &controls, 40, 1.0), Err(SpdError::Unsupported(_))));
|
||
assert!(out.ops.is_empty(), "nothing is emitted before the refusal");
|
||
// At the edge: 31 controls and the target make 32 qubits, refused.
|
||
assert!(matches!(controlled_phase(&mut out, &(0..31).collect::<Vec<u32>>(), 31, 1.0), Err(SpdError::Unsupported(_))));
|
||
// Through the adapter, for every multi-control gate it expands.
|
||
for base in [BaseGate::X, BaseGate::Y, BaseGate::Z, BaseGate::T, BaseGate::P] {
|
||
let mut c = Circuit::new(41);
|
||
let param = if base == BaseGate::P { 5 } else { 0 };
|
||
c.ops.push(Gate { base, controls: (0..40).collect(), target: 40, param });
|
||
c.validate().unwrap();
|
||
assert!(matches!(from_dyadic(&c), Err(SpdError::Unsupported(_))), "{base:?}");
|
||
}
|
||
// A small gate still expands: two controls, 2^3 − 1 rotations.
|
||
let mut out = RotCircuit::new(3);
|
||
controlled_phase(&mut out, &[0, 1], 2, 1.0).unwrap();
|
||
assert_eq!(out.ops.len(), 7);
|
||
}
|
||
|
||
#[cfg(feature = "quantum")]
|
||
#[test]
|
||
fn every_dyadic_algorithm_agrees_with_the_byte_exact_core() {
|
||
use crate::quantum::Circuit;
|
||
// Random dyadic circuits with every base gate and up to three controls:
|
||
// ⟨Z_q⟩ from the byte-exact simulator's probabilities and from the
|
||
// rotation form agree.
|
||
let mut rng = Rng(31);
|
||
for trial in 0..150 {
|
||
let n = 3 + rng.below(4) as u8;
|
||
let mut c = Circuit::new(n);
|
||
for _ in 0..25 {
|
||
let t = rng.below(u64::from(n)) as u8;
|
||
let mut ctrls = Vec::new();
|
||
for _ in 0..rng.below(4) {
|
||
let q = rng.below(u64::from(n)) as u8;
|
||
if q != t && !ctrls.contains(&q) {
|
||
ctrls.push(q);
|
||
}
|
||
}
|
||
let (base, param) = match rng.below(9) {
|
||
0 => (crate::quantum::BaseGate::H, 0),
|
||
1 => (crate::quantum::BaseGate::X, 0),
|
||
2 => (crate::quantum::BaseGate::Y, 0),
|
||
3 => (crate::quantum::BaseGate::Z, 0),
|
||
4 => (crate::quantum::BaseGate::S, 0),
|
||
5 => (crate::quantum::BaseGate::Sdg, 0),
|
||
6 => (crate::quantum::BaseGate::T, 0),
|
||
7 => (crate::quantum::BaseGate::Tdg, 0),
|
||
_ => (crate::quantum::BaseGate::P, 1 + rng.below(6) as u16),
|
||
};
|
||
if base == crate::quantum::BaseGate::H {
|
||
ctrls.clear();
|
||
}
|
||
c.ops.push(crate::quantum::Gate { base, controls: ctrls, target: t, param });
|
||
}
|
||
let sv = c.simulate().unwrap();
|
||
let w = sv.prob_weights();
|
||
let total: f64 = w.iter().map(|&x| x as f64).sum();
|
||
let rc = from_dyadic(&c).unwrap();
|
||
for q in 0..n {
|
||
let want: f64 = w.iter().enumerate().map(|(j, &x)| if (j >> q) & 1 == 0 { x as f64 } else { -(x as f64) }).sum::<f64>() / total;
|
||
let got = expect_z(&rc, u32::from(q), &SpdConfig { threshold: 0.0, max_weight: None }).unwrap();
|
||
assert!((got.value - want).abs() < 1e-6, "trial {trial} q {q}: {} vs {want}", got.value);
|
||
}
|
||
}
|
||
}
|
||
}
|