prism-q 0.32.0

//! Sparse state-vector simulation backend.
//!
//! Stores only non-zero amplitudes in a map keyed by basis-state index, giving
//! O(k) memory where k is the number of non-zero basis states. Entries whose
//! squared amplitude is at or below a pruning threshold (default 1e-16,
//! raised via `SparseBackend::set_prune_epsilon`) are dropped after gates
//! that can shrink or cancel amplitudes, and the kept entries are rescaled
//! so the state keeps its norm.
//!
//! Every gate walks the map, so the map hashes basis-state indices with the
//! crate's multiply-xor hasher rather than the stdlib default.
//!
//! # Memory layout
//!
//! - `FxHashMap<usize, Complex64>` keyed by basis-state index, non-zero
//!   amplitudes only.
//! - A retained swap-buffer map lets amplitude-moving gates rebuild the state
//!   without reallocating.
//!
//! # Gate support
//!
//! The full gate set, including MCU and every fused or batched variant.
//! Diagonal gates scale amplitudes in place without moving keys. Shot
//! sampling and Pauli expectations run natively on the map.
//!
//! # When to prefer this backend
//!
//! - States with few non-zero amplitudes (computational basis states, limited superposition).
//! - Large qubit counts where the state stays sparse throughout the circuit.
//! - Classical-like circuits with limited branching.
//!
//! # When NOT to use this backend
//!
//! - After a layer of Hadamard gates (state becomes maximally dense).
//! - Small qubit counts where dense statevector is faster due to HashMap overhead.
//! - Circuits wider than `usize::BITS` qubits, which the basis index cannot
//!   address. [`Backend::init`] rejects them and `Auto` routes them to MPS.

use num_complex::Complex64;
use rand::RngExt;
use rand::SeedableRng;
use rand_chacha::ChaCha8Rng;

#[cfg(feature = "parallel")]
use rayon::prelude::*;

#[cfg(feature = "parallel")]
const MIN_STATES_FOR_PAR: usize = 4096;

#[cfg(feature = "parallel")]
const MIN_SHOTS_FOR_PAR: usize = 32;

/// Shots drawn per ChaCha8 substream. One stream per shot block rather than
/// per shot: cipher setup and the first keystream block cost more than an
/// entire draw when the CDF is short, so the block grain amortizes them.
const SHOTS_PER_STREAM: usize = 256;

use crate::backend::{
    Backend, BasisSamples, dense_probability_len, dense_statevector_len, is_phase_one, overlap,
    reduced_density, reserve_dense_output, schmidt,
};
use crate::circuit::Instruction;
use crate::error::Result;
use crate::gates::{Gate, diag_entries_phase, is_antidiagonal_2x2, is_diagonal_2x2};
use crate::hash::FxHashMap;
use crate::sim::unified_pauli::PauliTerm;

const DEFAULT_EPSILON: f64 = 1e-16;

/// Widest circuit the basis-index representation addresses.
///
/// Every gate kernel masks the global basis index with `1usize << qubit`, so a
/// qubit index must stay below `usize::BITS`.
pub(crate) const MAX_SPARSE_INDEX_QUBITS: usize = usize::BITS as usize;

/// Sparse state-vector backend, O(k) where k is the number of non-zero amplitudes.
pub struct SparseBackend {
    num_qubits: usize,
    state: FxHashMap<usize, Complex64>,
    swap_buf: FxHashMap<usize, Complex64>,
    classical_bits: Vec<bool>,
    rng: ChaCha8Rng,
    epsilon: f64,
    pruned_weight: f64,
    /// Whether the threshold has exceeded the default since the last `init`,
    /// so lowering it mid-run cannot hide weight already discarded.
    raised: bool,
    /// [`crate::backend::max_sparse_entries`] read once at construction, so the
    /// per-gate growth check is a field compare rather than an atomic load.
    entry_cap: usize,
}

impl SparseBackend {
    pub fn new(seed: u64) -> Self {
        Self {
            num_qubits: 0,
            state: FxHashMap::default(),
            swap_buf: FxHashMap::default(),
            classical_bits: Vec::new(),
            rng: ChaCha8Rng::seed_from_u64(seed),
            epsilon: DEFAULT_EPSILON,
            pruned_weight: 0.0,
            raised: false,
            entry_cap: crate::backend::max_sparse_entries(),
        }
    }

    /// Number of nonzero amplitudes currently stored, the `k` in the O(k) cost
    /// of every gate walk.
    pub fn entry_count(&self) -> usize {
        self.state.len()
    }

    /// Set the pruning threshold on squared amplitude magnitude.
    ///
    /// Entries with `norm_sqr` at or below `epsilon` are dropped after gates
    /// that can shrink or cancel amplitudes; 0 drops exact zeros only. The
    /// construction default of 1e-16 removes only numerical dust and the run
    /// reports as exact; a larger threshold trades state weight for a smaller
    /// map, the run reports `Approximate`, and the dropped weight feeds the
    /// metadata's `fidelity_lower_bound` as a first-order estimate. Lowering
    /// the threshold again keeps the run reporting `Approximate` until the
    /// next [`Backend::init`]; the threshold itself survives `init`, the
    /// accumulated weight does not. Every prune that discards weight rescales
    /// the kept entries to the norm the state had before it, so probabilities,
    /// samples, and expectations all see a unit total and the discarded mass
    /// shows up only in the bound.
    ///
    /// # Panics
    /// Panics unless `0 <= epsilon < 1`.
    #[cfg(test)]
    pub(crate) fn set_prune_epsilon(&mut self, epsilon: f64) {
        assert!(
            (0.0..1.0).contains(&epsilon),
            "prune epsilon must lie in [0, 1)"
        );
        self.epsilon = epsilon;
        self.raised = self.raised || epsilon > DEFAULT_EPSILON;
    }

    /// The rescale pass runs only when the dropped mass moves the scale off
    /// 1.0 in f64, so the default threshold's dust never pays for it.
    #[inline(always)]
    fn prune(&mut self) {
        let eps = self.epsilon;
        let mut dropped = 0.0;
        let mut kept = 0.0;
        self.state.retain(|_, amp| {
            let weight = amp.norm_sqr();
            if weight > eps {
                kept += weight;
                true
            } else {
                dropped += weight;
                false
            }
        });
        self.pruned_weight += dropped;
        if dropped > 0.0 && kept > 0.0 {
            let scale = ((kept + dropped) / kept).sqrt();
            if scale != 1.0 {
                for amp in self.state.values_mut() {
                    *amp *= scale;
                }
            }
        }
    }

    /// Reject a gate whose worst-case fan-out would grow the map past the
    /// entry budget. `factor` is the per-source-entry fan-out of the caller
    /// (2 for a branching 1q gate, 4 for a dense 2q), so rejection can fire
    /// one gate early on a state that would have deduplicated below the cap.
    #[inline(always)]
    fn check_entry_growth(&self, factor: usize) -> Result<()> {
        let projected = self.state.len().saturating_mul(factor);
        if projected > self.entry_cap {
            return Err(self.entry_growth_error(projected));
        }
        Ok(())
    }

