wai_quantum/

quantum_toolchain.rs

//! WAI Quantum Toolchain — the compute core of a classical-hardware quantum IDE
//! (extensions/quantum-sim § Toolchain, `wai.quantum.toolchain`).
//!
//! A framework-neutral toolchain over the byte-exact fixed-point simulator
//! ([`crate::quantum`]) that runs entirely on the hardware you already have, with
//! no cloud service or account — and whose every run is reproducible to the bit
//! and sealable into a signed receipt.
//!
//! Four layers sit on the simulator:
//! - **Algorithm library** — the canonical algorithms (Bell, GHZ, QFT, Grover,
//!   Bernstein–Vazirani, Deutsch–Jozsa, teleportation, superdense coding) as
//!   circuit generators you can run and inspect, not read about.
//! - **State analysis** — the readouts an IDE shows about a reconstructed state:
//!   probabilities, Pauli expectations, per-qubit Bloch vectors, single-qubit
//!   entanglement entropy and purity. (Display-layer `f64`; the byte-exact artifact
//!   remains the statevector hash.)
//! - **Backend selection** — detect whether a circuit is Clifford and recommend the
//!   cheapest sufficient classical method (stabilizer tableau → dense statevector →
//!   tensor network), the "classical referee" cascade the analysis calls for.
//! - **OpenQASM interop** — import/export the lingua franca, so circuits move in and
//!   out of any other toolchain without a rewrite.
//!
//! Honest boundary: the simulator's gate set is the **dyadic** family (Clifford +
//! T + `P(k)=diag(1,e^{2πi/2^k})` with arbitrary controls) — every phase exactly
//! representable, which is *why* the result is byte-identical on every machine.
//! Continuous rotations (`RX(θ)` for irrational θ) are out of the exact set by
//! construction, not omission.
//!
//! Gatesets: the algorithms and the OpenQASM importer build circuits from
//! nothing, so in [`crate::quantum::Gateset::V2`], `wai.quantum.circuit2`. A
//! circuit derived from another (a noisy trajectory, a routed circuit) keeps
//! the other's gateset.

use crate::quantum::{circuit_hash, Amp, BaseGate, Circuit, Gate, QuantumError, StateVector, FRAC};

// ===========================================================================
// Algorithm library
// ===========================================================================

/// A multi-controlled-Z on all `n` qubits (controls `0..n-1`, target `n-1`) —
/// the phase-flip primitive Grover's oracle and diffusion are built from.
fn mcz(c: &mut Circuit, n: u8) {
    if n == 1 {
        c.z(0);
        return;
    }
    c.ops.push(Gate { base: BaseGate::Z, controls: (0..n - 1).collect(), target: n - 1, param: 0 });
}

/// `|Φ⁺⟩ = (|00⟩+|11⟩)/√2` — the canonical entangled pair.
pub fn bell() -> Circuit {
    let mut c = Circuit::new(2);
    c.h(0);
    c.cx(0, 1);
    c
}

/// The `n`-qubit GHZ state `(|0…0⟩+|1…1⟩)/√2` (a linear CX chain).
pub fn ghz(n: u8) -> Circuit {
    let mut c = Circuit::new(n.max(1));
    c.h(0);
    for i in 1..n {
        c.cx(i - 1, i);
    }
    c
}

/// **Grover search** over `n` qubits for a single marked basis state, run for the
/// optimal `⌊(π/4)√(2ⁿ)⌋` iterations — amplitude amplification made visible.
pub fn grover(n: u8, marked: u32) -> Circuit {
    let n = n.clamp(1, 12);
    let mut c = Circuit::new(n);
    for q in 0..n {
        c.h(q);
    }
    let dim = 1u64 << n;
    let iters = ((std::f64::consts::PI / 4.0) * (dim as f64).sqrt()).floor() as u32;
    for _ in 0..iters.max(1) {
        // oracle: phase-flip |marked>
        for q in 0..n {
            if (marked >> q) & 1 == 0 {
                c.x(q);
            }
        }
        mcz(&mut c, n);
        for q in 0..n {
            if (marked >> q) & 1 == 0 {
                c.x(q);
            }
        }
        // diffusion about the mean
        for q in 0..n {
            c.h(q);
        }
        for q in 0..n {
            c.x(q);
        }
        mcz(&mut c, n);
        for q in 0..n {
            c.x(q);
        }
        for q in 0..n {
            c.h(q);
        }
    }
    c
}

/// **Bernstein–Vazirani** — recover the hidden `nbits`-bit `secret` in a *single*
/// oracle query. Qubits `0..nbits` carry the answer; qubit `nbits` is the phase
/// ancilla.
pub fn bernstein_vazirani(secret: u32, nbits: u8) -> Circuit {
    let nbits = nbits.clamp(1, 20);
    let anc = nbits;
    let mut c = Circuit::new(nbits + 1);
    c.x(anc);
    for q in 0..=nbits {
        c.h(q);
    }
    for q in 0..nbits {
        if (secret >> q) & 1 == 1 {
            c.cx(q, anc);
        }
    }
    for q in 0..nbits {
        c.h(q);
    }
    c
}

/// **Deutsch–Jozsa** — one query decides whether a function is constant or
/// balanced. A balanced oracle CXs every input into the ancilla; measuring the
/// inputs gives all-zeros iff the function is constant.
pub fn deutsch_jozsa(balanced: bool, nbits: u8) -> Circuit {
    let nbits = nbits.clamp(1, 20);
    let anc = nbits;
    let mut c = Circuit::new(nbits + 1);
    c.x(anc);
    for q in 0..=nbits {
        c.h(q);
    }
    if balanced {
        for q in 0..nbits {
            c.cx(q, anc);
        }
    }
    for q in 0..nbits {
        c.h(q);
    }
    c
}

/// **Quantum teleportation** (deferred-measurement form): qubit 0's prepared state
/// is transported to qubit 2 through a shared Bell pair and unitary corrections.
/// The final reduced state of qubit 2 equals qubit 0's initial state.
pub fn teleportation() -> Circuit {
    let mut c = Circuit::new(3);
    // prepare a nontrivial state on qubit 0
    c.h(0);
    c.t(0);
    // Bell pair on 1,2
    c.h(1);
    c.cx(1, 2);
    // Bell-basis change on 0,1
    c.cx(0, 1);
    c.h(0);
    // deferred corrections
    c.cx(1, 2);
    c.cz(0, 2);
    c
}

/// **Superdense coding** — send two classical bits `bits ∈ 0..4` over one qubit of
/// a Bell pair. Decoding and measuring recovers `bits`.
pub fn superdense(bits: u8) -> Circuit {
    let mut c = Circuit::new(2);
    c.h(0);
    c.cx(0, 1);
    // encode so the decoded 2-bit index equals `bits`: Z shows up as q0 (value 1),
    // X as q1 (value 2).
    if bits & 1 != 0 {
        c.z(0);
    }
    if bits & 2 != 0 {
        c.x(0);
    }
    c.cx(0, 1);
    c.h(0);
    c
}

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

/// **Trotterized transverse-field Ising** dynamics — the many-body-dynamics circuit
/// class the tensor-network and sparse-Pauli backends were built for. Each step
/// applies a ZZ coupling on every line edge (`CX · P(k) · CX`) and a transverse X
/// field (`H · P(k) · H`), all with dyadic (byte-exact) rotation angles.
pub fn trotter_ising(n: u8, steps: u8) -> Circuit {
    let n = n.clamp(2, 12);
    let k = 3u16; // dyadic angle 2π/8 = π/4
    let mut c = Circuit::new(n);
    for _ in 0..steps.clamp(1, 20) {
        for q in 0..n - 1 {
            c.cx(q, q + 1);
            c.p(k, q + 1);
            c.cx(q, q + 1);
        }
        for q in 0..n {
            c.h(q);
            c.p(k, q);
            c.h(q);
        }
    }
    c
}

