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

//! Wide stabilizer simulation — `wai.quantum.tableau`.
//!
//! A Clifford circuit on `n` qubits, held as the Aaronson–Gottesman tableau
//! (`n` destabilizers and `n` stabilizers, each a signed Pauli string), packed
//! qubit by qubit: the X bits of one qubit across all `2n` rows are one run of
//! words, and so are its Z bits. A gate then touches one or two such runs, a
//! few word operations per 64 rows, so a circuit costs `O(gates · n / 64)`.
//!
//! Two readouts, both exact:
//!
//! - **Outcomes** ([`Tableau::support`]). The computational-basis outcomes of a
//!   stabilizer state are uniform over an affine subspace `x₀ + span(v₁…v_k)`
//!   of `GF(2)^n`. Gaussian elimination on the stabilizers (with their signs
//!   tracked through every product) leaves `k` generators with independent X
//!   parts and `n − k` that are `±Z`-strings; each of the latter is one linear
//!   constraint `z·x = r` on the outcomes, and the constraints' solutions are
//!   the subspace. It is found once, in `O(n³/64)`; a shot is then `x₀` plus a
//!   random combination of the `v`, `O(k · m / 64)` for `m` measured qubits
//!   ([`Sampler`]).
//! - **Expectation values** ([`Tableau::expectation`]). For a Pauli string `P`,
//!   `⟨P⟩` is `0` if `P` anticommutes with a stabilizer; otherwise `P` is, up to
//!   sign, the product of the stabilizers whose destabilizers anticommute with
//!   it, and the sign of that product is `⟨P⟩`. `O(n²/64)`.
//!
//! Signs follow Aaronson and Gottesman's `g` (the power of `i` when two Pauli
//! letters multiply), evaluated a word at a time.

/// The most qubits a tableau holds: the tableau is `4n²` bits, so 32768
/// qubits take 512 MiB, and the outcome subspace about as much again.
pub const MAX_QUBITS: u32 = 1 << 15;

/// One Clifford gate.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum Clifford {
    H(u32),
    S(u32),
    Sdg(u32),
    X(u32),
    Y(u32),
    Z(u32),
    /// √X.
    SX(u32),
    /// √X†.
    SXdg(u32),
    /// `CX(control, target)`.
    CX(u32, u32),
    /// `CY(control, target)`.
    CY(u32, u32),
    CZ(u32, u32),
    Swap(u32, u32),
}

#[derive(Clone, Debug, PartialEq, Eq)]
pub enum TableauError {
    TooManyQubits(u32),
    QubitOutOfRange { qubit: u32, n: u32 },
    /// A two-qubit gate or a Pauli string naming one qubit twice.
    RepeatedQubit(u32),
    /// A Pauli factor that is not `X`, `Y`, `Z` or `I`.
    BadPauli(char),
}

impl std::fmt::Display for TableauError {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            TableauError::TooManyQubits(n) => write!(f, "{n} qubits; a tableau holds at most {MAX_QUBITS}"),
            TableauError::QubitOutOfRange { qubit, n } => write!(f, "qubit {qubit} of a {n}-qubit tableau"),
            TableauError::RepeatedQubit(q) => write!(f, "qubit {q} named twice"),
            TableauError::BadPauli(c) => write!(f, "{c:?} is not a Pauli factor"),
        }
    }
}

impl std::error::Error for TableauError {}

/// The tableau of a stabilizer state: rows `0..n` the destabilizers, rows
/// `n..2n` the stabilizers, packed qubit by qubit.
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct Tableau {
    n: usize,
    /// Words per qubit's run of `2n` row bits.
    rw: usize,
    xs: Vec<u64>,
    zs: Vec<u64>,
    /// The rows' signs, `2n` bits.
    rs: Vec<u64>,
}

fn words(bits: usize) -> usize {
    bits.div_ceil(64)
}

fn get(v: &[u64], i: usize) -> bool {
    v[i / 64] >> (i % 64) & 1 == 1
}

fn set(v: &mut [u64], i: usize) {
    v[i / 64] |= 1 << (i % 64);
}

fn splitmix(s: &mut u64) -> u64 {
    *s = s.wrapping_add(0x9E37_79B9_7F4A_7C15);
    let mut z = *s;
    z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
    z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
    z ^ (z >> 31)
}

