open-standards/wai/wai-quantum/src/quantum_spd.rs

//! 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.
//! - **Threaded, with the same bits.** A large sum is updated by a pool of
//!   workers, started once per run, each over pieces of the table. A run of
//!   consecutive Clifford gates is applied in one pass, each term through
//!   the gates in order. A rotation's new terms have distinct strings, so each
//!   existing coefficient receives at most one addition. The rest are
//!   appended in table order, and dropped magnitudes are summed in table
//!   order. The thread count therefore changes the time and never a bit. On a
//!   16-core desktop processor the 127-qubit kicked-Ising run (20 steps,
//!   threshold 10⁻⁴, 2.2 million terms at its peak) takes 9.4 s, against 49 s
//!   on one thread.

use std::fmt;
use std::sync::atomic::AtomicU32;
use std::sync::atomic::Ordering::Relaxed;

use crate::pool::Pool;

pub use crate::repro::sin_cos;

/// The propagation with its terms on a device, with the CPU's bits
/// ([`device::expectation_on`]).
pub mod device;

// ---------------------------------------------------------------------------
// 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.
///
/// Large tables are updated by a pool of workers, each over a contiguous
/// share of the terms. Where a term sits in the hash index depends on which
/// thread inserted it first, but the index only finds terms; the order of the
/// terms and every coefficient's arithmetic are the sequential ones, so the
/// result is the same bits at any thread count.
struct Terms<const W: usize> {
    keys: Vec<Key<W>>,
    coef: Vec<f64>,
    index: Vec<AtomicU32>,
    stale: bool,
    /// Slots whose coefficient is `0.0`.
    dead: usize,
}

const EMPTY: u32 = u32::MAX;

/// Below this many terms a table is updated on one thread. Tests share far
/// smaller tables, so that every test runs the threaded paths too.
const PARALLEL_TERMS: usize = if cfg!(test) { 64 } else { 1 << 15 };

/// The coefficients, written by several threads at once in [`rotate`]'s
/// merge. Each new term's string is distinct from every other's, so each
/// coefficient is written by one thread at most.
#[derive(Clone, Copy)]
struct Coefs(*mut f64);

// SAFETY: see `Coefs`: the writes go to distinct indices.
unsafe impl Send for Coefs {}
unsafe impl Sync for Coefs {}

impl Coefs {
    // Taken through a method, so a closure captures the whole `Coefs` (which
    // is `Sync`) rather than its pointer (which is not).
    fn ptr(self) -> *mut f64 {
        self.0
    }
}

impl<const W: usize> Terms<W> {
    fn new() -> Self {
        Terms { keys: Vec::new(), coef: Vec::new(), index: (0..16).map(|_| AtomicU32::new(EMPTY)).collect(), stale: false, dead: 0 }
    }

    /// Shares for a pass over this table: one below [`PARALLEL_TERMS`],
    /// otherwise four a worker (so a busy core holds up less of the pass),
    /// each at least a quarter of that size.
    fn parts(&self, pool: &Pool) -> usize {
        if self.keys.len() < PARALLEL_TERMS || pool.threads() == 1 { 1 } else { (self.keys.len() / (PARALLEL_TERMS / 4)).clamp(1, 4 * pool.threads()) }
    }

    /// The slot of `k`, if present. The index must be current.
    fn find(&self, k: &Key<W>) -> Option<usize> {
        let cap = self.index.len();
        let mut slot = (k.hash() & (cap as u64 - 1)) as usize;
        loop {
            let i = self.index[slot].load(Relaxed);
            if i == EMPTY {
                return None;
            }
            if self.keys[i as usize] == *k {
                return Some(i as usize);
            }
            slot = (slot + 1) & (cap - 1);
        }
    }

    /// Enter term `i` in an index that does not hold it, by compare-and-swap
    /// on the first empty slot, so threads can enter different terms at once.
    fn enter(index: &[AtomicU32], k: &Key<W>, i: u32) {
        let cap = index.len();
        let mut slot = (k.hash() & (cap as u64 - 1)) as usize;
        while index[slot].compare_exchange(EMPTY, i, Relaxed, Relaxed).is_err() {
            slot = (slot + 1) & (cap - 1);
        }
    }