/// A random **Clifford** circuit (H / S / CX) — the stabilizer backend's home turf.
pub fn random_clifford(n: u8, depth: u16, seed: u64) -> Circuit {
    let n = n.clamp(1, 12);
    let mut st = seed.wrapping_mul(0x2545_F491).wrapping_add(1);
    let mut c = Circuit::new(n);
    for _ in 0..depth {
        match splitmix(&mut st) % 3 {
            0 => {
                c.h((splitmix(&mut st) % n as u64) as u8);
            }
            1 => {
                c.s((splitmix(&mut st) % n as u64) as u8);
            }
            _ if n >= 2 => {
                let a = (splitmix(&mut st) % n as u64) as u8;
                let mut b = (splitmix(&mut st) % n as u64) as u8;
                if b == a {
                    b = (b + 1) % n;
                }
                c.cx(a, b);
            }
            _ => {
                c.h(0);
            }
        }
    }
    c
}

/// Build a library algorithm by id. `param` carries the marked/secret/bits value.
pub fn build_algorithm(id: &str, n: u8, param: u32) -> Option<Circuit> {
    Some(match id {
        "bell" => bell(),
        "ghz" => ghz(n.clamp(2, 12)),
        "qft" => Circuit::qft(n.clamp(1, 12)),
        "grover" => grover(n.clamp(2, 10), param),
        "bv" => bernstein_vazirani(param, n.clamp(1, 10)),
        "dj" => deutsch_jozsa(param != 0, n.clamp(1, 10)),
        "teleport" => teleportation(),
        "superdense" => superdense((param & 3) as u8),
        "trotter" => trotter_ising(n.clamp(2, 10), param.clamp(1, 12) as u8),
        "rclifford" => random_clifford(n.clamp(2, 10), 40, param as u64),
        _ => return None,
    })
}

/// `(id, display name, one-line description)` for the whole library.
pub fn catalog() -> Vec<(&'static str, &'static str, &'static str)> {
    vec![
        ("bell", "Bell pair", "The canonical entangled pair (|00⟩+|11⟩)/√2 — entanglement in two gates."),
        ("ghz", "GHZ state", "n-qubit maximal entanglement (|0…0⟩+|1…1⟩)/√2."),
        ("qft", "Quantum Fourier Transform", "The quantum FFT — Hadamard + controlled-phase ladder."),
        ("grover", "Grover search", "Amplitude amplification finds a marked item in √N steps."),
        ("bv", "Bernstein–Vazirani", "Recover a hidden bitstring in one query — interference at work."),
        ("dj", "Deutsch–Jozsa", "One query decides constant vs balanced."),
        ("teleport", "Teleportation", "Move a state across a Bell pair via corrections."),
        ("superdense", "Superdense coding", "Two classical bits down one qubit."),
        ("trotter", "Trotter Ising", "Transverse-field Ising dynamics — the tensor-network / sparse-Pauli home turf."),
        ("rclifford", "Random Clifford", "A random Clifford circuit — the stabilizer backend's home turf."),
    ]
}

// ===========================================================================
// State analysis  (display-layer f64 — NOT the byte-exact artifact)
// ===========================================================================

const SCALE: f64 = (1u64 << FRAC) as f64;

#[inline]
fn re(a: Amp) -> f64 {
    a.re as f64 / SCALE
}
#[inline]
fn im(a: Amp) -> f64 {
    a.im as f64 / SCALE
}

fn total_weight(sv: &StateVector) -> f64 {
    sv.prob_weights().iter().map(|&w| w as f64).sum::<f64>()
}

/// Normalized per-basis-state probabilities.
pub fn probabilities(sv: &StateVector) -> Vec<f64> {
    let t = total_weight(sv);
    if t <= 0.0 {
        return vec![0.0; sv.amps.len()];
    }
    sv.prob_weights().iter().map(|&w| w as f64 / t).collect()
}

/// `⟨Z_q⟩` — the computational-basis magnetization of qubit `q`.
pub fn expect_z(sv: &StateVector, q: u8) -> f64 {
    let bit = 1usize << q;
    let t = total_weight(sv);
    if t <= 0.0 {
        return 0.0;
    }
    let mut acc = 0.0;
    for (i, w) in sv.prob_weights().iter().enumerate() {
        let s = if i & bit == 0 { 1.0 } else { -1.0 };
        acc += s * (*w as f64);
    }
    acc / t
}

/// `⟨X_q⟩` and `⟨Y_q⟩` — the off-diagonal (coherence) expectations.
pub fn expect_xy(sv: &StateVector, q: u8) -> (f64, f64) {
    let bit = 1usize << q;
    let t = total_weight(sv);
    if t <= 0.0 {
        return (0.0, 0.0);
    }
    let mut ex = 0.0;
    let mut ey = 0.0;
    for i in 0..sv.amps.len() {
        if i & bit == 0 {
            let a = sv.amps[i];
            let b = sv.amps[i | bit];
            // conj(a)·b
            ex += 2.0 * (re(a) * re(b) + im(a) * im(b));
            ey += 2.0 * (re(a) * im(b) - im(a) * re(b));
        }
    }
    (ex * SCALE * SCALE / t, ey * SCALE * SCALE / t)
}

/// The Bloch vector `(x, y, z)` of qubit `q` (its reduced single-qubit state).
pub fn bloch(sv: &StateVector, q: u8) -> [f64; 3] {
    let (x, y) = expect_xy(sv, q);
    [x, y, expect_z(sv, q)]
}

/// Single-qubit **entanglement entropy** of qubit `q` with the rest, in bits:
/// `0` when separable, `1` when maximally entangled. From the Bloch-vector length
/// `r`: eigenvalues of the reduced density matrix are `(1±r)/2`.
pub fn entanglement_entropy(sv: &StateVector, q: u8) -> f64 {
    let b = bloch(sv, q);
    let r = (b[0] * b[0] + b[1] * b[1] + b[2] * b[2]).sqrt().min(1.0);
    let lam = (1.0 + r) / 2.0;
    binary_entropy(lam)
}

/// Purity `Tr(ρ_q²) = (1+r²)/2` of qubit `q` (`1` = pure, `0.5` = maximally mixed).
pub fn purity(sv: &StateVector, q: u8) -> f64 {
    let b = bloch(sv, q);
    let r2 = (b[0] * b[0] + b[1] * b[1] + b[2] * b[2]).min(1.0);
    (1.0 + r2) / 2.0
}

fn binary_entropy(p: f64) -> f64 {
    if p <= 0.0 || p >= 1.0 {
        return 0.0;
    }
    -p * p.log2() - (1.0 - p) * (1.0 - p).log2()
}

/// The `top` most-probable basis states as `(index, probability)`, descending.
pub fn top_outcomes(sv: &StateVector, top: usize) -> Vec<(usize, f64)> {
    let mut v: Vec<(usize, f64)> = probabilities(sv).into_iter().enumerate().collect();
    v.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal));
    v.truncate(top);
    v
}

// ===========================================================================
// Noise simulation — depolarizing channel via Monte-Carlo trajectories
// ===========================================================================