/// Rows held one after another, `wn` words of X then of Z each, with signs.
struct Rows {
    wn: usize,
    x: Vec<u64>,
    z: Vec<u64>,
    r: Vec<bool>,
}

impl Rows {
    fn new(count: usize, n: usize) -> Rows {
        let wn = words(n);
        Rows { wn, x: vec![0; count * wn], z: vec![0; count * wn], r: vec![false; count] }
    }

    fn xr(&self, i: usize) -> &[u64] {
        &self.x[i * self.wn..(i + 1) * self.wn]
    }

    fn zr(&self, i: usize) -> &[u64] {
        &self.z[i * self.wn..(i + 1) * self.wn]
    }

    fn swap(&mut self, a: usize, b: usize) {
        if a == b {
            return;
        }
        for w in 0..self.wn {
            self.x.swap(a * self.wn + w, b * self.wn + w);
            self.z.swap(a * self.wn + w, b * self.wn + w);
        }
        self.r.swap(a, b);
    }

    /// Row `h` ← row `i` · row `h`, sign included (the rows commute).
    fn mul_into(&mut self, h: usize, i: usize) {
        let (wn, mut acc) = (self.wn, 0i64);
        for w in 0..wn {
            let (x1, z1) = (self.x[i * wn + w], self.z[i * wn + w]);
            let (x2, z2) = (self.x[h * wn + w], self.z[h * wn + w]);
            acc += phase_word(x1, z1, x2, z2);
            self.x[h * wn + w] = x1 ^ x2;
            self.z[h * wn + w] = z1 ^ z2;
        }
        let total = (2 * self.r[h] as i64 + 2 * self.r[i] as i64 + acc).rem_euclid(4);
        self.r[h] = total == 2;
    }
}

/// `Σ g(x1, z1, x2, z2)` over one word's 64 qubits: the power of `i` the
/// product of the two Paulis' letters carries.
fn phase_word(x1: u64, z1: u64, x2: u64, z2: u64) -> i64 {
    let y1 = x1 & z1;
    let xo = x1 & !z1;
    let zo = !x1 & z1;
    let plus = (y1 & z2 & !x2) | (xo & z2 & x2) | (zo & x2 & !z2);
    let minus = (y1 & x2 & !z2) | (xo & z2 & !x2) | (zo & x2 & z2);
    plus.count_ones() as i64 - minus.count_ones() as i64
}

impl Tableau {
    /// `|0…0⟩` on `n` qubits.
    pub fn new(n: u32) -> Result<Tableau, TableauError> {
        if n > MAX_QUBITS {
            return Err(TableauError::TooManyQubits(n));
        }
        let n = n as usize;
        let rw = words(2 * n);
        let mut t = Tableau { n, rw, xs: vec![0; n * rw], zs: vec![0; n * rw], rs: vec![0; rw] };
        for q in 0..n {
            // Destabilizer q is X_q; stabilizer q is Z_q.
            set(&mut t.xs[q * rw..(q + 1) * rw], q);
            set(&mut t.zs[q * rw..(q + 1) * rw], n + q);
        }
        Ok(t)
    }

    /// The circuit `ops` on `n` qubits from `|0…0⟩`.
    pub fn run(n: u32, ops: &[Clifford]) -> Result<Tableau, TableauError> {
        let mut t = Tableau::new(n)?;
        for &op in ops {
            t.apply(op)?;
        }
        Ok(t)
    }

    pub fn qubits(&self) -> u32 {
        self.n as u32
    }

    fn check(&self, q: u32) -> Result<usize, TableauError> {
        if (q as usize) < self.n { Ok(q as usize) } else { Err(TableauError::QubitOutOfRange { qubit: q, n: self.n as u32 }) }
    }

    fn pair(&self, a: u32, b: u32) -> Result<(usize, usize), TableauError> {
        let (a, b) = (self.check(a)?, self.check(b)?);
        if a == b {
            return Err(TableauError::RepeatedQubit(a as u32));
        }
        Ok((a, b))
    }

    fn h(&mut self, a: usize) {
        let rw = self.rw;
        for w in 0..rw {
            let (x, z) = (self.xs[a * rw + w], self.zs[a * rw + w]);
            self.rs[w] ^= x & z;
            self.xs[a * rw + w] = z;
            self.zs[a * rw + w] = x;
        }
    }

    fn s(&mut self, a: usize) {
        let rw = self.rw;
        for w in 0..rw {
            let (x, z) = (self.xs[a * rw + w], self.zs[a * rw + w]);
            self.rs[w] ^= x & z;
            self.zs[a * rw + w] = z ^ x;
        }
    }