    /// 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;
    }

    /// Rebuild the index for at least `extra` more terms than the table holds.
    /// An index of the right size is cleared and reused; one far too large
    /// for the table, or too small, is replaced.
    fn reindex_for(&mut self, extra: usize, pool: &Pool) {
        let mut cap = 16usize;
        while cap < 2 * (self.keys.len() + extra) + 2 {
            cap *= 2;
        }
        let parts = self.parts(pool);
        if self.index.len() < cap || self.index.len() > 8 * cap {
            self.index = (0..cap).map(|_| AtomicU32::new(EMPTY)).collect();
        } else {
            let index = &self.index;
            pool.split(index.len(), parts, &|r| {
                for slot in &index[r] {
                    slot.store(EMPTY, Relaxed);
                }
            });
        }
        let (keys, index) = (&self.keys, &self.index);
        pool.split(keys.len(), parts, &|r| {
            for i in r {
                Self::enter(index, &keys[i], i as u32);
            }
        });
        self.stale = false;
    }

    fn reindex(&mut self, pool: &Pool) {
        self.reindex_for(0, pool);
    }

    /// Add `c` to the coefficient of `k`, inserting it at the end if absent.
    fn add(&mut self, k: Key<W>, c: f64, pool: &Pool) {
        debug_assert!(!self.stale);
        if 2 * (self.keys.len() + 1) > self.index.len() {
            self.reindex(pool);
        }
        let cap = self.index.len();
        let mut slot = (k.hash() & (cap as u64 - 1)) as usize;
        loop {
            let i = self.index[slot].load(Relaxed);
            if i == EMPTY {
                self.index[slot].store(self.keys.len() as u32, Relaxed);
                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|`, added in table order.
    fn truncate(&mut self, threshold: f64, max_weight: Option<u32>, pool: &Pool) -> f64 {
        let parts = self.parts(pool);
        let per = self.keys.len().div_ceil(parts).max(1);
        let keys = &self.keys;
        // Each share zeroes its dropped terms and keeps their magnitudes in
        // order; the magnitudes are then summed in table order.
        let drop_share = |at: usize, coef: &mut [f64]| -> Vec<f64> {
            let mut out = Vec::new();
            for (j, c) in coef.iter_mut().enumerate() {
                if *c == 0.0 {
                    continue;
                }
                if c.abs() < threshold || max_weight.is_some_and(|m| keys[at + j].weight() > m) {
                    out.push(c.abs());
                    *c = 0.0;
                }
            }
            out
        };
        let shares = pool.chunks(&mut self.coef, per, &drop_share);
        let mut dropped = 0.0;
        for v in shares.iter().flatten() {
            dropped += v;
            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),
    }
}

/// A gate that maps each string to one signed string: a fixed Clifford gate,
/// or a quarter-turn rotation about a string.
enum Clifford<'a, const W: usize> {
    Gate(&'a Op),
    Quarter(Key<W>, u8),
}

/// Conjugate one term by a fixed Clifford gate.
#[inline]
fn clifford_term<const W: usize>(op: &Op, k: &mut Key<W>, c: &mut f64) {
    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 {
        *c = -*c;
    }
}

/// Conjugate one term 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`.
#[inline]
fn quarter_term<const W: usize>(p: &Key<W>, k: u8, q: &mut Key<W>, c: &mut f64) {
    if q.commutes(p) {
        return;
    }
    if k == 2 {
        *c = -*c;
        return;
    }
    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;
    }
    *q = r;
    if neg {
        *c = -*c;
    }
}

/// Apply consecutive Clifford gates to every term in one pass: each term
/// goes through the gates in order. Conjugation by a unitary is a bijection on
/// Pauli strings, so no two terms collide, and each term's result is what the
/// gates one pass at a time would give it.
fn cliffords<const W: usize>(t: &mut Terms<W>, gates: &[Clifford<'_, W>], pool: &Pool) {
    if t.dead > 0 {
        t.compact();
    }
    let parts = t.parts(pool);
    map_terms(&mut t.keys, &mut t.coef, pool, parts, &|k, c| {
        for g in gates {
            match g {
                Clifford::Gate(op) => clifford_term(op, k, c),
                Clifford::Quarter(p, q) => quarter_term(p, *q, k, c),
            }
        }
    });
    t.stale = true;
}

/// Apply `f` to every term, in `parts` shares of the table.
fn map_terms<const W: usize>(keys: &mut [Key<W>], coef: &mut [f64], pool: &Pool, parts: usize, f: &(dyn Fn(&mut Key<W>, &mut f64) + Sync)) {
    let per = keys.len().div_ceil(parts.max(1)).max(1);
    pool.chunks2(keys, coef, per, &|ks, cs| {
        for (k, c) in ks.iter_mut().zip(cs.iter_mut()) {
            f(k, c);
        }
    });
}

/// 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.
///
/// On a large table the work is shared: each thread scans its share and
/// lists its new terms in order, and the lists are joined in share order.
/// Multiplying by `P` maps distinct strings to distinct strings, so no two new
/// terms share a string. Each existing coefficient therefore receives at most
/// one addition, and the strings not yet in the table are appended in list
/// order. The result is the sequential merge's, bit for bit.
fn rotate<const W: usize>(t: &mut Terms<W>, p: &Key<W>, sin: f64, cos: f64, pending: &mut Vec<(Key<W>, f64)>, pool: &Pool) {
    if t.stale {
        t.reindex(pool);
    }
    pending.clear();
    let parts = t.parts(pool);
    if parts == 1 {
        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, pool);
        }
        return;
    }
    // The scan: each share scales its anticommuting terms and lists the new ones.
    let per = t.keys.len().div_ceil(parts);
    let keys = &t.keys;
    let shares: Vec<Vec<(Key<W>, f64)>> = pool.chunks(&mut t.coef, per, &|at, coef| {
        let mut out = Vec::new();
        for (j, c) in coef.iter_mut().enumerate() {
            if *c == 0.0 {
                continue;
            }
            let q = keys[at + j];
            if q.commutes(p) {
                continue;
            }
            let (r, e) = q.mul(p);
            let v = *c * sin;
            out.push((r, if (e + 3) & 3 == 2 { -v } else { v }));
            *c *= cos;
        }
        out
    });
    // The merge, share by share: add each new term to its string's
    // coefficient where the string is in the table, and list the rest in order.
    let coefs = Coefs(t.coef.as_mut_ptr());
    let table: &Terms<W> = t;
    let len = table.keys.len();
    let merged: Vec<(isize, Vec<usize>)> = pool.split(shares.len(), shares.len(), &|r| {
        let (mut revived_less_killed, mut absent) = (0isize, Vec::new());
        let share = &shares[r.start];
        for (j, &(k, v)) in share.iter().enumerate() {
            match table.find(&k) {
                Some(i) => {
                    debug_assert!(i < len);
                    // SAFETY: `i` is a slot of the table, and no other new
                    // term has this string, so no other thread writes slot `i`.
                    let slot = unsafe { &mut *coefs.ptr().add(i) };
                    let before = *slot;
                    let after = before + v;
                    *slot = after;
                    match (before == 0.0, after == 0.0) {
                        (true, false) => revived_less_killed += 1,
                        (false, true) => revived_less_killed -= 1,
                        _ => {}
                    }
                }
                None => absent.push(j),
            }
        }
        (revived_less_killed, absent)
    });
    let mut fresh = 0;
    for (d, absent) in &merged {
        t.dead = (t.dead as isize - d) as usize;
        fresh += absent.len();
    }
    if 2 * (t.keys.len() + fresh) > t.index.len() {
        t.reindex_for(fresh, pool);
    }
    let start = t.keys.len();
    t.keys.reserve(fresh);
    t.coef.reserve(fresh);
    for (share, (_, absent)) in shares.iter().zip(&merged) {
        for &j in absent {
            let (k, v) = share[j];
            t.keys.push(k);
            t.coef.push(v);
            if v == 0.0 {
                t.dead += 1;
            }
        }
    }
    let (keys, index) = (&t.keys, &t.index);
    pool.split(keys.len() - start, parts, &|r| {
        for i in r {
            Terms::<W>::enter(index, &keys[start + i], (start + i) as u32);
        }
    });
}

// ---------------------------------------------------------------------------
// 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.
/// Large sums are updated on every available thread; the result is the
/// same bits at any thread count ([`expectation_threaded`]).
pub fn expectation(
    circuit: &RotCircuit,
    observable: &[(f64, PauliString)],
    ones: &[u32],
    cfg: &SpdConfig,
) -> Result<SpdResult, SpdError> {
    #[cfg(not(target_arch = "wasm32"))]
    let threads = std::thread::available_parallelism().map_or(1, |t| t.get());
    #[cfg(target_arch = "wasm32")]
    let threads = 1;
    expectation_threaded(circuit, observable, ones, cfg, threads)
}

/// [`expectation`] on at most `threads` threads (wasm32 runs one). The thread
/// count changes the time, never a bit of the result.
pub fn expectation_threaded(
    circuit: &RotCircuit,
    observable: &[(f64, PauliString)],
    ones: &[u32],
    cfg: &SpdConfig,
    threads: usize,
) -> Result<SpdResult, SpdError> {
    expectation_capped(circuit, observable, ones, cfg, threads, usize::MAX).map(|(r, _)| r)
}

/// [`expectation_threaded`] holding at most `max_terms` terms: whenever a
/// rotation leaves more, the threshold is doubled and the sum truncated again
/// until it fits, and every coefficient dropped that way is added to the
/// bound like any other. A run that never reaches the cap is the uncapped
/// run, bit for bit. Returns the result and the threshold the run ended at.
pub fn expectation_capped(
    circuit: &RotCircuit,
    observable: &[(f64, PauliString)],
    ones: &[u32],
    cfg: &SpdConfig,
    threads: usize,
    max_terms: usize,
) -> Result<(SpdResult, f64), 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, threads, max_terms),
        2 => run::<2>(circuit, observable, ones, cfg, threads, max_terms),
        3 | 4 => run::<4>(circuit, observable, ones, cfg, threads, max_terms),
        5..=8 => run::<8>(circuit, observable, ones, cfg, threads, max_terms),
        _ => run::<16>(circuit, observable, ones, cfg, threads, max_terms),
    })
}