/// `⟨Z_q⟩` under a per-gate **depolarizing** channel, by averaging over
/// `trajectories` random error insertions: after each gate, each qubit it touched
/// suffers a uniformly random `X`/`Y`/`Z` with probability `depol_permille`/1000.
/// Each trajectory reuses the byte-exact dense simulator, in the circuit's own
/// gateset; deterministic in `seed`.
/// `⟨Z_q⟩` decays toward 0 as noise or depth grows — decoherence, made measurable.
pub fn noisy_expect_z(
    circuit: &Circuit,
    q: u8,
    depol_permille: u32,
    trajectories: u32,
    seed: u64,
) -> f64 {
    let mut st = seed.wrapping_mul(0x9E37_79B9).wrapping_add(1);
    let mut acc = 0.0;
    let traj = trajectories.max(1);
    for _ in 0..traj {
        let mut noisy = circuit.empty_like();
        for g in &circuit.ops {
            noisy.ops.push(g.clone());
            for &qb in g.controls.iter().chain(std::iter::once(&g.target)) {
                if (splitmix(&mut st) % 1000) < depol_permille as u64 {
                    let base = match splitmix(&mut st) % 3 {
                        0 => BaseGate::X,
                        1 => BaseGate::Y,
                        _ => BaseGate::Z,
                    };
                    noisy.ops.push(Gate { base, controls: vec![], target: qb, param: 0 });
                }
            }
        }
        if let Ok(sv) = noisy.simulate() {
            acc += expect_z(&sv, q);
        }
    }
    acc / traj as f64
}

// ===========================================================================
// Backend selection — the cheapest-sufficient classical cascade
// ===========================================================================

#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum Backend {
    /// Stabilizer tableau — Clifford circuits in O(n²), scales past statevector.
    Stabilizer,
    /// Exact dense statevector — 2ⁿ amplitudes, the workhorse.
    Statevector,
    /// Tensor-network contraction — the classical referee for large, low-
    /// entanglement, non-Clifford circuits.
    TensorNetwork,
}

/// Is every gate a Clifford operation? (Determines whether the polynomial-time
/// stabilizer path applies.)
pub fn is_clifford(c: &Circuit) -> bool {
    c.ops.iter().all(|g| match g.controls.len() {
        0 => matches!(
            g.base,
            BaseGate::I | BaseGate::X | BaseGate::Y | BaseGate::Z | BaseGate::H | BaseGate::S | BaseGate::Sdg
        ) || (g.base == BaseGate::P && (g.param == 1 || g.param == 2)),
        1 => matches!(g.base, BaseGate::X | BaseGate::Y | BaseGate::Z)
            || (g.base == BaseGate::P && g.param == 1),
        _ => false,
    })
}

/// Recommend the cheapest sufficient classical backend for a circuit, with a
/// one-line rationale.
pub fn recommend_backend(c: &Circuit) -> (Backend, &'static str) {
    if is_clifford(c) {
        (Backend::Stabilizer, "all-Clifford → stabilizer tableau, O(n²); scales far past statevector")
    } else if c.n_qubits <= 22 {
        (Backend::Statevector, "non-Clifford, ≤22 qubits → exact dense statevector")
    } else {
        (Backend::TensorNetwork, "non-Clifford, >22 qubits → tensor-network contraction (the referee)")
    }
}

/// Portable-exact classical work of a dense simulation: `n_ops · 2ⁿ` amplitude
/// updates — recomputable from `(gates, qubits)` alone, exponential in qubits.
pub fn statevector_work(c: &Circuit) -> u64 {
    (c.ops.len() as u64).saturating_mul(1u64 << c.n_qubits.min(63))
}

// ===========================================================================
// Transpiler — route onto nearest-neighbour (line) hardware, verified
// ===========================================================================

/// A circuit routed onto a linear nearest-neighbour coupling.
pub struct Routed {
    pub circuit: Circuit,
    /// `final_map[logical] = physical`.
    pub final_map: Vec<u8>,
    pub swaps: u32,
    pub source_2q: u32,
    pub source_nonlocal: u32,
}

/// Route `src` onto a line topology `0-1-…-(n-1)`, inserting SWAPs (as 3·CX) so
/// every one-control gate acts on adjacent physical qubits. Gates with ≥2 controls
/// are placed on the mapped qubits directly (still correct, not made local).
/// The routed circuit keeps `src`'s gateset.
pub fn route_line(src: &Circuit) -> Routed {
    let n = src.n_qubits;
    let mut phys: Vec<u8> = (0..n).collect();
    let mut logq: Vec<u8> = (0..n).collect();
    let mut out = src.empty_like();
    let mut swaps = 0u32;
    let (mut s2q, mut nonlocal) = (0u32, 0u32);
    let do_swap = |out: &mut Circuit, phys: &mut [u8], logq: &mut [u8], a: u8, b: u8| {
        out.swap(a, b);
        let (la, lb) = (logq[a as usize], logq[b as usize]);
        logq.swap(a as usize, b as usize);
        phys[la as usize] = b;
        phys[lb as usize] = a;
    };
    for g in &src.ops {
        match g.controls.len() {
            0 => {
                out.ops.push(Gate {
                    base: g.base,
                    controls: vec![],
                    target: phys[g.target as usize],
                    param: g.param,
                });
            }
            1 => {
                s2q += 1;
                let (la, lb) = (g.controls[0], g.target);
                let mut pa = phys[la as usize];
                let pb = phys[lb as usize];
                if (pa as i32 - pb as i32).abs() > 1 {
                    nonlocal += 1;
                    // shuffle the control toward the target until adjacent
                    while (pa as i32 - pb as i32).abs() > 1 {
                        let step = if pa < pb { pa + 1 } else { pa - 1 };
                        do_swap(&mut out, &mut phys, &mut logq, pa, step);
                        swaps += 1;
                        pa = step;
                    }
                }
                out.ops.push(Gate {
                    base: g.base,
                    controls: vec![phys[la as usize]],
                    target: phys[lb as usize],
                    param: g.param,
                });
            }
            _ => {
                s2q += 1;
                out.ops.push(Gate {
                    base: g.base,
                    controls: g.controls.iter().map(|&c| phys[c as usize]).collect(),
                    target: phys[g.target as usize],
                    param: g.param,
                });
            }
        }
    }
    Routed { circuit: out, final_map: phys, swaps, source_2q: s2q, source_nonlocal: nonlocal }
}

/// SWAP network restoring the identity mapping from `final_map`.
fn restore_swaps(final_map: &[u8], n: u8) -> Vec<(u8, u8)> {
    let mut logq = vec![0u8; n as usize];
    for (l, &p) in final_map.iter().enumerate() {
        logq[p as usize] = l as u8;
    }
    let mut phys: Vec<u8> = final_map.to_vec();
    let mut out = Vec::new();
    for p in 0..n as usize {
        while logq[p] != p as u8 {
            let cur = phys[p] as usize; // where logical p currently sits
            out.push((p as u8, cur as u8));
            let (lp, lc) = (logq[p], logq[cur]);
            logq[p] = lc;
            logq[cur] = lp;
            phys[lp as usize] = cur as u8;
            phys[lc as usize] = p as u8;
        }
    }
    out
}

/// Is the routed circuit equivalent to the source? Checks that `routed ++
/// (swaps restoring identity)` reconstructs the same byte-exact statevector.
pub fn routing_equivalent(src: &Circuit, routed: &Routed) -> bool {
    let mut full = routed.circuit.clone();
    for (a, b) in restore_swaps(&routed.final_map, src.n_qubits) {
        full.swap(a, b);
    }
    match (full.simulate(), src.simulate()) {
        (Ok(a), Ok(b)) => a.statevector_hash() == b.statevector_hash(),
        _ => false,
    }
}

// ===========================================================================
// OpenQASM interop — the framework-neutral lingua franca
// ===========================================================================