    fn sdg(&mut self, a: usize) {
        let rw = self.rw;
        for w in 0..rw {
            let (x, z) = (self.xs[a * rw + w], self.zs[a * rw + w]);
            self.rs[w] ^= x & !z;
            self.zs[a * rw + w] = z ^ x;
        }
    }

    fn cx(&mut self, a: usize, b: usize) {
        let rw = self.rw;
        for w in 0..rw {
            let (xa, za, xb, zb) = (self.xs[a * rw + w], self.zs[a * rw + w], self.xs[b * rw + w], self.zs[b * rw + w]);
            self.rs[w] ^= xa & zb & !(xb ^ za);
            self.xs[b * rw + w] = xb ^ xa;
            self.zs[a * rw + w] = za ^ zb;
        }
    }

    /// Applies one gate.
    pub fn apply(&mut self, op: Clifford) -> Result<(), TableauError> {
        let rw = self.rw;
        match op {
            Clifford::H(a) => {
                let a = self.check(a)?;
                self.h(a);
            }
            Clifford::S(a) => {
                let a = self.check(a)?;
                self.s(a);
            }
            Clifford::Sdg(a) => {
                let a = self.check(a)?;
                self.sdg(a);
            }
            Clifford::X(a) => {
                let a = self.check(a)?;
                for w in 0..rw {
                    self.rs[w] ^= self.zs[a * rw + w];
                }
            }
            Clifford::Z(a) => {
                let a = self.check(a)?;
                for w in 0..rw {
                    self.rs[w] ^= self.xs[a * rw + w];
                }
            }
            Clifford::Y(a) => {
                let a = self.check(a)?;
                for w in 0..rw {
                    self.rs[w] ^= self.xs[a * rw + w] ^ self.zs[a * rw + w];
                }
            }
            Clifford::SX(a) => {
                let a = self.check(a)?;
                self.h(a);
                self.s(a);
                self.h(a);
            }
            Clifford::SXdg(a) => {
                let a = self.check(a)?;
                self.h(a);
                self.sdg(a);
                self.h(a);
            }
            Clifford::CX(a, b) => {
                let (a, b) = self.pair(a, b)?;
                self.cx(a, b);
            }
            Clifford::CZ(a, b) => {
                let (a, b) = self.pair(a, b)?;
                self.h(b);
                self.cx(a, b);
                self.h(b);
            }
            Clifford::CY(a, b) => {
                let (a, b) = self.pair(a, b)?;
                self.sdg(b);
                self.cx(a, b);
                self.s(b);
            }
            Clifford::Swap(a, b) => {
                let (a, b) = self.pair(a, b)?;
                for w in 0..rw {
                    self.xs.swap(a * rw + w, b * rw + w);
                    self.zs.swap(a * rw + w, b * rw + w);
                }
            }
        }
        Ok(())
    }

    /// Rows `from..to` of the tableau, held row by row.
    fn rows(&self, from: usize, to: usize) -> Rows {
        let (n, rw) = (self.n, self.rw);
        let mut out = Rows::new(to - from, n);
        let wn = out.wn;
        for q in 0..n {
            for (bits, dst) in [(&self.xs, &mut out.x), (&self.zs, &mut out.z)] {
                let col = &bits[q * rw..(q + 1) * rw];
                for (w, &bits_w) in col.iter().enumerate().take(to.div_ceil(64)).skip(from / 64) {
                    let mut word = bits_w;
                    while word != 0 {
                        let row = w * 64 + word.trailing_zeros() as usize;
                        word &= word - 1;
                        if (from..to).contains(&row) {
                            dst[(row - from) * wn + q / 64] |= 1 << (q % 64);
                        }
                    }
                }
            }
        }
        for row in from..to {
            out.r[row - from] = get(&self.rs, row);
        }
        out
    }

