wai-quantum 0.5.0

//! The program a fabric member runs for a circuit, and the CPU executor of
//! exactly that program. A member is correct when its runs equal [`execute`]'s;
//! `execute` is certified against the native simulator for every
//! configuration.
//!
//! A [`Program`] holds the circuit's gates as [`GateWords`], in op-log order,
//! and the steps that apply them. A [`Step::Gate`] is one gate over the whole
//! state. A [`Step::Tile`] is a run of consecutive gates applied tile by tile:
//! the state is cut into tiles of `size` consecutive amplitudes, and each tile
//! goes through every gate of the run, in op-log order, on its own, so that a
//! member can apply a whole run in one dispatch, one group of invocations per
//! tile. A gate joins a run when it acts inside a tile: a gate whose matrix is
//! diagonal (`m01 = m10 = 0`) acts on each amplitude alone, whatever its
//! target, and any other gate acts on pairs that lie inside one tile when its
//! target's bit is below the tile's size. Every other gate is a step of its
//! own. [`Config`] says how a
//! circuit is planned: [`Config::PER_GATE`], one step per gate, is what
//! [`compile`] makes; [`Config::TILED`] makes tile runs of up to 1024
//! amplitudes.
//!
//! The arithmetic is [`lanes::apply_words`]'s in every configuration: the
//! stored lanes are the low 32 bits of the exact values, and a run reports the
//! first op whose exact result left 32 bits, which no valid circuit of at most
//! [`lanes::horizon`] gates can do.
//!
//! # Why every configuration gives the same bits
//!
//! A planner may only move data exactly; it never changes what is computed.
//!
//! - **A tile run is its gates, one after another.** Each gate of a run acts
//!   on amplitudes of one tile only, so the tiles evolve independently, and
//!   within a tile the gates apply in op-log order to the values the previous
//!   gate left. Every pair or amplitude a gate updates is computed from the
//!   same inputs by the same exact rows as in [`lanes::apply_words`], and the
//!   same pairs are updated (a control is never the target, so both members of
//!   a pair pass or fail the control test together). Splitting a run in two is
//!   exact for the same reason.
//! - **A diagonal gate acts per amplitude.** `apply_words` computes
//!   `a_i' = m00·a_i + m01·a_j` with `m01 = 0`, and the exact `fx(0, x)` is 0
//!   for every `x`, so the value is that of `m00·a_i + 0·0`; likewise
//!   `a_j' = 0·0 + m11·a_j`. The range check covers the same four parts.
//! - **The identity is skipped** only when its words are exactly
//!   `[2^30, 0, 0, 0, 0, 0, 2^30, 0]`: `fx(2^30, x) = x` and `fx(0, x) = 0`
//!   for every 32-bit `x`, so the gate changes nothing and cannot leave 32
//!   bits. Circuit/1's `P(19..=32)` is `(2^30 − 1, 0)`, not the identity, and
//!   is never skipped.
//! - **The report is the lowest op index** whose exact result left 32 bits.
//!   The lanes after a gate are the same in every configuration, so the set of
//!   gates whose exact result leaves 32 bits is the same, and so is its least
//!   element.
//!
//! What no configuration does, because each changes bits: negate a rounded
//! product (`fx(−m, x)` differs from `−fx(m, x)` at a tie), multiply two
//! gates' matrices into one, reorder gates (even ones on disjoint qubits), or
//! skip a gate that is not exactly the identity.

use super::lanes::{self, GateWords, Lane};
use crate::quantum::{Circuit, MAX_QUBITS, ONE, QuantumError};

/// Bumped whenever the meaning of a [`Program`] changes. An implementation id
/// that runs more than [`Step::Gate`] includes it
/// ([`super::wgsl::tile_implementation_id`]); the per-gate WGSL kernel's id
/// ([`super::wgsl::implementation_id`]) versions the one step it runs with
/// [`super::wgsl::HOST_LAYOUT`] instead. Revision 2 added [`Step::Tile`].
pub const PLAN_REVISION: u32 = 2;

/// The largest tile, as a power of two: 1024 amplitudes, four for each
/// invocation of a tile kernel's workgroup of 256.
pub const TILE_BITS_MAX: u8 = 10;

/// The identity's words: `m00 = m11 = 2^30`, every other entry zero.
const IDENTITY: [u32; 8] = [ONE as u32, 0, 0, 0, 0, 0, ONE as u32, 0];

/// How a circuit is planned.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub struct Config {
    /// Amplitudes per tile: `2^min(n, tile_bits)`. Zero plans no tiles (one
    /// step per gate). At most [`TILE_BITS_MAX`].
    pub tile_bits: u8,
    /// Gates per tile run, at most; further capped at `2^30 / 2^n` (at least
    /// one), so that one run's work stays bounded on the widest states.
    pub run_max: u32,
    /// Leave out every gate whose words are exactly the identity's.
    pub skip_identity: bool,
}