/// Emit OpenQASM 3.0 for a circuit (the dyadic gate set maps to standard gates;
/// `P(k)` becomes `p(2*pi/2^k)`).
///
/// The angles are the ideal ones, whatever the circuit's gateset. A
/// [`crate::quantum::Gateset::V1`] circuit exports angles its own table does
/// not simulate: circuit/1's `P(k)` deviates from `2π/2^k` for `k ≥ 8` and
/// applies no phase for `k ≥ 19`.
pub fn to_qasm(c: &Circuit) -> String {
    let mut s = String::from("OPENQASM 3.0;\ninclude \"stdgates.inc\";\n");
    s.push_str(&format!("qubit[{}] q;\n", c.n_qubits));
    for g in &c.ops {
        let t = g.target;
        let line = match (g.base, g.controls.as_slice()) {
            (BaseGate::I, []) => continue,
            (BaseGate::X, []) => format!("x q[{t}];"),
            (BaseGate::Y, []) => format!("y q[{t}];"),
            (BaseGate::Z, []) => format!("z q[{t}];"),
            (BaseGate::H, []) => format!("h q[{t}];"),
            (BaseGate::S, []) => format!("s q[{t}];"),
            (BaseGate::Sdg, []) => format!("sdg q[{t}];"),
            (BaseGate::T, []) => format!("t q[{t}];"),
            (BaseGate::Tdg, []) => format!("tdg q[{t}];"),
            (BaseGate::P, []) => format!("p(2*pi/{}) q[{t}];", 1u64 << g.param),
            (BaseGate::X, [c0]) => format!("cx q[{c0}],q[{t}];"),
            (BaseGate::Y, [c0]) => format!("cy q[{c0}],q[{t}];"),
            (BaseGate::Z, [c0]) => format!("cz q[{c0}],q[{t}];"),
            (BaseGate::P, [c0]) => format!("cp(2*pi/{}) q[{c0}],q[{t}];", 1u64 << g.param),
            (BaseGate::X, [a, b]) => format!("ccx q[{a}],q[{b}],q[{t}];"),
            (BaseGate::Z, [a, b]) => format!("ccz q[{a}],q[{b}],q[{t}];"),
            (base, ctrls) => format!("// unmapped: {base:?} controls={ctrls:?} target={t}"),
        };
        s.push_str(&line);
        s.push('\n');
    }
    s
}

/// Parse a practical OpenQASM 2/3 subset into a circuit.
///
/// This targets what the mainstream toolchains actually export, not a tidy subset
/// of the grammar: both register forms (`qubit[n] q;` and `qreg q[n];`), both
/// spellings of measurement (`measure q[0] -> c[0];` and the OpenQASM 3 assignment
/// `c[0] = measure q[0];`), the Clifford+T names, the rotations `rx`/`ry`/`rz`,
/// `sx`/`sxdg`, the universal `u`/`u2`/`u3`, and `p`/`cp`/`crz`.
/// `barrier`/`reset`/includes are ignored.
///
/// # Angles
///
/// Accepted as exact rational multiples of `π` — `pi`, `pi/4`, `3*pi/4`, `pi*0.5`,
/// `0`, and negatives — provided the reduced denominator is a power of two. Such an
/// angle needs no approximation: `m·2π/2^k` is emitted as one `P` gate per set bit
/// of `m`, so it costs at most `DYADIC_MAX` gates and is byte-exact.
///
/// Anything else is **refused rather than rounded**. A simulator that silently
/// approximated `π/3` would still return a hash, and that hash would be a lie.
///
/// `rz`, `rx`, `ry`, `sx` and `u*` are exact up to a global phase, which no
/// measurement can observe; `crz` is exact outright, because its leftover factor
/// lands on the control where it *would* be observable.
///
/// # Gateset
///
/// The circuit is built in [`crate::quantum::Gateset::V2`]. A gate that names
/// one qubit twice (`ccx q[0],q[0],q[1]`, `cx q[1],q[1]`) is refused as an
/// invalid circuit when it is read, never built, so the circuit stays valid
/// if its gateset is then changed.
pub fn from_qasm(src: &str) -> Result<Circuit, String> {
    let mut n: Option<u8> = None;
    // first pass: qubit count
    for raw in src.lines() {
        let line = strip(raw);
        if let Some(rest) = line.strip_prefix("qubit[")
            && let Some(end) = rest.find(']')
        {
            n = rest[..end].trim().parse().ok();
        } else if let Some(rest) = line.strip_prefix("qreg ")
            && let (Some(a), Some(b)) = (rest.find('['), rest.find(']'))
        {
            n = rest[a + 1..b].trim().parse().ok();
        }
    }
    let n = n.ok_or("no qubit register declaration (qubit[n] q; or qreg q[n];)")?;
    let mut c = Circuit::new(n);
    for raw in src.lines() {
        let line = strip(raw);
        if line.is_empty() {
            continue;
        }
        // Measurement is recorded, not simulated here, and the exporters disagree
        // about how to spell it: OpenQASM 2 writes `measure q[0] -> c[0];` while
        // Qiskit's OpenQASM 3 exporter writes `c[0] = measure q[0];`, whose first
        // token is the *classical* register. Skip either shape before parsing a
        // gate name out of the line.
        if line.contains("measure") {
            continue;
        }
        let head: String = line.chars().take_while(|c| c.is_ascii_alphanumeric() || *c == '_').collect();
        let head = head.as_str();
        match head {
            "OPENQASM" | "include" | "qubit" | "qreg" | "creg" | "bit" | "measure" | "barrier"
            | "reset" | "gate" => continue,
            _ => {}
        }
        let qs = parse_qubits(&line);
        let push = |c: &mut Circuit, base: BaseGate, controls: Vec<u8>, target: u8, param: u16| -> Result<(), String> {
            distinct_qubits(&controls, target)?;
            c.ops.push(Gate { base, controls, target, param });
            Ok(())
        };
        match head {
            "id" | "i" => {}
            "x" => require1(&qs).and_then(|t| push(&mut c, BaseGate::X, vec![], t, 0))?,
            "y" => require1(&qs).and_then(|t| push(&mut c, BaseGate::Y, vec![], t, 0))?,
            "z" => require1(&qs).and_then(|t| push(&mut c, BaseGate::Z, vec![], t, 0))?,
            "h" => require1(&qs).and_then(|t| push(&mut c, BaseGate::H, vec![], t, 0))?,
            "s" => require1(&qs).and_then(|t| push(&mut c, BaseGate::S, vec![], t, 0))?,
            "sdg" => require1(&qs).and_then(|t| push(&mut c, BaseGate::Sdg, vec![], t, 0))?,
            "t" => require1(&qs).and_then(|t| push(&mut c, BaseGate::T, vec![], t, 0))?,
            "tdg" => require1(&qs).and_then(|t| push(&mut c, BaseGate::Tdg, vec![], t, 0))?,
            "cx" | "cnot" => require2(&qs).and_then(|(a, t)| push(&mut c, BaseGate::X, vec![a], t, 0))?,
            "cy" => require2(&qs).and_then(|(a, t)| push(&mut c, BaseGate::Y, vec![a], t, 0))?,
            "cz" => require2(&qs).and_then(|(a, t)| push(&mut c, BaseGate::Z, vec![a], t, 0))?,
            "swap" => require2(&qs).and_then(|(a, b)| {
                distinct_qubits(&[a], b)?;
                c.swap(a, b);
                Ok(())
            })?,
            "ccx" | "toffoli" => require3(&qs).and_then(|(a, b, t)| push(&mut c, BaseGate::X, vec![a, b], t, 0))?,
            "ccz" => require3(&qs).and_then(|(a, b, t)| push(&mut c, BaseGate::Z, vec![a, b], t, 0))?,
            "p" | "u1" | "phase" => {
                let t = parse_angle(&arg_n(&line, 0)?)?;
                let q = require1(&qs)?;
                push_turn(&mut c, &[], q, t)?;
            }
            "cp" | "cu1" | "cphase" => {
                let t = parse_angle(&arg_n(&line, 0)?)?;
                let (a, q) = require2(&qs)?;
                push_turn(&mut c, &[a], q, t)?;
            }
            // Rotations. These are exact for dyadic angles but carry an
            // unobservable global phase relative to the textbook `r*` matrices:
            // `rz(θ)` is `P(θ)` up to `e^{-iθ/2}`. Measurement statistics are
            // identical, which is what a circuit means.
            "rz" => {
                let t = parse_angle(&arg_n(&line, 0)?)?;
                let q = require1(&qs)?;
                push_turn(&mut c, &[], q, t)?;
            }
            "rx" => {
                let t = parse_angle(&arg_n(&line, 0)?)?;
                let q = require1(&qs)?;
                push_rx(&mut c, q, t)?;
            }
            "ry" => {
                let t = parse_angle(&arg_n(&line, 0)?)?;
                let q = require1(&qs)?;
                push_ry(&mut c, q, t)?;
            }
            // `sx` is `rx(π/2)` — a quarter turn, so `P(2)` between two `H`s.
            "sx" | "sxdg" => {
                let q = require1(&qs)?;
                let t = if head == "sx" {
                    Turn { m: 1, k: 2 }
                } else {
                    Turn { m: 3, k: 2 }
                };
                push_rx(&mut c, q, t)?;
            }
            // `u(θ,φ,λ) = rz(φ)·ry(θ)·rz(λ)`; the rightmost applies first.
            "u" | "u3" | "u2" => {
                let (theta, phi, lam) = if head == "u2" {
                    // u2(φ,λ) = u3(π/2, φ, λ)
                    (Turn { m: 1, k: 2 }, parse_angle(&arg_n(&line, 0)?)?, parse_angle(&arg_n(&line, 1)?)?)
                } else {
                    (
                        parse_angle(&arg_n(&line, 0)?)?,
                        parse_angle(&arg_n(&line, 1)?)?,
                        parse_angle(&arg_n(&line, 2)?)?,
                    )
                };
                let q = require1(&qs)?;
                push_turn(&mut c, &[], q, lam)?;
                push_ry(&mut c, q, theta)?;
                push_turn(&mut c, &[], q, phi)?;
            }
            // Controlled-rz is exact with no phase caveat: `CP(θ)` preceded by
            // `P(−θ/2)` on the control reproduces `diag(1,1,e^{−iθ/2},e^{iθ/2})`.
            "crz" => {
                // Halve the SIGNED angle, then wrap. Wrapping first would pick the
                // wrong branch, because crz has period 4π and its half does not.
                let (p, den) = angle_turns(&arg_n(&line, 0)?)?;
                if den * 2 > (1i128 << crate::quantum::DYADIC_MAX) {
                    return Err("crz angle is finer than the dyadic floor".into());
                }
                let (a, q) = require2(&qs)?;
                distinct_qubits(&[a], q)?;
                push_turn(&mut c, &[], a, to_turn(-p, den * 2))?;
                push_turn(&mut c, &[a], q, to_turn(p, den))?;
            }
            other => {
                return Err(format!(
                    "unsupported gate: '{other}' (supported: id x y z h s sdg t tdg sx sxdg \
                     rx ry rz p/u1/phase u/u2/u3 cx cy cz cp/cu1 crz swap ccx ccz)"
                ));
            }
        }
    }
    c.validate().map_err(|e| e.to_string())?;
    Ok(c)
}