    /// `⟨ψ|P|ψ⟩` for the Pauli string `P` (`(qubit, 'X'|'Y'|'Z'|'I')`):
    /// `1`, `−1`, or `0`.
    pub fn expectation(&self, pauli: &[(u32, char)]) -> Result<i8, TableauError> {
        let n = self.n;
        let wn = words(n);
        let (mut px, mut pz) = (vec![0u64; wn], vec![0u64; wn]);
        let mut seen = vec![false; n];
        for &(q, c) in pauli {
            let q = self.check(q)?;
            if seen[q] {
                return Err(TableauError::RepeatedQubit(q as u32));
            }
            seen[q] = true;
            match c.to_ascii_uppercase() {
                'X' => set(&mut px, q),
                'Z' => set(&mut pz, q),
                'Y' => {
                    set(&mut px, q);
                    set(&mut pz, q);
                }
                'I' => {}
                other => return Err(TableauError::BadPauli(other)),
            }
        }
        let rows = self.rows(0, 2 * n);
        let anticommutes = |i: usize| -> bool {
            let (rx, rz) = (rows.xr(i), rows.zr(i));
            let mut c = 0u32;
            for w in 0..wn {
                c += ((px[w] & rz[w]) ^ (pz[w] & rx[w])).count_ones();
            }
            c % 2 == 1
        };
        if (n..2 * n).any(anticommutes) {
            return Ok(0);
        }
        // P is ± the product of the stabilizers whose destabilizers it
        // anticommutes with.
        let mut q = Rows::new(1, n);
        let mut scratch = Rows::new(2, n);
        for i in 0..n {
            if anticommutes(i) {
                scratch.x[..wn].copy_from_slice(&q.x[..wn]);
                scratch.z[..wn].copy_from_slice(&q.z[..wn]);
                scratch.r[0] = q.r[0];
                scratch.x[wn..].copy_from_slice(rows.xr(n + i));
                scratch.z[wn..].copy_from_slice(rows.zr(n + i));
                scratch.r[1] = rows.r[n + i];
                scratch.mul_into(0, 1);
                q.x.copy_from_slice(&scratch.x[..wn]);
                q.z.copy_from_slice(&scratch.z[..wn]);
                q.r[0] = scratch.r[0];
            }
        }
        debug_assert!(q.x == px && q.z == pz, "the product of stabilizers is P up to sign");
        Ok(if q.r[0] { -1 } else { 1 })
    }

    /// The computational-basis outcomes, as the affine subspace they are
    /// uniform over.
    pub fn support(&self) -> Support {
        let n = self.n;
        let mut rows = self.rows(n, 2 * n);
        let wn = rows.wn;
        // Eliminate the X parts: pivot rows first, ±Z-strings after.
        let mut k = 0;
        for q in 0..n {
            let Some(p) = (k..n).find(|&i| get(rows.xr(i), q)) else { continue };
            rows.swap(k, p);
            for i in 0..n {
                if i != k && get(rows.xr(i), q) {
                    rows.mul_into(i, k);
                }
            }
            k += 1;
        }
        // Each ±Z-string is a constraint z·x = r; reduce them to echelon form.
        let mut eqs: Vec<(Vec<u64>, bool)> = (k..n).map(|i| (rows.zr(i).to_vec(), rows.r[i])).collect();
        let mut pivots: Vec<usize> = Vec::new();
        let mut rank = 0;
        for q in 0..n {
            let Some(p) = (rank..eqs.len()).find(|&i| get(&eqs[i].0, q)) else { continue };
            eqs.swap(rank, p);
            let (pz, pr) = eqs[rank].clone();
            for (i, e) in eqs.iter_mut().enumerate() {
                if i != rank && get(&e.0, q) {
                    for (a, b) in e.0.iter_mut().zip(&pz) {
                        *a ^= b;
                    }
                    e.1 ^= pr;
                }
            }
            pivots.push(q);
            rank += 1;
        }
        let mut is_pivot = vec![false; n];
        for &c in &pivots {
            is_pivot[c] = true;
        }
        // x₀: the free variables 0, each pivot variable its equation's sign.
        let mut x0 = vec![0u64; wn];
        for (i, &c) in pivots.iter().enumerate() {
            if eqs[i].1 {
                set(&mut x0, c);
            }
        }
        // One basis vector per free variable.
        let mut basis = Vec::with_capacity(n - rank);
        for f in (0..n).filter(|&c| !is_pivot[c]) {
            let mut v = vec![0u64; wn];
            set(&mut v, f);
            for (i, &c) in pivots.iter().enumerate() {
                if get(&eqs[i].0, f) {
                    set(&mut v, c);
                }
            }
            basis.push(v);
        }
        eqs.truncate(rank);
        Support { n, x0, basis, constraints: eqs }
    }
}

