wai_quantum/

quantum.rs

//! Deterministic quantum-circuit **simulation** as transportable, receiptable
//! content — `wai.quantum.circuit2` and the frozen `wai.quantum.circuit`
//! (extensions/quantum-sim).
//!
//! The framing first: WAI's core property is *determinism* — byte-identical
//! reconstruction on every machine, "conformance is a hash". Real quantum
//! HARDWARE is the natural enemy of that (NISQ noise + stochastic measurement are
//! not reproducible across devices). So quantum does NOT enter WAI as a hardware
//! dependency and gives NO quantum speedup here. It enters exactly as worlds,
//! films and scores do (lever 3, *instructions-at-the-sink*): the wire carries a
//! compact **circuit** (a gate op-log), the sink CLASSICALLY simulates it, and the
//! reconstructed statevector hashes identically everywhere.
//!
//! The reconstruction runs in a fixed-point complex floor (`wai.det.amp64`, the
//! arithmetic a circuit/1 contract calls `wai.det.fixed64`): every amplitude is
//! an [`Amp`] — real and imaginary parts as `i64` at scale `2^FRAC` — and every
//! gate is an integer fixed-point matrix–vector update with a pinned rounding
//! rule and a *stated error bound* versus the ideal unitary (the same contract
//! as the WAI integer transform floor). No float anywhere on the decode, so:
//!
//! - **statevector-equivalence** — the whole amplitude vector is byte-identical on
//!   every machine (a portable statevector hash), and
//! - **shot-histogram-equivalence** — measurement is sampled by the pinned
//!   `splitmix64` PRNG (the same fixed-point sampler the WAI deterministic generator uses), so
//!   a shot histogram at a pinned `(seed, shots)` is byte-identical too.
//!
//! The transportable artifact is the *circuit* (kilobytes, linear in gate count);
//! the reconstruction is the *statevector* (`2^n` amplitudes — exponential). That
//! gap is the point: a tiny circuit reconstructs to a gigabyte statevector, and
//! because the reconstruction cost is exponential the joules-accounted receipt
//! (§ extension) is a genuinely novel artifact — an energy-metered, byte-exact,
//! signable record of a quantum computation.
//!
//! Gate set: a *controlled-1-qubit-gate* IR over Clifford+T + dyadic phase
//! (`I,X,Y,Z,H,S,S†,T,T†,P(k)` with any number of controls). Clifford+T is dense
//! in `SU(2^n)`, so this is universal. The dyadic phase `P(k) = e^{2πi/2^k}` is
//! read from the circuit's [`Gateset`] table:
//!
//! - [`Gateset::V2`] (`wai.quantum.circuit2`, magic `WQC2`) holds the integers
//!   nearest to each phase, so its controlled-phase ladder gives the **QFT**
//!   ("quantum FFT"), transported as a circuit and reconstructed byte-exact.
//!   Circuits built in code are V2 unless built otherwise.
//! - [`Gateset::V1`] (`wai.quantum.circuit`, magic `WQC1`) is frozen so that
//!   every circuit/1 file, hash and receipt reproduces. Its `P(k)` deviates from
//!   the ideal phase for `k ≥ 8` and applies none for `k ≥ 19`.
//!
//! A circuit decoded from a container takes the container's gateset, and its
//! [`circuit_hash`] names the capability. Pure i64, no float, no `ort`;
//! compiles to wasm32.
//!
//! Out of scope (and honestly so): variational-quantum-circuit / "quantum AI"
//! priors whose measurement is stochastic — those are at best a float/behavioral
//! capability that can never back `intent=replicate`, so they do not belong in
//! this exact-reconstruction extension. [`StateVector::fidelity_fx`] is provided as the
//! quantum-information-theoretic conformance metric such a future path would use.

use crate::ijson::{self, Json, RANGE_LIMIT};

/// Fractional bits of the fixed-point amplitude floor. Amplitudes live in
/// `[-1, 1]`, so `2^FRAC` (≈ `1.07e9` at 30) is the unit and the i128 products in
/// [`Amp::mul`] stay far inside range.
pub const FRAC: u32 = 30;
/// The fixed-point unit `1.0`.
pub const ONE: i64 = 1 << FRAC;
const ROUND: i128 = 1 << (FRAC - 1);
/// Largest `k` for which `P(k) = e^{2πi/2^k}` is tabulated: [`PHASE_V1`] and
/// [`PHASE_V2`] each hold `P(1..=DYADIC_MAX)`.
pub const DYADIC_MAX: usize = 32;

/// Domain separators for the conformance hashes (the WAI discipline).
const DOMAIN_STATEVECTOR: &[u8] = b"wai:quantum-statevector\x01";
const DOMAIN_CIRCUIT: &[u8] = b"wai:quantum-circuit\x01";
const DOMAIN_HISTOGRAM: &[u8] = b"wai:quantum-histogram\x01";

// ---------------------------------------------------------------------------
// Fixed-point complex amplitude
// ---------------------------------------------------------------------------

/// A complex amplitude in the fixed-point floor: `re, im` are `i64` at scale
/// `2^FRAC`. Byte-exact by construction — `Eq` is total (no float), so two
/// statevectors compare and hash identically iff every amplitude matches.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub struct Amp {
    pub re: i64,
    pub im: i64,
}

/// Fixed-point multiply with a pinned round-half-up rule. `i128` intermediate so
/// the `i64·i64` product never overflows; the `+ROUND` then arithmetic shift is
/// identical on every machine (the determinism the whole floor rests on).
#[inline]
pub fn fxmul(a: i64, b: i64) -> i64 {
    ((a as i128 * b as i128 + ROUND) >> FRAC) as i64
}

/// Fixed-point square root: `v` is at `2^FRAC`, result at `2^FRAC`
/// (`sqrt(v/ONE)·ONE = isqrt(v·ONE)`). Integer `isqrt`, no float; negatives
/// (never produced by the half-angle recurrence) clamp to 0.
#[inline]
pub fn sqrt_fx(v: i64) -> i64 {
    if v <= 0 {
        return 0;
    }
    (((v as u128) << FRAC).isqrt()) as i64
}

impl Amp {
    /// `0 + 0i`.
    pub const ZERO: Amp = Amp { re: 0, im: 0 };
    /// `1 + 0i` (the fixed-point unit).
    pub const ONE: Amp = Amp { re: ONE, im: 0 };

    #[inline]
    pub fn add(self, o: Amp) -> Amp {
        Amp { re: self.re + o.re, im: self.im + o.im }
    }

    /// Complex fixed-point multiply: `(a+bi)(c+di) = (ac−bd) + (ad+bc)i`.
    #[inline]
    pub fn mul(self, o: Amp) -> Amp {
        Amp {
            re: fxmul(self.re, o.re) - fxmul(self.im, o.im),
            im: fxmul(self.re, o.im) + fxmul(self.im, o.re),
        }
    }

    /// Complex conjugate.
    #[inline]
    pub fn conj(self) -> Amp {
        Amp { re: self.re, im: -self.im }
    }

    /// `|amp|^2` at scale `2^(2·FRAC)`, in `i128` (the un-shifted probability
    /// weight — summed and sampled without a divide).
    #[inline]
    pub fn norm2(self) -> i128 {
        self.re as i128 * self.re as i128 + self.im as i128 * self.im as i128
    }
}

// ---------------------------------------------------------------------------
// Dyadic phase tables — P(k) = e^{2πi / 2^k} at scale 2^FRAC, one per gateset.
// Literal constants: no float and no recurrence runs on the decode path.
// ---------------------------------------------------------------------------

/// The `wai.quantum.circuit` phase table, frozen bit for bit.
///
/// Index `k` holds `P(k) = e^{2πi/2^k}` as `(re, im)` at scale `2^FRAC`, for
/// `1 ≤ k ≤ DYADIC_MAX`; index 0 is unused and holds [`Amp::ZERO`].
///
/// These are the values of the integer half-angle recurrence every circuit/1
/// artefact was made with: `P(1) = (−2^30, 0)`, `P(2) = (0, 2^30)`, and for
/// `k ≥ 3`, with `c = P(k−1).re`,
/// `re = isqrt(⌊(2^30 + c)/2⌋·2^30)` and `im = isqrt(⌊(2^30 − c)/2⌋·2^30)`,
/// where `isqrt` is the floor square root. The floors compound: from `P(8)` on
/// the imaginary part is too large, by up to 7032 units at `P(18)`, and
/// `P(19..=32)` apply no phase at all. The table is kept exactly as it is so that
/// every circuit/1 file, hash and receipt reproduces. [`PHASE_V2`] is the
/// correctly rounded table. The extension states this table in its §3.1 and
/// Appendix B.1, rendered from this constant.
pub const PHASE_V1: [Amp; 33] = [
    Amp::ZERO,
    Amp { re: -1073741824, im: 0 },         // k = 1
    Amp { re: 0, im: 1073741824 },          // k = 2
    Amp { re: 759250124, im: 759250124 },   // k = 3
    Amp { re: 992008094, im: 410903207 },   // k = 4
    Amp { re: 1053110175, im: 209476638 },  // k = 5
    Amp { re: 1068571463, im: 105245102 },  // k = 6
    Amp { re: 1072448454, im: 52686012 },   // k = 7
    Amp { re: 1073418432, im: 26350953 },   // k = 8
    Amp { re: 1073660972, im: 13176485 },   // k = 9
    Amp { re: 1073721610, im: 6588405 },    // k = 10
    Amp { re: 1073736770, im: 3294284 },    // k = 11
    Amp { re: 1073740560, im: 1647223 },    // k = 12
    Amp { re: 1073741507, im: 823774 },     // k = 13
    Amp { re: 1073741744, im: 411887 },     // k = 14
    Amp { re: 1073741803, im: 207243 },     // k = 15
    Amp { re: 1073741818, im: 103621 },     // k = 16
    Amp { re: 1073741822, im: 56755 },      // k = 17
    Amp { re: 1073741823, im: 32768 },      // k = 18
    Amp { re: 1073741823, im: 0 },          // k = 19
    Amp { re: 1073741823, im: 0 },          // k = 20
    Amp { re: 1073741823, im: 0 },          // k = 21
    Amp { re: 1073741823, im: 0 },          // k = 22
    Amp { re: 1073741823, im: 0 },          // k = 23
    Amp { re: 1073741823, im: 0 },          // k = 24
    Amp { re: 1073741823, im: 0 },          // k = 25
    Amp { re: 1073741823, im: 0 },          // k = 26
    Amp { re: 1073741823, im: 0 },          // k = 27
    Amp { re: 1073741823, im: 0 },          // k = 28
    Amp { re: 1073741823, im: 0 },          // k = 29
    Amp { re: 1073741823, im: 0 },          // k = 30
    Amp { re: 1073741823, im: 0 },          // k = 31
    Amp { re: 1073741823, im: 0 },          // k = 32
];

/// The `wai.quantum.circuit2` phase table: each entry is the integer nearest to
/// `2^30·cos(2π/2^k)` and to `2^30·sin(2π/2^k)`.
///
/// Index `k` holds `P(k)` as `(re, im)` at scale `2^FRAC`, for
/// `1 ≤ k ≤ DYADIC_MAX`; index 0 is unused and holds [`Amp::ZERO`]. `P(1)` and
/// `P(2)` are exact. For `k ≥ 3` every true value is irrational, so none is a
/// half-integer and the nearest integer is unique; the closest any comes to a
/// tie is 0.00456 unit, at the imaginary part of `P(22)`. Every component is
/// within ½ unit of the ideal, which the crate's tests prove with integer
/// interval arithmetic. The literal values are normative: the extension
/// states them in its §3.2 and Appendix B.2, rendered from this constant.
///
/// `re P(3) == im P(3)`, so `T` keeps `re == im`. The resolution of `2^-30` is a
/// floor no rounding removes: `P(26..=32)` carry imaginary parts of 101 down to
/// 2, so the smallest angles are correct to the nearest unit and no further.
pub const PHASE_V2: [Amp; 33] = [
    Amp::ZERO,
    Amp { re: -1073741824, im: 0 },         // k = 1
    Amp { re: 0, im: 1073741824 },          // k = 2
    Amp { re: 759250125, im: 759250125 },   // k = 3
    Amp { re: 992008094, im: 410903207 },   // k = 4
    Amp { re: 1053110176, im: 209476638 },  // k = 5
    Amp { re: 1068571464, im: 105245103 },  // k = 6
    Amp { re: 1072448455, im: 52686014 },   // k = 7
    Amp { re: 1073418433, im: 26350943 },   // k = 8
    Amp { re: 1073660973, im: 13176464 },   // k = 9
    Amp { re: 1073721611, im: 6588356 },    // k = 10
    Amp { re: 1073736771, im: 3294193 },    // k = 11
    Amp { re: 1073740561, im: 1647099 },    // k = 12
    Amp { re: 1073741508, im: 823550 },     // k = 13
    Amp { re: 1073741745, im: 411775 },     // k = 14
    Amp { re: 1073741804, im: 205887 },     // k = 15
    Amp { re: 1073741819, im: 102944 },     // k = 16
    Amp { re: 1073741823, im: 51472 },      // k = 17
    Amp { re: 1073741824, im: 25736 },      // k = 18
    Amp { re: 1073741824, im: 12868 },      // k = 19
    Amp { re: 1073741824, im: 6434 },       // k = 20
    Amp { re: 1073741824, im: 3217 },       // k = 21
    Amp { re: 1073741824, im: 1608 },       // k = 22
    Amp { re: 1073741824, im: 804 },        // k = 23
    Amp { re: 1073741824, im: 402 },        // k = 24
    Amp { re: 1073741824, im: 201 },        // k = 25
    Amp { re: 1073741824, im: 101 },        // k = 26
    Amp { re: 1073741824, im: 50 },         // k = 27
    Amp { re: 1073741824, im: 25 },         // k = 28
    Amp { re: 1073741824, im: 13 },         // k = 29
    Amp { re: 1073741824, im: 6 },          // k = 30
    Amp { re: 1073741824, im: 3 },          // k = 31
    Amp { re: 1073741824, im: 2 },          // k = 32
];

/// Which phase table a circuit is reconstructed with, and so which capability
/// carries it.
///
/// Both gatesets share the arithmetic, the gate matrices and the container
/// layout. They differ only in the phase table, which `H`, `T`, `T†` and `P(k)`
/// read (`H` uses `re P(3)` as `1/√2`).
///
/// - [`Gateset::V1`] is `wai.quantum.circuit` (magic `WQC1`), frozen so that
///   existing content reproduces bit for bit.
/// - [`Gateset::V2`] is `wai.quantum.circuit2` (magic `WQC2`), with the correctly
///   rounded table. It is the default.
#[derive(Clone, Copy, Debug, Default, PartialEq, Eq, Hash)]
#[non_exhaustive]
pub enum Gateset {
    /// `wai.quantum.circuit`, with the frozen [`PHASE_V1`] table.
    V1,
    /// `wai.quantum.circuit2`, with the correctly rounded [`PHASE_V2`] table.
    #[default]
    V2,
}

impl Gateset {
    /// Every gateset, oldest first.
    pub const ALL: [Gateset; 2] = [Gateset::V1, Gateset::V2];

    /// The capability string: `wai.quantum.circuit` or `wai.quantum.circuit2`.
    pub const fn capability(self) -> &'static str {
        match self {
            Gateset::V1 => "wai.quantum.circuit",
            Gateset::V2 => "wai.quantum.circuit2",
        }
    }

    /// The container magic: `WQC1` or `WQC2`.
    pub const fn magic(self) -> [u8; 4] {
        match self {
            Gateset::V1 => *b"WQC1",
            Gateset::V2 => *b"WQC2",
        }
    }

    /// The contract's `ext`: `<capability>/<contract revision>`, so
    /// `wai.quantum.circuit2/1` is revision 1 of the circuit2 capability.
    pub const fn ext(self) -> &'static str {
        match self {
            Gateset::V1 => "wai.quantum.circuit/1",
            Gateset::V2 => "wai.quantum.circuit2/1",
        }
    }

    /// The contract's `numeric` identifier.
    pub const fn numeric(self) -> &'static str {
        match self {
            Gateset::V1 => "wai.det.fixed64",
            Gateset::V2 => "wai.det.amp64",
        }
    }

    /// The contract's `gateset` string, which names the reconstruction's gate
    /// algebra and so differs wherever its constants differ.
    pub const fn name(self) -> &'static str {
        match self {
            Gateset::V1 => "cliffordT+dyadicP",
            Gateset::V2 => "cliffordT+dyadicP.rn",
        }
    }

    /// The phase table: [`PHASE_V1`] or [`PHASE_V2`].
    pub const fn table(self) -> &'static [Amp; 33] {
        match self {
            Gateset::V1 => &PHASE_V1,
            Gateset::V2 => &PHASE_V2,
        }
    }

    /// `P(k) = e^{2πi/2^k}` from this gateset's table. `validate` guarantees
    /// `1 ≤ k ≤ DYADIC_MAX` for every `P` in a circuit.
    #[inline]
    pub fn phase(self, k: usize) -> Amp {
        self.table()[k]
    }

    /// `1/√2` in this gateset: `re P(3)`, the single source of the one irrational
    /// the Clifford+T set needs. 759250124 in V1, 759250125 in V2.
    #[inline]
    pub fn inv_sqrt2(self) -> i64 {
        self.table()[3].re
    }

    /// The canonical contract bytes for an `n_qubits` circuit: no whitespace,
    /// members in a fixed order, rendered directly rather than through a JSON
    /// serializer so that no serializer setting can change them.
    pub fn contract_bytes(self, n_qubits: u8) -> Vec<u8> {
        format!(
            "{{\"ext\":\"{}\",\"n_qubits\":{},\"numeric\":\"{}\",\"frac\":{},\"gateset\":\"{}\"}}",
            self.ext(),
            n_qubits,
            self.numeric(),
            FRAC,
            self.name(),
        )
        .into_bytes()
    }

    /// BLAKE3 of the table: `i64le(re) ‖ i64le(im)` for `k = 1..=DYADIC_MAX`,
    /// 512 bytes, with no domain prefix.
    pub fn table_digest(self) -> [u8; 32] {
        let mut h = blake3::Hasher::new();
        for a in &self.table()[1..] {
            h.update(&a.re.to_le_bytes());
            h.update(&a.im.to_le_bytes());
        }
        *h.finalize().as_bytes()
    }

    /// The gateset a container magic names, or `None` for any other bytes.
    pub fn from_magic(m: &[u8]) -> Option<Gateset> {
        Gateset::ALL.into_iter().find(|g| m == g.magic().as_slice())
    }

    /// The gateset a capability string names, or `None` for any other string.
    pub fn from_capability(s: &str) -> Option<Gateset> {
        Gateset::ALL.into_iter().find(|g| s == g.capability())
    }
}