/// Refuse a gate that names one qubit twice, a control repeated or equal to
/// the target, as an invalid circuit in every gateset.
fn distinct_qubits(controls: &[u8], target: u8) -> Result<(), String> {
    for (i, &q) in controls.iter().enumerate() {
        if q == target {
            return Err(QuantumError::Invalid(format!("control {q} equals target")).to_string());
        }
        if controls[..i].contains(&q) {
            return Err(QuantumError::Invalid(format!("control {q} repeated")).to_string());
        }
    }
    Ok(())
}

fn strip(line: &str) -> String {
    let no_comment = line.split("//").next().unwrap_or("");
    no_comment.trim().trim_end_matches(';').trim().to_string()
}

fn parse_qubits(line: &str) -> Vec<u8> {
    let mut out = Vec::new();
    let bytes = line.as_bytes();
    let mut i = 0;
    while i < bytes.len() {
        if bytes[i] == b'[' {
            let mut j = i + 1;
            while j < bytes.len() && bytes[j] != b']' {
                j += 1;
            }
            if let Ok(v) = line[i + 1..j].trim().parse::<u8>() {
                out.push(v);
            }
            i = j + 1;
        } else {
            i += 1;
        }
    }
    out
}

fn require1(qs: &[u8]) -> Result<u8, String> {
    qs.first().copied().ok_or_else(|| "expected 1 qubit".into())
}
fn require2(qs: &[u8]) -> Result<(u8, u8), String> {
    if qs.len() >= 2 {
        Ok((qs[0], qs[1]))
    } else {
        Err("expected 2 qubits".into())
    }
}
fn require3(qs: &[u8]) -> Result<(u8, u8, u8), String> {
    if qs.len() >= 3 {
        Ok((qs[0], qs[1], qs[2]))
    } else {
        Err("expected 3 qubits".into())
    }
}

/// Parse a `p(...)`/`cp(...)` argument into a dyadic `k` where the phase is
/// `2π/2^k`. Accepts `2*pi/2^k`, `2*pi/N` (N a power of two), `pi/N` (→ k=log2(N)+1),
/// and bare `pi` (→ k=1).
/// An angle held exactly, as `m / 2^k` of a full turn with `0 ≤ m < 2^k`.
///
/// The simulator's only parameterised gate is `P(k) = diag(1, e^{2πi/2^k})`, so
/// a single gate reaches just the angles `2π/2^k`. Every *dyadic* angle is still
/// exactly reachable, because the binary expansion of `m` splits it: each set bit
/// `i` contributes `2^i/2^k` of a turn, which is `P(k − i)`. That bounds any angle
/// at `DYADIC_MAX` gates rather than the `m` repetitions a naive reading suggests.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
struct Turn {
    m: u64,
    k: u32,
}

/// Parse a decimal literal as an exact rational. `"0.5"` is `1/2` — no float
/// rounding enters, so `pi*0.5` stays exact.
fn decimal_ratio(tok: &str) -> Option<(i128, i128)> {
    if tok.is_empty() {
        return None;
    }
    match tok.split_once('.') {
        None => tok.parse::<i128>().ok().map(|v| (v, 1)),
        Some((int, frac)) => {
            if frac.is_empty() || !frac.bytes().all(|b| b.is_ascii_digit()) {
                return None;
            }
            let scale = 10i128.checked_pow(frac.len() as u32)?;
            let whole: i128 = if int.is_empty() { 0 } else { int.parse().ok()? };
            let f: i128 = frac.parse().ok()?;
            Some((whole * scale + f, scale))
        }
    }
}

/// Parse an OpenQASM angle expression into an exact fraction of a turn.
///
/// Handles what the mainstream exporters actually emit: `pi`, `pi/4`, `3*pi/4`,
/// `2*pi/8`, `pi*0.5`, `0.5*pi`, `0`, and any of them negated. The angle must be a
/// rational multiple of `π` whose denominator is a power of two once reduced —
/// anything else is refused rather than approximated, because approximating it
/// would quietly cost the byte-exactness the whole crate rests on.
///
/// A bare number with no `pi` is read as radians. That is the one place input is
/// *interpreted* rather than read: it is accepted only when it lands within `1e-9`
/// of an exact dyadic angle, which is the case when an exporter has printed a
/// rounded decimal for an angle it meant exactly.
/// Wrap a signed fraction of a turn into the `[0, 1)` representative the `P`
/// ladder emits.
fn to_turn(p: i128, q: i128) -> Turn {
    let k = q.trailing_zeros();
    let m = ((p % q) + q) % q;
    Turn { m: m as u64, k }
}