    #[cold]
    fn entry_growth_error(&self, projected: usize) -> crate::error::PrismError {
        crate::error::PrismError::IncompatibleBackend {
            backend: "sparse".to_string(),
            reason: format!(
                "state holds {} entries and this gate can reach {projected}, \
                 exceeding the cap of {} entries on this machine \
                 (set PRISM_MAX_SPARSE_QUBITS to override)",
                self.state.len(),
                self.entry_cap
            ),
        }
    }

    #[inline(always)]
    fn apply_single_qubit(&mut self, target: usize, mat: [[Complex64; 2]; 2]) -> Result<()> {
        if is_diagonal_2x2(&mat) {
            self.apply_diagonal_1q(target, mat[0][0], mat[1][1]);
            return Ok(());
        }
        if is_antidiagonal_2x2(&mat) {
            self.apply_antidiagonal_1q(target, mat[0][1], mat[1][0]);
            return Ok(());
        }
        self.check_entry_growth(2)?;

        let mask = 1usize << target;
        let zero = Complex64::new(0.0, 0.0);
        self.swap_buf.clear();
        self.swap_buf.reserve(self.state.len() * 2);

        for (&idx, &amp) in &self.state {
            let bit = (idx >> target) & 1;
            let partner = idx ^ mask;

            *self.swap_buf.entry(idx).or_insert(zero) += mat[bit][bit] * amp;
            *self.swap_buf.entry(partner).or_insert(zero) += mat[1 - bit][bit] * amp;
        }

        std::mem::swap(&mut self.state, &mut self.swap_buf);
        self.prune();
        Ok(())
    }

    /// A diagonal 2x2 scales amplitudes in place: no partner entries, no map rebuild.
    /// A sub-unit diagonal (a Kraus operator) can shrink amplitudes below epsilon,
    /// so only that case pays the prune pass.
    #[inline(always)]
    fn apply_diagonal_1q(&mut self, target: usize, d0: Complex64, d1: Complex64) {
        for (idx, amp) in self.state.iter_mut() {
            *amp *= if (*idx >> target) & 1 == 1 { d1 } else { d0 };
        }
        if d0.norm_sqr() < 1.0 - 1e-12 || d1.norm_sqr() < 1.0 - 1e-12 {
            self.prune();
        }
    }

    /// An antidiagonal 2x2 maps each basis state to exactly one partner, the
    /// 1q case of the invariant `apply_cx` states, so the map is a 1:1 remap
    /// with no accumulation. A sub-unit scale (a normalized Kraus branch
    /// routed through a fused payload) can still shrink amplitudes below
    /// epsilon, so only that case pays the prune pass, as `apply_monomial_2q`
    /// does.
    #[inline(always)]
    fn apply_antidiagonal_1q(&mut self, target: usize, a01: Complex64, a10: Complex64) {
        let mask = 1usize << target;
        self.swap_buf.clear();
        self.swap_buf.reserve(self.state.len());
        self.swap_buf.extend(self.state.drain().map(|(idx, amp)| {
            let scale = if idx & mask == 0 { a10 } else { a01 };
            (idx ^ mask, amp * scale)
        }));
        std::mem::swap(&mut self.state, &mut self.swap_buf);
        if a01.norm_sqr() < 1.0 - 1e-12 || a10.norm_sqr() < 1.0 - 1e-12 {
            self.prune();
        }
    }

    /// CX is a deterministic 1:1 index mapping. No near-zero amplitudes are created.
    #[inline(always)]
    fn apply_cx(&mut self, control: usize, target: usize) {
        let ctrl_mask = 1usize << control;
        let tgt_mask = 1usize << target;
        self.swap_buf.clear();
        self.swap_buf.reserve(self.state.len());
        self.swap_buf.extend(self.state.drain().map(|(idx, amp)| {
            if idx & ctrl_mask != 0 {
                (idx ^ tgt_mask, amp)
            } else {
                (idx, amp)
            }
        }));
        std::mem::swap(&mut self.state, &mut self.swap_buf);
    }

    #[inline(always)]
    fn apply_cz(&mut self, q0: usize, q1: usize) {
        let mask0 = 1usize << q0;
        let mask1 = 1usize << q1;
        for (&idx, amp) in self.state.iter_mut() {
            if idx & mask0 != 0 && idx & mask1 != 0 {
                *amp = -*amp;
            }
        }
    }

    #[inline(always)]
    fn apply_swap(&mut self, q0: usize, q1: usize) {
        let m0 = 1usize << q0;
        let m1 = 1usize << q1;
        self.swap_buf.clear();
        self.swap_buf.reserve(self.state.len());
        self.swap_buf.extend(self.state.drain().map(|(idx, amp)| {
            let bit0 = (idx >> q0) & 1;
            let bit1 = (idx >> q1) & 1;
            if bit0 != bit1 {
                (idx ^ m0 ^ m1, amp)
            } else {
                (idx, amp)
            }
        }));
        std::mem::swap(&mut self.state, &mut self.swap_buf);
    }

    #[inline(always)]
    fn apply_cu(&mut self, control: usize, target: usize, mat: [[Complex64; 2]; 2]) -> Result<()> {
        self.check_entry_growth(2)?;
        let ctrl_mask = 1usize << control;
        let tgt_mask = 1usize << target;
        let zero = Complex64::new(0.0, 0.0);
        self.swap_buf.clear();
        self.swap_buf.reserve(self.state.len() * 2);

        for (&idx, &amp) in &self.state {
            if idx & ctrl_mask == 0 {
                *self.swap_buf.entry(idx).or_insert(zero) += amp;
            } else {
                let bit = (idx >> target) & 1;
                let partner = idx ^ tgt_mask;
                *self.swap_buf.entry(idx).or_insert(zero) += mat[bit][bit] * amp;
                *self.swap_buf.entry(partner).or_insert(zero) += mat[1 - bit][bit] * amp;
            }
        }

        std::mem::swap(&mut self.state, &mut self.swap_buf);
        self.prune();
        Ok(())
    }

    #[inline(always)]
    fn apply_mcu(
        &mut self,
        controls: &[usize],
        target: usize,
        mat: [[Complex64; 2]; 2],
    ) -> Result<()> {
        self.check_entry_growth(2)?;
        let ctrl_mask: usize = controls.iter().map(|&q| 1usize << q).fold(0, |a, b| a | b);
        let tgt_mask = 1usize << target;
        let zero = Complex64::new(0.0, 0.0);
        self.swap_buf.clear();
        self.swap_buf.reserve(self.state.len() * 2);

        for (&idx, &amp) in &self.state {
            if idx & ctrl_mask != ctrl_mask {
                *self.swap_buf.entry(idx).or_insert(zero) += amp;
            } else {
                let bit = (idx >> target) & 1;
                let partner = idx ^ tgt_mask;
                *self.swap_buf.entry(idx).or_insert(zero) += mat[bit][bit] * amp;
                *self.swap_buf.entry(partner).or_insert(zero) += mat[1 - bit][bit] * amp;
            }
        }

        std::mem::swap(&mut self.state, &mut self.swap_buf);
        self.prune();
        Ok(())
    }