/// The computational-basis outcomes of a stabilizer state: uniform over
/// `x₀ + span(basis)`. Bit `q` of an outcome is qubit `q`.
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct Support {
    n: usize,
    x0: Vec<u64>,
    basis: Vec<Vec<u64>>,
    /// The constraints `z·x = r` the outcomes satisfy, in echelon form.
    constraints: Vec<(Vec<u64>, bool)>,
}

impl Support {
    /// `k`: the outcomes are `2^k` equally likely bit strings.
    pub fn dimension(&self) -> usize {
        self.basis.len()
    }

    /// Whether `bits` (bit `q` = qubit `q`, `n` bits) is an outcome: every
    /// constraint `z·x = r` holds.
    pub fn contains(&self, bits: &[u64]) -> bool {
        self.constraints.iter().all(|(z, r)| {
            let parity = z.iter().zip(bits.iter().chain(std::iter::repeat(&0))).map(|(a, b)| (a & b).count_ones()).sum::<u32>() % 2 == 1;
            parity == *r
        })
    }

    /// A sampler for the qubits `measured`, in that order.
    pub fn sampler(&self, measured: &[u32]) -> Sampler {
        let m = measured.len();
        let wm = words(m);
        let project = |v: &[u64]| -> Vec<u64> {
            let mut out = vec![0u64; wm];
            for (j, &q) in measured.iter().enumerate() {
                if (q as usize) < self.n && get(v, q as usize) {
                    set(&mut out, j);
                }
            }
            out
        };
        let x0 = project(&self.x0);
        let basis: Vec<Vec<u64>> = self.basis.iter().map(|v| project(v)).filter(|v| v.iter().any(|&w| w != 0)).collect();
        Sampler { m, x0, basis }
    }
}

/// Shots of the measured qubits: `x₀` plus a uniformly random combination of
/// the basis, projected onto them.
#[derive(Clone, Debug)]
pub struct Sampler {
    m: usize,
    x0: Vec<u64>,
    basis: Vec<Vec<u64>>,
}

impl Sampler {
    /// The measured qubits' count.
    pub fn width(&self) -> usize {
        self.m
    }

    /// The independent directions the measured outcomes vary in (at most the
    /// support's dimension).
    pub fn dimension(&self) -> usize {
        self.basis.len()
    }

    /// One shot into `out` (bit `j` = the `j`-th measured qubit), drawing from
    /// `rng`.
    pub fn shot(&self, rng: &mut u64, out: &mut Vec<u64>) {
        out.clear();
        out.extend_from_slice(&self.x0);
        let mut bits = 0u64;
        for (j, v) in self.basis.iter().enumerate() {
            if j % 64 == 0 {
                bits = splitmix(rng);
            }
            if bits >> (j % 64) & 1 == 1 {
                for (o, w) in out.iter_mut().zip(v) {
                    *o ^= w;
                }
            }
        }
    }

