wai-quantum 0.3.18

//! WAI Quantum Toolchain — the compute core of a classical-hardware quantum IDE
//! (extensions/quantum-sim § Toolchain, `wai.quantum.toolchain`).
//!
//! Everyone teaches quantum circuits-first, inside one vendor's SDK, attached to a
//! cloud QPU, behind a wall of linear algebra. This is the opposite: 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, no account, no vendor lock-in — 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.

use crate::quantum::{circuit_hash, Amp, BaseGate, Circuit, Gate, 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; 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::new(circuit.n_qubits);
        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).
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 = Circuit::new(n);
    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)`).
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. Supports `qubit[n] q;` /
/// `qreg q[n];`, the standard named gates, and `p(...)`/`cp(...)` with a
/// power-of-two dyadic phase. `measure`/`barrier`/includes are ignored.
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;
        }
        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| {
            c.ops.push(Gate { base, controls, target, param });
        };
        match head {
            "id" | "i" => {}
            "x" => require1(&qs).map(|t| push(&mut c, BaseGate::X, vec![], t, 0))?,
            "y" => require1(&qs).map(|t| push(&mut c, BaseGate::Y, vec![], t, 0))?,
            "z" => require1(&qs).map(|t| push(&mut c, BaseGate::Z, vec![], t, 0))?,
            "h" => require1(&qs).map(|t| push(&mut c, BaseGate::H, vec![], t, 0))?,
            "s" => require1(&qs).map(|t| push(&mut c, BaseGate::S, vec![], t, 0))?,
            "sdg" => require1(&qs).map(|t| push(&mut c, BaseGate::Sdg, vec![], t, 0))?,
            "t" => require1(&qs).map(|t| push(&mut c, BaseGate::T, vec![], t, 0))?,
            "tdg" => require1(&qs).map(|t| push(&mut c, BaseGate::Tdg, vec![], t, 0))?,
            "cx" | "cnot" => require2(&qs).map(|(a, t)| push(&mut c, BaseGate::X, vec![a], t, 0))?,
            "cy" => require2(&qs).map(|(a, t)| push(&mut c, BaseGate::Y, vec![a], t, 0))?,
            "cz" => require2(&qs).map(|(a, t)| push(&mut c, BaseGate::Z, vec![a], t, 0))?,
            "swap" => require2(&qs).map(|(a, b)| {
                c.swap(a, b);
            })?,
            "ccx" | "toffoli" => require3(&qs).map(|(a, b, t)| push(&mut c, BaseGate::X, vec![a, b], t, 0))?,
            "ccz" => require3(&qs).map(|(a, b, t)| push(&mut c, BaseGate::Z, vec![a, b], t, 0))?,
            "p" | "u1" | "phase" => {
                let k = dyadic_k(&line)?;
                require1(&qs).map(|t| push(&mut c, BaseGate::P, vec![], t, k))?
            }
            "cp" | "cu1" | "cphase" => {
                let k = dyadic_k(&line)?;
                require2(&qs).map(|(a, t)| push(&mut c, BaseGate::P, vec![a], t, k))?
            }
            other => return Err(format!("unsupported gate: '{other}'")),
        }
    }
    c.validate().map_err(|e| e.to_string())?;
    Ok(c)
}

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).
fn dyadic_k(line: &str) -> Result<u16, String> {
    let a = line.find('(').ok_or("phase gate needs an argument")?;
    let b = line[a..].find(')').ok_or("unterminated argument")? + a;
    let arg = line[a + 1..b].replace(' ', "").to_lowercase();
    let pow2 = |m: u64| -> Option<u16> {
        if m.is_power_of_two() {
            Some(m.trailing_zeros() as u16)
        } else {
            None
        }
    };
    if arg == "pi" {
        return Ok(1);
    }
    if let Some(rest) = arg.strip_prefix("2*pi/") {
        let den: u64 = rest.parse().map_err(|_| "bad denominator")?;
        return pow2(den).ok_or_else(|| format!("non-dyadic phase 2*pi/{den}"));
    }
    if let Some(rest) = arg.strip_prefix("pi/") {
        let den: u64 = rest.parse().map_err(|_| "bad denominator")?;
        return pow2(den).map(|k| k + 1).ok_or_else(|| format!("non-dyadic phase pi/{den}"));
    }
    Err(format!("unsupported (non-dyadic) phase argument: '{arg}'"))
}

/// The content hash of a circuit — the identity a run receipt binds.
pub fn circuit_id(c: &Circuit) -> [u8; 32] {
    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());
    }

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