// ---------------------------------------------------------------------------
// Gates — a controlled-1-qubit-gate IR
// ---------------------------------------------------------------------------

/// The 1-qubit base gates. `P` carries a dyadic index in [`Gate::param`].
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum BaseGate {
    I,
    X,
    Y,
    Z,
    H,
    S,
    Sdg,
    T,
    Tdg,
    /// `P(k) = diag(1, e^{2πi/2^k})`.
    P,
}

impl BaseGate {
    pub(crate) fn opcode(self) -> u8 {
        match self {
            BaseGate::I => 0,
            BaseGate::X => 1,
            BaseGate::Y => 2,
            BaseGate::Z => 3,
            BaseGate::H => 4,
            BaseGate::S => 5,
            BaseGate::Sdg => 6,
            BaseGate::T => 7,
            BaseGate::Tdg => 8,
            BaseGate::P => 9,
        }
    }
    fn from_opcode(b: u8) -> Option<BaseGate> {
        Some(match b {
            0 => BaseGate::I,
            1 => BaseGate::X,
            2 => BaseGate::Y,
            3 => BaseGate::Z,
            4 => BaseGate::H,
            5 => BaseGate::S,
            6 => BaseGate::Sdg,
            7 => BaseGate::T,
            8 => BaseGate::Tdg,
            9 => BaseGate::P,
            _ => return None,
        })
    }

    /// The 2×2 fixed-point matrix `[[m00,m01],[m10,m11]]` for this gate in the
    /// gateset `gs`, whose table `H` (`re P(3)` as `1/√2`), `T`, `T†` and `P`
    /// read. `param` must be in `1..=DYADIC_MAX` for `P`, as `validate` ensures.
    pub(crate) fn matrix(self, param: u16, gs: Gateset) -> [[Amp; 2]; 2] {
        let z = Amp::ZERO;
        let one = Amp::ONE;
        let i = Amp { re: 0, im: ONE };
        let neg_i = Amp { re: 0, im: -ONE };
        let s = Amp { re: gs.inv_sqrt2(), im: 0 };
        let neg_s = Amp { re: -gs.inv_sqrt2(), im: 0 };
        let neg_one = Amp { re: -ONE, im: 0 };
        match self {
            BaseGate::I => [[one, z], [z, one]],
            BaseGate::X => [[z, one], [one, z]],
            BaseGate::Y => [[z, neg_i], [i, z]],
            BaseGate::Z => [[one, z], [z, neg_one]],
            BaseGate::H => [[s, s], [s, neg_s]],
            BaseGate::S => [[one, z], [z, i]],
            BaseGate::Sdg => [[one, z], [z, neg_i]],
            BaseGate::T => [[one, z], [z, gs.phase(3)]],
            BaseGate::Tdg => [[one, z], [z, gs.phase(3).conj()]],
            BaseGate::P => [[one, z], [z, gs.phase(param as usize)]],
        }
    }
}

/// One instruction: a base gate on `target`, applied only when every qubit in
/// `controls` is set. `controls = []` is a plain 1-qubit gate; one control gives
/// CX/CZ/CP; two gives Toffoli-class; etc.
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct Gate {
    pub base: BaseGate,
    pub controls: Vec<u8>,
    pub target: u8,
    pub param: u16,
}

// ---------------------------------------------------------------------------
// Circuit
// ---------------------------------------------------------------------------

/// A quantum circuit: `n_qubits`, an ordered op-log, and the [`Gateset`] it is
/// reconstructed with. This is the transportable object — compact, linear in
/// gate count.
///
/// The gateset is part of the circuit: two circuits with the same ops and
/// different gatesets are different circuits, and compare unequal.
/// [`Circuit::new`] gives [`Gateset::V2`]; a circuit decoded from a container
/// takes the container's; a circuit derived from another keeps the other's
/// ([`Circuit::empty_like`]). The fields stay public, so a circuit/1 build is
/// `let mut c = Circuit::qft(n); c.gateset = Gateset::V1;`.
#[derive(Clone, Debug, PartialEq, Eq)]
#[non_exhaustive]
pub struct Circuit {
    pub n_qubits: u8,
    pub ops: Vec<Gate>,
    pub gateset: Gateset,
}

/// Errors from parsing a `WQC` container or building an out-of-range circuit.
#[derive(Debug, PartialEq, Eq)]
#[non_exhaustive]
pub enum QuantumError {
    /// A container whose structure is broken: a bad magic, a section out of
    /// bounds, repeated or missing, an `n_qubits` that is not an integer in
    /// `0..=255`, a malformed op-log, or a measurement that is not I-JSON. In
    /// `WQC2`, also bytes after the last op or a measurement that is not the
    /// canonical bytes.
    Malformed(String),
    /// A circuit its gateset does not allow: a qubit index ≥ `n_qubits`, a
    /// control equal to the target, more than 255 controls, more than
    /// `u32::MAX` ops, a `P(k)` with `k` outside `1..=DYADIC_MAX`, or, in
    /// [`Gateset::V2`], a repeated control or a parameter on a gate other than
    /// `P`. Also a circuit too large for a container's `u32` section sizes.
    Invalid(String),
    /// A container's contract section does not declare the capability its
    /// magic names: it is not a JSON object carrying that capability's values,
    /// or, in `WQC2`, it is not the canonical bytes.
    Contract(String),
    /// `n_qubits` exceeds what this build will materialize (`2^n` amplitudes).
    TooManyQubits(u8),
}

impl QuantumError {
    /// The error's class, as the conformance corpus names it: `"malformed"`,
    /// `"invalid"`, `"contract"` or `"too-many-qubits"`.
    pub fn kind(&self) -> &'static str {
        match self {
            QuantumError::Malformed(_) => "malformed",
            QuantumError::Invalid(_) => "invalid",
            QuantumError::Contract(_) => "contract",
            QuantumError::TooManyQubits(_) => "too-many-qubits",
        }
    }
}

impl std::fmt::Display for QuantumError {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            QuantumError::Malformed(e) => write!(f, "malformed WQC: {e}"),
            QuantumError::Invalid(e) => write!(f, "invalid circuit: {e}"),
            QuantumError::Contract(e) => write!(f, "contract mismatch: {e}"),
            QuantumError::TooManyQubits(n) => write!(f, "too many qubits: {n}"),
        }
    }
}
impl std::error::Error for QuantumError {}

/// The largest circuit this build will simulate (a state vector is `2^n`
/// amplitudes × 16 bytes; 26 qubits ≈ 1 GiB). The *format* supports more — this
/// is a materialization guard, not a spec limit.
pub const MAX_QUBITS: u8 = 26;

impl Circuit {
    /// An empty circuit on `n_qubits` in [`Gateset::V2`], `wai.quantum.circuit2`.
    pub fn new(n_qubits: u8) -> Self {
        Circuit::with_gateset(n_qubits, Gateset::V2)
    }

    /// An empty circuit on `n_qubits` in the gateset `gateset`.
    pub fn with_gateset(n_qubits: u8, gateset: Gateset) -> Self {
        Circuit { n_qubits, ops: Vec::new(), gateset }
    }

    /// An empty circuit with this one's qubit count and gateset, for building a
    /// circuit derived from it without resetting its gateset.
    pub fn empty_like(&self) -> Self {
        Circuit::with_gateset(self.n_qubits, self.gateset)
    }

    fn push(&mut self, base: BaseGate, controls: Vec<u8>, target: u8, param: u16) -> &mut Self {
        self.ops.push(Gate { base, controls, target, param });
        self
    }

    // Ergonomic builders (used by the QFT constructor and the tests/vectors).
    pub fn x(&mut self, q: u8) -> &mut Self { self.push(BaseGate::X, vec![], q, 0) }
    pub fn y(&mut self, q: u8) -> &mut Self { self.push(BaseGate::Y, vec![], q, 0) }
    pub fn z(&mut self, q: u8) -> &mut Self { self.push(BaseGate::Z, vec![], q, 0) }
    pub fn h(&mut self, q: u8) -> &mut Self { self.push(BaseGate::H, vec![], q, 0) }
    pub fn s(&mut self, q: u8) -> &mut Self { self.push(BaseGate::S, vec![], q, 0) }
    pub fn t(&mut self, q: u8) -> &mut Self { self.push(BaseGate::T, vec![], q, 0) }
    /// `P(k) = diag(1, e^{2πi/2^k})` on qubit `q`.
    pub fn p(&mut self, k: u16, q: u8) -> &mut Self { self.push(BaseGate::P, vec![], q, k) }
    pub fn cx(&mut self, c: u8, t: u8) -> &mut Self { self.push(BaseGate::X, vec![c], t, 0) }
    pub fn cz(&mut self, c: u8, t: u8) -> &mut Self { self.push(BaseGate::Z, vec![c], t, 0) }
    /// Controlled `P(k)` — the QFT's rotation gate.
    pub fn cp(&mut self, k: u16, c: u8, t: u8) -> &mut Self { self.push(BaseGate::P, vec![c], t, k) }
    /// Toffoli (`CCX`).
    pub fn ccx(&mut self, c0: u8, c1: u8, t: u8) -> &mut Self {
        self.push(BaseGate::X, vec![c0, c1], t, 0)
    }
    /// `SWAP(a,b)` as three CX (kept in the op-log as gates, not a new opcode).
    pub fn swap(&mut self, a: u8, b: u8) -> &mut Self {
        self.cx(a, b).cx(b, a).cx(a, b)
    }

    /// The **Quantum Fourier Transform** over all `n` qubits — the "quantum FFT",
    /// transported as a circuit (Hadamard + controlled-phase ladder + bit-reversal
    /// swaps) and reconstructed byte-exact. Built with [`Circuit::new`], so in
    /// [`Gateset::V2`], whose table gives every phase the ladder needs.
    pub fn qft(n: u8) -> Circuit {
        let mut c = Circuit::new(n);
        for j in 0..n {
            c.h(j);
            for l in (j + 1)..n {
                // controlled R_m between control l and target j, m = l−j+1
                let m = (l - j + 1) as u16;
                c.cp(m, l, j);
            }
        }
        for j in 0..(n / 2) {
            c.swap(j, n - 1 - j);
        }
        c
    }

    /// Check the circuit against its gateset's rules.
    ///
    /// In both gatesets every target and control is a qubit of the circuit, no
    /// control equals its target, a gate has at most 255 controls (the op-log
    /// counts them in a `u8`), and `P(k)` has `1 ≤ k ≤ DYADIC_MAX`. The circuit
    /// has at most `u32::MAX` ops, since the op-log counts them in a `u32`.
    /// [`Gateset::V2`] also requires a gate's controls to be pairwise distinct,
    /// so at most 254, and the parameter of every gate other than `P` to be 0.
    pub fn validate(&self) -> Result<(), QuantumError> {
        let n = self.n_qubits;
        let strict = match self.gateset {
            Gateset::V1 => false,
            Gateset::V2 => true,
        };
        if u32::try_from(self.ops.len()).is_err() {
            return Err(QuantumError::Invalid(format!("{} ops, more than u32::MAX", self.ops.len())));
        }
        for g in &self.ops {
            if g.target >= n {
                return Err(QuantumError::Invalid(format!("target {} ≥ n_qubits {}", g.target, n)));
            }
            if g.controls.len() > usize::from(u8::MAX) {
                return Err(QuantumError::Invalid(format!("{} controls, more than 255", g.controls.len())));
            }
            // One bit per qubit index, for the distinct-controls rule.
            let mut seen = [0u64; 4];
            for &c in &g.controls {
                if c >= n {
                    return Err(QuantumError::Invalid(format!("control {c} ≥ n_qubits {n}")));
                }
                if c == g.target {
                    return Err(QuantumError::Invalid(format!("control {c} equals target")));
                }
                if strict {
                    let (word, bit) = (usize::from(c >> 6), 1u64 << (c & 63));
                    if seen[word] & bit != 0 {
                        return Err(QuantumError::Invalid(format!("control {c} repeated")));
                    }
                    seen[word] |= bit;
                }
            }
            if g.base == BaseGate::P {
                let k = g.param as usize;
                if k < 1 || k > DYADIC_MAX {
                    return Err(QuantumError::Invalid(format!("P(k) with k={k} out of 1..={DYADIC_MAX}")));
                }
            } else if strict && g.param != 0 {
                return Err(QuantumError::Invalid(format!("{:?} with param {}; only P takes one", g.base, g.param)));
            }
        }
        Ok(())
    }

    /// Classically simulate the circuit from `|0…0⟩` with its gateset's table,
    /// returning the fixed-point state vector. Pure i64: byte-identical on every
    /// machine.
    pub fn simulate(&self) -> Result<StateVector, QuantumError> {
        self.simulate_from(0)
    }

    /// Simulate starting from computational basis state `|start⟩`.
    pub fn simulate_from(&self, start: usize) -> Result<StateVector, QuantumError> {
        self.validate()?;
        if self.n_qubits > MAX_QUBITS {
            return Err(QuantumError::TooManyQubits(self.n_qubits));
        }
        let dim = 1usize << self.n_qubits;
        let mut amps = vec![Amp::ZERO; dim];
        amps[start % dim] = Amp::ONE;
        for g in &self.ops {
            apply(&mut amps, self.n_qubits, g, self.gateset);
        }
        Ok(StateVector { n_qubits: self.n_qubits, amps })
    }
}