/// As [`parse_angle`], but keeping the **signed** angle rather than its residue
/// modulo a full turn.
///
/// The difference matters for half-angle gates. `crz` has period `4π`, not `2π`,
/// so `θ` and `θ ± 2π` have halves that differ by `π`. Wrapping first and halving
/// second silently picks the wrong branch. For `rz`/`rx`/`ry` the same leftover is
/// a *global* phase and unobservable, which is why only `crz` needs this.
fn angle_turns(arg: &str) -> Result<(i128, i128), String> {
    let s = arg.replace(' ', "").to_lowercase();
    if s.is_empty() {
        return Err("empty angle".into());
    }
    let (sign, body) = match s.strip_prefix('-') {
        Some(r) => (-1i128, r),
        None => (1i128, s.strip_prefix('+').unwrap_or(&s)),
    };

    // Split into factors, remembering whether each is multiplied or divided in.
    let mut parts: Vec<(char, String)> = Vec::new();
    let (mut cur, mut op) = (String::new(), '*');
    for ch in body.chars() {
        if ch == '*' || ch == '/' {
            parts.push((op, std::mem::take(&mut cur)));
            op = ch;
        } else {
            cur.push(ch);
        }
    }
    parts.push((op, cur));

    let (mut num, mut den, mut pi_pow) = (sign, 1i128, 0i32);
    for (op, tok) in &parts {
        if tok == "pi" {
            pi_pow += if *op == '*' { 1 } else { -1 };
            continue;
        }
        let (a, b) = decimal_ratio(tok)
            .ok_or_else(|| format!("cannot read '{tok}' in angle '{arg}'"))?;
        if a == 0 && *op == '/' {
            return Err("division by zero in angle".into());
        }
        if *op == '*' {
            num *= a;
            den *= b;
        } else {
            num *= b;
            den *= a;
        }
    }

    if num == 0 {
        return Ok((0, 2));
    }
    if pi_pow == 0 {
        // Bare radians. Accept only if it is a dyadic angle to within 1e-9.
        let v = (num as f64 / den as f64) / std::f64::consts::TAU;
        for k in 1..=crate::quantum::DYADIC_MAX as u32 {
            let scaled = v * (1u64 << k) as f64;
            let r = scaled.round();
            if (scaled - r).abs() < 1e-9 {
                return Ok((r as i128, 1i128 << k));
            }
        }
        return Err(format!(
            "angle '{arg}' is not a dyadic multiple of 2π; this simulator is exact, \
             so it refuses angles it cannot represent rather than rounding them"
        ));
    }
    if pi_pow != 1 {
        return Err(format!("angle '{arg}' is not linear in pi"));
    }

    // angle = (num/den)·π, so in turns that is num/(2·den).
    let mut p = num;
    let mut q = den * 2;
    if q < 0 {
        p = -p;
        q = -q;
    }
    let g = {
        let (mut a, mut b) = (p.unsigned_abs(), q.unsigned_abs());
        while b != 0 {
            let t = a % b;
            a = b;
            b = t;
        }
        a.max(1) as i128
    };
    p /= g;
    q /= g;
    if q <= 0 || (q as u128).count_ones() != 1 {
        return Err(format!(
            "angle '{arg}' is not a dyadic multiple of 2π (reduced denominator {q}); \
             this simulator is exact, so it refuses angles it cannot represent"
        ));
    }
    let k = q.trailing_zeros();
    if k > crate::quantum::DYADIC_MAX as u32 {
        return Err(format!("angle '{arg}' is finer than the dyadic floor (k={k})"));
    }
    Ok((p, q))
}

/// Parse an OpenQASM angle expression into an exact fraction of a turn.
///
/// Handles what the mainstream exporters actually emit: `pi`, `pi/4`, `3*pi/4`,
/// `2*pi/8`, `pi*0.5`, `0.5*pi`, `0`, and any of them negated. The angle must be a
/// rational multiple of `π` whose denominator is a power of two once reduced —
/// anything else is refused rather than approximated, because approximating it
/// would quietly cost the byte-exactness the whole crate rests on.
fn parse_angle(arg: &str) -> Result<Turn, String> {
    let (p, q) = angle_turns(arg)?;
    Ok(to_turn(p, q))
}

/// Read the `i`-th parenthesised argument of a gate line.
fn arg_n(line: &str, i: usize) -> Result<String, String> {
    let a = line.find('(').ok_or("gate needs an argument")?;
    let b = line[a..].find(')').ok_or("unterminated argument")? + a;
    line[a + 1..b]
        .split(',')
        .nth(i)
        .map(|s| s.trim().to_string())
        .ok_or_else(|| format!("expected at least {} argument(s)", i + 1))
}

/// Emit an exact dyadic phase as `P` gates, one per set bit of `m`. A gate
/// naming one qubit twice is refused, even when the phase is a whole turn and
/// emits nothing.
fn push_turn(c: &mut Circuit, controls: &[u8], target: u8, t: Turn) -> Result<(), String> {
    distinct_qubits(controls, target)?;
    for i in 0..=t.k {
        if (t.m >> i) & 1 == 1 {
            let k = t.k - i;
            // k == 0 is a whole turn: the identity, so it contributes nothing.
            if k >= 1 {
                c.ops.push(Gate {
                    base: BaseGate::P,
                    controls: controls.to_vec(),
                    target,
                    param: k as u16,
                });
            }
        }
    }
    Ok(())
}

/// `rx(θ) = H · rz(θ) · H`, exactly, in the dyadic basis.
fn push_rx(c: &mut Circuit, target: u8, t: Turn) -> Result<(), String> {
    c.h(target);
    push_turn(c, &[], target, t)?;
    c.h(target);
    Ok(())
}

/// `ry(θ) = S · rx(θ) · S†`. Circuit order applies the rightmost factor first.
fn push_ry(c: &mut Circuit, target: u8, t: Turn) -> Result<(), String> {
    let one = |b: BaseGate| Gate { base: b, controls: vec![], target, param: 0 };
    c.ops.push(one(BaseGate::Sdg));
    push_rx(c, target, t)?;
    c.ops.push(one(BaseGate::S));
    Ok(())
}

/// The content hash of a circuit — the identity a run receipt binds. It is
/// [`circuit_hash`], so it covers the circuit's gateset, and fails as that
/// does on a circuit its gateset does not allow.
pub fn circuit_id(c: &Circuit) -> Result<[u8; 32], QuantumError> {
    circuit_hash(c)
}

#[cfg(test)]
mod tests {
    #[test]
    #[ignore]
    fn probe_sv_perf() {
        use std::time::Instant;
        use crate::quantum::Circuit;
        println!("\n  statevector: time per gate vs width");
        for n in [8u8, 12, 16, 18, 20, 22] {
            let mut c = Circuit::new(n);
            for _ in 0..3 { for q in 0..n { c.h(q); } for q in 0..n - 1 { c.cx(q, q + 1); } }
            let gates = c.ops.len();
            let t0 = Instant::now();
            let sv = c.simulate().unwrap();
            let el = t0.elapsed();
            println!("   n={n:<3} dim={:<9} {gates:>4} gates  {:>9.2?} total  {:>8.1} ns/gate  {:>7.2} ns/amp-gate",
                1usize << n, el, el.as_nanos() as f64 / gates as f64,
                el.as_nanos() as f64 / (gates as f64 * (1usize << n) as f64));
            std::hint::black_box(sv.statevector_hash());
        }
        println!("\n  gate-kind cost at n=20 (200 gates each)");
        for kind in ["h", "t", "z", "cx", "cz"] {
            let n = 20u8;
            let mut c = Circuit::new(n);
            for i in 0..200 {
                let q = (i % n as usize) as u8;
                match kind {
                    "h" => { c.h(q); }
                    "t" => { c.t(q); }
                    "z" => { c.z(q); }
                    "cx" => { c.cx(q, (q + 1) % n); }
                    _ => { c.cz(q, (q + 1) % n); }
                }
            }
            let t0 = Instant::now();
            let sv = c.simulate().unwrap();
            let el = t0.elapsed();
            println!("   {kind:<3} {:>8.1} ns/gate", el.as_nanos() as f64 / 200.0);
            std::hint::black_box(sv.statevector_hash());
        }
    }