    /// `shots` shots from `seed`, counted: `(outcome, count)`, most frequent
    /// first, then by outcome.
    pub fn counts(&self, shots: u64, seed: u64) -> Vec<(Vec<u64>, u64)> {
        let mut rng = seed ^ 0x7AB1_EA00_5EED_0001;
        let mut map: std::collections::HashMap<Vec<u64>, u64> = std::collections::HashMap::new();
        let mut buf = Vec::new();
        for _ in 0..shots {
            self.shot(&mut rng, &mut buf);
            *map.entry(buf.clone()).or_insert(0) += 1;
        }
        let mut v: Vec<(Vec<u64>, u64)> = map.into_iter().collect();
        v.sort_by(|a, b| b.1.cmp(&a.1).then(a.0.cmp(&b.0)));
        v
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    fn next(s: &mut u64) -> u64 {
        splitmix(s)
    }

    fn random_clifford(n: u32, ops: usize, seed: u64) -> Vec<Clifford> {
        let mut s = seed;
        (0..ops)
            .map(|_| {
                let a = (next(&mut s) % n as u64) as u32;
                let mut b = (next(&mut s) % n as u64) as u32;
                if b == a {
                    b = (a + 1) % n;
                }
                match next(&mut s) % 12 {
                    0 => Clifford::H(a),
                    1 => Clifford::S(a),
                    2 => Clifford::Sdg(a),
                    3 => Clifford::X(a),
                    4 => Clifford::Y(a),
                    5 => Clifford::Z(a),
                    6 => Clifford::SX(a),
                    7 => Clifford::SXdg(a),
                    8 => Clifford::CX(a, b),
                    9 => Clifford::CY(a, b),
                    10 => Clifford::CZ(a, b),
                    _ => Clifford::Swap(a, b),
                }
            })
            .collect()
    }

    #[test]
    fn the_zero_state_and_a_wide_ghz_state_read_exactly() {
        let t = Tableau::new(5).unwrap();
        assert_eq!(t.support().dimension(), 0);
        assert_eq!(t.expectation(&[(3, 'Z')]).unwrap(), 1);
        assert_eq!(t.expectation(&[(3, 'X')]).unwrap(), 0);
        let n = 2000u32;
        let mut ops = vec![Clifford::H(0)];
        ops.extend((0..n - 1).map(|q| Clifford::CX(q, q + 1)));
        let t = Tableau::run(n, &ops).unwrap();
        let s = t.support();
        assert_eq!(s.dimension(), 1);
        assert_eq!(t.expectation(&[(0, 'Z'), (n - 1, 'Z')]).unwrap(), 1);
        assert_eq!(t.expectation(&[(7, 'Z')]).unwrap(), 0);
        let all_x: Vec<(u32, char)> = (0..n).map(|q| (q, 'X')).collect();
        assert_eq!(t.expectation(&all_x).unwrap(), 1);
        // Y on every qubit of an even GHZ: ⟨Y^{⊗n}⟩ = Re(i^n) = 1 for n ≡ 0 mod 4.
        let all_y: Vec<(u32, char)> = (0..n).map(|q| (q, 'Y')).collect();
        assert_eq!(t.expectation(&all_y).unwrap(), 1);
        let measured: Vec<u32> = (0..n).collect();
        let c = s.sampler(&measured).counts(1000, 3);
        assert_eq!(c.len(), 2);
        assert!(c.iter().all(|(o, k)| *k > 400 && (o.iter().all(|&w| w == 0) || o.iter().map(|w| w.count_ones()).sum::<u32>() == n)));
    }

    #[test]
    fn a_gate_and_its_inverse_leave_the_tableau_as_it_was() {
        let base = Tableau::run(6, &random_clifford(6, 40, 9)).unwrap();
        for a in [
            vec![Clifford::H(2), Clifford::H(2)],
            vec![Clifford::S(1), Clifford::Sdg(1)],
            vec![Clifford::SX(4), Clifford::SXdg(4)],
            vec![Clifford::CX(0, 3), Clifford::CX(0, 3)],
            vec![Clifford::CY(5, 1), Clifford::CY(5, 1)],
            vec![Clifford::CZ(2, 4), Clifford::CZ(4, 2)],
            vec![Clifford::Swap(0, 5), Clifford::Swap(5, 0)],
            vec![Clifford::S(3), Clifford::S(3), Clifford::S(3), Clifford::S(3)],
        ] {
            let mut t = base.clone();
            for op in a.clone() {
                t.apply(op).unwrap();
            }
            assert_eq!(t, base, "{a:?}");
        }
    }

    #[test]
    fn refusals_are_typed() {
        assert_eq!(Tableau::new(MAX_QUBITS + 1).unwrap_err(), TableauError::TooManyQubits(MAX_QUBITS + 1));
        let mut t = Tableau::new(3).unwrap();
        assert_eq!(t.apply(Clifford::H(3)).unwrap_err(), TableauError::QubitOutOfRange { qubit: 3, n: 3 });
        assert_eq!(t.apply(Clifford::CX(1, 1)).unwrap_err(), TableauError::RepeatedQubit(1));
        assert_eq!(t.expectation(&[(0, 'Q')]).unwrap_err(), TableauError::BadPauli('Q'));
        assert_eq!(t.expectation(&[(0, 'Z'), (0, 'X')]).unwrap_err(), TableauError::RepeatedQubit(0));
    }

    /// Times at scale; run with `--ignored --nocapture`.
    #[test]
    #[ignore]
    fn timings_at_scale() {
        for n in [1000u32, 2000, 5000, 10000] {
            let ops = random_clifford(n, 10 * n as usize, n as u64);
            let t0 = std::time::Instant::now();
            let t = Tableau::run(n, &ops).unwrap();
            let gates = t0.elapsed();
            let t1 = std::time::Instant::now();
            let s = t.support();
            let support = t1.elapsed();
            let measured: Vec<u32> = (0..n).collect();
            let t2 = std::time::Instant::now();
            let c = s.sampler(&measured).counts(1000, 1);
            let shots = t2.elapsed();
            let t3 = std::time::Instant::now();
            let e = t.expectation(&[(0, 'Z'), (n / 2, 'X')]).unwrap();
            let expect = t3.elapsed();
            println!("n={n}: {} gates {gates:?}; support dim {} in {support:?}; 1000 shots ({} distinct) {shots:?}; one expectation ({e}) {expect:?}", ops.len(), s.dimension(), c.len());
        }
    }

    #[cfg(feature = "quantum_sv")]
    mod against_the_state_vector {
        use super::*;
        use crate::gates::{Gate, complex};
        use crate::quantum_sv::StateVector;

        fn matrix(op: Clifford) -> Gate {
            let c = |re: f64, im: f64| complex(re, im);
            let o = c(0.0, 0.0);
            let one = c(1.0, 0.0);
            let ctrl = |a: u32, b: u32, u: [crate::gates::C; 4]| Gate { qubits: vec![a, b], matrix: vec![one, o, o, o, o, one, o, o, o, o, u[0], u[1], o, o, u[2], u[3]] };
            match op {
                Clifford::H(a) => Gate::h(a),
                Clifford::S(a) => Gate::s(a),
                Clifford::Sdg(a) => Gate::sdg(a),
                Clifford::X(a) => Gate::x(a),
                Clifford::Y(a) => Gate::y(a),
                Clifford::Z(a) => Gate::z(a),
                Clifford::SX(a) => Gate { qubits: vec![a], matrix: vec![c(0.5, 0.5), c(0.5, -0.5), c(0.5, -0.5), c(0.5, 0.5)] },
                Clifford::SXdg(a) => Gate { qubits: vec![a], matrix: vec![c(0.5, -0.5), c(0.5, 0.5), c(0.5, 0.5), c(0.5, -0.5)] },
                Clifford::CX(a, b) => Gate::cx(a, b),
                Clifford::CY(a, b) => ctrl(a, b, [o, c(0.0, -1.0), c(0.0, 1.0), o]),
                Clifford::CZ(a, b) => Gate::cz(a, b),
                Clifford::Swap(a, b) => Gate::swap(a, b),
            }
        }

        #[test]
        fn outcomes_and_expectations_are_the_state_vectors_on_random_circuits() {
            for seed in 0..40 {
                let n = 2 + (seed % 9) as u32;
                let ops = random_clifford(n, 60, seed);
                let t = Tableau::run(n, &ops).unwrap();
                let mut sv = StateVector::zero(n).unwrap();
                sv.run(&ops.iter().map(|&o| matrix(o)).collect::<Vec<_>>(), 1).unwrap();
                let s = t.support();
                let p = 1.0 / (1u64 << s.dimension()) as f64;
                let mut support_size = 0;
                for x in 0..1usize << n {
                    let prob = sv.re[x] * sv.re[x] + sv.im[x] * sv.im[x];
                    // The state vector's index has qubit 0 as its most significant bit.
                    let mut bits = vec![0u64; 1];
                    for q in 0..n as usize {
                        if x >> (n as usize - 1 - q) & 1 == 1 {
                            bits[0] |= 1 << q;
                        }
                    }
                    let inside = s.contains(&bits);
                    if prob > 1e-9 {
                        support_size += 1;
                        assert!(inside, "seed {seed}: outcome {x:b} has probability {prob} but is outside the support");
                        assert!((prob - p).abs() < 1e-9, "seed {seed}: probability {prob}, expected {p}");
                    } else {
                        assert!(!inside, "seed {seed}: outcome {x:b} is in the support with probability {prob}");
                    }
                }
                assert_eq!(support_size, 1 << s.dimension(), "seed {seed}");
                let mut r = seed ^ 0xabc;
                for _ in 0..20 {
                    let string: Vec<(u32, char)> = (0..n).filter_map(|q| match next(&mut r) % 4 { 0 => None, 1 => Some((q, 'X')), 2 => Some((q, 'Y')), _ => Some((q, 'Z')) }).collect();
                    let want = sv.expectation(&[(1.0, string.clone())]);
                    let got = t.expectation(&string).unwrap() as f64;
                    assert!((got - want).abs() < 1e-9, "seed {seed}: {string:?}: {got} vs {want}");
                }
            }
        }
    }
}