/// Apply one controlled-1-qubit gate in place. Iterates the amplitude pairs that
/// differ only in the target bit, updating those whose control bits are all set.
/// Apply one gate in place.
///
/// The obvious implementation walks every index and multiplies a 2x2 matrix into
/// each pair. Most gates do not need that. `fxmul(ONE, x) == x` and
/// `fxmul(0, x) == 0` hold exactly in this fixed point — `ROUND` is `1 << (FRAC-1)`,
/// so the rounded product of an integer with one is that integer — which makes the
/// shortcuts below **bit-identical to the general path by construction**, not merely
/// close. Byte-exactness is the contract here; a faster simulator that changed a
/// single amplitude would be a broken one.
///
/// - identity does nothing at all;
/// - a **diagonal** gate (`Z S S† T T† P`, and their controlled forms) only phases
///   the `|1⟩` half, so it touches half the amplitudes with one multiply instead of
///   four plus two adds — and a phase of exactly `-1` is a negation, with none;
/// - an **antidiagonal** gate is two multiplies, and `X` (hence `CX`, the most
///   common two-qubit gate there is) is a pure swap with no arithmetic whatsoever;
/// - an all-**real** matrix (`H`, and any real rotation) halves the multiplies,
///   since the imaginary cross terms are multiplications by zero.
///
/// Measured at 20 qubits against the previous implementation, interleaved in one
/// process: `Z` 8.9x, `X` 7.1x, `CZ` 4.6x, `CX` 4.3x, `T` 3.4x, `S` 3.1x, `Y` 1.9x,
/// `H` 1.8x.
///
/// The loop shape is deliberately left alone. Blocking it, and indexing the pairs
/// by bit-insertion, were both tried and both came out slower than the plain scan:
/// the arithmetic is what costs here, not the iteration.
///
/// The shortcuts read only the matrix, so they hold for either gateset's table:
/// circuit2's `P(k ≥ 18)` has `re == ONE` and `im ≠ 0`, which is the diagonal
/// path, and only an exact `(−ONE, 0)` takes the negation.
fn apply(amps: &mut [Amp], n: u8, g: &Gate, gs: Gateset) {
    if matches!(g.base, BaseGate::I) {
        return;
    }
    let m = g.base.matrix(g.param, gs);
    let tbit = 1usize << g.target;
    let ctrl_mask: usize = g.controls.iter().fold(0usize, |acc, &c| acc | (1usize << c));
    let dim = 1usize << n;

    let is_zero = |a: Amp| a.re == 0 && a.im == 0;
    let is_one = |a: Amp| a.re == ONE && a.im == 0;

    // diag(1, phase): the |0> half is untouched, so walk only the |1> half and do
    // one multiply instead of four and two adds. A phase of exactly -1 negates.
    if is_one(m[0][0]) && is_zero(m[0][1]) && is_zero(m[1][0]) {
        let ph = m[1][1];
        let negate = ph.re == -ONE && ph.im == 0;
        let mut i = 0usize;
        while i < dim {
            if i & tbit != 0 && (i & ctrl_mask) == ctrl_mask {
                amps[i] = if negate {
                    Amp { re: -amps[i].re, im: -amps[i].im }
                } else {
                    ph.mul(amps[i])
                };
            }
            i += 1;
        }
        return;
    }

    // antidiagonal: the halves exchange. X (so CX) is an exact swap, no arithmetic.
    if is_zero(m[0][0]) && is_zero(m[1][1]) {
        let (a01, a10) = (m[0][1], m[1][0]);
        let plain_swap = is_one(a01) && is_one(a10);
        let mut i = 0usize;
        while i < dim {
            if i & tbit == 0 && (i & ctrl_mask) == ctrl_mask {
                let j = i | tbit;
                if plain_swap {
                    amps.swap(i, j);
                } else {
                    let (a0, a1) = (amps[i], amps[j]);
                    amps[i] = a01.mul(a1);
                    amps[j] = a10.mul(a0);
                }
            }
            i += 1;
        }
        return;
    }

    // All four entries real — H, and any real rotation. A complex multiply by a
    // real scalar is two fxmuls rather than four plus a subtract, because the
    // imaginary cross terms are `fxmul(0, x)`, which is exactly zero. Same
    // products, same order, half the multiplies.
    if m[0][0].im == 0 && m[0][1].im == 0 && m[1][0].im == 0 && m[1][1].im == 0 {
        let (p00, p01) = (m[0][0].re, m[0][1].re);
        let (p10, p11) = (m[1][0].re, m[1][1].re);
        let mut i = 0usize;
        while i < dim {
            if i & tbit == 0 && (i & ctrl_mask) == ctrl_mask {
                let j = i | tbit;
                let a0 = amps[i];
                let a1 = amps[j];
                amps[i] = Amp {
                    re: fxmul(p00, a0.re) + fxmul(p01, a1.re),
                    im: fxmul(p00, a0.im) + fxmul(p01, a1.im),
                };
                amps[j] = Amp {
                    re: fxmul(p10, a0.re) + fxmul(p11, a1.re),
                    im: fxmul(p10, a0.im) + fxmul(p11, a1.im),
                };
            }
            i += 1;
        }
        return;
    }

    let mut i = 0usize;
    while i < dim {
        if i & tbit == 0 && (i & ctrl_mask) == ctrl_mask {
            let j = i | tbit;
            let a0 = amps[i];
            let a1 = amps[j];
            amps[i] = m[0][0].mul(a0).add(m[0][1].mul(a1));
            amps[j] = m[1][0].mul(a0).add(m[1][1].mul(a1));
        }
        i += 1;
    }
}

// ---------------------------------------------------------------------------
// State vector + conformance
// ---------------------------------------------------------------------------

/// A reconstructed fixed-point state vector — the hashed quantity.
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct StateVector {
    pub n_qubits: u8,
    pub amps: Vec<Amp>,
}

/// `splitmix64` — the pinned integer PRNG shared with the WAI deterministic generator; the
/// fixed-point sampler the shot-histogram determinism requires.
fn splitmix64(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)
}

impl StateVector {
    /// The `|k⟩` computational basis state on `n` qubits.
    pub fn basis(n: u8, k: usize) -> StateVector {
        let dim = 1usize << n;
        let mut amps = vec![Amp::ZERO; dim];
        amps[k % dim] = Amp::ONE;
        StateVector { n_qubits: n, amps }
    }

    /// Canonical serialization of the amplitudes: `re_le(i64) ‖ im_le(i64)` per
    /// amplitude, in index order — the bytes the statevector hash covers.
    pub fn canonical_bytes(&self) -> Vec<u8> {
        let mut out = Vec::with_capacity(self.amps.len() * 16);
        for a in &self.amps {
            out.extend_from_slice(&a.re.to_le_bytes());
            out.extend_from_slice(&a.im.to_le_bytes());
        }
        out
    }

    /// **statevector-equivalence** — the portable BLAKE3 identity of the whole
    /// reconstructed vector. Two conforming sinks agree on this iff every amplitude
    /// is byte-identical, on every machine, with no tolerance parameter.
    pub fn statevector_hash(&self) -> [u8; 32] {
        let mut h = blake3::Hasher::new();
        h.update(DOMAIN_STATEVECTOR);
        h.update(&[self.n_qubits]);
        h.update(&FRAC.to_le_bytes());
        h.update(&self.canonical_bytes());
        *h.finalize().as_bytes()
    }

    /// Per-basis-state probability weights `|amp|^2` at scale `2^(2·FRAC)`, `i128`
    /// (un-normalized; summing gives ≈ `2^(2·FRAC)`). No divide — the sampler walks
    /// these directly.
    pub fn prob_weights(&self) -> Vec<i128> {
        self.amps.iter().map(|a| a.norm2()).collect()
    }

    /// Sample `shots` computational-basis measurements with the pinned
    /// `splitmix64` seeded by `seed`; returns a per-basis-state count vector.
    /// Deterministic on every machine.
    pub fn sample_shots(&self, seed: u64, shots: u64) -> Vec<u64> {
        let weights = self.prob_weights();
        // prefix sums in i128; total ≈ 2^(2·FRAC)
        let mut cum = Vec::with_capacity(weights.len());
        let mut total: i128 = 0;
        for w in &weights {
            total += *w;
            cum.push(total);
        }
        let mut counts = vec![0u64; self.amps.len()];
        if total <= 0 {
            return counts;
        }
        let mut state = seed;
        for _ in 0..shots {
            let r = (splitmix64(&mut state) as u128 % total as u128) as i128;
            // first index whose prefix sum strictly exceeds r
            let idx = match cum.binary_search_by(|c| {
                if *c <= r { std::cmp::Ordering::Less } else { std::cmp::Ordering::Greater }
            }) {
                Ok(i) | Err(i) => i,
            };
            let slot = idx.min(counts.len() - 1);
            counts[slot] += 1;
        }
        counts
    }

    /// **shot-histogram-equivalence** — the portable BLAKE3 identity of a shot
    /// histogram at a pinned `(seed, shots)`.
    pub fn histogram_hash(&self, seed: u64, shots: u64) -> [u8; 32] {
        let counts = self.sample_shots(seed, shots);
        let mut h = blake3::Hasher::new();
        h.update(DOMAIN_HISTOGRAM);
        h.update(&[self.n_qubits]);
        h.update(&seed.to_le_bytes());
        h.update(&shots.to_le_bytes());
        for c in &counts {
            h.update(&c.to_le_bytes());
        }
        *h.finalize().as_bytes()
    }

    /// State **fidelity** `|⟨self|other⟩|^2` in the fixed-point floor (result at
    /// `2^FRAC`) — the quantum-information-theoretic similarity metric. For two
    /// exactly-reconstructed vectors it is a *diagnostic* (identical states give
    /// ≈ `ONE`); it is the metric a future float/behavioral quantum prior — which
    /// cannot back `intent=replicate` — would conform against. `i128` accumulation
    /// so the overlap sum never overflows.
    pub fn fidelity_fx(&self, other: &StateVector) -> i64 {
        assert_eq!(self.amps.len(), other.amps.len(), "fidelity needs equal dimension");
        let mut re: i128 = 0;
        let mut im: i128 = 0;
        for (a, b) in self.amps.iter().zip(&other.amps) {
            // ⟨a|b⟩ = Σ conj(a)·b
            let p = a.conj().mul(*b);
            re += p.re as i128;
            im += p.im as i128;
        }
        // |overlap|^2 at 2^FRAC: (re^2 + im^2) >> FRAC, re/im already at 2^FRAC
        ((re * re + im * im) >> FRAC) as i64
    }
}

// ---------------------------------------------------------------------------
// WQC container (magic "WQC1" or "WQC2"): the on-wire circuit object
// ---------------------------------------------------------------------------

/// A pinned measurement request carried in a `WQC` (section `0x03`), so a file can
/// declare "conform on the shot histogram at this `(seed, shots)`".
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub struct Measure {
    pub seed: u64,
    pub shots: u64,
}

const SECT_CONTRACT: u8 = 0x01;
const SECT_OPLOG: u8 = 0x02;
const SECT_MEASURE: u8 = 0x03;

/// The canonical measurement bytes, `{"seed":S,"shots":N,"basis":"computational"}`,
/// rendered directly, as the contract is.
fn measure_bytes(m: Measure) -> Vec<u8> {
    format!("{{\"seed\":{},\"shots\":{},\"basis\":\"computational\"}}", m.seed, m.shots).into_bytes()
}

/// The op-log of a circuit that passed `validate`: `u32le n_ops`, then per op
/// `u8 opcode | u8 n_ctrl | u8 target | u16le param | n_ctrl × u8`. `validate`
/// bounds both counts, so neither conversion can fail.
fn oplog_bytes(c: &Circuit) -> Vec<u8> {
    let n_ops = u32::try_from(c.ops.len()).expect("validate bounds the op count by u32::MAX");
    let mut out = Vec::new();
    out.extend_from_slice(&n_ops.to_le_bytes());
    for g in &c.ops {
        out.push(g.base.opcode());
        out.push(u8::try_from(g.controls.len()).expect("validate bounds a gate's controls by 255"));
        out.push(g.target);
        out.extend_from_slice(&g.param.to_le_bytes());
        out.extend_from_slice(&g.controls);
    }
    out
}

/// How many ops to reserve for an op-log of `len` bytes that claims `n_ops`.
/// The count is the file's claim, so it never sizes an allocation alone: every
/// op takes at least 5 bytes, so the section holds at most `(len − 4) / 5`.
fn oplog_capacity(n_ops: u32, len: usize) -> usize {
    usize::try_from(n_ops).unwrap_or(usize::MAX).min(len.saturating_sub(4) / 5)
}

/// Read an op-log: the ops, and how many bytes they took. Every op is
/// bounds-checked, and an opcode above 9 is malformed. Bytes after the last op
/// are the caller's to judge.
fn parse_oplog(b: &[u8]) -> Result<(Vec<Gate>, usize), QuantumError> {
    if b.len() < 4 {
        return Err(QuantumError::Malformed("oplog header".into()));
    }
    let n_ops = u32::from_le_bytes([b[0], b[1], b[2], b[3]]);
    let mut ops = Vec::with_capacity(oplog_capacity(n_ops, b.len()));
    let mut pos = 4;
    for _ in 0..n_ops {
        if b.len() - pos < 5 {
            return Err(QuantumError::Malformed("op header".into()));
        }
        let base = BaseGate::from_opcode(b[pos])
            .ok_or_else(|| QuantumError::Malformed(format!("opcode {}", b[pos])))?;
        let n_ctrl = usize::from(b[pos + 1]);
        let target = b[pos + 2];
        let param = u16::from_le_bytes([b[pos + 3], b[pos + 4]]);
        pos += 5;
        if b.len() - pos < n_ctrl {
            return Err(QuantumError::Malformed("op controls".into()));
        }
        let controls = b[pos..pos + n_ctrl].to_vec();
        pos += n_ctrl;
        ops.push(Gate { base, controls, target, param });
    }
    Ok((ops, pos))
}

/// The byte range of a section, `blob_start + off` for `len` bytes, or `None`
/// if it runs past `total` or its arithmetic overflows. Checked, so that it
/// cannot wrap where `usize` is 32 bits.
fn section_range(blob_start: usize, off: u32, len: u32, total: usize) -> Option<std::ops::Range<usize>> {
    let start = blob_start.checked_add(usize::try_from(off).ok()?)?;
    let end = start.checked_add(usize::try_from(len).ok()?)?;
    (end <= total).then_some(start..end)
}

/// How a member of a contract or measurement object reads in an error: its
/// JSON, cut at 64 characters, since it comes from an untrusted file.
fn shown(v: Option<&Json>) -> String {
    let Some(v) = v else {
        return "missing".to_string();
    };
    let s = v.to_string();
    match s.char_indices().nth(64) {
        Some((cut, _)) => format!("{}…", &s[..cut]),
        None => s,
    }
}

/// Read a contract section under the gateset its container's magic names, and
/// return its `n_qubits`.
///
/// The section must be an I-JSON object (extensions/quantum-sim §2, read by
/// [`ijson::parse`]), or it is a contract mismatch, and `n_qubits` must be an
/// unsigned JSON integer of at most 255, or the container is malformed. circuit/1 requires `ext`, `numeric`, `frac` (the integer 30)
/// and `gateset` to carry its values, and leaves member order, whitespace and
/// other members free. circuit2 requires the canonical bytes exactly.
fn read_contract(gateset: Gateset, raw: &[u8]) -> Result<u8, QuantumError> {
    let obj = ijson::parse(raw).map_err(|e| QuantumError::Contract(format!("the contract is not I-JSON: {e}")))?;
    if !matches!(obj, Json::Object(_)) {
        return Err(QuantumError::Contract(format!("the contract is not a JSON object: {}", shown(Some(&obj)))));
    }
    let n_qubits = || {
        obj.get("n_qubits")
            .and_then(Json::as_u64)
            .and_then(|n| u8::try_from(n).ok())
            .ok_or_else(|| {
                QuantumError::Malformed(format!("n_qubits {} is not an integer in 0..=255", shown(obj.get("n_qubits"))))
            })
    };
    match gateset {
        Gateset::V1 => {
            let text = |key: &str, value: &str| {
                if obj.get(key).and_then(Json::as_str) == Some(value) {
                    Ok(())
                } else {
                    Err(QuantumError::Contract(format!("{key} is {}, not {value:?}", shown(obj.get(key)))))
                }
            };
            text("ext", gateset.ext())?;
            text("numeric", gateset.numeric())?;
            if obj.get("frac").and_then(Json::as_u64) != Some(u64::from(FRAC)) {
                return Err(QuantumError::Contract(format!("frac is {}, not {FRAC}", shown(obj.get("frac")))));
            }
            text("gateset", gateset.name())?;
            n_qubits()
        }
        Gateset::V2 => {
            let n = n_qubits()?;
            if raw != gateset.contract_bytes(n).as_slice() {
                return Err(QuantumError::Contract(format!(
                    "not the canonical {} contract for {n} qubits",
                    gateset.capability()
                )));
            }
            Ok(n)
        }
    }
}

/// Read a measurement section. It must be I-JSON ([`ijson::parse`]). circuit/1 then
/// reads it as it always has: a missing or non-integer `seed` or `shots` reads
/// as 0.
/// circuit2 requires a `u64` `seed` and `shots` and the canonical bytes for
/// them, which fixes the basis as `computational`.
fn read_measure(gateset: Gateset, raw: &[u8]) -> Result<Measure, QuantumError> {
    let mv = ijson::parse(raw).map_err(|e| QuantumError::Malformed(format!("the measurement is not I-JSON: {e}")))?;
    let field = |key: &str| mv.get(key).and_then(Json::as_u64);
    match gateset {
        Gateset::V1 => Ok(Measure { seed: field("seed").unwrap_or(0), shots: field("shots").unwrap_or(0) }),
        Gateset::V2 => {
            let u64_field = |key: &str| {
                field(key).ok_or_else(|| QuantumError::Malformed(format!("measure {key} is {}, not a u64", shown(mv.get(key)))))
            };
            let m = Measure { seed: u64_field("seed")?, shots: u64_field("shots")? };
            if raw != measure_bytes(m).as_slice() {
                return Err(QuantumError::Malformed("not the canonical measure bytes".into()));
            }
            Ok(m)
        }
    }
}

/// Serialize a circuit, and optionally a measurement, to the container its
/// gateset names: `WQC1` for [`Gateset::V1`], `WQC2` for [`Gateset::V2`].
///
/// The layout is the same for both: the magic, `u16le` section count, a table
/// of `(u8 kind, u32le offset, u32le len)` with offsets relative to the end of
/// the table, then the contract (`0x01`), the op-log (`0x02`) and the
/// measurement (`0x03`). A circuit/1 circuit is written byte for byte as every
/// earlier release wrote it.
///
/// Fails, with the error `validate` gives, on a circuit its gateset does not
/// allow. A circuit that has simulated cannot fail here, unless its op-log runs
/// past what a `u32` section length can address.
pub fn to_wqc(c: &Circuit, measure: Option<Measure>) -> Result<Vec<u8>, QuantumError> {
    c.validate()?;
    let mut sections = vec![(SECT_CONTRACT, c.gateset.contract_bytes(c.n_qubits)), (SECT_OPLOG, oplog_bytes(c))];
    if let Some(m) = measure {
        sections.push((SECT_MEASURE, measure_bytes(m)));
    }

    let mut out = Vec::new();
    out.extend_from_slice(&c.gateset.magic());
    // At most three sections.
    out.extend_from_slice(&(sections.len() as u16).to_le_bytes());
    let mut off: u64 = 0;
    for (kind, data) in &sections {
        let (Ok(off32), Ok(len32)) = (u32::try_from(off), u32::try_from(data.len())) else {
            return Err(QuantumError::Invalid(format!(
                "section {kind:#04x}, {} bytes at offset {off}, does not fit the container's u32 offsets",
                data.len()
            )));
        };
        out.push(*kind);
        out.extend_from_slice(&off32.to_le_bytes());
        out.extend_from_slice(&len32.to_le_bytes());
        off += u64::from(len32);
    }
    for (_, data) in &sections {
        out.extend_from_slice(data);
    }
    Ok(out)
}