impl Config {
    /// One [`Step::Gate`] per op, nothing skipped: what [`compile`] makes.
    pub const PER_GATE: Config = Config { tile_bits: 0, run_max: 1, skip_identity: false };
    /// Tiles of up to 1024 amplitudes, runs of up to 1024 gates, the identity
    /// skipped.
    pub const TILED: Config = Config { tile_bits: TILE_BITS_MAX, run_max: 1024, skip_identity: true };

    /// Amplitudes per tile on `n_qubits`: `2^min(n, tile_bits)`, or 0 when
    /// the configuration plans no tiles or the state has no qubits.
    pub fn tile_size(&self, n_qubits: u8) -> u32 {
        if self.tile_bits == 0 || n_qubits == 0 { 0 } else { 1 << n_qubits.min(self.tile_bits).min(TILE_BITS_MAX) }
    }

    /// Gates per tile run on `n_qubits`, at most: `run_max`, capped at
    /// `2^30 / 2^n` and at least one. A member may split a run into several
    /// dispatches to bound each one's work (splitting a run is exact).
    pub fn run_cap(&self, n_qubits: u8) -> u32 {
        let budget = if n_qubits >= 30 { 1 } else { (1u32 << 30) >> n_qubits };
        self.run_max.min(budget).max(1)
    }
}

/// One gate of a program.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub struct PlannedGate {
    pub words: GateWords,
    /// The gate's position in the op-log (what a range report names).
    pub index: u32,
    /// `m01 = m10 = (0, 0)`: the gate acts on each amplitude alone.
    pub diag: bool,
}

impl PlannedGate {
    /// The gate with these words at op-log position `index`; `diag` is read
    /// from the words.
    pub fn new(words: GateWords, index: u32) -> PlannedGate {
        PlannedGate { words, index, diag: words.m[2..6] == [0; 4] }
    }

    /// Whether a member may run the gate: `diag` says what the words say (a
    /// tile kernel trusts it to act per amplitude), and the index is not
    /// `u32::MAX`, which the tile kernel's report keeps for "no op". Every
    /// gate [`PlannedGate::new`] and [`compile`] make is well formed:
    /// `validate` holds an op-log's length, and so every index, below
    /// `u32::MAX`.
    pub fn is_well_formed(&self) -> bool {
        self.diag == (self.words.m[2..6] == [0; 4]) && self.index != u32::MAX
    }

    /// Whether the words are exactly the identity's (whatever the target and
    /// controls): the gate changes nothing and cannot leave 32 bits.
    pub fn is_identity(&self) -> bool {
        self.words.m == IDENTITY
    }

    /// Whether the gate acts inside tiles of `size` amplitudes: diagonal, or
    /// its target's bit below `size`.
    fn fits_tile(&self, size: u32) -> bool {
        self.diag || self.words.tbit < size
    }
}

/// One step of a program.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
#[non_exhaustive]
pub enum Step {
    /// `gates[gate]` over the whole state: one dispatch.
    Gate { gate: u32 },
    /// `gates[first..first + count]`, in order, on each tile of `size`
    /// consecutive amplitudes (a power of two): one dispatch.
    Tile { first: u32, count: u32, size: u32 },
}

/// What a member runs for one circuit from one basis state.
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct Program {
    pub n_qubits: u8,
    /// The basis state the run starts from, already reduced modulo `2^n`.
    pub start: usize,
    pub gates: Vec<PlannedGate>,
    pub steps: Vec<Step>,
}

/// Validates exactly as [`Circuit::simulate_from`] does, with the same errors
/// in the same order (`validate()`, then [`QuantumError::TooManyQubits`]),
/// then plans one [`Step::Gate`] per op, in op-log order:
/// [`compile_with`]`(c, start, &Config::PER_GATE)`.
pub fn compile(c: &Circuit, start: usize) -> Result<Program, QuantumError> {
    compile_with(c, start, &Config::PER_GATE)
}

/// Validates as [`compile`] does, then plans the ops under `cfg`
/// ([`schedule`]).
///
/// Panics if `cfg.tile_bits` is above [`TILE_BITS_MAX`] or `cfg.run_max` is
/// zero.
pub fn compile_with(c: &Circuit, start: usize, cfg: &Config) -> Result<Program, QuantumError> {
    c.validate()?;
    if c.n_qubits > MAX_QUBITS {
        return Err(QuantumError::TooManyQubits(c.n_qubits));
    }
    // `validate` holds the op count to `u32`.
    let gates = c.ops.iter().enumerate().map(|(k, g)| PlannedGate::new(lanes::gate_words(g, c.gateset), k as u32)).collect();
    Ok(schedule(c.n_qubits, start, gates, cfg))
}