    #[inline(always)]
    fn apply_cu_phase(&mut self, control: usize, target: usize, phase: Complex64) {
        let ctrl_mask = 1usize << control;
        let tgt_mask = 1usize << target;
        for (&idx, amp) in self.state.iter_mut() {
            if idx & ctrl_mask != 0 && idx & tgt_mask != 0 {
                *amp *= phase;
            }
        }
    }

    #[inline(always)]
    fn apply_mcu_phase(&mut self, controls: &[usize], target: usize, phase: Complex64) {
        let ctrl_mask: usize = controls.iter().map(|&q| 1usize << q).fold(0, |a, b| a | b);
        let tgt_mask = 1usize << target;
        for (&idx, amp) in self.state.iter_mut() {
            if idx & ctrl_mask == ctrl_mask && idx & tgt_mask != 0 {
                *amp *= phase;
            }
        }
    }

    #[inline(always)]
    fn apply_rzz(&mut self, q0: usize, q1: usize, theta: f64) {
        let phase_same = Complex64::from_polar(1.0, -theta / 2.0);
        let phase_diff = Complex64::from_polar(1.0, theta / 2.0);
        for (idx, amp) in self.state.iter_mut() {
            let parity = ((*idx >> q0) ^ (*idx >> q1)) & 1;
            *amp *= if parity == 0 { phase_same } else { phase_diff };
        }
    }

    /// One walk over the map with the per-edge phase pair precomputed, instead of
    /// one walk (and two `from_polar`) per edge.
    fn apply_batch_rzz(&mut self, edges: &[(usize, usize, f64)]) {
        let phases: Vec<(usize, usize, [Complex64; 2])> = edges
            .iter()
            .map(|&(q0, q1, theta)| {
                let same = Complex64::from_polar(1.0, -theta / 2.0);
                let diff = Complex64::from_polar(1.0, theta / 2.0);
                (q0, q1, [same, diff])
            })
            .collect();

        for (idx, amp) in self.state.iter_mut() {
            let mut combined = Complex64::new(1.0, 0.0);
            for &(q0, q1, pair) in &phases {
                combined *= pair[((*idx >> q0) ^ (*idx >> q1)) & 1];
            }
            *amp *= combined;
        }
    }

    fn apply_batch_phase(&mut self, control: usize, phases: &[(usize, Complex64)]) {
        let ctrl_mask = 1usize << control;
        let one = Complex64::new(1.0, 0.0);
        for (&idx, amp) in self.state.iter_mut() {
            if idx & ctrl_mask == 0 {
                continue;
            }
            let mut combined = one;
            for &(target, phase) in phases {
                if idx & (1usize << target) != 0 {
                    combined *= phase;
                }
            }
            if !is_phase_one(combined) {
                *amp *= combined;
            }
        }
    }

    /// Exact structural test for a monomial 4x4: one nonzero per source basis
    /// state and per destination. Structural zeros only, no tolerance, so a
    /// matrix product carrying float dust off its pattern stays on the general
    /// path.
    fn monomial_4x4(mat: &[[Complex64; 4]; 4]) -> Option<([usize; 4], [Complex64; 4])> {
        let zero = Complex64::new(0.0, 0.0);
        let mut dest = [usize::MAX; 4];
        let mut scale = [zero; 4];
        let mut used = 0u8;
        for (col, mat_row) in mat.iter().enumerate() {
            for (row, &coeff) in mat_row.iter().enumerate() {
                if coeff != zero {
                    if dest[row] != usize::MAX {
                        return None;
                    }
                    dest[row] = col;
                    scale[row] = coeff;
                }
            }
        }
        for &d in &dest {
            if d == usize::MAX || used & (1 << d) != 0 {
                return None;
            }
            used |= 1 << d;
        }
        Some((dest, scale))
    }

    /// A monomial 4x4 maps each occupied basis state to exactly one
    /// destination, so no cancellation can arise, the invariant `apply_cx`
    /// states. A sub-unit scale (a normalized Kraus branch routed through a
    /// fused payload) can still shrink amplitudes below epsilon, so only that
    /// case pays the prune pass, as `apply_diagonal_1q` does.
    fn apply_monomial_2q(&mut self, q0: usize, q1: usize, dest: [usize; 4], scale: [Complex64; 4]) {
        if dest == [0, 1, 2, 3] {
            for (idx, amp) in self.state.iter_mut() {
                let row = ((*idx >> q0) & 1) * 2 + ((*idx >> q1) & 1);
                *amp *= scale[row];
            }
        } else {
            let mask0 = 1usize << q0;
            let mask1 = 1usize << q1;
            self.swap_buf.clear();
            self.swap_buf.reserve(self.state.len());
            self.swap_buf.extend(self.state.drain().map(|(idx, amp)| {
                let row = ((idx >> q0) & 1) * 2 + ((idx >> q1) & 1);
                let col = dest[row];
                let new_idx = idx & !(mask0 | mask1) | (((col >> 1) & 1) << q0) | ((col & 1) << q1);
                (new_idx, amp * scale[row])
            }));
            std::mem::swap(&mut self.state, &mut self.swap_buf);
        }
        if scale.iter().any(|c| c.norm_sqr() < 1.0 - 1e-12) {
            self.prune();
        }
    }

    fn apply_fused_2q(&mut self, q0: usize, q1: usize, mat: &[[Complex64; 4]; 4]) -> Result<()> {
        if let Some((dest, scale)) = Self::monomial_4x4(mat) {
            self.apply_monomial_2q(q0, q1, dest, scale);
            return Ok(());
        }
        self.check_entry_growth(4)?;
        let mask0 = 1usize << q0;
        let mask1 = 1usize << q1;
        let zero = Complex64::new(0.0, 0.0);
        self.swap_buf.clear();
        self.swap_buf.reserve(self.state.len() * 2);

        for (&idx, &amp) in &self.state {
            let bit0 = (idx >> q0) & 1;
            let bit1 = (idx >> q1) & 1;
            let row = bit0 * 2 + bit1;
            let base = idx & !(mask0 | mask1);

            for (col, mat_row) in mat.iter().enumerate() {
                let coeff = mat_row[row];
                if coeff == zero {
                    continue;
                }
                let col_bit0 = (col >> 1) & 1;
                let col_bit1 = col & 1;
                let dest = base | (col_bit0 << q0) | (col_bit1 << q1);
                *self.swap_buf.entry(dest).or_insert(zero) += coeff * amp;
            }
        }

        std::mem::swap(&mut self.state, &mut self.swap_buf);
        self.prune();
        Ok(())
    }

    fn masked_prob(&self, mask: usize, bit_set: bool) -> f64 {
        #[cfg(feature = "parallel")]
        if self.state.len() >= MIN_STATES_FOR_PAR {
            return self
                .state
                .par_iter()
                .filter(|&(&idx, _)| (idx & mask != 0) == bit_set)
                .map(|(_, amp)| amp.norm_sqr())
                .sum();
        }
        self.state
            .iter()
            .filter(|&(&idx, _)| (idx & mask != 0) == bit_set)
            .map(|(_, amp)| amp.norm_sqr())
            .sum()
    }