/// Parse a `WQC1` or `WQC2` container into a circuit, in the gateset its magic
/// names, and its optional measurement.
///
/// The checks run in this order, and the first to fail is the error:
///
/// 1. The magic is `WQC1` or `WQC2`, else malformed.
/// 2. Every section lies inside the file, with arithmetic that cannot wrap; no
///    known section (`0x01`, `0x02`, `0x03`) appears twice; the contract and
///    op-log are present. Unknown sections are ignored. Else malformed.
/// 3. The contract declares the capability the magic names. In `WQC1` it must
///    be an I-JSON object carrying circuit/1's `ext`, `numeric`, `frac` and
///    `gateset`, in any formatting; in `WQC2` it must be the canonical bytes.
///    Else a contract mismatch, except that an `n_qubits` that is not an
///    unsigned JSON integer in `0..=255` is malformed.
/// 4. The op-log parses within its section: every op lies inside it and has
///    an opcode of at most 9. In `WQC2` it must end exactly after its last op;
///    in `WQC1` any bytes after it are ignored. Else malformed.
/// 5. The circuit passes `validate` under its gateset, else invalid.
/// 6. The measurement, if present, is I-JSON, and in `WQC2` it is the
///    canonical bytes. Else malformed.
///
/// I-JSON is extensions/quantum-sim §2's reading of a section: UTF-8, no
/// member name twice in one object, no surrogate or noncharacter, no number of
/// magnitude 2^1024 − 2^970 or more, and at most 127 nested arrays and
/// objects, each decided from the section's text. [`reject_corpus`] holds a
/// payload for each refusal, and [`accept_corpus`] payloads at the edges of
/// the rules that must be read.
pub fn from_wqc(bytes: &[u8]) -> Result<(Circuit, Option<Measure>), QuantumError> {
    if bytes.len() < 6 {
        return Err(QuantumError::Malformed("magic".into()));
    }
    let gateset = Gateset::from_magic(&bytes[0..4]).ok_or_else(|| QuantumError::Malformed("magic".into()))?;
    let strict = match gateset {
        Gateset::V1 => false,
        Gateset::V2 => true,
    };

    let n_sections = usize::from(u16::from_le_bytes([bytes[4], bytes[5]]));
    let table_start = 6;
    let blob_start = n_sections
        .checked_mul(9)
        .and_then(|len| len.checked_add(table_start))
        .filter(|&end| end <= bytes.len())
        .ok_or_else(|| QuantumError::Malformed("section table".into()))?;
    let u32_at = |p: usize| u32::from_le_bytes([bytes[p], bytes[p + 1], bytes[p + 2], bytes[p + 3]]);
    let mut contract: Option<&[u8]> = None;
    let mut oplog: Option<&[u8]> = None;
    let mut measure_raw: Option<&[u8]> = None;
    for s in 0..n_sections {
        // Inside the table, which ends at blob_start ≤ bytes.len().
        let p = table_start + s * 9;
        let kind = bytes[p];
        let range = section_range(blob_start, u32_at(p + 1), u32_at(p + 5), bytes.len())
            .ok_or_else(|| QuantumError::Malformed(format!("section {kind:#04x} bounds")))?;
        let slot = match kind {
            SECT_CONTRACT => &mut contract,
            SECT_OPLOG => &mut oplog,
            SECT_MEASURE => &mut measure_raw,
            _ => continue, // unknown section: ignore (forward-compatible)
        };
        if slot.replace(&bytes[range]).is_some() {
            return Err(QuantumError::Malformed(format!("duplicate section {kind:#04x}")));
        }
    }
    let contract = contract.ok_or_else(|| QuantumError::Malformed("missing contract".into()))?;
    let oplog = oplog.ok_or_else(|| QuantumError::Malformed("missing oplog".into()))?;

    let n_qubits = read_contract(gateset, contract)?;

    let (ops, used) = parse_oplog(oplog)?;
    if strict && used != oplog.len() {
        return Err(QuantumError::Malformed(format!("{} bytes after the last op", oplog.len() - used)));
    }

    let circuit = Circuit { n_qubits, ops, gateset };
    circuit.validate()?;

    let measure = measure_raw.map(|raw| read_measure(gateset, raw)).transpose()?;
    Ok((circuit, measure))
}

/// **circuit-equivalence** — the portable BLAKE3 identity of a circuit:
/// `BLAKE3("wai:quantum-circuit\x01" ‖ contract ‖ op-log)`, over the canonical
/// contract bytes of its gateset and its op-log. The contract names the
/// capability, so the hash binds the gateset and its table: the same ops in
/// circuit/1 and circuit2 are different circuits. The receipt binds this
/// alongside the statevector hash.
///
/// Fails, as [`to_wqc`] does, on a circuit its gateset does not allow; a
/// circuit that has simulated cannot fail.
pub fn circuit_hash(c: &Circuit) -> Result<[u8; 32], QuantumError> {
    c.validate()?;
    let mut h = blake3::Hasher::new();
    h.update(DOMAIN_CIRCUIT);
    h.update(&c.gateset.contract_bytes(c.n_qubits));
    h.update(&oplog_bytes(c));
    Ok(*h.finalize().as_bytes())
}

// ---- The reject corpus ---------------------------------------------------------

/// A container from raw section-table entries `(kind, offset, len)` and a blob.
fn wqc_table(magic: &[u8; 4], table: &[(u8, u32, u32)], blob: &[u8]) -> Vec<u8> {
    let mut out = magic.to_vec();
    out.extend_from_slice(&u16::try_from(table.len()).expect("a short table").to_le_bytes());
    for &(kind, off, len) in table {
        out.push(kind);
        out.extend_from_slice(&off.to_le_bytes());
        out.extend_from_slice(&len.to_le_bytes());
    }
    out.extend_from_slice(blob);
    out
}

/// A container holding `sections` in order, each table entry pointing at its
/// own bytes.
fn wqc_raw(magic: &[u8; 4], sections: &[(u8, &[u8])]) -> Vec<u8> {
    let (mut table, mut blob) = (Vec::new(), Vec::new());
    for &(kind, data) in sections {
        let at = |n: usize| u32::try_from(n).expect("a small section");
        table.push((kind, at(blob.len()), at(data.len())));
        blob.extend_from_slice(data);
    }
    wqc_table(magic, &table, &blob)
}

/// One op-log entry, as the wire has it.
fn op(opcode: u8, controls: &[u8], target: u8, param: u16) -> Vec<u8> {
    let mut out = vec![opcode, u8::try_from(controls.len()).expect("a few controls"), target];
    out.extend_from_slice(&param.to_le_bytes());
    out.extend_from_slice(controls);
    out
}

/// An op-log holding `ops`.
fn oplog_of(ops: &[Vec<u8>]) -> Vec<u8> {
    let mut out = u32::try_from(ops.len()).expect("a few ops").to_le_bytes().to_vec();
    for o in ops {
        out.extend_from_slice(o);
    }
    out
}


/// `contract`, an object, with `members` added at its end.
fn with_members(contract: &[u8], members: &[u8]) -> Vec<u8> {
    let (body, close) = contract.split_at(contract.len() - 1);
    assert_eq!(close, b"}", "a contract is an object");
    [body, b",", members, b"}"].concat()
}

/// `bytes` with the first `from` replaced by `to`; `from` must be there.
fn swap(bytes: &[u8], from: &str, to: &str) -> Vec<u8> {
    let s = std::str::from_utf8(bytes).expect("a contract is text");
    assert!(s.contains(from), "{from:?} is not in {s}");
    s.replacen(from, to, 1).into_bytes()
}