/// `⟨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,
    threads: usize,
    max_terms: usize,
) -> (SpdResult, f64) {
    crate::pool::scoped(threads, |pool| {
        let (r, t, _) = propagate::<W>(circuit, observable, ones, cfg, pool, max_terms);
        (r, t)
    })
}

/// [`expectation_capped`] on one thread, with the table it ended with: every
/// slot's string (in [`device::key_of`]'s layout) and coefficient, in order,
/// dead slots included. What a device's certificate is held to.
#[cfg(any(test, feature = "quantum_fabric"))]
pub(crate) fn expectation_table(circuit: &RotCircuit, observable: &[(f64, PauliString)], ones: &[u32], cfg: &SpdConfig, max_terms: usize) -> Result<(SpdResult, f64, Vec<u64>, Vec<f64>), SpdError> {
    circuit.validate()?;
    for (_, s) in observable {
        check_string(s, circuit.n)?;
    }
    if let Some(&q) = ones.iter().find(|&&q| q >= circuit.n) {
        return Err(SpdError::QubitOutOfRange { qubit: q, n: circuit.n });
    }
    fn table<const W: usize>(circuit: &RotCircuit, observable: &[(f64, PauliString)], ones: &[u32], cfg: &SpdConfig, max_terms: usize) -> (SpdResult, f64, Vec<u64>, Vec<f64>) {
        crate::pool::scoped(1, |pool| {
            let (r, threshold, t) = propagate::<W>(circuit, observable, ones, cfg, pool, max_terms);
            let keys = t.keys.iter().flat_map(|k| k.x.iter().chain(k.z.iter()).copied()).collect();
            (r, threshold, keys, t.coef)
        })
    }
    Ok(match circuit.n.div_ceil(64) {
        1 => table::<1>(circuit, observable, ones, cfg, max_terms),
        2 => table::<2>(circuit, observable, ones, cfg, max_terms),
        3 | 4 => table::<4>(circuit, observable, ones, cfg, max_terms),
        5..=8 => table::<8>(circuit, observable, ones, cfg, max_terms),
        _ => table::<16>(circuit, observable, ones, cfg, max_terms),
    })
}