/// Plans `gates`, in op-log order, on a state of `n_qubits` from the basis
/// state `start % 2^n`: with `cfg.skip_identity`, every gate whose words are
/// exactly the identity's (and whose target and controls fit the state) is
/// left out; then each maximal run of consecutive gates that act inside a
/// tile of [`Config::tile_size`] amplitudes, split every [`Config::run_cap`]
/// gates, is a [`Step::Tile`], and every other gate a [`Step::Gate`].
///
/// Panics if `n_qubits` is above [`MAX_QUBITS`], `cfg.tile_bits` above
/// [`TILE_BITS_MAX`], or `cfg.run_max` zero.
pub fn schedule(n_qubits: u8, start: usize, gates: Vec<PlannedGate>, cfg: &Config) -> Program {
    assert!(n_qubits <= MAX_QUBITS, "a program holds at most {MAX_QUBITS} qubits, not {n_qubits}");
    assert!(cfg.tile_bits <= TILE_BITS_MAX, "tiles hold at most 2^{TILE_BITS_MAX} amplitudes, not 2^{}", cfg.tile_bits);
    assert!(cfg.run_max >= 1, "a tile run holds at least one gate");
    // An identity whose target or controls a kernel refuses stays, to be refused.
    let len = 1usize << n_qubits;
    let gates: Vec<PlannedGate> = if cfg.skip_identity { gates.into_iter().filter(|g| !(g.is_identity() && tile_shape(&g.words, len))).collect() } else { gates };
    let (size, cap) = (cfg.tile_size(n_qubits), cfg.run_cap(n_qubits) as usize);
    let mut steps = Vec::new();
    let mut k = 0;
    while k < gates.len() {
        if size != 0 && gates[k].fits_tile(size) {
            let first = k;
            while k < gates.len() && k - first < cap && gates[k].fits_tile(size) {
                k += 1;
            }
            // `validate` holds the op count, and so `first` and `k`, to `u32`.
            steps.push(Step::Tile { first: first as u32, count: (k - first) as u32, size });
        } else {
            steps.push(Step::Gate { gate: k as u32 });
            k += 1;
        }
    }
    Program { n_qubits, start: start % len, gates, steps }
}

/// A run in lanes.
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct LaneRun {
    pub amps: Vec<Lane>,
    /// The lowest op-log index whose exact result left 32 bits; `None` if
    /// none did. When it is set, `amps` holds only the low 32 bits of the
    /// reference's values from that op on.
    pub out_of_range: Option<u32>,
}

/// The program on the CPU: `amps[start] = (2^30, 0)`, then each step as
/// [`execute_from`] runs it. Panics only if the program is malformed, which
/// [`compile`], [`compile_with`] and [`schedule`] never make.
pub fn execute(p: &Program) -> LaneRun {
    let mut amps = vec![[0i32; 2]; 1usize << p.n_qubits];
    amps[p.start] = [ONE as i32, 0];
    execute_from(p, &amps).expect("a planned program's words fit its state")
}

/// The program's steps on the CPU from `state` (which replaces the basis
/// state `p.start`): each [`Step::Gate`] with [`lanes::apply_words`], each
/// [`Step::Tile`] tile by tile with a tile kernel's own index arithmetic.
/// `None` if `state` is not `2^n` lanes or the program is malformed: a gate
/// that is not [`PlannedGate::is_well_formed`], a step naming gates it does
/// not hold, a tile size that is not a power of two of at most
/// `2^TILE_BITS_MAX` and the state's length, an empty run, or words a kernel
/// refuses (those [`lanes::apply_words`] refuses, and in a run a gate that
/// neither is diagonal nor has its target's bit below the tile's size). A
/// member refuses the same programs.
pub fn execute_from(p: &Program, state: &[Lane]) -> Option<LaneRun> {
    if p.n_qubits > MAX_QUBITS || state.len() != 1usize << p.n_qubits || !p.gates.iter().all(PlannedGate::is_well_formed) {
        return None;
    }
    let mut amps = state.to_vec();
    let mut out_of_range: Option<u32> = None;
    let mut flag = |index: u32| out_of_range = Some(out_of_range.map_or(index, |k| k.min(index)));
    for step in &p.steps {
        match *step {
            Step::Gate { gate } => {
                let g = p.gates.get(gate as usize)?;
                if !lanes::apply_words(&mut amps, &g.words)? {
                    flag(g.index);
                }
            }
            Step::Tile { first, count, size } => {
                let run = p.gates.get(first as usize..(first as usize).checked_add(count as usize)?)?;
                let size = size as usize;
                if run.is_empty() || !size.is_power_of_two() || size > 1 << TILE_BITS_MAX || size > amps.len() {
                    return None;
                }
                for g in run {
                    if !tile_shape(&g.words, amps.len()) || !(g.diag || (g.words.tbit as usize) < size) {
                        return None;
                    }
                }
                for base in (0..amps.len()).step_by(size) {
                    for g in run {
                        if !tile_gate(&mut amps[base..base + size], base, g) {
                            flag(g.index);
                        }
                    }
                }
            }
        }
    }
    Some(LaneRun { amps, out_of_range })
}