/// The reject corpus of extensions/quantum-sim §2, as `(name, payload, class)`:
/// payloads a reader must refuse, each with the class [`QuantumError::kind`]
/// gives it. Each breaks one rule and keeps every rule checked before it, so
/// the check order of [`from_wqc`] decides the class; most are a two-qubit Bell
/// container with that one thing broken. The conformance corpus's `reject/`
/// tree is written from this list.
pub fn reject_corpus() -> Vec<(&'static str, Vec<u8>, &'static str)> {
    let bell = oplog_of(&[op(4, &[], 0, 0), op(1, &[0], 1, 0)]);
    let canon1 = Gateset::V1.contract_bytes(2);
    let canon2 = Gateset::V2.contract_bytes(2);
    let v1 = |contract: &[u8]| wqc_raw(b"WQC1", &[(SECT_CONTRACT, contract), (SECT_OPLOG, &bell)]);
    let v1_oplog = |oplog: &[u8]| wqc_raw(b"WQC1", &[(SECT_CONTRACT, &canon1), (SECT_OPLOG, oplog)]);
    let v2 = |contract: &[u8], oplog: &[u8], measure: Option<&[u8]>| {
        let mut sections = vec![(SECT_CONTRACT, contract), (SECT_OPLOG, oplog)];
        sections.extend(measure.map(|m| (SECT_MEASURE, m)));
        wqc_raw(b"WQC2", &sections)
    };
    let v2_contract = |contract: &[u8]| v2(contract, &bell, None);
    let v2_oplog = |oplog: &[u8]| v2(&canon2, oplog, None);
    let v2_measure = |measure: &[u8]| v2(&canon2, &bell, Some(measure));
    let mut magic_wqc3 = v2_contract(&canon2);
    magic_wqc3[0..4].copy_from_slice(b"WQC3");
    // An op-log of 9 bytes that claims u32::MAX ops.
    let n_ops_u32max = [u32::MAX.to_le_bytes().as_slice(), &op(4, &[], 0, 0)].concat();
    // An op-log section whose offset plus length wraps a u32.
    let wrap = wqc_table(
        b"WQC1",
        &[(SECT_CONTRACT, 0, u32::try_from(canon1.len()).expect("a short contract")), (SECT_OPLOG, 0xFFFF_FFF0, 0x20)],
        &[canon1.as_slice(), &bell].concat(),
    );
    // The Bell container cut off inside its section table: two entries claimed, one present.
    let mut table_past_end = v1(&canon1);
    table_past_end.truncate(6 + 9);
    let at = |n: usize| u32::try_from(n).expect("a short section");
    let (c1, b1) = (at(canon1.len()), at(bell.len()));
    let both = [canon1.as_slice(), &bell].concat();
    // A contract whose last member's string holds the byte 0xFF.
    let not_utf8 = [&canon1[..canon1.len() - 1], b",\"note\":\"\xff\"}"].concat();
    // An op that names two controls where the section holds one.
    let truncated_controls = [1u32.to_le_bytes().as_slice(), &[1, 2, 1, 0, 0, 0]].concat();
    let measure = b"{\"seed\":1,\"shots\":10,\"basis\":\"computational\"}";
    // A circuit/1 contract nesting 127 arrays inside its object: 128 levels,
    // one past what §2 allows.
    let deep = [&canon1[..canon1.len() - 1], b",\"x\":", &b"[".repeat(127), &b"]".repeat(127), b"}"].concat();
    let reordered2 = format!(
        r#"{{"n_qubits":2,"ext":"{}","numeric":"{}","frac":30,"gateset":"{}"}}"#,
        Gateset::V2.ext(),
        Gateset::V2.numeric(),
        Gateset::V2.name()
    );
    vec![
        ("magic_wqc3", magic_wqc3, "malformed"),
        ("v1_ext_other", v1(&swap(&canon1, "circuit/1", "circuit/2")), "contract"),
        ("v1_numeric_other", v1(&swap(&canon1, "wai.det.fixed64", "wai.det.amp64")), "contract"),
        ("v1_frac_31", v1(&swap(&canon1, "\"frac\":30", "\"frac\":31")), "contract"),
        ("v1_frac_float", v1(&swap(&canon1, "\"frac\":30", "\"frac\":30.0")), "contract"),
        ("v1_gateset_rn", v1(&swap(&canon1, "dyadicP\"", "dyadicP.rn\"")), "contract"),
        ("v1_contract2_in_wqc1", v1(&canon2), "contract"),
        ("v1_fields_missing", v1(br#"{"n_qubits":3}"#), "contract"),
        ("v1_contract_not_object", v1(&[b"[".as_slice(), &canon1, b"]"].concat()), "contract"),
        ("v1_n_qubits_256", v1(&swap(&canon1, "\"n_qubits\":2", "\"n_qubits\":256")), "malformed"),
        ("v1_n_qubits_258", v1(&swap(&canon1, "\"n_qubits\":2", "\"n_qubits\":258")), "malformed"),
        (
            "v1_contract_twice",
            wqc_raw(b"WQC1", &[(SECT_CONTRACT, &canon1), (SECT_OPLOG, &bell), (SECT_CONTRACT, &Gateset::V1.contract_bytes(3))]),
            "malformed",
        ),
        (
            "v1_oplog_twice",
            wqc_raw(b"WQC1", &[(SECT_CONTRACT, &canon1), (SECT_OPLOG, &bell), (SECT_OPLOG, &oplog_of(&[op(1, &[], 0, 0)]))]),
            "malformed",
        ),
        // §2.1: the contract is I-JSON.
        ("v1_contract_duplicate_member", v1(&swap(&canon1, "\"frac\":30", "\"frac\":30,\"frac\":30")), "contract"),
        ("v1_contract_lone_surrogate", v1(&swap(&canon1, "}", ",\"note\":\"\\ud800\"}")), "contract"),
        ("v1_contract_too_deep", v1(&deep), "contract"),
        ("v1_n_qubits_minus0", v1(&swap(&canon1, "\"n_qubits\":2", "\"n_qubits\":-0")), "malformed"),
        ("v1_contract_lone_low_surrogate", v1(&with_members(&canon1, br#""note":"\udc00""#)), "contract"),
        ("v1_contract_noncharacter", v1(&with_members(&canon1, br#""note":"\ufdd0""#)), "contract"),
        ("v1_contract_raw_noncharacter", v1(&with_members(&canon1, "\"note\":\"\u{ffff}\"".as_bytes())), "contract"),
        (
            "v1_contract_integer_rounds_to_infinity",
            v1(&with_members(&canon1, format!("\"x\":{RANGE_LIMIT}").as_bytes())),
            "contract",
        ),
        ("v1_contract_number_rounds_to_infinity", v1(&with_members(&canon1, br#""x":1.7976931348623159e308"#)), "contract"),
        // 10^309, written with a run of zeros its exponent offsets.
        (
            "v1_contract_long_literal",
            v1(&with_members(&canon1, format!("\"x\":0.{}1e655669", "0".repeat(655_359)).as_bytes())),
            "contract",
        ),
        ("v1_contract_bom", v1(&[b"\xef\xbb\xbf".as_slice(), &canon1].concat()), "contract"),
        ("v1_contract_nan", v1(&with_members(&canon1, br#""x":NaN"#)), "contract"),
        ("v1_contract_control_character", v1(&with_members(&canon1, b"\"note\":\"a\x01b\"")), "contract"),
        // §2.3: so is the measurement.
        (
            "v1_measure_duplicate_member",
            wqc_raw(b"WQC1", &[(SECT_CONTRACT, &canon1), (SECT_OPLOG, &bell), (SECT_MEASURE, br#"{"seed":1,"seed":2,"shots":10}"#)]),
            "malformed",
        ),
        ("v1_contract_not_utf8", v1(&not_utf8), "contract"),
        ("v1_contract_number_out_of_range", v1(&swap(&canon1, "}", ",\"x\":1e400}")), "contract"),
        ("v1_n_qubits_fraction", v1(&swap(&canon1, "\"n_qubits\":2", "\"n_qubits\":2.0")), "malformed"),
        ("v1_n_qubits_missing", v1(&swap(&canon1, "\"n_qubits\":2,", "")), "malformed"),
        (
            "v1_measure_not_json",
            wqc_raw(b"WQC1", &[(SECT_CONTRACT, &canon1), (SECT_OPLOG, &bell), (SECT_MEASURE, br#"{"seed":1,"#)]),
            "malformed",
        ),
        // §2: the payload, its table and its sections.
        ("short_payload", b"WQC2\x00".to_vec(), "malformed"),
        ("table_past_end", table_past_end, "malformed"),
        ("section_past_end", wqc_table(b"WQC1", &[(SECT_CONTRACT, 0, c1), (SECT_OPLOG, c1, b1 + 1)], &both), "malformed"),
        (
            "unknown_section_past_end",
            wqc_table(b"WQC1", &[(SECT_CONTRACT, 0, c1), (SECT_OPLOG, c1, b1), (0x7F, c1 + b1, 1)], &both),
            "malformed",
        ),
        ("no_contract", wqc_raw(b"WQC1", &[(SECT_OPLOG, &bell)]), "malformed"),
        ("no_oplog", wqc_raw(b"WQC1", &[(SECT_CONTRACT, &canon1)]), "malformed"),
        // §2.2: the op-log parses, and the circuit keeps its rules.
        ("v1_oplog_short", v1_oplog(&[0, 0, 0]), "malformed"),
        ("v1_op_truncated_controls", v1_oplog(&truncated_controls), "malformed"),
        ("v1_opcode_10", v1_oplog(&oplog_of(&[op(10, &[], 0, 0)])), "malformed"),
        ("v1_target_out_of_range", v1_oplog(&oplog_of(&[op(4, &[], 2, 0)])), "invalid"),
        ("v1_control_out_of_range", v1_oplog(&oplog_of(&[op(1, &[2], 1, 0)])), "invalid"),
        ("v1_control_is_target", v1_oplog(&oplog_of(&[op(1, &[1], 1, 0)])), "invalid"),
        ("v1_phase_0", v1_oplog(&oplog_of(&[op(9, &[], 0, 0)])), "invalid"),
        ("v1_phase_33", v1_oplog(&oplog_of(&[op(9, &[], 0, 33)])), "invalid"),
        ("v1_n_ops_u32max", v1_oplog(&n_ops_u32max), "malformed"),
        ("v1_offset_wrap", wrap, "malformed"),
        ("v2_contract_whitespace", v2_contract(&swap(&canon2, "\"n_qubits\":2", "\"n_qubits\": 2")), "contract"),
        ("v2_contract_key_order", v2_contract(reordered2.as_bytes()), "contract"),
        ("v2_contract_extra_member", v2_contract(&swap(&canon2, "}", ",\"note\":\"x\"}")), "contract"),
        ("v2_contract1_in_wqc2", v2_contract(&canon1), "contract"),
        ("v2_n_qubits_256", v2_contract(&swap(&canon2, "\"n_qubits\":2", "\"n_qubits\":256")), "malformed"),
        ("v2_contract_not_object", v2_contract(&[b"[".as_slice(), &canon2, b"]"].concat()), "contract"),
        // The I-JSON rules come before n_qubits: out of range is a contract
        // refusal, in range but not an integer a malformed one.
        (
            "v2_n_qubits_rounds_to_infinity",
            v2_contract(&swap(&canon2, "\"n_qubits\":2", "\"n_qubits\":1.7976931348623159e308")),
            "contract",
        ),
        (
            "v2_n_qubits_largest_finite",
            v2_contract(&swap(&canon2, "\"n_qubits\":2", "\"n_qubits\":1.7976931348623158e308")),
            "malformed",
        ),
        ("v2_phase_0", v2_oplog(&oplog_of(&[op(9, &[], 0, 0)])), "invalid"),
        (
            "v2_measure_twice",
            wqc_raw(b"WQC2", &[(SECT_CONTRACT, &canon2), (SECT_OPLOG, &bell), (SECT_MEASURE, measure), (SECT_MEASURE, measure)]),
            "malformed",
        ),
        ("v2_measure_seed_signed", v2_measure(br#"{"seed":-1,"shots":10,"basis":"computational"}"#), "malformed"),
        (
            "v2_duplicate_controls",
            v2(&Gateset::V2.contract_bytes(3), &oplog_of(&[op(1, &[0, 0], 2, 0)]), None),
            "invalid",
        ),
        ("v2_param_on_h", v2_oplog(&oplog_of(&[op(4, &[], 0, 7), op(1, &[0], 1, 0)])), "invalid"),
        ("v2_oplog_trailing_byte", v2_oplog(&[bell.as_slice(), &[0]].concat()), "malformed"),
        ("v2_opcode_10", v2_oplog(&oplog_of(&[op(10, &[], 0, 0)])), "malformed"),
        ("v2_measure_spaced", v2_measure(br#"{"seed": 1, "shots": 10, "basis": "computational"}"#), "malformed"),
        ("v2_measure_no_shots", v2_measure(br#"{"seed":1,"basis":"computational"}"#), "malformed"),
        ("v2_measure_basis_x", v2_measure(br#"{"seed":1,"shots":10,"basis":"x"}"#), "malformed"),
    ]
}

/// The accept corpus of extensions/quantum-sim §2, as `(name, payload, circuit,
/// measurement)`: payloads a reader must accept, each with the circuit and the
/// measurement it reads as. They sit at the edges the reject corpus does not
/// reach from the other side: circuit/1's free formatting and further members,
/// every I-JSON rule at its limit, circuit/1's lax measurement, op-log bytes
/// after the last op, unknown sections, and the widest `n_qubits` and `u64`
/// values. The conformance corpus's `accept/` tree is written from this list.
pub fn accept_corpus() -> Vec<(&'static str, Vec<u8>, Circuit, Option<Measure>)> {
    let bell_ops = |gateset| Circuit {
        n_qubits: 2,
        ops: vec![
            Gate { base: BaseGate::H, controls: vec![], target: 0, param: 0 },
            Gate { base: BaseGate::X, controls: vec![0], target: 1, param: 0 },
        ],
        gateset,
    };
    let (bell1, bell2) = (bell_ops(Gateset::V1), bell_ops(Gateset::V2));
    let wide = |c: &Circuit| Circuit { n_qubits: 255, ..c.clone() };
    let bell = oplog_of(&[op(4, &[], 0, 0), op(1, &[0], 1, 0)]);
    let canon1 = Gateset::V1.contract_bytes(2);
    let canon2 = Gateset::V2.contract_bytes(2);
    let v1 = |contract: &[u8]| wqc_raw(b"WQC1", &[(SECT_CONTRACT, contract), (SECT_OPLOG, &bell)]);
    let v1_measure = |measure: &[u8]| wqc_raw(b"WQC1", &[(SECT_CONTRACT, &canon1), (SECT_OPLOG, &bell), (SECT_MEASURE, measure)]);
    let m = |seed, shots| Some(Measure { seed, shots });
    let reordered = format!(
        r#"{{"gateset":"{}","frac":30,"numeric":"{}","n_qubits":2,"ext":"{}"}}"#,
        Gateset::V1.name(),
        Gateset::V1.numeric(),
        Gateset::V1.ext()
    );
    let spaced = String::from_utf8(canon1.clone()).expect("a contract is text").replace(':', " : ").replace(',', " ,\n\t");
    let spaced = format!(" \r\n{spaced}\r\n ");
    let depth = format!(r#""x":{}{}"#, "[".repeat(126), "]".repeat(126));
    let largest_integer = format!("{}1", &RANGE_LIMIT[..RANGE_LIMIT.len() - 1]);
    let unknown = |magic: &[u8; 4], contract: &[u8]| {
        wqc_raw(magic, &[(0x7F, b"later"), (SECT_CONTRACT, contract), (SECT_OPLOG, &bell), (0x04, b"")])
    };
    let mut param_on_h = bell1.clone();
    param_on_h.ops[0].param = 7;
    let mut repeated_control = bell1.clone();
    repeated_control.ops[1].controls = vec![0, 0];
    vec![
        ("v1_members_reordered", v1(reordered.as_bytes()), bell1.clone(), None),
        ("v1_whitespace", v1(spaced.as_bytes()), bell1.clone(), None),
        (
            "v1_further_members",
            v1(&with_members(&canon1, br#""note":"x","list":[0,-1,2.5,-3e-2,1E+2,true,false,null,{},[]],"obj":{"a":{"b":[]}},"again":{"note":1}"#)),
            bell1.clone(),
            None,
        ),
        (
            "v1_escapes_and_characters",
            v1(&with_members(&canon1, "\"note\":\"\\\"\\\\\\/\\b\\f\\n\\r\\t\\u00e9\\ud83d\\ude00\\udbff\\udffd\u{e9}\u{fdcf}\u{fdf0}\u{fffd}\u{10fffd}\"".as_bytes())),
            bell1.clone(),
            None,
        ),
        ("v1_depth_127", v1(&with_members(&canon1, depth.as_bytes())), bell1.clone(), None),
        (
            "v1_numbers_largest_finite",
            v1(&with_members(&canon1, format!(r#""x":1.7976931348623158e308,"y":-1.7976931348623157e308,"z":{largest_integer},"w":-{largest_integer}"#).as_bytes())),
            bell1.clone(),
            None,
        ),
        // Exactly 1 and exactly 10, written long: a run of zeros the exponent
        // offsets, and an exponent with leading zeros.
        (
            "v1_numbers_long_literals",
            v1(&with_members(&canon1, format!("\"x\":1{}e-655360,\"y\":1e{}1", "0".repeat(655_360), "0".repeat(4300)).as_bytes())),
            bell1.clone(),
            None,
        ),
        (
            "v1_numbers_toward_zero",
            v1(&with_members(&canon1, br#""x":1e-400,"y":-0,"z":0e99999999999999999999,"w":-0.0E-0"#)),
            bell1.clone(),
            None,
        ),
        ("v1_n_qubits_255", v1(&Gateset::V1.contract_bytes(255)), wide(&bell1), None),
        ("v1_measure_lax", v1_measure(br#" { "basis" : "x" , "shots" : 10 , "seed" : 1 } "#), bell1.clone(), m(1, 10)),
        ("v1_measure_missing_shots", v1_measure(br#"{"seed":1,"basis":"computational"}"#), bell1.clone(), m(1, 0)),
        (
            "v1_measure_u64_max",
            v1_measure(br#"{"seed":18446744073709551615,"shots":18446744073709551615}"#),
            bell1.clone(),
            m(u64::MAX, u64::MAX),
        ),
        (
            "v1_measure_not_u64",
            v1_measure(br#"{"seed":18446744073709551616,"shots":-0,"x":2.5}"#),
            bell1.clone(),
            m(0, 0),
        ),
        ("v1_measure_other_types", v1_measure(br#"{"seed":"1","shots":1e1,"basis":null}"#), bell1.clone(), m(0, 0)),
        ("v1_measure_not_object", v1_measure(b"[1,2]"), bell1.clone(), m(0, 0)),
        (
            "v1_oplog_trailing_bytes",
            wqc_raw(b"WQC1", &[(SECT_CONTRACT, &canon1), (SECT_OPLOG, &[bell.as_slice(), &[0, 1, 2]].concat())]),
            bell1.clone(),
            None,
        ),
        (
            "v1_param_on_h",
            wqc_raw(b"WQC1", &[(SECT_CONTRACT, &canon1), (SECT_OPLOG, &oplog_of(&[op(4, &[], 0, 7), op(1, &[0], 1, 0)]))]),
            param_on_h,
            None,
        ),
        (
            "v1_repeated_control",
            wqc_raw(b"WQC1", &[(SECT_CONTRACT, &canon1), (SECT_OPLOG, &oplog_of(&[op(4, &[], 0, 0), op(1, &[0, 0], 1, 0)]))]),
            repeated_control,
            None,
        ),
        ("v1_unknown_sections", unknown(b"WQC1", &canon1), bell1.clone(), None),
        ("v2_unknown_sections", unknown(b"WQC2", &canon2), bell2.clone(), None),
        (
            "v2_n_qubits_255",
            wqc_raw(b"WQC2", &[(SECT_CONTRACT, &Gateset::V2.contract_bytes(255)), (SECT_OPLOG, &bell)]),
            wide(&bell2),
            None,
        ),
        (
            "v2_measure_u64_max",
            wqc_raw(
                b"WQC2",
                &[
                    (SECT_CONTRACT, &canon2),
                    (SECT_OPLOG, &bell),
                    (SECT_MEASURE, br#"{"seed":18446744073709551615,"shots":18446744073709551615,"basis":"computational"}"#),
                ],
            ),
            bell2,
            m(u64::MAX, u64::MAX),
        ),
    ]
}

#[cfg(test)]
mod tests {
    use super::*;
    /// The pre-optimisation `apply`, kept verbatim as an oracle. The fast paths
    /// are justified by an argument (`fxmul(ONE,x) == x`, `fxmul(0,x) == 0`), and
    /// an argument about fixed-point rounding is exactly the kind that is
    /// convincing and wrong, so it is checked against the code it replaced.
    fn apply_reference(amps: &mut [Amp], n: u8, g: &Gate, gs: Gateset) {
        let m = g.base.matrix(g.param, gs);
        let tbit = 1usize << g.target;
        let ctrl_mask: usize = g.controls.iter().fold(0usize, |acc, &c| acc | (1usize << c));
        let dim = 1usize << n;
        let mut i = 0usize;
        while i < dim {
            if i & tbit == 0 && (i & ctrl_mask) == ctrl_mask {
                let j = i | tbit;
                let a0 = amps[i];
                let a1 = amps[j];
                amps[i] = m[0][0].mul(a0).add(m[0][1].mul(a1));
                amps[j] = m[1][0].mul(a0).add(m[1][1].mul(a1));
            }
            i += 1;
        }
    }


    /// Optimised `apply` against the pre-optimisation oracle, in one process,
    /// interleaved, best-of-N. Comparing separate runs on a busy machine produced
    /// swings big enough to invent speedups and regressions that were not there —
    /// including an apparent slowdown in a code path that had not changed.
    #[test]
    #[ignore]
    fn probe_apply_speedup() {
        use std::time::Instant;
        let n = 20u8;
        let dim = 1usize << n;
        let reps = 5;
        let gates = 60;
        println!("\\n  n={n}  best of {reps}, interleaved");
        println!("  gate   reference     optimised     speedup");
        for kind in ["h", "t", "z", "s", "x", "cx", "cz", "y"] {
            let mut ops: Vec<Gate> = Vec::new();
            for i in 0..gates {
                let q = (i % n as usize) as u8;
                let o = (q + 1) % n;
                ops.push(match kind {
                    "h" => Gate { base: BaseGate::H, controls: vec![], target: q, param: 0 },
                    "t" => Gate { base: BaseGate::T, controls: vec![], target: q, param: 0 },
                    "z" => Gate { base: BaseGate::Z, controls: vec![], target: q, param: 0 },
                    "s" => Gate { base: BaseGate::S, controls: vec![], target: q, param: 0 },
                    "x" => Gate { base: BaseGate::X, controls: vec![], target: q, param: 0 },
                    "y" => Gate { base: BaseGate::Y, controls: vec![], target: q, param: 0 },
                    "cx" => Gate { base: BaseGate::X, controls: vec![o], target: q, param: 0 },
                    _ => Gate { base: BaseGate::Z, controls: vec![o], target: q, param: 0 },
                });
            }
            let (mut best_ref, mut best_opt) = (u128::MAX, u128::MAX);
            for _ in 0..reps {
                let mut a = vec![Amp::ZERO; dim];
                a[0] = Amp::ONE;
                let t0 = Instant::now();
                for g in &ops { apply_reference(&mut a, n, g, Gateset::V2); }
                best_ref = best_ref.min(t0.elapsed().as_nanos());
                std::hint::black_box(&a);

                let mut b = vec![Amp::ZERO; dim];
                b[0] = Amp::ONE;
                let t1 = Instant::now();
                for g in &ops { apply(&mut b, n, g, Gateset::V2); }
                best_opt = best_opt.min(t1.elapsed().as_nanos());
                std::hint::black_box(&b);
            }
            let (r, o) = (best_ref as f64 / gates as f64, best_opt as f64 / gates as f64);
            println!("  {kind:<5} {r:>10.0} ns {o:>10.0} ns {:>9.2}x", r / o);
        }
    }

    /// Every amplitude, not just the hash: random circuits over the whole gate set
    /// must come out of the optimised path bit-for-bit as they came out of the old
    /// one. A faster simulator that moved one low bit would have broken every
    /// receipt this crate has ever signed. Both tables, and every `P(k)`: in
    /// circuit2, `P(k ≥ 18)` has `re == ONE` with `im ≠ 0`, which must take the
    /// diagonal path and nothing else.
    #[test]
    fn optimised_apply_is_bit_identical_to_the_general_path() {
        let mut st = 0x1234_5678_9ABC_DEF0u64;
        let mut rnd = |m: usize| {
            st = st.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
            ((st >> 33) as usize) % m.max(1)
        };
        for (case, gs) in Gateset::ALL.into_iter().flat_map(|gs| std::iter::repeat_n(gs, 300)).enumerate() {
            let n = 3 + rnd(3) as u8;
            let dim = 1usize << n;
            let mut c = Circuit::with_gateset(n, gs);
            for _ in 0..24 {
                let t = rnd(n as usize) as u8;
                let mut other = rnd(n as usize) as u8;
                if other == t {
                    other = (t + 1) % n;
                }
                match rnd(13) {
                    0 => { c.x(t); }
                    1 => { c.y(t); }
                    2 => { c.z(t); }
                    3 => { c.h(t); }
                    4 => { c.s(t); }
                    5 => { c.t(t); }
                    6 => { c.p(1 + rnd(DYADIC_MAX) as u16, t); }
                    7 => { c.cx(other, t); }
                    8 => { c.cz(other, t); }
                    9 => { c.cp(1 + rnd(DYADIC_MAX) as u16, other, t); }
                    10 => c.ops.push(Gate { base: BaseGate::Sdg, controls: vec![], target: t, param: 0 }),
                    11 => c.ops.push(Gate { base: BaseGate::Tdg, controls: vec![], target: t, param: 0 }),
                    _ => c.ops.push(Gate { base: BaseGate::I, controls: vec![], target: t, param: 0 }),
                }
            }
            // controlled-Y and Toffoli too, which exercise the antidiagonal and
            // multi-control paths
            if n >= 3 {
                c.ops.push(Gate { base: BaseGate::Y, controls: vec![0], target: n - 1, param: 0 });
                c.ccx(0, 1, n - 1);
            }

            let mut want = vec![Amp::ZERO; dim];
            want[0] = Amp::ONE;
            for g in &c.ops {
                apply_reference(&mut want, n, g, gs);
            }
            let got = c.simulate().unwrap();
            assert_eq!(got.amps, want, "case {case}: optimised apply diverged on {} gates, n={n}, {gs:?}", c.ops.len());
        }
    }


    /// Tolerance for "≈ 1.0" checks — fixed-point rounding accumulates a little,
    /// but determinism (not accuracy) is the contract. ~1e-4 of ONE.
    const TOL: i64 = ONE / 10_000;

    fn approx(a: i64, b: i64, tol: i64) -> bool {
        (a - b).abs() <= tol
    }

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

    /// The conformance corpus's `dyadic_ladder`: X(0); then Y, Z, S, S†, T, T†
    /// and I on qubit 1; then CP(k; control 0, target 1) for k = 1..=32, twice.
    /// Every table entry multiplies a nonzero amplitude, so a change of one unit
    /// in any entry moves the statevector hash.
    fn dyadic_ladder(gs: Gateset) -> Circuit {
        let mut c = Circuit::with_gateset(2, gs);
        c.x(0).y(1).z(1).s(1);
        c.push(BaseGate::Sdg, vec![], 1, 0);
        c.t(1);
        c.push(BaseGate::Tdg, vec![], 1, 0);
        c.push(BaseGate::I, vec![], 1, 0);
        for _ in 0..2 {
            for k in 1..=DYADIC_MAX as u16 {
                c.cp(k, 0, 1);
            }
        }
        c
    }

    /// circuit/1's table is the recurrence its artefacts were made with, and
    /// stays that way. The defect is pinned on purpose, so that nobody repairs
    /// circuit/1 in place: a repaired table would stop every circuit/1 file and
    /// receipt from reproducing. The repair is circuit2.
    #[test]
    fn phase_table_v1_is_frozen() {
        // The recurrence: truncating halving, floor square root. Both halves are
        // non-negative, so truncating division is floor division.
        let mut t = [Amp::ZERO; DYADIC_MAX + 1];
        t[1] = Amp { re: -ONE, im: 0 };
        t[2] = Amp { re: 0, im: ONE };
        for k in 3..=DYADIC_MAX {
            let c = t[k - 1].re;
            t[k] = Amp { re: sqrt_fx((ONE + c) / 2), im: sqrt_fx((ONE - c) / 2) };
        }
        assert_eq!(t, PHASE_V1, "PHASE_V1 is not the circuit/1 recurrence");
        assert_eq!(
            hx(&Gateset::V1.table_digest()),
            "3fa9805ff96feb819a1cf371fedcd3b5228f93eddf933d0edcf401cbf8c28466"
        );
        // From k = 19 on, circuit/1 applies no phase at all.
        for k in 19..=DYADIC_MAX {
            assert_eq!(PHASE_V1[k], Amp { re: ONE - 1, im: 0 }, "PHASE_V1[{k}]");
        }
        assert_eq!(Gateset::V1.inv_sqrt2(), 759_250_124);
        // End to end, so the unit suite alone catches circuit/1 drift: the ladder
        // reproduces its circuit/1 corpus hash.
        let sv = dyadic_ladder(Gateset::V1).simulate().unwrap();
        assert_eq!(
            hx(&sv.statevector_hash()),
            "9754af24edd026601297a0f3334eb7171174978841a2938c9b372db942a3851d"
        );
    }

    /// Every circuit2 entry is the nearest integer, proved with integers alone.
    ///
    /// Enclose `2^62·cos(2π/2^k)` and `2^62·sin(2π/2^k)` between integer bounds,
    /// starting exactly at `k = 2`, and halve the angle with rounding chosen to
    /// keep each bound on its side: `cos(θ/2) = √((1 + cos θ)/2)` and
    /// `sin(θ/2) = sin θ / (2 cos(θ/2))`, a quotient, so nothing cancels. If both
    /// ends of an enclosure round to the same integer at `2^-30`, the true value
    /// does too, and that integer is within ½ unit of it. Every endpoint is also
    /// required to sit more than `2^20` ulps from a tie.
    #[test]
    fn phase_table_v2_is_certified_nearest() {
        const S: u128 = 1 << 62;
        let mul = |a: u128, b: u128| a.checked_mul(b).expect("the enclosure overflowed u128");
        let ceil_isqrt = |x: u128| {
            let r = x.isqrt();
            if r * r == x { r } else { r + 1 }
        };
        // Round half up from 2^-62 to 2^-30.
        let r = |x: u128| ((x + (1 << 31)) >> 32) as i64;
        let far_from_tie = |x: u128| ((x & 0xFFFF_FFFF) as i128 - (1i128 << 31)).abs() > 1 << 20;

        // (cos, sin)(π/2) = (0, 1), exactly.
        let (mut c_lo, mut c_hi, mut s_lo, mut s_hi) = (0u128, 0u128, S, S);
        for k in 3..=DYADIC_MAX {
            let c_lo2 = (mul(S + c_lo, S) >> 1).isqrt();
            let c_hi2 = ceil_isqrt(mul(S + c_hi, S).div_ceil(2));
            let s_lo2 = mul(s_lo, S) / (2 * c_hi2);
            let s_hi2 = mul(s_hi, S).div_ceil(2 * c_lo2);
            assert!(c_lo2 <= c_hi2 && s_lo2 <= s_hi2, "P({k}): an empty enclosure");
            let want = PHASE_V2[k];
            assert_eq!((r(c_lo2), r(c_hi2)), (want.re, want.re), "re P({k}) is not certified");
            assert_eq!((r(s_lo2), r(s_hi2)), (want.im, want.im), "im P({k}) is not certified");
            for x in [c_lo2, c_hi2, s_lo2, s_hi2] {
                assert!(far_from_tie(x), "P({k}): an endpoint lies within 2^20 ulps of a tie");
            }
            (c_lo, c_hi, s_lo, s_hi) = (c_lo2, c_hi2, s_lo2, s_hi2);
        }
        assert_eq!(PHASE_V2[0], Amp::ZERO);
        assert_eq!(PHASE_V2[1], Amp { re: -ONE, im: 0 });
        assert_eq!(PHASE_V2[2], Amp { re: 0, im: ONE });
        assert_eq!(
            hx(&Gateset::V2.table_digest()),
            "d34cb793715348e7c92db027c0bc3dddb6008f5730470bafaed978ea7fb23510"
        );
        assert_eq!(Gateset::V2.inv_sqrt2(), 759_250_125);
    }

    /// Every circuit2 component is within one unit of the ideal, against an
    /// independent floating-point reference, and the same check fails on
    /// circuit/1, so it can tell a good table from a bad one.
    ///
    /// The reference is the crate's own `sin_cos`, which gives the same bits on
    /// every IEEE platform and is within an ulp; it is compiled under these
    /// features, all of which `full` enables.
    #[cfg(any(
        feature = "quantum_spd",
        feature = "quantum_pauli",
        feature = "quantum_mps",
        feature = "quantum_vml",
        feature = "quantum_phasor",
        feature = "quantum_kernel",
        feature = "quantum_tn",
        feature = "quantum_sv",
        feature = "quantum_tdvp",
        feature = "quantum_nqs"
    ))]
    #[test]
    fn every_v2_phase_is_within_one_unit_of_ideal() {
        // ½ unit, plus slack for the reference's own error (about 1e-7 unit).
        const BOUND: f64 = 0.5 + 1e-6;
        // The largest deviation of each entry from 2^30·(cos, sin)(2π/2^k).
        let deviation = |t: &[Amp; 33]| -> Vec<(usize, f64)> {
            (1..=DYADIC_MAX)
                .map(|k| {
                    let (s, c) = crate::repro::sin_cos(std::f64::consts::TAU / (1u64 << k) as f64);
                    let one = ONE as f64;
                    let dre = (t[k].re as f64 - c * one).abs();
                    let dim = (t[k].im as f64 - s * one).abs();
                    (k, dre.max(dim))
                })
                .collect()
        };
        let mut worst = 0.0f64;
        for (k, d) in deviation(&PHASE_V2) {
            assert!(d <= BOUND, "PHASE_V2[{k}] is {d} units from the ideal");
            worst = worst.max(d);
        }
        assert!(worst < 1.0);
        // Negative control: circuit/1 fails the same check from k = 8 on (its
        // P(8) imaginary part is 9.52 units high).
        let v1_bad: Vec<usize> =
            deviation(&PHASE_V1).into_iter().filter(|&(_, d)| d > BOUND).map(|(k, _)| k).collect();
        assert!(v1_bad.contains(&8), "the check passed circuit/1's P(8): {v1_bad:?}");
        assert!((8..=DYADIC_MAX).all(|k| v1_bad.contains(&k)), "circuit/1 failures: {v1_bad:?}");
    }

    /// The identities each table must keep: `P(1) = −1` and `P(2) = i` exactly,
    /// `T = e^{iπ/4}` with `re == im == 1/√2`, and `|P(k)|² = 1` to within one
    /// unit in circuit2 and three in circuit/1.
    #[test]
    fn table_identities_per_gateset() {
        for gs in Gateset::ALL {
            assert_eq!(gs.phase(1), Amp { re: -ONE, im: 0 }, "{gs:?} P(1)");
            assert_eq!(gs.phase(2), Amp { re: 0, im: ONE }, "{gs:?} P(2)");
            assert_eq!(gs.phase(3).re, gs.phase(3).im, "{gs:?} T");
            assert_eq!(gs.phase(3).re, gs.inv_sqrt2(), "{gs:?} 1/√2");
            assert_eq!(gs.table()[0], Amp::ZERO, "{gs:?} index 0");
            let units: i128 = match gs {
                Gateset::V1 => 3,
                Gateset::V2 => 1,
            };
            for k in 1..=DYADIC_MAX {
                // |P|² − 1 at scale 2^60, exactly.
                let d = gs.phase(k).norm2() - (ONE as i128) * (ONE as i128);
                assert!(
                    d.abs() <= units << FRAC,
                    "{gs:?}: |P({k})|² − 1 = {} units",
                    d as f64 / ONE as f64
                );
            }
        }
    }

    /// The identifiers each gateset names map back to it, and its contract bytes
    /// are the canonical template. circuit/1's must equal the serializer render
    /// every existing circuit hash was taken over (`serde_contract_v1`).
    #[test]
    fn gateset_identifiers_are_consistent() {
        assert_eq!(Gateset::default(), Gateset::V2);
        for gs in Gateset::ALL {
            assert_eq!(Gateset::from_magic(&gs.magic()), Some(gs));
            assert_eq!(Gateset::from_capability(gs.capability()), Some(gs));
            assert_eq!(gs.ext(), format!("{}/1", gs.capability()));
        }
        assert_eq!(Gateset::from_magic(b"WQC3"), None);
        assert_eq!(Gateset::from_magic(b"WQC"), None);
        assert_eq!(Gateset::from_magic(b"WQC1\0"), None);
        assert_eq!(Gateset::from_capability("wai.quantum.circuit/1"), None);
        assert_eq!(Gateset::from_capability("wai.quantum.circuit3"), None);
        assert_eq!(
            Gateset::V1.contract_bytes(5),
            br#"{"ext":"wai.quantum.circuit/1","n_qubits":5,"numeric":"wai.det.fixed64","frac":30,"gateset":"cliffordT+dyadicP"}"#
        );
        assert_eq!(
            Gateset::V2.contract_bytes(5),
            br#"{"ext":"wai.quantum.circuit2/1","n_qubits":5,"numeric":"wai.det.amp64","frac":30,"gateset":"cliffordT+dyadicP.rn"}"#
        );
        assert_eq!((Gateset::V1.contract_bytes(5).len(), Gateset::V2.contract_bytes(5).len()), (112, 114));
        for n in 0..=u8::MAX {
            assert_eq!(Gateset::V1.contract_bytes(n), serde_contract_v1(n), "n_qubits {n}");
        }
    }

    /// circuit/1's contract as 0.3.x rendered it, through the serializer with
    /// member order preserved. Every circuit/1 hash and file was taken over
    /// these bytes.
    fn serde_contract_v1(n_qubits: u8) -> Vec<u8> {
        serde_json::to_vec(&serde_json::json!({
            "ext": "wai.quantum.circuit/1",
            "n_qubits": n_qubits,
            "numeric": "wai.det.fixed64",
            "frac": FRAC,
            "gateset": "cliffordT+dyadicP",
        }))
        .unwrap()
    }

    /// The measurement section as 0.3.x rendered it, through the serializer.
    fn serde_measure(m: Measure) -> Vec<u8> {
        serde_json::to_vec(&serde_json::json!({ "seed": m.seed, "shots": m.shots, "basis": "computational" })).unwrap()
    }

    #[test]
    fn simulation_is_deterministic_byte_for_byte() {
        // The core claim: same circuit → identical amplitudes, twice.
        let c = Circuit::qft(6);
        let a = c.simulate().unwrap();
        let b = c.simulate().unwrap();
        assert_eq!(a.amps, b.amps, "quantum sim must be byte-exact reproducible");
        assert_eq!(a.statevector_hash(), b.statevector_hash());
    }

    #[test]
    fn bell_state() {
        // H(0); CX(0,1) → (|00⟩ + |11⟩)/√2
        let mut c = Circuit::new(2);
        c.h(0).cx(0, 1);
        let sv = c.simulate().unwrap();
        let s = c.gateset.inv_sqrt2();
        assert!(approx(sv.amps[0b00].re, s, TOL));
        assert_eq!(sv.amps[0b01], Amp::ZERO);
        assert_eq!(sv.amps[0b10], Amp::ZERO);
        assert!(approx(sv.amps[0b11].re, s, TOL));
        // fidelity with itself ≈ 1
        assert!(approx(sv.fidelity_fx(&sv), ONE, TOL));
    }

    #[test]
    fn ghz_state() {
        // H(0); CX(0,1); CX(1,2) → (|000⟩ + |111⟩)/√2
        let mut c = Circuit::new(3);
        c.h(0).cx(0, 1).cx(1, 2);
        let sv = c.simulate().unwrap();
        let s = c.gateset.inv_sqrt2();
        assert!(approx(sv.amps[0b000].re, s, TOL));
        assert!(approx(sv.amps[0b111].re, s, TOL));
        for k in 1..7 {
            assert_eq!(sv.amps[k], Amp::ZERO, "index {k} should be empty");
        }
    }

    #[test]
    fn qft_of_zero_is_uniform_superposition() {
        // QFT|0…0⟩ = (1/√N) Σ|k⟩ — every amplitude has magnitude 1/√N.
        let n = 4u8;
        let dim = 1usize << n;
        let sv = Circuit::qft(n).simulate().unwrap();
        // |amp|^2 = 1/N; norm2() is at 2^(2·FRAC), so (norm2 >> FRAC) is at 2^FRAC,
        // i.e. the expected value is ONE/N.
        let want2 = ONE / dim as i64;
        for k in 0..dim {
            let mag2 = ((sv.amps[k].norm2()) >> FRAC) as i64;
            assert!(approx(mag2, want2, ONE / 1000), "amp {k} magnitude^2 {mag2} != {want2}");
        }
    }

    #[test]
    fn orthogonal_states_have_zero_fidelity() {
        let a = StateVector::basis(3, 0b000);
        let b = StateVector::basis(3, 0b111);
        assert_eq!(a.fidelity_fx(&b), 0);
        assert!(approx(a.fidelity_fx(&a), ONE, 1));
    }

    #[test]
    fn shot_histogram_is_deterministic() {
        // Bell state: shots concentrate on |00⟩ and |11⟩, reproducibly.
        let mut c = Circuit::new(2);
        c.h(0).cx(0, 1);
        let sv = c.simulate().unwrap();
        let h1 = sv.histogram_hash(0xC0FFEE, 10_000);
        let h2 = sv.histogram_hash(0xC0FFEE, 10_000);
        assert_eq!(h1, h2, "pinned-seed histogram must be reproducible");
        let counts = sv.sample_shots(0xC0FFEE, 10_000);
        assert_eq!(counts[0b01], 0);
        assert_eq!(counts[0b10], 0);
        assert!(counts[0b00] > 4000 && counts[0b00] < 6000, "≈ half on |00⟩: {}", counts[0b00]);
        assert_eq!(counts[0b00] + counts[0b11], 10_000);
    }

    /// The four circuits of the circuit/1 corpus, built in `gs`: bell, with
    /// the corpus measurement, then ghz3, qft3 and qft5.
    fn corpus_circuits(gs: Gateset) -> [(&'static str, Circuit, Option<Measure>); 4] {
        let mut bell = Circuit::with_gateset(2, gs);
        bell.h(0).cx(0, 1);
        let mut ghz3 = Circuit::with_gateset(3, gs);
        ghz3.h(0).cx(0, 1).cx(1, 2);
        let qft = |n| {
            let mut c = Circuit::qft(n);
            c.gateset = gs;
            c
        };
        [
            ("bell", bell, Some(Measure { seed: 12_648_430, shots: 100_000 })),
            ("ghz3", ghz3, None),
            ("qft3", qft(3), None),
            ("qft5", qft(5), None),
        ]
    }

    /// A circuit/1 container is written byte for byte as 0.3.x wrote it: the
    /// four corpus files reproduce, and so do their manifest circuit hashes.
    /// The measurement, now rendered directly, matches the serializer's render.
    #[test]
    fn circuit1_files_are_written_as_before() {
        // (file, bytes, BLAKE3 of the file, manifest circuit_hash)
        let want = [
            ("bell", 216, "97b294a85afe628060a51f12e7f9c2104ebc0b7d04ef36ff15b154ed1c54728e",
             "894d9a9582b19f27b8d08d06938357237fc00f7f30bf4bfc2be0c04fa5e47a4e"),
            ("ghz3", 157, "8f2d924416a19bc96407f36669ca3f55bcd6a57b3256bd16a669dc96bd3fb490",
             "34e44f7d68037cbb3aee677f56b894dc48df6d408aa068a16d3b18076e5e7c5c"),
            ("qft3", 191, "e1a7e49e2804130d802d8e927f63d6633b966c8e6b7cd3eb7278e01fd1f884c5",
             "7e8414be4410557a05b52bcc1084e1c5f93b6caa0f28fa45f4d2ab661efd2a60"),
            ("qft5", 261, "92b4e1497f52109922a0f2ea7ec8de2b1d9e92f155dafbf6e5f69eca4427878c",
             "f36a9af99a251cae6f42144b2be3c57523f15e1e2ab3492cf15cac709fc048be"),
        ];
        for ((name, c, m), (file, len, digest, hash)) in corpus_circuits(Gateset::V1).into_iter().zip(want) {
            assert_eq!(name, file);
            let bytes = to_wqc(&c, m).unwrap();
            assert_eq!(bytes.len(), len, "{name}.wqc");
            assert_eq!(hx(blake3::hash(&bytes).as_bytes()), digest, "{name}.wqc");
            assert_eq!(hx(&circuit_hash(&c).unwrap()), hash, "{name}");
        }
        let values = [0, 1, 42, 1024, 100_000, 12_648_430, u64::from(u32::MAX), u64::MAX];
        for seed in values {
            for shots in values {
                let m = Measure { seed, shots };
                assert_eq!(measure_bytes(m), serde_measure(m), "{m:?}");
            }
        }
    }

    /// Every hash of both conformance corpora, reproduced from code, so the
    /// unit suite alone catches a drift in either gateset. The circuit/1 rows
    /// for bell, ghz3, qft3 and qft5 are the manifest of the files 0.3.x
    /// wrote, nine hashes in all; the other rows were computed by an
    /// integer-only implementation that shares no code with this crate. None
    /// of these values may be edited to make this test pass.
    #[test]
    fn the_conformance_corpora_are_pinned() {
        // (file, bytes, n_ops, circuit_hash, statevector_hash)
        type Row = (&'static str, usize, usize, &'static str, &'static str);
        let v1: [Row; 7] = [
            ("bell", 216, 2, "894d9a9582b19f27b8d08d06938357237fc00f7f30bf4bfc2be0c04fa5e47a4e",
             "bd5df69662f4e61878d17baee1f3cf5705c3726d6740d1a4bf7be5133c099efb"),
            ("ghz3", 157, 3, "34e44f7d68037cbb3aee677f56b894dc48df6d408aa068a16d3b18076e5e7c5c",
             "08c5e6fb50b91a00d28795506cc5d99ce1e9a3a8b4c41e4bfd05f5d829af9542"),
            ("qft3", 191, 9, "7e8414be4410557a05b52bcc1084e1c5f93b6caa0f28fa45f4d2ab661efd2a60",
             "0fe1e5b8e2ac8e9afbfa926b3877967ebfead92efb7ed8e1ff6b117668926d6e"),
            ("qft5", 261, 21, "f36a9af99a251cae6f42144b2be3c57523f15e1e2ab3492cf15cac709fc048be",
             "c268d77f2ddf84fecb028e80ce989123912870f04e5b9a1e3d885c9ce235c27d"),
            ("qft5_x22", 276, 24, "40c847fe9f85b6d79d245b6e65420197e6bbedb9c0b1ab77baaf1cfe077ff4f5",
             "95a2d2fd41f6a13c19562729a4b553ee12db6adfaa361b5061c16d801354d0cd"),
            ("dyadic_ladder", 564, 72, "2fa5b02c952dcaf55ab4cfd2ea5b3fae999ac912a32a4e06dd2b6ee96bb58d4f",
             "9754af24edd026601297a0f3334eb7171174978841a2938c9b372db942a3851d"),
            ("qft8_ones", 460, 56, "5ca94d3a960b8ec65a1eb802699f1ed8e20f00e7fbd9297849f0c60047404a1a",
             "fd683498744552e9e2e8657a2bdc244e4f03d4ee4dd292c9705ba38b19fb709c"),
        ];
        let v2: [Row; 7] = [
            ("bell", 218, 2, "6f6199943059fc4d90554d7c33f1ac69723603f5128962a6e3cac9c765484d2b",
             "db68a50b08a8c2eabc89a881db9fcb130fff89dc36182d4c7f90da5bf4720b15"),
            ("ghz3", 159, 3, "c5911b898a899c999ce01730bb825c84678569493ac513760e15dcbccd84227f",
             "51850f28817200e9e03980b60f34879c7bd0147a07b0c1746528b801cae2e2f8"),
            ("qft3", 193, 9, "c23a1b1e31b87705999cbc54866f4c9e3e4c4d2b274d6260672121abb5ca89b8",
             "70fe25ec22008f86076328b2c0908ede778e6d8227e563e8c7dcd89f7e844bf2"),
            ("qft5", 263, 21, "8559c83769ac8c0cb3aabc64672233aa569cd4ce92aa6545f1449982b9aa8b8c",
             "50b01fab0bfee966e553a947aea6aa09cebf8f21ba201cd548d703048d66f2ba"),
            ("qft5_x22", 278, 24, "c126604030b16f4f1e88e6683894f9316262f49b3419be06ab74c7dbd2441651",
             "5def2fca1d29e9e2b9ca2020cd2fa2dc056d07ce2e21b6e101b6116a03efd61c"),
            ("dyadic_ladder", 566, 72, "8a6dcf0d41ecefd1a9761fef49b385d7dffe7543a9382efda02bfe4f7a0896fe",
             "28bdb7265c10d798012ac78b604327700acd8cdac0e1f88ed4045f8b0b78bc7c"),
            ("qft8_ones", 462, 56, "286f044750b19937f6e7bf6f534053aa74a21777516c773dbd7c254df7301b34",
             "e059d43a574780d819a872f7904bf96918ba5880e1c7c32b38159306acd559fe"),
        ];
        for (gs, rows) in [(Gateset::V1, v1), (Gateset::V2, v2)] {
            // The four originals, then the three table-sensitive vectors.
            let x_then_qft = |n: u8, xs: &[u8]| {
                let mut c = Circuit::with_gateset(n, gs);
                for &q in xs {
                    c.x(q);
                }
                c.ops.extend(Circuit::qft(n).ops);
                c
            };
            let mut vectors: Vec<(&str, Circuit, Option<Measure>)> = corpus_circuits(gs).into();
            vectors.push(("qft5_x22", x_then_qft(5, &[1, 2, 4]), None));
            vectors.push(("dyadic_ladder", dyadic_ladder(gs), None));
            vectors.push(("qft8_ones", x_then_qft(8, &[0, 1, 2, 3, 4, 5, 6, 7]), None));
            for ((name, c, m), (file, len, n_ops, circuit, statevector)) in vectors.into_iter().zip(rows) {
                assert_eq!(name, file);
                assert_eq!(c.gateset, gs, "{name}");
                assert_eq!(to_wqc(&c, m).unwrap().len(), len, "{gs:?} {name}.wqc");
                assert_eq!(c.ops.len(), n_ops, "{gs:?} {name}");
                assert_eq!(hx(&circuit_hash(&c).unwrap()), circuit, "{gs:?} {name}");
                let sv = c.simulate().unwrap();
                assert_eq!(hx(&sv.statevector_hash()), statevector, "{gs:?} {name}");
                // Bell is the measured vector. Its sampled outcomes do not
                // move between the tables, so neither does its histogram.
                if let Some(Measure { seed, shots }) = m {
                    assert_eq!(name, "bell");
                    assert_eq!(
                        hx(&sv.histogram_hash(seed, shots)),
                        "0893f48ef32ca9f2865f6617537f0a4a5334f67e3c944bba845beb2f377acdb3",
                        "{gs:?} {name}"
                    );
                }
            }
        }
    }

    /// Both containers round-trip: the magic names the gateset, and the
    /// decoded circuit equals the written one, gateset included, with the same
    /// hash and state. circuit2's corpus circuits have the plan's sizes and
    /// circuit hashes. circuit/1 keeps reading what it always read: its
    /// contract in any member order, whitespace or with extra members, hashed
    /// as the canonical one; bytes after the op-log; a measurement that is
    /// spaced or lacks `shots`. Unknown sections are ignored in both.
    #[test]
    fn wqc_round_trips() {
        for gs in Gateset::ALL {
            for (name, c, m) in corpus_circuits(gs) {
                for m in [None, m.or(Some(Measure { seed: 42, shots: 1024 }))] {
                    let bytes = to_wqc(&c, m).unwrap();
                    assert_eq!(bytes[0..4], gs.magic(), "{gs:?} {name}");
                    let (c2, m2) = from_wqc(&bytes).unwrap();
                    assert_eq!((&c2, m2), (&c, m), "{gs:?} {name}");
                    assert_eq!(circuit_hash(&c2).unwrap(), circuit_hash(&c).unwrap(), "{gs:?} {name}");
                    assert_eq!(c2.simulate().unwrap(), c.simulate().unwrap(), "{gs:?} {name}");
                }
            }
        }
        // circuit2's four: (file, bytes, circuit_hash). QFT-5 is 263 bytes
        // here and 261 in circuit/1.
        let want = [
            ("bell", 218, "6f6199943059fc4d90554d7c33f1ac69723603f5128962a6e3cac9c765484d2b"),
            ("ghz3", 159, "c5911b898a899c999ce01730bb825c84678569493ac513760e15dcbccd84227f"),
            ("qft3", 193, "c23a1b1e31b87705999cbc54866f4c9e3e4c4d2b274d6260672121abb5ca89b8"),
            ("qft5", 263, "8559c83769ac8c0cb3aabc64672233aa569cd4ce92aa6545f1449982b9aa8b8c"),
        ];
        for ((name, c, m), (file, len, hash)) in corpus_circuits(Gateset::V2).into_iter().zip(want) {
            assert_eq!(name, file);
            assert_eq!(to_wqc(&c, m).unwrap().len(), len, "circuit2 {name}.wqc");
            assert_eq!(hx(&circuit_hash(&c).unwrap()), hash, "circuit2 {name}");
        }

        // The accept corpus reads as it says, and its Bell circuit is the
        // vectors' bell.
        let [(_, bell, _), ..] = corpus_circuits(Gateset::V1);
        let accepted = accept_corpus();
        assert_eq!(accepted[0].2, bell);
        let names: std::collections::BTreeSet<_> = accepted.iter().map(|c| c.0).collect();
        assert_eq!(names.len(), accepted.len(), "case names repeat");
        for (name, file, circuit, measure) in accepted {
            let (back, m) = from_wqc(&file).unwrap_or_else(|e| panic!("{name}: {e}"));
            assert_eq!((&back, m), (&circuit, measure), "{name}");
        }
    }

    /// Every case of the reject corpus fails with its class; so do the
    /// structural edge cases; the rules apply in the order the spec gives.
    /// More than 255 controls is refused everywhere a circuit is checked, and
    /// the widest circuit2 gate round-trips. An op-log that claims `u32::MAX`
    /// ops reserves for the ops its bytes can hold, not for the claim, and no
    /// section range can wrap.
    #[test]
    fn container_rejects() {
        let bell = oplog_of(&[op(4, &[], 0, 0), op(1, &[0], 1, 0)]);
        let (canon1, canon2) = (Gateset::V1.contract_bytes(2), Gateset::V2.contract_bytes(2));
        let measure = measure_bytes(Measure { seed: 1, shots: 10 });
        // The bases the cases are built from are accepted.
        for file in [
            wqc_raw(b"WQC1", &[(SECT_CONTRACT, &canon1), (SECT_OPLOG, &bell)]),
            wqc_raw(b"WQC2", &[(SECT_CONTRACT, &canon2), (SECT_OPLOG, &bell)]),
            wqc_raw(b"WQC2", &[(SECT_CONTRACT, &canon2), (SECT_OPLOG, &bell), (SECT_MEASURE, &measure)]),
        ] {
            from_wqc(&file).unwrap();
        }

        let cases = reject_corpus();
        let names: std::collections::BTreeSet<_> = cases.iter().map(|c| c.0).collect();
        assert_eq!(names.len(), cases.len(), "case names repeat");
        let kinds: std::collections::BTreeSet<_> = cases.iter().map(|c| c.2).collect();
        assert_eq!(kinds, std::collections::BTreeSet::from(["contract", "invalid", "malformed"]));
        for (name, file, kind) in &cases {
            match from_wqc(file) {
                Ok(c) => panic!("{name} was accepted: {c:?}"),
                Err(e) => assert_eq!(e.kind(), *kind, "{name}: {e}"),
            }
        }
        // The I-JSON cases are refused by the rule their names give.
        for (name, rule) in [
            ("v1_contract_duplicate_member", "duplicate member"),
            ("v1_contract_lone_surrogate", "surrogate"),
            ("v1_contract_lone_low_surrogate", "surrogate"),
            ("v1_contract_too_deep", "depth"),
            ("v1_contract_not_utf8", "utf-8"),
            ("v1_contract_number_out_of_range", "number range"),
            ("v1_contract_integer_rounds_to_infinity", "number range"),
            ("v1_contract_number_rounds_to_infinity", "number range"),
            ("v1_contract_noncharacter", "noncharacter"),
            ("v1_contract_raw_noncharacter", "noncharacter"),
            ("v1_contract_bom", "syntax"),
            ("v1_contract_nan", "syntax"),
            ("v1_contract_control_character", "syntax"),
            ("v2_n_qubits_rounds_to_infinity", "number range"),
            ("v1_contract_long_literal", "number range"),
            ("v1_measure_duplicate_member", "duplicate member"),
            ("v1_measure_not_json", "syntax"),
        ] {
            let (_, file, _) = cases.iter().find(|c| c.0 == name).unwrap_or_else(|| panic!("no case {name}"));
            let e = from_wqc(file).unwrap_err().to_string();
            assert!(e.contains(&format!("I-JSON: {rule} at byte")), "{name}: {e}");
        }

        // Structure, beyond the corpus.
        let structural: [(&str, Vec<u8>, &str); 9] = [
            ("short", b"WQC1\x00".to_vec(), "malformed"),
            ("table_cut", wqc_raw(b"WQC1", &[(SECT_CONTRACT, &canon1), (SECT_OPLOG, &bell)])[..20].to_vec(), "malformed"),
            ("no_contract", wqc_raw(b"WQC2", &[(SECT_OPLOG, &bell)]), "malformed"),
            ("no_oplog", wqc_raw(b"WQC2", &[(SECT_CONTRACT, &canon2)]), "malformed"),
            ("past_end", wqc_table(b"WQC1", &[(SECT_CONTRACT, 0, 112), (SECT_OPLOG, 112, 16)], &[canon1.as_slice(), &bell].concat()), "malformed"),
            ("contract_not_json", wqc_raw(b"WQC1", &[(SECT_CONTRACT, b"{\"ext\":"), (SECT_OPLOG, &bell)]), "contract"),
            (
                "v2_n_qubits_256",
                wqc_raw(b"WQC2", &[(SECT_CONTRACT, &swap(&canon2, "\"n_qubits\":2", "\"n_qubits\":256")), (SECT_OPLOG, &bell)]),
                "malformed",
            ),
            (
                "v1_n_qubits_negative",
                wqc_raw(b"WQC1", &[(SECT_CONTRACT, &swap(&canon1, "\"n_qubits\":2", "\"n_qubits\":-2")), (SECT_OPLOG, &bell)]),
                "malformed",
            ),
            (
                "v1_measure_not_json",
                wqc_raw(b"WQC1", &[(SECT_CONTRACT, &canon1), (SECT_OPLOG, &bell), (SECT_MEASURE, b"seed=1")]),
                "malformed",
            ),
        ];
        for (name, file, kind) in &structural {
            assert_eq!(from_wqc(file).map(|_| ()).unwrap_err().kind(), *kind, "{name}");
        }
        // An error quotes at most 64 characters of what the file holds.
        let long = swap(&canon1, "wai.quantum.circuit/1", &"é".repeat(10_000));
        let e = from_wqc(&wqc_raw(b"WQC1", &[(SECT_CONTRACT, &long), (SECT_OPLOG, &bell)])).unwrap_err();
        assert_eq!(e.kind(), "contract");
        assert!(e.to_string().chars().count() < 160, "{e}");

        // The order: a repeated section before the contract, the contract
        // before the op-log, the op-log before validation, validation before
        // the measurement.
        let first = |sections: &[(u8, &[u8])]| from_wqc(&wqc_raw(b"WQC2", sections)).unwrap_err().kind();
        let canon3 = Gateset::V2.contract_bytes(3);
        let spaced = swap(&canon3, ":3", ": 3");
        let repeats = oplog_of(&[op(1, &[0, 0], 2, 0)]);
        assert_eq!(first(&[(SECT_CONTRACT, &spaced), (SECT_OPLOG, &repeats), (SECT_OPLOG, &repeats)]), "malformed");
        assert_eq!(first(&[(SECT_CONTRACT, &spaced), (SECT_OPLOG, &[0, 0])]), "contract");
        assert_eq!(first(&[(SECT_CONTRACT, &canon3), (SECT_OPLOG, &[repeats.as_slice(), &[0]].concat())]), "malformed");
        assert_eq!(first(&[(SECT_CONTRACT, &canon3), (SECT_OPLOG, &repeats), (SECT_MEASURE, b"[]")]), "invalid");

        // More than 255 controls: refused by validate, to_wqc and circuit_hash
        // alike, in both gatesets.
        for gs in Gateset::ALL {
            let mut c = Circuit::with_gateset(2, gs);
            c.ops.push(Gate { base: BaseGate::X, controls: vec![0; 256], target: 1, param: 0 });
            let e = c.validate().unwrap_err();
            assert_eq!(e.kind(), "invalid", "{gs:?}");
            assert_eq!(to_wqc(&c, None).unwrap_err(), e, "{gs:?}");
            assert_eq!(circuit_hash(&c).unwrap_err(), e, "{gs:?}");
        }
        // 254 distinct controls on 255 qubits is the widest circuit2 gate.
        let mut wide = Circuit::with_gateset(255, Gateset::V2);
        wide.ops.push(Gate { base: BaseGate::X, controls: (0..254).collect(), target: 254, param: 0 });
        let bytes = to_wqc(&wide, None).unwrap();
        let (back, _) = from_wqc(&bytes).unwrap();
        assert_eq!(back, wide);
        assert_eq!(circuit_hash(&back).unwrap(), circuit_hash(&wide).unwrap());

        // A 9-byte op-log claiming u32::MAX ops: the reservation is what 9
        // bytes can hold, one op, and the decode is malformed at the second.
        assert_eq!(oplog_capacity(u32::MAX, 9), 1);
        assert_eq!(oplog_capacity(u32::MAX, 3), 0);
        assert_eq!(oplog_capacity(3, 1 << 20), 3);
        let n_ops_u32max = [u32::MAX.to_le_bytes().as_slice(), &op(4, &[], 0, 0)].concat();
        assert_eq!(parse_oplog(&n_ops_u32max).unwrap_err(), QuantumError::Malformed("op header".into()));

        // Section ranges are checked, so an offset that would wrap a 32-bit
        // usize is out of bounds there, as it is here.
        assert_eq!(section_range(24, 0, 112, 136), Some(24..136));
        assert_eq!(section_range(24, 0, 113, 136), None);
        assert_eq!(section_range(usize::MAX - 8, 0xFFFF_FFF0, 0x20, usize::MAX), None);
        assert_eq!(section_range(usize::MAX - 0x10, 0x8, 0x20, usize::MAX), None);
        assert_eq!(section_range(usize::MAX - 0x40, 0x10, 0x20, usize::MAX), Some(usize::MAX - 0x30..usize::MAX - 0x10));
    }


    #[test]
    fn validation_rejects_bad_indices() {
        let mut c = Circuit::new(2);
        c.push(BaseGate::X, vec![], 5, 0); // target out of range
        assert!(matches!(c.validate(), Err(QuantumError::Invalid(_))));

        let mut c = Circuit::new(2);
        c.push(BaseGate::X, vec![1], 1, 0); // control == target
        assert!(matches!(c.validate(), Err(QuantumError::Invalid(_))));

        let mut c = Circuit::new(1);
        c.push(BaseGate::P, vec![], 0, 999); // P(k) out of dyadic range
        assert!(matches!(c.validate(), Err(QuantumError::Invalid(_))));

        // The op-log counts controls in a u8: 256 is too many in either
        // gateset, and circuit/1, which allows repeats, takes 255.
        for gs in Gateset::ALL {
            let mut c = Circuit::with_gateset(2, gs);
            c.push(BaseGate::X, vec![0; 256], 1, 0);
            assert!(matches!(c.validate(), Err(QuantumError::Invalid(_))), "{gs:?}");
        }
        let mut c = Circuit::with_gateset(2, Gateset::V1);
        c.push(BaseGate::X, vec![0; 255], 1, 0);
        assert_eq!(c.validate(), Ok(()));
    }

    /// A circuit carries its gateset, and a circuit compares unequal to the same
    /// ops in the other one.
    #[test]
    fn circuits_carry_their_gateset() {
        assert_eq!(Circuit::new(3).gateset, Gateset::V2);
        assert_eq!(Circuit::qft(3).gateset, Gateset::V2);
        for gs in Gateset::ALL {
            let mut c = Circuit::with_gateset(3, gs);
            assert_eq!((c.n_qubits, c.gateset), (3, gs));
            c.h(0).cx(0, 1);
            let e = c.empty_like();
            assert_eq!((e.n_qubits, e.gateset, e.ops.len()), (3, gs, 0));
        }
        let v2 = Circuit::qft(3);
        let mut v1 = v2.clone();
        v1.gateset = Gateset::V1;
        assert_ne!(v1, v2);
        assert_ne!(v1.simulate().unwrap(), v2.simulate().unwrap());
    }

    #[test]
    fn error_kinds_name_their_class() {
        let cases = [
            (QuantumError::Malformed("m".into()), "malformed", "malformed WQC: m"),
            (QuantumError::Invalid("i".into()), "invalid", "invalid circuit: i"),
            (QuantumError::Contract("c".into()), "contract", "contract mismatch: c"),
            (QuantumError::TooManyQubits(27), "too-many-qubits", "too many qubits: 27"),
        ];
        for (e, kind, shown) in cases {
            assert_eq!((e.kind(), e.to_string().as_str()), (kind, shown));
        }
    }

    /// `P(k)` on `|1⟩` lands on the table entry itself, bit for bit, in both
    /// gatesets, plain and controlled: `fxmul` by `ONE` is exact.
    #[test]
    fn phase_is_copied_exactly() {
        for gs in Gateset::ALL {
            for k in 1..=DYADIC_MAX as u16 {
                let mut c = Circuit::with_gateset(1, gs);
                c.x(0).p(k, 0);
                assert_eq!(c.simulate().unwrap().amps, [Amp::ZERO, gs.phase(k.into())], "{gs:?} P({k})");
                let mut c = Circuit::with_gateset(2, gs);
                c.x(0).x(1).cp(k, 0, 1);
                let amps = c.simulate().unwrap().amps;
                assert_eq!(amps[0b11], gs.phase(k.into()), "{gs:?} CP({k})");
                assert!(amps[..0b11].iter().all(|a| *a == Amp::ZERO), "{gs:?} CP({k})");
            }
        }
    }

    /// A small phase applied many times. `P(19)` applied 2^18 times to `|1⟩` is
    /// a rotation by π. circuit2 gets there. circuit/1's `P(19)` is
    /// `(2^30 − 1, 0)`: it never rotates, and each application loses one unit.
    #[test]
    fn repeated_small_phase() {
        let reps = 1usize << 18;
        let run = |gs| {
            let mut c = Circuit::with_gateset(1, gs);
            c.x(0);
            c.ops.extend(std::iter::repeat_n(Gate { base: BaseGate::P, controls: vec![], target: 0, param: 19 }, reps));
            let amps = c.simulate().unwrap().amps;
            assert_eq!(amps[0], Amp::ZERO);
            amps[1]
        };
        assert_eq!(run(Gateset::V1), Amp { re: 1_073_479_680, im: 0 });
        assert_eq!(run(Gateset::V1).re, ONE - reps as i64);
        let v2 = run(Gateset::V2);
        let angle = (v2.im as f64).atan2(v2.re as f64);
        assert!((angle.abs() - std::f64::consts::PI).abs() < 1e-4, "circuit2 turned {angle} rad, not π");
        // The norm drifts as slowly as the rounding allows (1.0000189 here).
        let norm = (v2.norm2() as f64).sqrt() / ONE as f64;
        assert!((norm - 1.0).abs() < 1e-4, "|amp| = {norm}");
    }

    /// circuit2 is strict where circuit/1 is not: a repeated control and a
    /// parameter on a gate that takes none are circuit/1 circuits, and are not
    /// circuit2 ones.
    #[test]
    fn v2_rejects_what_v1_tolerates() {
        let use_gate = |gs, n, g: Gate| {
            let mut c = Circuit::with_gateset(n, gs);
            c.ops.push(g);
            c
        };
        let mut lax = vec![Gate { base: BaseGate::X, controls: vec![0, 0], target: 1, param: 0 }];
        for base in [BaseGate::I, BaseGate::X, BaseGate::Y, BaseGate::Z, BaseGate::H, BaseGate::S, BaseGate::Sdg, BaseGate::T, BaseGate::Tdg] {
            lax.push(Gate { base, controls: vec![], target: 0, param: 1 });
        }
        lax.push(Gate { base: BaseGate::H, controls: vec![], target: 0, param: u16::MAX });
        for g in lax {
            let v1 = use_gate(Gateset::V1, 2, g.clone());
            assert_eq!(v1.validate(), Ok(()), "{g:?}");
            assert!(v1.simulate().is_ok(), "{g:?}");
            let v2 = use_gate(Gateset::V2, 2, g.clone());
            let e = v2.validate().unwrap_err();
            assert_eq!(e.kind(), "invalid", "{g:?}");
            assert_eq!(v2.simulate().unwrap_err(), e, "{g:?}");
        }
        // The parameter rule spares P, whose parameter is its k.
        let mut p = Circuit::new(1);
        p.p(5, 0);
        assert_eq!(p.validate(), Ok(()));
        // Distinct controls: 254 on 255 qubits is the most circuit2 can hold.
        let all: Vec<u8> = (0..254).collect();
        assert_eq!(use_gate(Gateset::V2, 255, Gate { base: BaseGate::X, controls: all.clone(), target: 254, param: 0 }).validate(), Ok(()));
        let mut repeat = all.clone();
        repeat.push(253);
        let g = Gate { base: BaseGate::X, controls: repeat, target: 254, param: 0 };
        assert_eq!(use_gate(Gateset::V1, 255, g.clone()).validate(), Ok(()));
        assert!(matches!(use_gate(Gateset::V2, 255, g).validate(), Err(QuantumError::Invalid(_))));
    }

    /// Accuracy against the ideal, by floating-point references whose phases
    /// come from the crate's own `sin_cos`, compiled under these features, all
    /// of which `full` enables.
    #[cfg(any(
        feature = "quantum_spd",
        feature = "quantum_pauli",
        feature = "quantum_mps",
        feature = "quantum_vml",
        feature = "quantum_phasor",
        feature = "quantum_kernel",
        feature = "quantum_tn",
        feature = "quantum_sv",
        feature = "quantum_tdvp",
        feature = "quantum_nqs"
    ))]
    mod against_ideal {
        use super::*;

        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)
        }

        fn cadd(a: C, b: C) -> C {
            (a.0 + b.0, a.1 + b.1)
        }

        /// `e^{2πi/2^k}`.
        fn root(k: u16) -> C {
            let (s, c) = crate::repro::sin_cos(std::f64::consts::TAU / (1u64 << k) as f64);
            (c, s)
        }

        /// The ideal matrix of a base gate.
        fn matrix(base: BaseGate, param: u16) -> [[C; 2]; 2] {
            let (z, one, i, neg_i) = ((0.0, 0.0), (1.0, 0.0), (0.0, 1.0), (0.0, -1.0));
            let h = std::f64::consts::FRAC_1_SQRT_2;
            let t = root(3);
            match base {
                BaseGate::I => [[one, z], [z, one]],
                BaseGate::X => [[z, one], [one, z]],
                BaseGate::Y => [[z, neg_i], [i, z]],
                BaseGate::Z => [[one, z], [z, (-1.0, 0.0)]],
                BaseGate::H => [[(h, 0.0), (h, 0.0)], [(h, 0.0), (-h, 0.0)]],
                BaseGate::S => [[one, z], [z, i]],
                BaseGate::Sdg => [[one, z], [z, neg_i]],
                BaseGate::T => [[one, z], [z, t]],
                BaseGate::Tdg => [[one, z], [z, (t.0, -t.1)]],
                BaseGate::P => [[one, z], [z, root(param)]],
            }
        }

        /// The circuit's ideal state from `|0…0⟩`, gate by gate.
        fn simulate(c: &Circuit) -> Vec<C> {
            let dim = 1usize << c.n_qubits;
            let mut a = vec![(0.0, 0.0); dim];
            a[0] = (1.0, 0.0);
            for g in &c.ops {
                let m = matrix(g.base, g.param);
                let tbit = 1usize << g.target;
                let mask = g.controls.iter().fold(0usize, |acc, &q| acc | (1usize << q));
                for i in 0..dim {
                    if i & tbit == 0 && i & mask == mask {
                        let j = i | tbit;
                        let (a0, a1) = (a[i], a[j]);
                        a[i] = cadd(cmul(m[0][0], a0), cmul(m[0][1], a1));
                        a[j] = cadd(cmul(m[1][0], a0), cmul(m[1][1], a1));
                    }
                }
            }
            a
        }

        /// QFT-14 of `|1…1⟩`, whose ladder multiplies nonzero amplitudes by every
        /// `P(2..=14)`. circuit2 stays within 3 units of the ideal amplitude
        /// (2.0 measured); circuit/1 is off by more than 8 (14.4 measured).
        #[test]
        fn qft14_accuracy() {
            let n = 14u8;
            let mut c = Circuit::new(n);
            for q in 0..n {
                c.x(q);
            }
            c.ops.extend(Circuit::qft(n).ops);
            let want = simulate(&c);
            let error = |gs| {
                let mut c = c.clone();
                c.gateset = gs;
                let one = ONE as f64;
                c.simulate()
                    .unwrap()
                    .amps
                    .iter()
                    .zip(&want)
                    .map(|(a, w)| {
                        let (dre, dim) = (a.re as f64 - w.0 * one, a.im as f64 - w.1 * one);
                        (dre * dre + dim * dim).sqrt()
                    })
                    .fold(0.0f64, f64::max)
            };
            let (v2, v1) = (error(Gateset::V2), error(Gateset::V1));
            assert!(v2 <= 3.0, "circuit2 QFT-14 is {v2} units from the ideal");
            assert!(v1 > 8.0, "circuit/1 QFT-14 is only {v1} units from the ideal");
        }

        /// One circuit2 gate on a normalised pair moves each component at most
        /// 2 units from the ideal (1.43 measured): ½ unit per table entry and ½
        /// per rounded product. `H`, `T`, `T†` and every `P(k)`, 400 seeded pairs
        /// each, through the simulator's own `apply`. circuit/1 fails the same
        /// check.
        #[test]
        fn one_v2_gate_is_within_two_units() {
            let mut gates = vec![(BaseGate::H, 0u16), (BaseGate::T, 0), (BaseGate::Tdg, 0)];
            gates.extend((1..=DYADIC_MAX as u16).map(|k| (BaseGate::P, k)));
            let worst = |gs| {
                let mut s = 0x5EED_0000_0000_0008u64;
                let mut uniform = || (splitmix64(&mut s) >> 11) as f64 / (1u64 << 53) as f64 * 2.0 - 1.0;
                let one = ONE as f64;
                let mut worst = 0.0f64;
                for &(base, param) in &gates {
                    let m = matrix(base, param);
                    for _ in 0..400 {
                        let z = [uniform(), uniform(), uniform(), uniform()];
                        let norm = z.iter().map(|x| x * x).sum::<f64>().sqrt();
                        let fixed = |x: f64| (x / norm * one) as i64;
                        let a = [Amp { re: fixed(z[0]), im: fixed(z[1]) }, Amp { re: fixed(z[2]), im: fixed(z[3]) }];
                        let mut got = a;
                        apply(&mut got, 1, &Gate { base, controls: vec![], target: 0, param }, gs);
                        let af = a.map(|x| (x.re as f64 / one, x.im as f64 / one));
                        for r in 0..2 {
                            let w = cadd(cmul(m[r][0], af[0]), cmul(m[r][1], af[1]));
                            worst = worst.max((got[r].re as f64 - w.0 * one).abs()).max((got[r].im as f64 - w.1 * one).abs());
                        }
                    }
                }
                worst
            };
            assert_eq!(gates.len() * 400, 14_000);
            let v2 = worst(Gateset::V2);
            assert!(v2 <= 2.0, "a circuit2 gate is {v2} units from the ideal");
            let v1 = worst(Gateset::V1);
            assert!(v1 > 2.0, "the check passed circuit/1: {v1}");
        }
    }
}