fn propagate<const W: usize>(circuit: &RotCircuit, observable: &[(f64, PauliString)], ones: &[u32], cfg: &SpdConfig, pool: &Pool, max_terms: usize) -> (SpdResult, f64, Terms<W>) {
    let mut t = Terms::<W>::new();
    for (c, s) in observable {
        t.add(Key::from_string(s), *c, pool);
    }
    let mut bound = t.truncate(0.0, None, pool);
    let mut peak = t.live();
    let mut work = 0u64;
    let mut pending = Vec::new();
    // Consecutive Clifford gates wait in `batch` and are applied in one pass,
    // before the next rotation and at the end. A Clifford gate keeps the
    // number of terms, so the work count and the peak are as gate by gate.
    let mut batch: Vec<Clifford<'_, W>> = Vec::new();
    let mut threshold = cfg.threshold;
    for op in circuit.ops.iter().rev() {
        work += t.live() as u64;
        match op {
            Op::Rot { axis, theta } => {
                if !batch.is_empty() {
                    cliffords(&mut t, &batch, pool);
                    batch.clear();
                }
                let (s, c) = sin_cos(*theta);
                rotate(&mut t, &Key::from_string(axis), s, c, &mut pending, pool);
                bound += t.truncate(threshold, cfg.max_weight, pool);
                // Over the cap: a doubled threshold until the sum fits.
                while t.live() > max_terms {
                    threshold = if threshold > 0.0 { threshold * 2.0 } else { f64::MIN_POSITIVE };
                    bound += t.truncate(threshold, cfg.max_weight, pool);
                }
            }
            Op::QuarterRot { axis, k } => batch.push(Clifford::Quarter(Key::from_string(axis), *k)),
            other => batch.push(Clifford::Gate(other)),
        }
        peak = peak.max(t.live());
    }
    if !batch.is_empty() {
        cliffords(&mut t, &batch, pool);
    }
    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] };
        }
    }
    let live = t.live();
    (SpdResult { value, truncation_bound: bound, peak_terms: peak, final_terms: live, term_updates: work }, threshold, t)
}