/// The predicate both kernels refuse a record by, for a state of `len` lanes:
/// `tbit` a power of two below `len`, `ctrl` clear of `tbit` and of every bit
/// at or above `len`.
fn tile_shape(w: &GateWords, len: usize) -> bool {
    let (tbit, ctrl) = (w.tbit as usize, w.ctrl as usize);
    tbit.is_power_of_two() && tbit < len && ctrl & tbit == 0 && ctrl & !(len - 1) == 0
}

/// The low 32 bits of an exact value, and whether it is a 32-bit one.
fn lane(v: [i64; 2]) -> (Lane, bool) {
    // `as` keeps the low 32 bits.
    ([v[0] as i32, v[1] as i32], i32::try_from(v[0]).is_ok() && i32::try_from(v[1]).is_ok())
}

/// One gate of a run on the tile starting at amplitude `base`, as a tile
/// kernel computes it: a diagonal gate per amplitude, row 0 with a zero
/// partner where the target's bit is clear and row 1 where it is set; any
/// other gate per pair, the pair `q` being `q` with a zero inserted at the
/// target's bit and that index with the bit set. The control test reads the
/// state-wide index. Returns whether every exact result fits in 32 bits.
fn tile_gate(tile: &mut [Lane], base: usize, g: &PlannedGate) -> bool {
    let (t, c) = (g.words.tbit as usize, g.words.ctrl as usize);
    let mut fits = true;
    if g.diag {
        for (e, a) in tile.iter_mut().enumerate() {
            let gi = base + e;
            if gi & c == c {
                let row = if gi & t == 0 { lanes::exact_row(&g.words, 0, *a, [0, 0]) } else { lanes::exact_row(&g.words, 1, [0, 0], *a) };
                let (v, ok) = lane(row);
                *a = v;
                fits &= ok;
            }
        }
    } else {
        for q in 0..tile.len() / 2 {
            let low = q & (t - 1);
            let li = ((q ^ low) << 1) | low;
            let lj = li | t;
            if (base + li) & c == c {
                let (a0, a1) = (tile[li], tile[lj]);
                let (v0, ok0) = lane(lanes::exact_row(&g.words, 0, a0, a1));
                let (v1, ok1) = lane(lanes::exact_row(&g.words, 1, a0, a1));
                (tile[li], tile[lj]) = (v0, v1);
                fits &= ok0 && ok1;
            }
        }
    }
    fits
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::quantum::{BaseGate, Gate, Gateset};
    use crate::quantum_fabric::tests::configurations;
    use crate::quantum_fabric::{Certification, Corpus, FABRIC_DEFAULT, Parity, apply_kernel, wgsl};

    fn widened(lanes: &[Lane]) -> Vec<[i64; 2]> {
        lanes.iter().map(|a| [i64::from(a[0]), i64::from(a[1])]).collect()
    }

    fn native(c: &Circuit, start: usize) -> Vec<[i64; 2]> {
        c.simulate_from(start).unwrap().amps.iter().map(|a| [a.re, a.im]).collect()
    }

    fn hex(b: &[u8]) -> String {
        b.iter().map(|x| format!("{x:02x}")).collect()
    }

    fn next(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)
    }

    #[test]
    fn execute_is_the_native_simulator() {
        for gs in Gateset::ALL {
            let random = Corpus::Random { seed: 7, cases: 300 }.circuits(gs);
            let fabric = FABRIC_DEFAULT.circuits(gs);
            for (name, circuits) in [("random-v1", random), ("fabric-v1", fabric)] {
                for c in &circuits {
                    for start in [0, 5] {
                        let p = compile(c, start).unwrap();
                        assert_eq!(p.start, start % (1 << c.n_qubits));
                        assert_eq!((p.gates.len(), p.steps.len()), (c.ops.len(), c.ops.len()));
                        let run = execute(&p);
                        assert_eq!(run.out_of_range, None, "{name} {gs:?}");
                        assert_eq!(widened(&run.amps), native(c, start), "{name} {gs:?} start {start}");
                    }
                }
            }
        }
        // The program is the op-log: one step per op, in order, each gate
        // naming its position; `diag` exactly when the off-diagonal is zero.
        let mut c = Circuit::with_gateset(3, Gateset::V1);
        c.h(0).t(1).cx(0, 2).p(19, 2).y(1).z(0);
        let p = compile(&c, 11).unwrap();
        assert_eq!(p.start, 3);
        assert_eq!(p.steps, (0..6).map(|gate| Step::Gate { gate }).collect::<Vec<_>>());
        assert_eq!(p.gates.iter().map(|g| g.index).collect::<Vec<_>>(), [0, 1, 2, 3, 4, 5]);
        assert_eq!(p.gates.iter().map(|g| g.diag).collect::<Vec<_>>(), [false, true, false, true, false, true]);
        assert_eq!(compile_with(&c, 11, &Config::PER_GATE).unwrap(), p);
        // No qubits: no gates, and the one amplitude.
        let empty = compile(&Circuit::with_gateset(0, Gateset::V2), 9).unwrap();
        assert_eq!((empty.start, empty.gates.len()), (0, 0));
        assert_eq!(execute(&empty), LaneRun { amps: vec![[ONE as i32, 0]], out_of_range: None });
        assert_eq!(compile_with(&Circuit::with_gateset(0, Gateset::V2), 9, &Config::TILED).unwrap(), empty);
        // It refuses what the native simulator refuses, with the same error,
        // in every configuration.
        for cfg in configurations() {
            for c in crate::quantum_fabric::tests::refused_circuits() {
                assert_eq!(compile_with(&c, 0, &cfg).unwrap_err(), c.simulate_from(0).unwrap_err(), "{c:?}");
            }
        }
        let mut wide = Circuit::new(MAX_QUBITS + 1);
        wide.ops.push(Gate { base: BaseGate::X, controls: vec![], target: 0, param: 0 });
        assert_eq!(compile(&wide, 0).unwrap_err(), QuantumError::TooManyQubits(MAX_QUBITS + 1));
        assert_eq!(compile_with(&wide, 0, &Config::TILED).unwrap_err(), QuantumError::TooManyQubits(MAX_QUBITS + 1));
        assert_eq!(compile(&Circuit::new(MAX_QUBITS + 1), 0).unwrap_err(), QuantumError::TooManyQubits(MAX_QUBITS + 1));
        // Validation comes first, as in the native simulator.
        let mut both = Circuit::new(MAX_QUBITS + 1);
        both.ops.push(Gate { base: BaseGate::P, controls: vec![], target: 0, param: 0 });
        assert_eq!(compile(&both, 0).unwrap_err(), both.simulate_from(0).unwrap_err());
        assert!(matches!(compile(&both, 0), Err(QuantumError::Invalid(_))));
    }

    /// The lanes32 digests, random-v1 and fabric-v1, per gateset.
    const PINS: [(Gateset, &str, &str); 2] = [
        (Gateset::V1, "d18454c639acfd9c7b66450f4ded1aac4ab46cc10a7bedbcf4958a96e204fb36", "41c726a0ac80f7defef92e15ecf893bdb7a3076b5e7f670fa0b7142e5c033df7"),
        (Gateset::V2, "3dbde276c1095bf827faec71d5980e050f8a55a5ba7d97cd3bbdaec3f43a13b2", "9a84b55af68f1a8b07390b234cf562dd8d85cef3185a697ab85690c8790a7809"),
    ];

    /// Every configuration's program equals the native simulator over
    /// random-v1 (300 circuits) and fabric-v1, both gatesets, and certifies
    /// with the lanes32 pins: a certificate of one configuration
    /// (implementation `lanes32-plan`) has the pinned corpus digest.
    #[test]
    fn every_configuration_is_the_native_simulator() {
        let configs = configurations();
        assert_eq!(configs.len(), 10);
        for (gs, random_pin, fabric_pin) in PINS {
            // fabric-v1 is random-v1's first 200 circuits, then the wide ones.
            let fabric = FABRIC_DEFAULT.circuits(gs);
            let more = Corpus::Random { seed: 7, cases: 300 }.circuits(gs);
            assert_eq!(fabric[..200], more[..200]);
            let circuits: Vec<&Circuit> = fabric.iter().chain(&more[200..]).collect();
            let reference: Vec<_> = circuits.iter().map(|c| lanes::reference_output(c).unwrap()).collect();
            for cfg in &configs {
                let (mut random_v1, mut fabric_v1) = (Certification::new(), Certification::new());
                for (k, (c, want)) in circuits.iter().zip(&reference).enumerate() {
                    let run = execute(&compile_with(c, 0, cfg).unwrap());
                    assert_eq!(run.out_of_range, None, "{cfg:?} {gs:?} case {k}");
                    let out = lanes::output(&run.amps);
                    assert_eq!(&out, want, "{cfg:?} {gs:?} case {k}");
                    let input = lanes::circuit_bytes(c);
                    if k < 200 {
                        random_v1.add(&input, want, &out);
                    }
                    if k < fabric.len() {
                        fabric_v1.add(&input, want, &out);
                    }
                }
                // A start other than zero, on random-v1's first circuits.
                for c in &circuits[..40] {
                    assert_eq!(widened(&execute(&compile_with(c, 5, cfg).unwrap()).amps), native(c, 5), "{cfg:?} {gs:?} start 5");
                }
                let random_v1 = random_v1.finish(apply_kernel(gs), "lanes32-plan", "local/test", &Corpus::Random { seed: 7, cases: 200 }.generator(gs));
                let fabric_v1 = fabric_v1.finish(apply_kernel(gs), "lanes32-plan", "local/test", &FABRIC_DEFAULT.generator(gs));
                for (cert, cases, pin) in [(random_v1, 200, random_pin), (fabric_v1, 224, fabric_pin)] {
                    assert_eq!((cert.parity, cert.cases), (Parity::BitExact, cases), "{cfg:?} {gs:?}");
                    assert_eq!(hex(&cert.corpus), pin, "{cfg:?} {gs:?}");
                }
            }
        }
    }

    /// A random circuit of `ops` gates on `n` qubits: every base gate, up to
    /// three distinct controls, `P(1..=32)`.
    fn random_circuit(n: u8, ops: usize, seed: u64, gs: Gateset) -> Circuit {
        let bases = [BaseGate::I, BaseGate::X, BaseGate::Y, BaseGate::Z, BaseGate::H, BaseGate::S, BaseGate::Sdg, BaseGate::T, BaseGate::Tdg, BaseGate::P];
        let mut s = seed;
        let mut c = Circuit::with_gateset(n, gs);
        for _ in 0..ops {
            let base = bases[(next(&mut s) % bases.len() as u64) as usize];
            let target = (next(&mut s) % u64::from(n)) as u8;
            let mut controls = Vec::new();
            for _ in 0..next(&mut s) % 4 {
                let q = (next(&mut s) % u64::from(n)) as u8;
                if q != target && !controls.contains(&q) {
                    controls.push(q);
                }
            }
            let param = if base == BaseGate::P { 1 + (next(&mut s) % 32) as u16 } else { 0 };
            c.ops.push(Gate { base, controls, target, param });
        }
        c
    }

    /// Long circuits: the tiled program equals the per-gate one and the
    /// native simulator, with runs long enough to be split at their cap.
    #[test]
    fn long_circuits_tile_as_they_run_gate_by_gate() {
        let short_runs = Config { tile_bits: 3, run_max: 7, skip_identity: true };
        for gs in Gateset::ALL {
            for (n, ops, seed) in [(4u8, 3000usize, 0x6c6f_6e67u64), (11, 600, 0x7469_6c65), (12, 400, 0x7275_6e73)] {
                let c = random_circuit(n, ops, seed + n as u64, gs);
                assert!(c.ops.len() as u64 <= lanes::horizon(n, gs));
                let per_gate = execute(&compile(&c, 0).unwrap());
                assert_eq!((per_gate.out_of_range, widened(&per_gate.amps)), (None, native(&c, 0)), "{gs:?} n = {n}");
                for cfg in [Config::TILED, short_runs] {
                    let p = compile_with(&c, 0, &cfg).unwrap();
                    let longest = p.steps.iter().filter_map(|s| match *s { Step::Tile { count, .. } => Some(count), _ => None }).max().unwrap();
                    assert!(longest > 1, "{cfg:?} {gs:?} n = {n}");
                    assert!(longest <= cfg.run_cap(n));
                    assert_eq!(execute(&p), per_gate, "{cfg:?} {gs:?} n = {n}");
                }
            }
            // A 4-qubit circuit is one tile, and every gate fits it: the
            // runs are cut only at the cap.
            let c = random_circuit(4, 3000, 3, gs);
            let p = compile_with(&c, 0, &Config::TILED).unwrap();
            let counts: Vec<u32> = p.steps.iter().map(|s| match *s { Step::Tile { count, size, .. } => { assert_eq!(size, 16); count } other => panic!("{other:?}") }).collect();
            assert_eq!(counts.iter().sum::<u32>() as usize, p.gates.len());
            assert!(counts[..counts.len() - 1].iter().all(|&k| k == 1024), "{counts:?}");
        }
    }

    /// Any words, whatever their range: the tiled program equals the gates
    /// one by one ([`lanes::run_words`]), lanes and report, with diagonal
    /// gates on every target, controls above and below the tile's bits, and
    /// runs cut by gates that do not fit a tile.
    #[test]
    fn tiles_are_exact_on_any_words() {
        let mut s = 0x7469_6c65_7321u64;
        let (mut flagged, mut clean, mut tiled) = (0, 0, 0);
        for case in 0..1500 {
            let n = 1 + (next(&mut s) % 12) as u8;
            let full = case % 2 == 0;
            let value = |s: &mut u64| -> i32 { if full { next(s) as i32 } else { (next(s) % (1 << 23)) as i32 - (1 << 22) } };
            let state: Vec<Lane> = (0..1usize << n).map(|_| [value(&mut s), value(&mut s)]).collect();
            let gates: Vec<GateWords> = (0..1 + next(&mut s) % 24)
                .map(|_| {
                    let t = (next(&mut s) % u64::from(n)) as u32;
                    let ctrl = (0..u32::from(n)).filter(|&b| b != t && next(&mut s) % 4 == 0).fold(0u32, |acc, b| acc | 1 << b);
                    let mut m = [0u32; 8];
                    let entry = |s: &mut u64| -> u32 { if full { next(s) as u32 } else { ((next(s) % (1 << 30)) as i64 - (1 << 29)) as i32 as u32 } };
                    m.iter_mut().for_each(|x| *x = entry(&mut s));
                    match next(&mut s) % 4 {
                        0 | 1 => m[2..6].fill(0),
                        2 => m = IDENTITY,
                        _ => {}
                    }
                    GateWords { m, tbit: 1 << t, ctrl }
                })
                .collect();
            let want = lanes::run_words(&state, &gates).expect("valid shapes");
            let planned: Vec<PlannedGate> = gates.iter().enumerate().map(|(k, w)| PlannedGate::new(*w, k as u32)).collect();
            for cfg in [Config::TILED, Config { tile_bits: 2, run_max: 3, skip_identity: false }] {
                let p = schedule(n, 0, planned.clone(), &cfg);
                if p.steps.iter().any(|s| matches!(s, Step::Tile { count, .. } if *count > 1)) {
                    tiled += 1;
                }
                assert_eq!(execute_from(&p, &state), Some(want.clone()), "case {case}, n = {n}, {cfg:?}");
            }
            if want.out_of_range.is_some() {
                flagged += 1;
            } else {
                clean += 1;
            }
        }
        assert!(flagged >= 100 && clean >= 100 && tiled >= 1000, "{flagged} flagged, {clean} clean, {tiled} with a run of two gates or more");
    }

    #[test]
    fn tiles_cover_every_path() {
        for gs in Gateset::ALL {
            let mut count = std::collections::BTreeMap::<&str, usize>::new();
            let mut add = |what: &'static str, k: usize| *count.entry(what).or_default() += k;
            for c in FABRIC_DEFAULT.circuits(gs) {
                let p = compile_with(&c, 0, &Config::TILED).unwrap();
                add("identity skipped", c.ops.len() - p.gates.len());
                for step in &p.steps {
                    match *step {
                        Step::Gate { .. } => add("per-gate step", 1),
                        Step::Tile { first, count, size } => {
                            add("tile run of two gates or more", usize::from(count > 1));
                            add("tile of 1024 amplitudes", usize::from(size == 1024));
                            add("tile below 1024 amplitudes", usize::from(size < 1024));
                            let tiles = (1u32 << c.n_qubits) / size;
                            add("tile grid of more than one row", usize::from(wgsl::tile_grid(tiles, wgsl::TILE_X_GROUPS_MAX).y > 1));
                            for g in &p.gates[first as usize..(first + count) as usize] {
                                add("diagonal gate on a target outside the tile", usize::from(g.diag && g.words.tbit >= size));
                                add("non-diagonal gate in a run", usize::from(!g.diag));
                                add("control above the tile's bits", usize::from(g.words.ctrl >= size));
                                add("control inside the tile", usize::from(g.words.ctrl & (size - 1) != 0));
                            }
                        }
                    }
                }
            }
            assert_eq!(count.len(), 10, "{gs:?}: {count:?}");
            for (what, k) in &count {
                assert!(*k > 0, "{gs:?}: no {what}");
            }
        }
        // The cap: 1024 gates per run up to 20 qubits, then half as many for
        // every qubit more, 16 at 26.
        for (n, cap) in [(1, 1024), (20, 1024), (21, 512), (24, 64), (26, 16)] {
            assert_eq!(Config::TILED.run_cap(n), cap, "n = {n}");
        }
        assert_eq!((Config::TILED.tile_size(0), Config::TILED.tile_size(3), Config::TILED.tile_size(10), Config::TILED.tile_size(26)), (0, 8, 1024, 1024));
        assert_eq!((Config::PER_GATE.tile_size(12), Config::PER_GATE.run_cap(12)), (0, 1));
    }

    #[test]
    fn the_v1_tail_phases_are_not_skipped() {
        // P(19..=32) of circuit/1 is (2^30 − 1, 0): not the identity.
        let p19 = lanes::gate_words(&Gate { base: BaseGate::P, controls: vec![], target: 0, param: 19 }, Gateset::V1);
        assert_eq!(p19.m, [ONE as u32, 0, 0, 0, 0, 0, (ONE - 1) as u32, 0]);
        assert!(!PlannedGate::new(p19, 0).is_identity() && PlannedGate::new(p19, 0).diag);
        let mut c = Circuit::with_gateset(1, Gateset::V1);
        c.h(0).p(19, 0);
        let p = compile_with(&c, 0, &Config::TILED).unwrap();
        assert_eq!(p.gates.len(), 2);
        assert_eq!(p.steps, [Step::Tile { first: 0, count: 2, size: 2 }]);
        let run = execute(&p);
        assert_eq!(widened(&run.amps), native(&c, 0));
        let h = Gateset::V1.inv_sqrt2() as i32;
        assert_eq!(run.amps, [[h, 0], [h - 1, 0]]);
        // It takes 2^29 + 1 to 2^29.
        let tail = schedule(1, 0, vec![PlannedGate::new(p19, 0)], &Config::TILED);
        assert_eq!(tail.steps, [Step::Tile { first: 0, count: 1, size: 2 }]);
        let run = execute_from(&tail, &[[0, 0], [(1 << 29) + 1, 0]]).unwrap();
        assert_eq!(run.amps, [[0, 0], [1 << 29, 0]]);
        // The identity itself is skipped, and an identity left in changes
        // nothing.
        let mut c = Circuit::with_gateset(2, Gateset::V1);
        c.h(0);
        c.ops.push(Gate { base: BaseGate::I, controls: vec![1], target: 0, param: 0 });
        c.h(1);
        let skipped = compile_with(&c, 0, &Config::TILED).unwrap();
        assert_eq!(skipped.gates.iter().map(|g| g.index).collect::<Vec<_>>(), [0, 2]);
        let kept = compile_with(&c, 0, &Config { skip_identity: false, ..Config::TILED }).unwrap();
        assert_eq!(kept.gates.len(), 3);
        assert_eq!(execute(&skipped), execute(&kept));
    }

    /// Programs a tile kernel refuses are refused here too.
    #[test]
    fn malformed_programs_are_refused() {
        let x2 = PlannedGate::new(lanes::gate_words(&Gate { base: BaseGate::X, controls: vec![], target: 2, param: 0 }, Gateset::V2), 0);
        let z0 = PlannedGate::new(lanes::gate_words(&Gate { base: BaseGate::Z, controls: vec![], target: 0, param: 0 }, Gateset::V2), 0);
        let state = vec![[0, 0]; 8];
        let program = |gates: Vec<PlannedGate>, steps: Vec<Step>| Program { n_qubits: 3, start: 0, gates, steps };
        // A non-diagonal gate whose target is the tile's size: its pairs leave the tile.
        assert_eq!(execute_from(&program(vec![x2], vec![Step::Tile { first: 0, count: 1, size: 4 }]), &state), None);
        assert!(execute_from(&program(vec![x2], vec![Step::Tile { first: 0, count: 1, size: 8 }]), &state).is_some());
        // Sizes and ranges a run cannot have.
        for (first, count, size) in [(0, 0, 4), (0, 2, 4), (1, 1, 4), (0, 1, 3), (0, 1, 16), (0, 1, 0)] {
            assert_eq!(execute_from(&program(vec![z0], vec![Step::Tile { first, count, size }]), &state), None, "{first} {count} {size}");
        }
        assert_eq!(execute_from(&program(vec![z0], vec![Step::Gate { gate: 1 }]), &state), None);
        // Words a kernel refuses, inside a run.
        for w in [GateWords { tbit: 3, ..z0.words }, GateWords { tbit: 8, ..z0.words }, GateWords { ctrl: 1, ..z0.words }, GateWords { ctrl: 8, ..z0.words }] {
            assert_eq!(execute_from(&program(vec![PlannedGate { words: w, ..z0 }], vec![Step::Tile { first: 0, count: 1, size: 8 }]), &state), None, "{w:?}");
        }
        // An identity a kernel would refuse is not skipped, so it is refused.
        let refused_identity = PlannedGate::new(GateWords { m: IDENTITY, tbit: 3, ctrl: 0 }, 0);
        let p = schedule(3, 0, vec![refused_identity], &Config::TILED);
        assert_eq!(p.gates.len(), 1);
        assert_eq!(execute_from(&p, &state), None);
        // A gate whose `diag` disagrees with its words: X planned as diagonal
        // would act per amplitude, and wrongly; a diagonal gate planned as not
        // diagonal is refused alike. In a run and alone.
        let x0 = PlannedGate::new(lanes::gate_words(&Gate { base: BaseGate::X, controls: vec![], target: 0, param: 0 }, Gateset::V2), 0);
        for g in [PlannedGate { diag: true, ..x0 }, PlannedGate { diag: false, ..z0 }] {
            assert!(!g.is_well_formed());
            for step in [Step::Tile { first: 0, count: 1, size: 8 }, Step::Gate { gate: 0 }] {
                assert_eq!(execute_from(&program(vec![g], vec![step]), &state), None, "{g:?} {step:?}");
            }
        }
        assert!(x0.is_well_formed() && z0.is_well_formed());
        // The index the tile kernel keeps for "no op".
        let last = PlannedGate { index: u32::MAX, ..z0 };
        assert!(!last.is_well_formed());
        for step in [Step::Tile { first: 0, count: 1, size: 8 }, Step::Gate { gate: 0 }] {
            assert_eq!(execute_from(&program(vec![last], vec![step]), &state), None);
        }
        assert!(execute_from(&program(vec![PlannedGate { index: u32::MAX - 1, ..z0 }], vec![Step::Tile { first: 0, count: 1, size: 8 }]), &state).is_some());
        // A state of the wrong length.
        assert_eq!(execute_from(&program(vec![z0], vec![Step::Gate { gate: 0 }]), &state[..4]), None);
    }
}