    use super::*;

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

    #[test]
    fn bell_is_maximally_entangled() {
        let sv = bell().simulate().unwrap();
        let p = probabilities(&sv);
        assert!(approx(p[0b00], 0.5, 1e-4) && approx(p[0b11], 0.5, 1e-4));
        assert!(approx(p[0b01], 0.0, 1e-4) && approx(p[0b10], 0.0, 1e-4));
        // each qubit is maximally entangled with the other
        assert!(approx(entanglement_entropy(&sv, 0), 1.0, 1e-3), "S={}", entanglement_entropy(&sv, 0));
        assert!(approx(purity(&sv, 0), 0.5, 1e-3));
    }

    #[test]
    fn ghz_all_or_nothing() {
        let sv = ghz(4).simulate().unwrap();
        let p = probabilities(&sv);
        assert!(approx(p[0], 0.5, 1e-4) && approx(p[15], 0.5, 1e-4));
        assert!(approx(entanglement_entropy(&sv, 2), 1.0, 1e-3));
    }

    #[test]
    fn plus_state_bloch() {
        // H|0> = |+>  → Bloch (1,0,0), zero entanglement (product state)
        let mut c = Circuit::new(1);
        c.h(0);
        let sv = c.simulate().unwrap();
        let b = bloch(&sv, 0);
        assert!(approx(b[0], 1.0, 1e-3) && approx(b[1], 0.0, 1e-3) && approx(b[2], 0.0, 1e-3), "bloch={b:?}");
        assert!(approx(entanglement_entropy(&sv, 0), 0.0, 1e-3));
    }

    #[test]
    fn grover_amplifies_the_marked_state() {
        for &(n, marked) in &[(3u8, 5u32), (4, 11), (5, 20)] {
            let sv = grover(n, marked).simulate().unwrap();
            let top = top_outcomes(&sv, 1)[0];
            assert_eq!(top.0, marked as usize, "grover n={n} should peak at {marked}");
            assert!(top.1 > 0.6, "marked prob {} too low", top.1);
        }
    }

    #[test]
    fn bernstein_vazirani_recovers_secret() {
        let secret = 0b10110u32;
        let sv = bernstein_vazirani(secret, 5).simulate().unwrap();
        // marginal over the 5 input qubits should concentrate on `secret`
        let top = top_outcomes(&sv, 1)[0];
        assert_eq!(top.0 & 0b11111, secret as usize, "BV should reveal the secret");
    }

    #[test]
    fn deutsch_jozsa_distinguishes() {
        let sv_c = deutsch_jozsa(false, 4).simulate().unwrap();
        let sv_b = deutsch_jozsa(true, 4).simulate().unwrap();
        // constant → inputs all-zero; balanced → inputs nonzero
        let inputs = |i: usize| i & 0b1111;
        assert_eq!(inputs(top_outcomes(&sv_c, 1)[0].0), 0, "constant → all-zero inputs");
        assert_ne!(inputs(top_outcomes(&sv_b, 1)[0].0), 0, "balanced → nonzero inputs");
    }

    #[test]
    fn superdense_roundtrips_all_four() {
        for bits in 0u8..4 {
            let sv = superdense(bits).simulate().unwrap();
            let top = top_outcomes(&sv, 1)[0];
            assert_eq!(top.0, bits as usize, "superdense {bits} decoded wrong");
            assert!(top.1 > 0.99);
        }
    }

    #[test]
    fn clifford_detection() {
        assert!(is_clifford(&bell()));
        assert!(is_clifford(&ghz(5)));
        assert!(!is_clifford(&teleportation())); // has T
        assert!(!is_clifford(&Circuit::qft(4))); // has P(k>2)
        assert_eq!(recommend_backend(&bell()).0, Backend::Stabilizer);
        assert_eq!(recommend_backend(&Circuit::qft(6)).0, Backend::Statevector);
    }

    #[test]
    fn qasm_round_trip() {
        for c in [bell(), ghz(3), Circuit::qft(4), grover(3, 5)] {
            let qasm = to_qasm(&c);
            let back = from_qasm(&qasm).expect("reparse");
            // same reconstructed state ⇒ same circuit semantics
            assert_eq!(
                c.simulate().unwrap().statevector_hash(),
                back.simulate().unwrap().statevector_hash(),
                "QASM round-trip changed the state"
            );
        }
    }

    #[test]
    fn qasm_imports_foreign_circuit() {
        let src = "OPENQASM 3.0;\ninclude \"stdgates.inc\";\nqubit[2] q;\nh q[0];\ncx q[0], q[1];\n";
        let c = from_qasm(src).unwrap();
        assert_eq!(c.simulate().unwrap().statevector_hash(), bell().simulate().unwrap().statevector_hash());
    }

    /// Equality of circuits *as physics*: the same state up to the global phase
    /// that no measurement can see. `rz` and `P` differ by exactly that factor, so
    /// hash equality would be the wrong test here.
    fn same_state(a: &str, b: &Circuit) {
        let got = from_qasm(a).unwrap_or_else(|e| panic!("parse {a:?}: {e}"));
        let f = got.simulate().unwrap().fidelity_fx(&b.simulate().unwrap());
        assert!(
            (f - crate::quantum::ONE).abs() < 1 << 12,
            "{a:?} gave fidelity {f} against the reference, want 1"
        );
    }

    /// Qiskit's OpenQASM 3 exporter writes measurement as an assignment, so the
    /// first token of the line is the classical register rather than a gate. This
    /// is its default output shape, so failing it means failing the format the
    /// most widely used toolchain emits.
    #[test]
    fn qasm3_assignment_measure_parses() {
        let src = "OPENQASM 3.0;\ninclude \"stdgates.inc\";\nqubit[2] q;\nbit[2] c;\n\
                   h q[0];\ncx q[0], q[1];\nc[0] = measure q[0];\nc[1] = measure q[1];\n";
        same_state(src, &bell());
    }

    /// The rotations the mainstream exporters emit, each checked against the gate
    /// it must equal rather than merely accepted by the parser.
    #[test]
    fn rotations_decompose_to_the_dyadic_basis() {
        let mut x = Circuit::new(1);
        x.x(0);
        let mut y = Circuit::new(1);
        y.y(0);
        let mut z = Circuit::new(1);
        z.z(0);
        let mut t = Circuit::new(1);
        t.t(0);

        same_state("qreg q[1];\nrx(pi) q[0];", &x);
        same_state("qreg q[1];\nry(pi) q[0];", &y);
        same_state("qreg q[1];\nrz(pi) q[0];", &z);
        same_state("qreg q[1];\np(pi/4) q[0];", &t);
        // u3(π,0,π) is X — the canonical single-qubit-gate spelling
        same_state("qreg q[1];\nu3(pi,0,pi) q[0];", &x);
        same_state("qreg q[1];\nu(pi,0,pi) q[0];", &x);
        // two √X make an X
        same_state("qreg q[1];\nsx q[0];\nsx q[0];", &x);
    }

    /// A general dyadic angle is not one `P` gate — `3π/4` is `P(3)` then `P(2)`,
    /// found from the binary expansion of the numerator. Check it against the
    /// explicit pair, and check the bit-splitting really happened.
    #[test]
    fn a_general_dyadic_angle_splits_across_gates() {
        let mut want = Circuit::new(1);
        want.h(0).p(3, 0).p(2, 0);
        let got = from_qasm("qreg q[1];\nh q[0];\np(3*pi/4) q[0];").unwrap();
        assert_eq!(
            got.simulate().unwrap().statevector_hash(),
            want.simulate().unwrap().statevector_hash()
        );
        assert_eq!(got.ops.len(), 3, "3π/4 should split into two P gates after the H");
    }