// ---------------------------------------------------------------------------
// 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 threads_do_not_change_a_bit() {
        let mut rng = Rng(77);
        for trial in 0..6 {
            let n = [9u32, 40, 70, 130, 200, 300][trial];
            let c = random_circuit(&mut rng, n, 120);
            let obs = vec![(0.7, rng.string(n, 4)), (-0.4, rng.string(n, 3)), (0.25, rng.string(n, 2))];
            // Exactly (threshold 0) only where the sum stays small.
            for threshold in if n <= 9 { vec![0.0, 1e-3] } else { vec![1e-5, 1e-3] } {
                let cfg = SpdConfig { threshold, max_weight: None };
                let one = expectation_threaded(&c, &obs, &[1, 2], &cfg, 1).unwrap();
                for threads in [2usize, 3, 8] {
                    let many = expectation_threaded(&c, &obs, &[1, 2], &cfg, threads).unwrap();
                    assert_eq!((many.value.to_bits(), many.truncation_bound.to_bits(), many.peak_terms, many.final_terms, many.term_updates), (one.value.to_bits(), one.truncation_bound.to_bits(), one.peak_terms, one.final_terms, one.term_updates), "trial {trial}, threshold {threshold}, {threads} threads");
                }
            }
        }
        // The utility circuit, past the size at which shares go to threads.
        let ki = kicked_ising(127, &heavy_hex_127(), 0.7, 6);
        let cfg = SpdConfig { threshold: 1e-4, max_weight: None };
        let one = expectation_threaded(&ki, &[(1.0, vec![(62, Pauli::Z)])], &[], &cfg, 1).unwrap();
        let many = expectation_threaded(&ki, &[(1.0, vec![(62, Pauli::Z)])], &[], &cfg, 6).unwrap();
        assert!(one.peak_terms > 64);
        assert_eq!((one.value.to_bits(), one.truncation_bound.to_bits(), one.peak_terms), (many.value.to_bits(), many.truncation_bound.to_bits(), many.peak_terms));
    }

    #[test]
    fn a_term_cap_raises_the_threshold_and_the_bound_still_holds() {
        let ki = kicked_ising(127, &heavy_hex_127(), 0.7, 5);
        let obs = [(1.0, vec![(62, Pauli::Z)])];
        let cfg = SpdConfig { threshold: 1e-5, max_weight: None };
        let (free, t0) = expectation_capped(&ki, &obs, &[], &cfg, 4, usize::MAX).unwrap();
        assert_eq!(t0, 1e-5);
        // A cap above the peak changes nothing.
        let (same, t1) = expectation_capped(&ki, &obs, &[], &cfg, 4, free.peak_terms + 1).unwrap();
        assert_eq!((same.value.to_bits(), same.truncation_bound.to_bits(), t1), (free.value.to_bits(), free.truncation_bound.to_bits(), 1e-5));
        // A cap a tenth of the peak raises the threshold; the capped value
        // stays within its own bound of the far more exact one.
        let (capped, t2) = expectation_capped(&ki, &obs, &[], &cfg, 4, free.peak_terms / 10).unwrap();
        assert!(t2 > 1e-5 && capped.peak_terms <= free.peak_terms);
        assert!((capped.value - free.value).abs() <= capped.truncation_bound + free.truncation_bound + 1e-12);
        // A small register, exact reference: the bound covers the true error.
        let mut rng = Rng(5);
        let c = random_circuit(&mut rng, 6, 60);
        let o = vec![(1.0, vec![(0u32, Pauli::Z), (3, Pauli::X)])];
        let exact = dense_expect(&c, &o[0].1, &[]);
        let (r, _) = expectation_capped(&c, &o, &[], &SpdConfig { threshold: 0.0, max_weight: None }, 2, 40).unwrap();
        assert!((r.value - exact).abs() <= r.truncation_bound + 1e-12, "{} vs {exact}, bound {}", r.value, r.truncation_bound);
    }

    #[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);
            }
        }
    }
}