wai-quantum 0.3.26

//! Deterministic quantum-circuit **simulation** as transportable, receiptable
//! content — `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 (the `wai.det.fixed64`
//! discipline applied to amplitudes): 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}` gives
//! the controlled-phase ladder the **QFT** ("quantum FFT") needs, transported as a
//! circuit and reconstructed byte-exact. 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 std::sync::OnceLock;

/// 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 — the QFT reaches this
/// many qubits. The table is built by an integer half-angle recurrence (no float).
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";

/// The capability string this module reconstructs.
pub const CAPABILITY: &str = "wai.quantum.circuit";

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

// ---------------------------------------------------------------------------
// Pinned dyadic phase table — e^{2πi / 2^k}, built by integer half-angle
// recurrence from (cos π/2, sin π/2) = (0, 1). No float ever runs.
// ---------------------------------------------------------------------------

fn dyadic_table() -> &'static [Amp] {
    static TABLE: OnceLock<Vec<Amp>> = OnceLock::new();
    TABLE.get_or_init(|| {
        let mut t = vec![Amp::ZERO; DYADIC_MAX + 1];
        // k = 1: e^{iπ}   = -1 + 0i
        t[1] = Amp { re: -ONE, im: 0 };
        // k = 2: e^{iπ/2} =  0 + 1i
        t[2] = Amp { re: 0, im: ONE };
        // k ≥ 3: half the angle each step. For θ ∈ (0, π/2], both cos and sin are
        // non-negative, so the positive integer sqrt is the right branch:
        //   cos(θ/2) = sqrt((1+cosθ)/2),  sin(θ/2) = sqrt((1−cosθ)/2)
        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),
            };
        }
        t
    })
}

/// `1/√2` in the fixed-point floor — derived from the dyadic table (`cos π/4`), so
/// there is a *single* source for the one irrational the Clifford+T set needs.
pub fn inv_sqrt2() -> i64 {
    dyadic_table()[3].re
}

/// `e^{2πi / 2^k}` in the fixed-point floor, for `1 ≤ k ≤ DYADIC_MAX`.
pub fn phase(k: usize) -> Amp {
    dyadic_table()[k]
}