    fn apply_reset(&mut self, qubit: usize) {
        let mask = 1usize << qubit;
        let prob_one = self.masked_prob(mask, true);
        let outcome = self.rng.random::<f64>() < prob_one;
        let inv_norm = crate::backend::measurement_inv_norm(outcome, prob_one);

        if !outcome {
            self.state.retain(|&idx, amp| {
                if idx & mask == 0 {
                    *amp *= inv_norm;
                    true
                } else {
                    false
                }
            });
            return;
        }

        self.swap_buf.clear();
        self.swap_buf.extend(
            self.state
                .drain()
                .filter(|(idx, _)| idx & mask != 0)
                .map(|(idx, amp)| (idx ^ mask, amp * inv_norm)),
        );
        std::mem::swap(&mut self.state, &mut self.swap_buf);
    }

    fn apply_measure(&mut self, qubit: usize, classical_bit: usize) {
        let mask = 1usize << qubit;
        let prob_one = self.masked_prob(mask, true);

        let outcome = self.rng.random::<f64>() < prob_one;
        self.classical_bits[classical_bit] = outcome;

        let inv_norm = crate::backend::measurement_inv_norm(outcome, prob_one);

        self.state.retain(|&idx, amp| {
            let matches = (idx & mask != 0) == outcome;
            if matches {
                *amp *= inv_norm;
            }
            matches
        });
    }

    fn dispatch_gate(&mut self, gate: &Gate, targets: &[usize]) -> Result<()> {
        match gate {
            Gate::Rzz(theta) => {
                self.apply_rzz(targets[0], targets[1], *theta);
            }
            Gate::Cx => {
                self.apply_cx(targets[0], targets[1]);
            }
            Gate::Cz => {
                self.apply_cz(targets[0], targets[1]);
            }
            Gate::Swap => {
                self.apply_swap(targets[0], targets[1]);
            }
            Gate::Cu(mat) => {
                if let Some(phase) = gate.controlled_phase() {
                    self.apply_cu_phase(targets[0], targets[1], phase);
                } else {
                    self.apply_cu(targets[0], targets[1], **mat)?;
                }
            }
            Gate::Mcu(data) => {
                let num_ctrl = data.num_controls as usize;
                if let Some(phase) = gate.controlled_phase() {
                    self.apply_mcu_phase(&targets[..num_ctrl], targets[num_ctrl], phase);
                } else {
                    self.apply_mcu(&targets[..num_ctrl], targets[num_ctrl], data.mat)?;
                }
            }
            Gate::BatchPhase(data) => {
                self.apply_batch_phase(targets[0], &data.phases);
            }
            Gate::BatchRzz(data) => {
                self.apply_batch_rzz(&data.edges);
            }
            Gate::DiagonalBatch(data) => {
                for (idx, amp) in self.state.iter_mut() {
                    *amp *= diag_entries_phase(*idx, &data.entries);
                }
            }
            Gate::MultiFused(data) => {
                for &(target, mat) in &data.gates {
                    self.apply_single_qubit(target, mat)?;
                }
            }
            Gate::Fused2q(mat) => {
                self.apply_fused_2q(targets[0], targets[1], mat)?;
            }
            Gate::Multi2q(data) => {
                for &(q0, q1, ref mat) in &data.gates {
                    self.apply_fused_2q(q0, q1, mat)?;
                }
            }
            other => {
                debug_assert!(
                    targets.len() == 1,
                    "sparse dispatch_gate: unexpected multi-qubit gate {:?}",
                    other
                );
                let mat = other.matrix_2x2();
                self.apply_single_qubit(targets[0], mat)?;
            }
        }
        Ok(())
    }
}