    /// Negative and zero angles, which exporters emit constantly and the old
    /// parser refused outright. A rotation and its inverse must cancel exactly.
    #[test]
    fn negative_and_zero_angles_are_exact() {
        let plain = Circuit::new(1);
        same_state("qreg q[1];\nrz(0) q[0];", &plain);
        same_state("qreg q[1];\np(0) q[0];", &plain);
        same_state("qreg q[1];\nrz(pi/2) q[0];\nrz(-pi/2) q[0];", &plain);
        same_state("qreg q[1];\nsx q[0];\nsxdg q[0];", &plain);
        // `pi*0.5` is the decimal spelling; 0.5 is an exact binary fraction
        let mut s = Circuit::new(1);
        s.s(0);
        same_state("qreg q[1];\np(pi*0.5) q[0];", &s);
        same_state("qreg q[1];\np(0.5*pi) q[0];", &s);
    }

    /// Controlled-`rz` is exact with no phase caveat at all, because the leftover
    /// factor lands on the control where it is observable — and representable.
    #[test]
    fn controlled_rz_is_exact() {
        // A rotation and its inverse must cancel to nothing — checked on a
        // superposed control, where crz and cp are observably different, so this
        // cannot pass by accidentally implementing the wrong one.
        let mut hh = Circuit::new(2);
        hh.h(0).h(1);
        same_state(
            "qreg q[2];\nh q[0];\nh q[1];\ncrz(pi/2) q[0],q[1];\ncrz(-pi/2) q[0],q[1];",
            &hh,
        );
        let crz = from_qasm("qreg q[2];\nh q[0];\nh q[1];\ncrz(pi/2) q[0],q[1];").unwrap();
        let mut cp = Circuit::new(2);
        cp.h(0).h(1).cp(2, 0, 1);
        assert_ne!(
            crz.simulate().unwrap().statevector_hash(),
            cp.simulate().unwrap().statevector_hash(),
            "crz must not silently be cp — they differ by a phase on the control"
        );
    }

    /// The refusal matters as much as the acceptance: an angle this simulator
    /// cannot hold exactly is rejected rather than rounded, because rounding it
    /// would quietly cost the byte-exactness everything else depends on.
    #[test]
    fn non_dyadic_angles_are_refused_not_rounded() {
        for src in ["qreg q[1];\nrz(pi/3) q[0];", "qreg q[1];\np(pi/5) q[0];"] {
            let e = from_qasm(src).expect_err("a non-dyadic angle must not be accepted");
            assert!(e.contains("dyadic"), "error should say why: {e}");
        }
    }

    #[test]
    fn catalog_all_build_and_simulate() {
        for (id, _, _) in catalog() {
            let c = build_algorithm(id, 4, 5).expect(id);
            assert!(c.simulate().is_ok(), "{id} failed to simulate");
        }
    }

    #[test]
    fn routing_preserves_the_state() {
        // route several circuits onto a line and prove equivalence via statevector.
        for c in [bell(), ghz(5), Circuit::qft(5), trotter_ising(6, 3)] {
            let r = route_line(&c);
            assert!(routing_equivalent(&c, &r), "routed circuit must reproduce the state");
            // every routed one-control gate is now nearest-neighbour on the line
            for g in &r.circuit.ops {
                if g.controls.len() == 1 {
                    let d = (g.controls[0] as i32 - g.target as i32).abs();
                    assert_eq!(d, 1, "routed 2q gate must be local: {} {}", g.controls[0], g.target);
                }
            }
        }
    }

    #[test]
    fn noise_decoheres_the_signal() {
        // a qubit prepared in |1> reads ⟨Z⟩ = -1 with no noise; depolarizing noise
        // pulls it toward 0.
        let mut c = Circuit::new(2);
        c.x(0);
        for _ in 0..6 {
            c.cx(0, 1);
            c.cx(0, 1); // identity pair, but each gate is a noise opportunity
        }
        let e0 = noisy_expect_z(&c, 0, 0, 300, 1); // no noise → exact
        let e1 = noisy_expect_z(&c, 0, 40, 500, 1); // light noise
        let e2 = noisy_expect_z(&c, 0, 150, 500, 1); // heavier noise
        assert!(approx(e0, -1.0, 1e-9), "no noise ⇒ exact: {e0}");
        // ⟨Z⟩ decoheres monotonically from -1 toward 0 as noise grows
        assert!(e1 > e0 && e1 < 0.0, "light noise decoheres: {e1}");
        assert!(e2 > e1, "more noise ⇒ more decoherence: {e2} > {e1}");
    }

    #[test]
    fn showcase_circuits_target_the_right_backends() {
        // Trotter Ising is non-Clifford (has P(3)=T-family) → routed to statevector/
        // tensor-network, its home turf.
        let ti = trotter_ising(6, 4);
        assert!(!is_clifford(&ti), "Trotter Ising should be non-Clifford");
        assert!(ti.simulate().is_ok());
        // Random Clifford → stabilizer tier.
        let rc = random_clifford(6, 60, 7);
        assert!(is_clifford(&rc), "random Clifford must be Clifford");
        assert_eq!(random_clifford(6, 60, 7), random_clifford(6, 60, 7), "deterministic in seed");
    }

    /// A circuit derived from another keeps its gateset: a noisy trajectory
    /// and a routed circuit run under the source's table, not the default.
    #[test]
    fn derived_circuits_keep_their_gateset() {
        use crate::quantum::Gateset;
        // H·T·H reads ⟨Z⟩ through the table, so the two gatesets give
        // different expectations.
        let htz = |gs| {
            let mut c = Circuit::with_gateset(1, gs);
            c.h(0).t(0).h(0);
            c
        };
        let exact = |gs| expect_z(&htz(gs).simulate().unwrap(), 0);
        assert_ne!(exact(Gateset::V1).to_bits(), exact(Gateset::V2).to_bits());
        for gs in Gateset::ALL {
            // No noise: each trajectory is the circuit itself.
            assert_eq!(noisy_expect_z(&htz(gs), 0, 0, 3, 1).to_bits(), exact(gs).to_bits(), "{gs:?}");
            let mut src = Circuit::qft(5);
            src.gateset = gs;
            let r = route_line(&src);
            assert_eq!(r.circuit.gateset, gs);
            assert!(routing_equivalent(&src, &r));
        }
    }

    /// A QASM gate naming one qubit twice is refused as an invalid circuit,
    /// never built, whatever gateset the circuit is later given.
    #[test]
    fn qasm_gates_that_repeat_a_qubit_are_invalid() {
        for (gate, want) in [
            ("ccx q[0],q[0],q[1];", "control 0 repeated"),
            ("ccz q[2],q[2],q[0];", "control 2 repeated"),
            ("cx q[1],q[1];", "control 1 equals target"),
            ("cz q[0],q[0];", "control 0 equals target"),
            ("ccx q[0],q[1],q[1];", "control 1 equals target"),
            ("swap q[1],q[1];", "control 1 equals target"),
            ("cp(pi/2) q[0],q[0];", "control 0 equals target"),
            // A whole turn emits no gate, and is still a gate naming q[0] twice.
            ("cp(0) q[0],q[0];", "control 0 equals target"),
            ("crz(pi/2) q[2],q[2];", "control 2 equals target"),
        ] {
            let e = from_qasm(&format!("qreg q[3];\n{gate}")).expect_err(gate);
            assert_eq!(e, format!("invalid circuit: {want}"), "{gate}");
        }
        let ok = from_qasm("qreg q[3];\nccx q[0],q[1],q[2];\ncp(pi/2) q[0],q[2];").unwrap();
        assert_eq!(ok.gateset, crate::quantum::Gateset::V2);
        assert_eq!(ok.ops.len(), 2);
    }
}