// ---------------------------------------------------------------------------
// 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 {
    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.
    fn matrix(self, param: u16) -> [[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: inv_sqrt2(), im: 0 };
        let neg_s = Amp { re: -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, phase(3)]],
            BaseGate::Tdg => [[one, z], [z, phase(3).conj()]],
            BaseGate::P => [[one, z], [z, 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` and an ordered op-log. This is the transportable
/// object — compact, linear in gate count.
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct Circuit {
    pub n_qubits: u8,
    pub ops: Vec<Gate>,
}

/// Errors from parsing a `WQC` container or building an out-of-range circuit.
#[derive(Debug, PartialEq, Eq)]
pub enum QuantumError {
    /// Bad magic, truncated section, or malformed op-log.
    Malformed(String),
    /// A qubit index ≥ `n_qubits`, a control equal to the target, or a `P(k)` with
    /// `k` outside `1..=DYADIC_MAX`.
    Invalid(String),
    /// `n_qubits` exceeds what this build will materialize (`2^n` amplitudes).
    TooManyQubits(u8),
}

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::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 {
    pub fn new(n_qubits: u8) -> Self {
        Circuit { n_qubits, ops: Vec::new() }
    }

    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. Uses the pinned dyadic phases.
    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
    }

    /// Validate qubit indices, control/target distinctness, and `P(k)` range.
    pub fn validate(&self) -> Result<(), QuantumError> {
        let n = self.n_qubits;
        for g in &self.ops {
            if g.target >= n {
                return Err(QuantumError::Invalid(format!("target {} ≥ n_qubits {}", g.target, n)));
            }
            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 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}")));
                }
            }
        }
        Ok(())
    }

    /// Classically simulate the circuit from `|0…0⟩`, 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);
        }
        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.
fn apply(amps: &mut [Amp], n: u8, g: &Gate) {
    if matches!(g.base, BaseGate::I) {
        return;
    }
    let m = g.base.matrix(g.param);
    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"): 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;

fn contract_json(n_qubits: u8) -> Vec<u8> {
    // canonical, fixed field order (serde_json preserve_order is on for the crate)
    let v = serde_json::json!({
        "ext": "wai.quantum.circuit/1",
        "n_qubits": n_qubits,
        "numeric": "wai.det.fixed64",
        "frac": FRAC,
        "gateset": "cliffordT+dyadicP",
    });
    serde_json::to_vec(&v).expect("contract json")
}

fn oplog_bytes(c: &Circuit) -> Vec<u8> {
    // u32 n_ops || per op: u8 opcode | u8 n_ctrl | u8 target | u16 param | ctrls...
    let mut out = Vec::new();
    out.extend_from_slice(&(c.ops.len() as u32).to_le_bytes());
    for g in &c.ops {
        out.push(g.base.opcode());
        out.push(g.controls.len() as u8);
        out.push(g.target);
        out.extend_from_slice(&g.param.to_le_bytes());
        out.extend_from_slice(&g.controls);
    }
    out
}

fn parse_oplog(n_qubits: u8, b: &[u8]) -> Result<Vec<Gate>, 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]]) as usize;
    let mut pos = 4;
    let mut ops = Vec::with_capacity(n_ops);
    for _ in 0..n_ops {
        if pos + 5 > b.len() {
            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 = b[pos + 1] as usize;
        let target = b[pos + 2];
        let param = u16::from_le_bytes([b[pos + 3], b[pos + 4]]);
        pos += 5;
        if pos + n_ctrl > b.len() {
            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 });
    }
    let _ = n_qubits;
    Ok(ops)
}

/// Serialize a circuit (and optional measurement) to a `WQC1` container.
pub fn to_wqc(c: &Circuit, measure: Option<Measure>) -> Vec<u8> {
    let contract = contract_json(c.n_qubits);
    let oplog = oplog_bytes(c);
    let measure_bytes = measure.map(|m| {
        let v = serde_json::json!({ "seed": m.seed, "shots": m.shots, "basis": "computational" });
        serde_json::to_vec(&v).expect("measure json")
    });

    let mut sections: Vec<(u8, Vec<u8>)> = vec![(SECT_CONTRACT, contract), (SECT_OPLOG, oplog)];
    if let Some(mb) = measure_bytes {
        sections.push((SECT_MEASURE, mb));
    }

    let mut out = Vec::new();
    out.extend_from_slice(b"WQC1");
    out.extend_from_slice(&(sections.len() as u16).to_le_bytes());
    // section table: u8 kind, u32 off (relative to start of the sections blob), u32 len
    let mut off: u32 = 0;
    for (kind, data) in &sections {
        out.push(*kind);
        out.extend_from_slice(&off.to_le_bytes());
        out.extend_from_slice(&(data.len() as u32).to_le_bytes());
        off += data.len() as u32;
    }
    for (_, data) in &sections {
        out.extend_from_slice(data);
    }
    out
}

/// Parse a `WQC1` container back into a circuit and optional measurement.
pub fn from_wqc(bytes: &[u8]) -> Result<(Circuit, Option<Measure>), QuantumError> {
    if bytes.len() < 6 || &bytes[0..4] != b"WQC1" {
        return Err(QuantumError::Malformed("magic".into()));
    }
    let n_sections = u16::from_le_bytes([bytes[4], bytes[5]]) as usize;
    let table_start = 6;
    let table_len = n_sections * 9;
    if bytes.len() < table_start + table_len {
        return Err(QuantumError::Malformed("section table".into()));
    }
    let blob_start = table_start + table_len;
    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 {
        let p = table_start + s * 9;
        let kind = bytes[p];
        let off = u32::from_le_bytes([bytes[p + 1], bytes[p + 2], bytes[p + 3], bytes[p + 4]]) as usize;
        let len = u32::from_le_bytes([bytes[p + 5], bytes[p + 6], bytes[p + 7], bytes[p + 8]]) as usize;
        let start = blob_start + off;
        if start + len > bytes.len() {
            return Err(QuantumError::Malformed("section bounds".into()));
        }
        let data = &bytes[start..start + len];
        match kind {
            SECT_CONTRACT => contract = Some(data),
            SECT_OPLOG => oplog = Some(data),
            SECT_MEASURE => measure_raw = Some(data),
            _ => {} // unknown section: ignore (forward-compatible)
        }
    }
    let contract = contract.ok_or_else(|| QuantumError::Malformed("missing contract".into()))?;
    let oplog = oplog.ok_or_else(|| QuantumError::Malformed("missing oplog".into()))?;

    let cv: serde_json::Value =
        serde_json::from_slice(contract).map_err(|e| QuantumError::Malformed(format!("contract json: {e}")))?;
    let n_qubits = cv["n_qubits"].as_u64().ok_or_else(|| QuantumError::Malformed("n_qubits".into()))? as u8;

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

    let measure = match measure_raw {
        Some(mb) => {
            let mv: serde_json::Value =
                serde_json::from_slice(mb).map_err(|e| QuantumError::Malformed(format!("measure json: {e}")))?;
            Some(Measure {
                seed: mv["seed"].as_u64().unwrap_or(0),
                shots: mv["shots"].as_u64().unwrap_or(0),
            })
        }
        None => None,
    };
    Ok((circuit, measure))
}

/// **circuit-equivalence** — the portable BLAKE3 identity of the circuit itself
/// (contract + op-log). The receipt binds this alongside the statevector hash.
pub fn circuit_hash(c: &Circuit) -> [u8; 32] {
    let mut h = blake3::Hasher::new();
    h.update(DOMAIN_CIRCUIT);
    h.update(&contract_json(c.n_qubits));
    h.update(&oplog_bytes(c));
    *h.finalize().as_bytes()
}

#[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) {
        let m = g.base.matrix(g.param);
        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); }
                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); }
                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.
    #[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 _ in 0..300 {
            let n = 3 + rnd(3) as u8;
            let dim = 1usize << n;
            let mut c = Circuit::new(n);
            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(5) as u16, t); }
                    7 => { c.cx(other, t); }
                    8 => { c.cz(other, t); }
                    9 => { c.cp(1 + rnd(4) 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);
            }
            let got = c.simulate().unwrap();
            assert_eq!(got.amps, want, "optimised apply diverged on {} gates, n={n}", 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
    }

    #[test]
    fn dyadic_table_is_unit_modulus() {
        // cos^2 + sin^2 ≈ 1 for every tabulated phase, and P(1)=Z, P(2)=S, P(3)=T.
        for k in 1..=DYADIC_MAX {
            let a = phase(k);
            let n2 = ((a.norm2()) >> FRAC) as i64; // at 2^FRAC
            assert!(approx(n2, ONE, ONE / 1000), "|e^(2πi/2^{k})|^2 = {n2} not ≈ ONE");
        }
        assert_eq!(phase(1), Amp { re: -ONE, im: 0 }); // e^{iπ}
        assert_eq!(phase(2), Amp { re: 0, im: ONE }); // e^{iπ/2}
        // T = e^{iπ/4}: re == im == 1/√2
        assert_eq!(phase(3).re, inv_sqrt2());
        assert_eq!(phase(3).im, inv_sqrt2());
    }

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

    #[test]
    fn wqc_round_trips() {
        let c = Circuit::qft(5);
        let m = Some(Measure { seed: 42, shots: 1024 });
        let bytes = to_wqc(&c, m);
        assert_eq!(&bytes[0..4], b"WQC1");
        let (c2, m2) = from_wqc(&bytes).unwrap();
        assert_eq!(c, c2, "circuit must round-trip through WQC");
        assert_eq!(m, m2);
        // and the reconstructed circuit hashes + simulates identically
        assert_eq!(circuit_hash(&c), circuit_hash(&c2));
        assert_eq!(c.simulate().unwrap().amps, c2.simulate().unwrap().amps);
    }


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