impl Backend for SparseBackend {
    fn name(&self) -> &'static str {
        "sparse"
    }

    fn as_any(&self) -> Option<&dyn std::any::Any> {
        Some(self)
    }

    fn resolved(&self) -> crate::sim::ResolvedBackend {
        crate::sim::ResolvedBackend::Sparse
    }

    /// At the default threshold pruning removes only numerical dust and the
    /// route is exact; a raised threshold can discard real weight, so the
    /// route reports `Approximate` whether or not this run pruned anything,
    /// the same convention the MPS backend uses for its bond cap.
    fn exactness(&self) -> crate::sim::Exactness {
        if self.raised {
            crate::sim::Exactness::Approximate {
                fidelity_lower_bound: Some((1.0 - self.pruned_weight).max(0.0)),
            }
        } else {
            crate::sim::Exactness::Exact
        }
    }

    fn init(&mut self, num_qubits: usize, num_classical_bits: usize) -> Result<()> {
        if num_qubits > MAX_SPARSE_INDEX_QUBITS {
            return Err(crate::error::PrismError::IncompatibleBackend {
                backend: "sparse".to_string(),
                reason: format!(
                    "a {num_qubits}-qubit circuit exceeds the {MAX_SPARSE_INDEX_QUBITS}-qubit \
                     basis-index width of the sparse backend"
                ),
            });
        }
        self.num_qubits = num_qubits;
        self.state.clear();
        self.state.insert(0, Complex64::new(1.0, 0.0));
        self.pruned_weight = 0.0;
        self.raised = self.epsilon > DEFAULT_EPSILON;
        crate::backend::init_classical_bits(&mut self.classical_bits, num_classical_bits);
        Ok(())
    }

    fn apply(&mut self, instruction: &Instruction) -> Result<()> {
        match instruction {
            Instruction::Gate { gate, targets } => self.dispatch_gate(gate, targets)?,
            Instruction::Measure {
                qubit,
                classical_bit,
            } => {
                self.apply_measure(*qubit, *classical_bit);
            }
            Instruction::Reset { qubit } => {
                self.apply_reset(*qubit);
            }
            Instruction::Barrier { .. } => {}
            Instruction::Conditional {
                condition,
                gate,
                targets,
            } => {
                if condition.evaluate(&self.classical_bits) {
                    self.dispatch_gate(gate, targets)?;
                }
            }
            Instruction::Region(region) => self.apply_region(region)?,
        }
        Ok(())
    }

    fn reset(&mut self, qubit: usize) -> Result<()> {
        self.apply_reset(qubit);
        Ok(())
    }

    fn apply_1q_matrix(&mut self, qubit: usize, matrix: &[[Complex64; 2]; 2]) -> Result<()> {
        self.apply_single_qubit(qubit, *matrix)
    }

    fn reduced_density_matrix_1q(&self, qubit: usize) -> Result<[[Complex64; 2]; 2]> {
        let mask = 1usize << qubit;
        let mut p0 = 0.0f64;
        let mut p1 = 0.0f64;
        let mut r = Complex64::new(0.0, 0.0);

        for (&idx, &amp) in &self.state {
            if idx & mask == 0 {
                p0 += amp.norm_sqr();
                if let Some(&amp_one) = self.state.get(&(idx | mask)) {
                    r += amp_one * amp.conj();
                }
            } else {
                p1 += amp.norm_sqr();
            }
        }

        Ok([
            [Complex64::new(p0, 0.0), r.conj()],
            [r, Complex64::new(p1, 0.0)],
        ])
    }

    /// Nonzeros sorted by their traced index, then an outer product within
    /// each run: `O(nnz log nnz)` for the sort plus `O(nnz * g)` products,
    /// `g` the largest run. The sort fixes the summation order, so the answer
    /// does not depend on map iteration order. Scaled to trace one.
    fn reduced_density_matrix(&mut self, subsystem: &[usize]) -> Result<Vec<Complex64>> {
        schmidt::validate_qubit_set(subsystem, self.num_qubits)?;
        let dim = reduced_density::reduced_density_side(self.name(), subsystem.len())?;
        let named = subsystem.iter().fold(0usize, |mask, &q| mask | (1 << q));
        let mut entries: Vec<(usize, usize, Complex64)> = self
            .state
            .iter()
            .map(|(&idx, &amp)| {
                let t = subsystem
                    .iter()
                    .enumerate()
                    .fold(0, |t, (i, &q)| t | (((idx >> q) & 1) << i));
                (idx & !named, t, amp)
            })
            .collect();
        entries.sort_unstable_by_key(|&(e, t, _)| (e, t));

        let mut rho = vec![Complex64::new(0.0, 0.0); dim * dim];
        for run in entries.chunk_by(|a, b| a.0 == b.0) {
            for &(_, t, a) in run {
                for &(_, tp, b) in run {
                    rho[t * dim + tp] += a * b.conj();
                }
            }
        }
        reduced_density::normalize_trace(&mut rho, dim);
        Ok(rho)
    }

    /// One pass over this map with a lookup in the other, `O(nnz)` and no
    /// width cap: only the basis states both hold carry weight into the inner
    /// product.
    fn overlap_sq(&self, other: &dyn Backend) -> Result<f64> {
        if let Some(sparse) = other
            .as_any()
            .and_then(|any| any.downcast_ref::<SparseBackend>())
        {
            if sparse.num_qubits == self.num_qubits {
                let mut inner = Complex64::new(0.0, 0.0);
                for (index, amp) in &self.state {
                    if let Some(theirs) = sparse.state.get(index) {
                        inner += amp.conj() * theirs;
                    }
                }
                let left: f64 = self.state.values().map(Complex64::norm_sqr).sum();
                let right: f64 = sparse.state.values().map(Complex64::norm_sqr).sum();
                return Ok(overlap::normalized(inner.norm_sqr(), left, right));
            }
        }
        overlap::export_overlap_sq(
            self.name(),
            self.num_qubits(),
            || self.export_statevector(),
            other,
        )
    }

    fn classical_results(&self) -> &[bool] {
        &self.classical_bits
    }

    fn probabilities(&self) -> Result<Vec<f64>> {
        let dim = dense_probability_len(self.name(), self.num_qubits)?;
        let mut probs = Vec::new();
        reserve_dense_output(&mut probs, dim, self.name(), "probabilities")?;
        probs.resize(dim, 0.0f64);
        for (&idx, amp) in &self.state {
            probs[idx] = amp.norm_sqr();
        }
        Ok(probs)
    }

    fn num_qubits(&self) -> usize {
        self.num_qubits
    }

    fn supports_native_sampling(&self) -> bool {
        true
    }

    /// Samples from a CDF over the `k` stored amplitudes instead of `2^n`.
    ///
    /// Entries are sorted by basis index so the draw is a deterministic
    /// function of the state, not of map iteration order. Each block of
    /// `SHOTS_PER_STREAM` shots draws sequentially from a ChaCha8
    /// substream keyed on the seed and the block index (streams 2 and up;
    /// the MPS sampler keys per shot, the block grain here amortizes cipher
    /// setup over draws cheaper than the setup). Blocks fill in parallel
    /// only when the CDF holds `MIN_STATES_FOR_PAR` entries; below that,
    /// fork-join dispatch costs more than the whole sequential pass. Either
    /// path walks the same partition, so the words are identical at any
    /// thread count; the drawn bitstrings differ from the dense route's
    /// single-stream draws by design.
    fn sample_basis_states(&mut self, num_shots: usize, seed: u64) -> Result<BasisSamples> {
        let mut indices: Vec<usize> = self.state.keys().copied().collect();
        indices.sort_unstable();
        let probs: Vec<f64> = indices
            .iter()
            .map(|idx| self.state[idx].norm_sqr())
            .collect();
        let cdf = crate::sim::shots::build_cdf(&probs);

        let mut samples = BasisSamples::new(num_shots, self.num_qubits);
        let words_per_shot = samples.words_per_shot();
        let sample_block = |block: usize, block_words: &mut [u64]| {
            let mut rng = ChaCha8Rng::seed_from_u64(seed);
            rng.set_stream(block as u64 + 2);
            for shot_words in block_words.chunks_mut(words_per_shot) {
                let r: f64 = rng.random();
                shot_words[0] = indices[crate::sim::shots::sample_from_cdf(&cdf, r)] as u64;
            }
        };

        #[cfg(feature = "parallel")]
        if num_shots >= MIN_SHOTS_FOR_PAR && cdf.len() >= MIN_STATES_FOR_PAR {
            samples
                .words_mut()
                .par_chunks_mut(SHOTS_PER_STREAM * words_per_shot)
                .enumerate()
                .for_each(|(block, block_words)| sample_block(block, block_words));
            return Ok(samples);
        }

        for (block, block_words) in samples
            .words_mut()
            .chunks_mut(SHOTS_PER_STREAM * words_per_shot)
            .enumerate()
        {
            sample_block(block, block_words);
        }
        Ok(samples)
    }

    fn supports_pauli_expectation(&self) -> bool {
        true
    }

    /// Every observable is accumulated in one walk of the stored amplitudes.
    /// A Z-only observable has `xmask == 0`, so its partner is the entry
    /// itself and it costs a popcount instead of a hash lookup; the two
    /// families are accumulated separately for that reason.
    fn pauli_expectations(&self, observables: &[Vec<PauliTerm>]) -> Result<Vec<f64>> {
        let masks = observables
            .iter()
            .map(|observable| crate::sim::pauli_masks(observable, self.num_qubits))
            .collect::<Result<Vec<_>>>()?;

        let norm: f64 = self.state.values().map(|amp| amp.norm_sqr()).sum();
        if norm == 0.0 {
            return Ok(vec![0.0; masks.len()]);
        }

        let z_only: Vec<usize> = masks
            .iter()
            .filter(|&&(xmask, _, _)| xmask == 0)
            .map(|&(_, zmask, _)| zmask)
            .collect();
        let general: Vec<(usize, usize)> = masks
            .iter()
            .filter(|&&(xmask, _, _)| xmask != 0)
            .map(|&(xmask, zmask, _)| (xmask, zmask))
            .collect();

        let mut z_sum = vec![0.0f64; z_only.len()];
        let mut g_sum = vec![Complex64::new(0.0, 0.0); general.len()];
        for (&idx, &amp) in &self.state {
            let parity_sign = |mask: usize| {
                if (idx & mask).count_ones() & 1 == 1 {
                    -1.0
                } else {
                    1.0
                }
            };
            let norm_sqr = amp.norm_sqr();
            for (slot, &zmask) in z_sum.iter_mut().zip(&z_only) {
                *slot += norm_sqr * parity_sign(zmask);
            }
            for (slot, &(xmask, zmask)) in g_sum.iter_mut().zip(&general) {
                let Some(&partner) = self.state.get(&(idx ^ xmask)) else {
                    continue;
                };
                *slot += partner.conj() * amp * parity_sign(zmask);
            }
        }

        Ok(crate::sim::finish_expectations(
            &masks, &z_sum, &g_sum, norm,
        ))
    }

    fn export_statevector(&self) -> Result<Vec<Complex64>> {
        let dim = dense_statevector_len(self.name(), "statevector export", self.num_qubits)?;
        let mut sv = Vec::new();
        reserve_dense_output(&mut sv, dim, self.name(), "statevector export")?;
        sv.resize(dim, Complex64::new(0.0, 0.0));
        for (&idx, &amp) in &self.state {
            sv[idx] = amp;
        }
        Ok(sv)
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::circuit::Circuit;
    use crate::sim;

    const EPS: f64 = 1e-12;

    fn run_sparse(circuit: &Circuit) -> SparseBackend {
        let mut b = SparseBackend::new(42);
        sim::run_on(&mut b, circuit).unwrap();
        b
    }

    fn run_sparse_probs(circuit: &Circuit) -> Vec<f64> {
        let b = run_sparse(circuit);
        b.probabilities().unwrap()
    }

    #[test]
    fn test_init_zero_state() {
        let mut b = SparseBackend::new(42);
        b.init(3, 0).unwrap();
        assert_eq!(b.state.len(), 1);
        assert!((b.state[&0].re - 1.0).abs() < EPS);
    }

    #[test]
    fn test_x_gate() {
        let mut c = Circuit::new(1, 0);
        c.add_gate(Gate::X, &[0]);
        let b = run_sparse(&c);
        assert_eq!(b.state.len(), 1);
        assert!(b.state.contains_key(&1));
        assert!((b.state[&1].norm() - 1.0).abs() < EPS);
    }

    #[test]
    fn test_h_creates_superposition() {
        let mut c = Circuit::new(1, 0);
        c.add_gate(Gate::H, &[0]);
        let b = run_sparse(&c);
        assert_eq!(b.state.len(), 2);
        assert!((b.state[&0].norm_sqr() - 0.5).abs() < EPS);
        assert!((b.state[&1].norm_sqr() - 0.5).abs() < EPS);
    }

    #[test]
    fn test_hh_is_identity() {
        let mut c = Circuit::new(1, 0);
        c.add_gate(Gate::H, &[0]);
        c.add_gate(Gate::H, &[0]);
        let b = run_sparse(&c);
        assert_eq!(b.state.len(), 1);
        assert!((b.state[&0].re - 1.0).abs() < EPS);
    }

    #[test]
    fn test_cx_bell_state() {
        let mut c = Circuit::new(2, 0);
        c.add_gate(Gate::H, &[0]);
        c.add_gate(Gate::Cx, &[0, 1]);
        let b = run_sparse(&c);
        assert_eq!(b.state.len(), 2);
        assert!((b.state[&0].norm_sqr() - 0.5).abs() < EPS);
        assert!((b.state[&3].norm_sqr() - 0.5).abs() < EPS);
    }

    #[test]
    fn test_cz_phase() {
        let mut c = Circuit::new(2, 0);
        c.add_gate(Gate::X, &[0]);
        c.add_gate(Gate::X, &[1]);
        c.add_gate(Gate::Cz, &[0, 1]);
        let b = run_sparse(&c);
        assert_eq!(b.state.len(), 1);
        assert!((b.state[&3].re - (-1.0)).abs() < EPS);
    }

    #[test]
    fn test_swap() {
        let mut c = Circuit::new(2, 0);
        c.add_gate(Gate::X, &[1]);
        c.add_gate(Gate::Swap, &[0, 1]);
        let b = run_sparse(&c);
        assert_eq!(b.state.len(), 1);
        assert!(b.state.contains_key(&1));
    }

    #[test]
    fn test_rx_pi() {
        let mut c = Circuit::new(1, 0);
        c.add_gate(Gate::Rx(std::f64::consts::PI), &[0]);
        let probs = run_sparse_probs(&c);
        assert!(probs[0].abs() < EPS);
        assert!((probs[1] - 1.0).abs() < EPS);
    }

    #[test]
    fn test_rz_preserves_sparsity() {
        let mut c = Circuit::new(1, 0);
        c.add_gate(Gate::Rz(1.234), &[0]);
        let b = run_sparse(&c);
        assert_eq!(b.state.len(), 1);
        assert!((b.state[&0].norm() - 1.0).abs() < EPS);
    }

    #[test]
    fn test_measure_collapses() {
        let mut c = Circuit::new(1, 1);
        c.add_gate(Gate::H, &[0]);
        c.add_measure(0, 0);
        let b = run_sparse(&c);
        assert_eq!(b.state.len(), 1);
        let outcome = b.classical_results()[0];
        if outcome {
            assert!(b.state.contains_key(&1));
        } else {
            assert!(b.state.contains_key(&0));
        }
    }

    #[test]
    fn test_measure_deterministic() {
        let mut c = Circuit::new(1, 1);
        c.add_gate(Gate::H, &[0]);
        c.add_measure(0, 0);

        let b1 = run_sparse(&c);
        let b2 = run_sparse(&c);
        assert_eq!(b1.classical_results()[0], b2.classical_results()[0]);
    }

    #[test]
    fn test_probs_bell() {
        let mut c = Circuit::new(2, 0);
        c.add_gate(Gate::H, &[0]);
        c.add_gate(Gate::Cx, &[0, 1]);
        let probs = run_sparse_probs(&c);
        assert!((probs[0] - 0.5).abs() < EPS);
        assert!(probs[1].abs() < EPS);
        assert!(probs[2].abs() < EPS);
        assert!((probs[3] - 0.5).abs() < EPS);
    }

    #[test]
    fn test_probs_zero_state() {
        let c = Circuit::new(3, 0);
        let probs = run_sparse_probs(&c);
        assert!((probs[0] - 1.0).abs() < EPS);
        let rest: f64 = probs[1..].iter().sum();
        assert!(rest.abs() < EPS);
    }

    #[test]
    fn test_pruning() {
        let mut b = SparseBackend::new(42);
        b.init(1, 0).unwrap();
        b.state.insert(1, Complex64::new(1e-20, 0.0));
        assert_eq!(b.state.len(), 2);
        b.prune();
        assert_eq!(b.state.len(), 1);
        assert!(b.state.contains_key(&0));
    }

    #[test]
    fn test_fused_gate() {
        let h_mat = Gate::H.matrix_2x2();
        let t_mat = Gate::T.matrix_2x2();
        let zero = Complex64::new(0.0, 0.0);
        let mut fused = [[zero; 2]; 2];
        for i in 0..2 {
            for j in 0..2 {
                for k in 0..2 {
                    fused[i][j] += t_mat[i][k] * h_mat[k][j];
                }
            }
        }

        let mut c1 = Circuit::new(1, 0);
        c1.add_gate(Gate::H, &[0]);
        c1.add_gate(Gate::T, &[0]);
        let p1 = run_sparse_probs(&c1);

        let mut c2 = Circuit::new(1, 0);
        c2.add_gate(Gate::Fused(Box::new(fused)), &[0]);
        let p2 = run_sparse_probs(&c2);

        for (a, b) in p1.iter().zip(p2.iter()) {
            assert!((a - b).abs() < EPS);
        }
    }

    #[test]
    fn test_default_epsilon_reports_exact() {
        let mut c = Circuit::new(2, 0);
        c.add_gate(Gate::H, &[0]);
        c.add_gate(Gate::Cx, &[0, 1]);
        let b = run_sparse(&c);
        assert_eq!(b.epsilon, DEFAULT_EPSILON);
        assert_eq!(b.exactness(), crate::sim::Exactness::Exact);
    }

    #[test]
    fn test_zero_epsilon_still_drops_exact_zeros() {
        let mut b = SparseBackend::new(42);
        b.set_prune_epsilon(0.0);
        b.init(1, 0).unwrap();
        b.state.insert(1, Complex64::new(0.0, 0.0));
        let zero = Complex64::new(0.0, 0.0);
        let half = Complex64::new(0.5, 0.0);
        b.apply_1q_matrix(0, &[[half, zero], [zero, half]]).unwrap();
        assert_eq!(b.state.len(), 1);
        assert!(b.state.contains_key(&0));
    }

    #[test]
    fn test_lowering_epsilon_keeps_the_approximate_report() {
        let theta = 2.0 * (0.1_f64).sqrt().asin();
        let mut b = SparseBackend::new(42);
        b.set_prune_epsilon(0.2);
        b.init(1, 0).unwrap();
        b.apply_1q_matrix(0, &Gate::Ry(theta).matrix_2x2()).unwrap();
        b.set_prune_epsilon(DEFAULT_EPSILON);
        match b.exactness() {
            crate::sim::Exactness::Approximate {
                fidelity_lower_bound: Some(bound),
            } => assert!((bound - 0.9).abs() < EPS),
            other => panic!("lowering the threshold hid discarded weight: {other:?}"),
        }
    }

    #[test]
    fn test_raised_epsilon_drops_weight_and_reports_bound() {
        // Ry puts weight 0.1 on |1>; a 0.2 threshold drops it at the prune
        // after the general-path gate.
        let theta = 2.0 * (0.1_f64).sqrt().asin();
        let mut b = SparseBackend::new(42);
        b.set_prune_epsilon(0.2);
        b.init(1, 0).unwrap();
        b.apply(&Instruction::Gate {
            gate: Gate::Ry(theta),
            targets: crate::circuit::SmallVec::from_slice(&[0]),
        })
        .unwrap();
        assert_eq!(b.state.len(), 1);
        assert!(b.state.contains_key(&0));
        match b.exactness() {
            crate::sim::Exactness::Approximate {
                fidelity_lower_bound: Some(bound),
            } => assert!((bound - 0.9).abs() < EPS),
            other => panic!("expected a bounded Approximate, got {other:?}"),
        }
    }

    // H puts 0.5 on each of |00>, |01>; Ry on q1 then splits each into
    // 0.45 and 0.05. A 0.06 threshold drops both 0.05 entries, and the two
    // kept entries must come back to 0.5 each rather than 0.45, which a
    // unit-total CDF fold would otherwise hand to the highest kept state.
    #[test]
    fn test_prune_renormalizes_kept_entries() {
        let theta = 2.0 * (0.1_f64).sqrt().asin();
        let mut b = SparseBackend::new(42);
        b.set_prune_epsilon(0.06);
        b.init(2, 0).unwrap();
        b.apply_1q_matrix(0, &Gate::H.matrix_2x2()).unwrap();
        b.apply_1q_matrix(1, &Gate::Ry(theta).matrix_2x2()).unwrap();
        assert_eq!(b.state.len(), 2);

        let probs = b.probabilities().unwrap();
        let total: f64 = probs.iter().sum();
        assert!((total - 1.0).abs() < EPS, "probabilities sum to {total}");
        assert!((probs[0] - 0.5).abs() < EPS);
        assert!((probs[1] - 0.5).abs() < EPS);
        assert_eq!(probs[2], 0.0);
        assert_eq!(probs[3], 0.0);

        let z1 = b.pauli_expectations(&[vec![PauliTerm::z(1)]]).unwrap();
        assert!((z1[0] - 1.0).abs() < EPS);
        match b.exactness() {
            crate::sim::Exactness::Approximate {
                fidelity_lower_bound: Some(bound),
            } => assert!((bound - 0.9).abs() < EPS),
            other => panic!("expected a bounded Approximate, got {other:?}"),
        }
    }

    #[test]
    fn test_epsilon_survives_init_weight_does_not() {
        let theta = 2.0 * (0.1_f64).sqrt().asin();
        let mut b = SparseBackend::new(42);
        b.set_prune_epsilon(0.2);
        b.init(1, 0).unwrap();
        b.apply_1q_matrix(0, &Gate::Ry(theta).matrix_2x2()).unwrap();
        assert!(b.pruned_weight > 0.0);
        b.init(1, 0).unwrap();
        match b.exactness() {
            crate::sim::Exactness::Approximate {
                fidelity_lower_bound: Some(bound),
            } => assert_eq!(bound, 1.0),
            other => panic!("expected Approximate after init, got {other:?}"),
        }
    }

    #[test]
    #[should_panic(expected = "prune epsilon")]
    fn test_prune_epsilon_rejects_one() {
        SparseBackend::new(42).set_prune_epsilon(1.0);
    }

    #[test]
    fn test_y_gate_amplitude() {
        let mut c = Circuit::new(1, 0);
        c.add_gate(Gate::Y, &[0]);
        let b = run_sparse(&c);
        assert_eq!(b.state.len(), 1);
        assert!((b.state[&1] - Complex64::new(0.0, 1.0)).norm() < EPS);
    }

    // Distinct phases on the two antidiagonal coefficients, so a swapped
    // a01/a10 assignment fails on amplitudes, not just norms.
    #[test]
    fn test_antidiagonal_fused_amplitudes() {
        let t = Gate::T.matrix_2x2();
        let zero = Complex64::new(0.0, 0.0);
        let fused = [[zero, t[1][1]], [t[0][0], zero]];

        let mut c = Circuit::new(1, 0);
        c.add_gate(Gate::H, &[0]);
        c.add_gate(Gate::Fused(Box::new(fused)), &[0]);
        let b = run_sparse(&c);

        let h = 1.0 / 2.0_f64.sqrt();
        assert_eq!(b.state.len(), 2);
        assert!((b.state[&0] - t[1][1] * h).norm() < EPS);
        assert!((b.state[&1] - t[0][0] * h).norm() < EPS);
    }

    #[test]
    fn test_antidiagonal_subunit_prunes() {
        let mut b = SparseBackend::new(42);
        b.init(1, 0).unwrap();
        b.state.insert(1, Complex64::new(1.5e-8, 0.0));
        let zero = Complex64::new(0.0, 0.0);
        let half = Complex64::new(0.5, 0.0);
        b.apply_1q_matrix(0, &[[zero, half], [half, zero]]).unwrap();
        assert_eq!(b.state.len(), 1);
        assert!(b.state.contains_key(&1));
    }

    #[test]
    fn test_ghz_4_sparse() {
        let mut c = Circuit::new(4, 0);
        c.add_gate(Gate::H, &[0]);
        for i in 0..3 {
            c.add_gate(Gate::Cx, &[i, i + 1]);
        }
        let b = run_sparse(&c);
        assert_eq!(b.state.len(), 2);
        assert!((b.state[&0].norm_sqr() - 0.5).abs() < EPS);
        assert!((b.state[&15].norm_sqr() - 0.5).abs() < EPS);
    }

    #[test]
    fn test_cu_phase_applies_phase() {
        let mut c = Circuit::new(2, 0);
        c.add_gate(Gate::X, &[0]);
        c.add_gate(Gate::X, &[1]);
        c.add_gate(Gate::cphase(std::f64::consts::FRAC_PI_4), &[0, 1]);
        let b = run_sparse(&c);
        assert_eq!(b.state.len(), 1);
        let expected = Complex64::from_polar(1.0, std::f64::consts::FRAC_PI_4);
        assert!((b.state[&3] - expected).norm() < EPS);
    }

    #[test]
    fn test_cu_phase_no_action_control_zero() {
        let mut c = Circuit::new(2, 0);
        c.add_gate(Gate::H, &[1]);
        c.add_gate(Gate::cphase(1.0), &[0, 1]);
        let b = run_sparse(&c);
        let h = 1.0 / 2.0_f64.sqrt();
        assert!((b.state[&0].re - h).abs() < EPS);
        assert!((b.state[&2].re - h).abs() < EPS);
        assert!(!b.state.contains_key(&1));
        assert!(!b.state.contains_key(&3));
    }

    #[test]
    fn test_cu_phase_matches_cz() {
        let mut c1 = Circuit::new(2, 0);
        c1.add_gate(Gate::H, &[0]);
        c1.add_gate(Gate::H, &[1]);
        c1.add_gate(Gate::cphase(std::f64::consts::PI), &[0, 1]);

        let mut c2 = Circuit::new(2, 0);
        c2.add_gate(Gate::H, &[0]);
        c2.add_gate(Gate::H, &[1]);
        c2.add_gate(Gate::Cz, &[0, 1]);

        let b1 = run_sparse(&c1);
        let b2 = run_sparse(&c2);

        for (&idx, &amp1) in &b1.state {
            let amp2 = b2
                .state
                .get(&idx)
                .copied()
                .unwrap_or(Complex64::new(0.0, 0.0));
            assert!((amp1 - amp2).norm() < EPS, "mismatch at idx {idx}");
        }
    }

    #[test]
    fn test_diagonal_1q_in_place() {
        let mut c = Circuit::new(1, 0);
        c.add_gate(Gate::H, &[0]);
        c.add_gate(Gate::P(1.234), &[0]);
        let b = run_sparse(&c);
        assert_eq!(b.state.len(), 2);
        let h = 1.0 / 2.0_f64.sqrt();
        assert!((b.state[&0].re - h).abs() < EPS);
        assert!((b.state[&1] - Complex64::from_polar(h, 1.234)).norm() < EPS);
    }

    #[test]
    fn test_batch_rzz_matches_individual() {
        use crate::gates::BatchRzzData;

        let edges = vec![(0usize, 1usize, 0.7f64), (1, 2, 1.3)];

        let mut c1 = Circuit::new(3, 0);
        for q in 0..3 {
            c1.add_gate(Gate::H, &[q]);
        }
        for &(q0, q1, theta) in &edges {
            c1.add_gate(Gate::Rzz(theta), &[q0, q1]);
        }
        let b1 = run_sparse(&c1);

        let mut c2 = Circuit::new(3, 0);
        for q in 0..3 {
            c2.add_gate(Gate::H, &[q]);
        }
        c2.add_gate(Gate::BatchRzz(Box::new(BatchRzzData { edges })), &[0, 1, 2]);
        let b2 = run_sparse(&c2);

        for (&idx, &amp1) in &b1.state {
            let amp2 = b2
                .state
                .get(&idx)
                .copied()
                .unwrap_or(Complex64::new(0.0, 0.0));
            assert!((amp1 - amp2).norm() < EPS, "mismatch at idx {idx}");
        }
    }

    #[test]
    fn test_diagonal_batch_matches_individual() {
        use crate::gates::{DiagEntry, DiagonalBatchData};

        let mut c1 = Circuit::new(3, 0);
        for q in 0..3 {
            c1.add_gate(Gate::H, &[q]);
        }
        c1.add_gate(Gate::S, &[0]);
        c1.add_gate(Gate::Cz, &[0, 1]);
        c1.add_gate(Gate::Rzz(0.9), &[1, 2]);
        let b1 = run_sparse(&c1);

        let s_mat = Gate::S.matrix_2x2();
        let mut c2 = Circuit::new(3, 0);
        for q in 0..3 {
            c2.add_gate(Gate::H, &[q]);
        }
        c2.add_gate(
            Gate::DiagonalBatch(Box::new(DiagonalBatchData {
                entries: vec![
                    DiagEntry::Phase1q {
                        qubit: 0,
                        d0: s_mat[0][0],
                        d1: s_mat[1][1],
                    },
                    DiagEntry::Phase2q {
                        q0: 0,
                        q1: 1,
                        phase: Complex64::new(-1.0, 0.0),
                    },
                    DiagEntry::Parity2q {
                        q0: 1,
                        q1: 2,
                        same: Complex64::from_polar(1.0, -0.45),
                        diff: Complex64::from_polar(1.0, 0.45),
                    },
                ],
            })),
            &[0, 1, 2],
        );
        let b2 = run_sparse(&c2);

        for (&idx, &amp1) in &b1.state {
            let amp2 = b2
                .state
                .get(&idx)
                .copied()
                .unwrap_or(Complex64::new(0.0, 0.0));
            assert!((amp1 - amp2).norm() < EPS, "mismatch at idx {idx}");
        }
    }

    #[test]
    fn test_batch_phase_matches_individual() {
        use crate::gates::BatchPhaseData;
        use smallvec::smallvec;

        let phase1 = Complex64::from_polar(1.0, 0.5);
        let phase2 = Complex64::from_polar(1.0, 1.2);

        let mut c1 = Circuit::new(3, 0);
        c1.add_gate(Gate::H, &[0]);
        c1.add_gate(Gate::H, &[1]);
        c1.add_gate(Gate::H, &[2]);
        c1.add_gate(Gate::cphase(0.5), &[0, 1]);
        c1.add_gate(Gate::cphase(1.2), &[0, 2]);
        let p1 = run_sparse_probs(&c1);

        let mut c2 = Circuit::new(3, 0);
        c2.add_gate(Gate::H, &[0]);
        c2.add_gate(Gate::H, &[1]);
        c2.add_gate(Gate::H, &[2]);
        c2.add_gate(
            Gate::BatchPhase(Box::new(BatchPhaseData {
                phases: smallvec![(1, phase1), (2, phase2)],
            })),
            &[0, 1, 2],
        );
        let p2 = run_sparse_probs(&c2);

        for (a, b) in p1.iter().zip(p2.iter()) {
            assert!((a - b).abs() < EPS, "probs mismatch: {a} vs {b}");
        }
    }
}