prism-q 0.35.0

//! Backend selection and execution planning.
//!
//! Resolves a [`BackendKind`] against a circuit into a buildable plan,
//! including the CPU-vs-GPU accel verdict and the temporal-Clifford split.

use crate::backend::density_matrix::DensityMatrixBackend;
use crate::backend::mps::MpsBackend;
use crate::backend::product::ProductStateBackend;
use crate::backend::sparse::{MAX_SPARSE_INDEX_QUBITS, SparseBackend};
use crate::backend::stabilizer::StabilizerBackend;
use crate::backend::statevector::StatevectorBackend;
use crate::backend::tensornetwork::TensorNetworkBackend;
use crate::backend::{
    Backend, DM_QUBIT_CAP_ENV, check_state_allocation, max_density_matrix_qubits,
    max_statevector_qubits,
};
use crate::circuit::{Circuit, Instruction};
use crate::error::{PrismError, Result};
use crate::gates::Gate;
use crate::sim::unified_pauli::SpdTruncation;

#[cfg(any(feature = "gpu", feature = "distributed"))]
use std::sync::Arc;

#[cfg(feature = "gpu")]
use crate::gpu::GpuContext;

#[cfg(feature = "distributed")]
use crate::backend::distributed_statevector::DistributedStatevectorBackend;
#[cfg(feature = "distributed")]
use crate::distributed::DistributedContext;

use super::StartState;
use super::metadata::ResolvedBackend;
use super::{RunOutcome, try_backend_probabilities};

pub(super) const AUTO_MPS_BOND_DIM: usize = 256;

/// Fixed cost of one gate at one cut in [`mps_apply_cost`] units, standing for
/// the per-gate work (tensor reshapes, gauge moves, allocation) a small bond does
/// not amortize. Fit on 85 circuits at 32 to 64 qubits: about 4 us per gate-cut
/// against about 15 ns per unit of bond cubed.
const MPS_GATE_CUT_OVERHEAD: f64 = 256.0;

/// Estimate the work an MPS capped at `max_bond` spends applying `circuit`, in
/// units of one bond dimension cubed.
///
/// Each gate on two or more qubits multiplies the bond at every cut it spans by
/// its operator Schmidt rank (2 for controlled and `Rzz` gates, 4 otherwise),
/// capped at `min(2^min(left, right), max_bond)`, and costs the cube of each
/// resulting bond, the SVD it triggers there, plus [`MPS_GATE_CUT_OVERHEAD`].
/// Gates on one pair with nothing between them on either qubit or any spanned
/// cut form one block, whose rank stays at most 4, as fusion delivers them
/// merged. At 15 ns per unit, MPS time measured 0.08 to 3.8 times the estimate
/// (median 1.04) across the fitted circuits; the bond is an upper bound, and on
/// shallow circuits whose gates entangle less than their rank allows the
/// estimate runs far high.
pub(super) fn mps_apply_cost(circuit: &Circuit, max_bond: usize) -> f64 {
    let n = circuit.num_qubits;
    if n < 2 {
        return 0.0;
    }
    let limit: Vec<f64> = (0..n - 1)
        .map(|cut| {
            2f64.powi((cut + 1).min(n - cut - 1) as i32)
                .min(max_bond as f64)
        })
        .collect();
    let mut bond = vec![1.0f64; n - 1];
    let mut base = vec![1.0f64; n - 1];
    let mut block_rank = vec![1.0f64; n - 1];
    let mut qubit_block = vec![0usize; n];
    let mut cut_block = vec![0usize; n - 1];
    let mut next_block = 1usize;
    let mut cost = 0.0;
    for instruction in &circuit.instructions {
        let Instruction::Gate { gate, targets } = instruction else {
            continue;
        };
        if targets.len() < 2 {
            continue;
        }
        let lo = *targets.iter().min().expect("two or more targets");
        let hi = *targets.iter().max().expect("two or more targets");
        let rank = match gate {
            Gate::Cx | Gate::Cz | Gate::Rzz(_) | Gate::Cu(_) | Gate::Mcu(_) => 2.0,
            _ => 4.0,
        };
        let block = qubit_block[lo];
        let continues = targets.len() == 2
            && block != 0
            && qubit_block[hi] == block
            && cut_block[lo..hi].iter().all(|&b| b == block);
        let block = if continues {
            block
        } else {
            next_block += 1;
            next_block - 1
        };
        for cut in lo..hi {
            if continues {
                block_rank[cut] = (block_rank[cut] * rank).min(4.0);
            } else {
                base[cut] = bond[cut];
                block_rank[cut] = rank;
            }
            bond[cut] = (base[cut] * block_rank[cut]).min(limit[cut]);
            cut_block[cut] = block;
            cost += MPS_GATE_CUT_OVERHEAD + bond[cut].powi(3);
        }
        for &q in targets.iter() {
            qubit_block[q] = if targets.len() == 2 { block } else { 0 };
        }
    }
    cost
}

pub(super) const MAX_AUTO_T_COUNT_SHOTS: usize = 40;

pub(super) const MAX_STABILIZER_RANK_QUBITS: usize = 25;

pub(super) const MIN_QUBITS_FOR_SPD_AUTO: usize = 12;

/// Statevector amplitude updates one SPD term step is priced at when `Auto`
/// answers marginals. A term step measured 100 to 300 ns against about 1.4 ns
/// per amplitude update; 64 leans toward the Pauli route so a circuit it wins
/// stays on it.
pub(super) const AUTO_SPD_AMPLITUDES_PER_TERM: usize = 64;

pub(super) const MIN_FACTORED_STABILIZER_QUBITS: usize = 24;

pub(super) const MIN_BLOCK_FOR_FACTORED_STAB: usize = 2;

/// Entry bound under which `Auto` takes the sparse map over the MPS past the
/// statevector cap. A sparse gate costs about 8 ns per entry and an MPS gate at
/// bond 2 about 2 us, so on a circuit that barely entangles the two tie near 256.
const AUTO_SPARSE_MAX_SUPPORT: usize = 128;

#[inline]
pub(super) fn stabilizer_rank_budget(num_qubits: usize) -> usize {
    let log2n = if num_qubits >= 2 {
        (num_qubits as f64).log2().ceil() as usize * 2
    } else {
        0
    };
    num_qubits.saturating_sub(log2n)
}

/// SPD term steps `Auto` spends on a marginals query before handing it to the
/// statevector: the dense run's `instructions * 2^n` amplitude updates, priced
/// at [`AUTO_SPD_AMPLITUDES_PER_TERM`] each.
#[inline]
pub(super) fn auto_spd_work_budget(circuit: &Circuit) -> usize {
    circuit.instructions.len().saturating_mul(
        1usize
            .checked_shl(circuit.num_qubits as u32)
            .unwrap_or(usize::MAX),
    ) / AUTO_SPD_AMPLITUDES_PER_TERM
}

/// Backend selection for a simulation run.
///
/// `Auto` resolves per call from circuit shape. Two routes run before the
/// family tree: circuits that decompose into independent blocks run per block
/// (Clifford-only circuits at 24 qubits and above with a block of 2 or more use
/// FactoredStabilizer), and shot paths sample Clifford+T circuits whose T
/// count fits the size-derived stabilizer-rank budget through the StabilizerRank
/// expansion, to 40 T gates. A probability
/// query never takes that expansion under `Auto`, since it is slower than the
/// statevector at every width both can hold. The pruned expansion is reachable
/// only through [`run_stabilizer_rank_approx`], never from `Auto`; marginal queries
/// on Clifford+T circuits at 12 qubits and above answer via Sparse Pauli
/// Dynamics, which hands a circuit that fits the statevector back to the tree
/// once its Pauli sum costs more than the dense run would. The remaining tree:
///
/// 1. No entangling gates        → ProductState (O(n))
/// 2. All Clifford gates         → Stabilizer (O(n²))
/// 3. Above the statevector memory cap:
///    a. Sparse-friendly         → Sparse (O(k) where k = non-zero amplitudes)
///    b. Otherwise               → MPS (bounded bond dimension)
/// 4. Partial independence       → Factored (per-group dense sub-states)
/// 5. Otherwise                  → Statevector (exact, general-purpose)
///
/// [`run_stabilizer_rank_approx`]: crate::run_stabilizer_rank_approx
#[derive(Debug, Clone)]
#[non_exhaustive]
pub enum BackendKind {
    Auto,
    Statevector,
    /// Tableau simulation for Clifford-only circuits.
    Stabilizer,
    /// Sparse state vector holding only nonzero amplitudes.
    Sparse,
    /// Matrix Product State simulation with a bounded bond dimension.
    Mps {
        max_bond_dim: usize,
    },
    /// Per-qubit product state for circuits without entangling gates.
    ProductState,
    /// Deferred-contraction tensor network for low-treewidth circuits.
    TensorNetwork,
    /// The tensor network with bond truncation at the peak cap.
    ///
    /// An intermediate over the tensor-network peak cap is factored and the
    /// new bond kept only as far as discarding at most `tolerance` of that
    /// cut's squared weight allows; index slicing covers whatever stays above
    /// the cap. `tolerance` must be finite and above zero, since zero is
    /// [`BackendKind::TensorNetwork`]. The run is approximate, so
    /// [`Simulate::require_exact`] rejects it.
    ///
    /// [`Simulate::require_exact`]: crate::sim::Simulate::require_exact
    TensorNetworkBounded {
        tolerance: f64,
    },
    /// Dynamic split-state simulation for sparse-entanglement circuits.
    Factored,
    /// Clifford+T decomposition into weighted stabilizer branches.
    StabilizerRank,
    /// Independent Clifford blocks, each on its own tableau.
    FactoredStabilizer,
    /// Exact density-matrix backend for mixed-state evolution.
    ///
    /// Explicit-dispatch only: [`BackendKind::Auto`] never selects it. Stores
    /// `4^n` `Complex64` amplitudes, so the qubit ceiling is roughly half the
    /// statevector cap (14 on a 16 GiB host); `PRISM_MAX_DM_QUBITS` moves it
    /// within that bound. Fused payloads are accepted, with every fusion floor
    /// gated on the `2n`-qubit buffer the backend sweeps rather than on the
    /// circuit width.
    ///
    /// Under an attached noise model this is the exact route: the mixture is
    /// evolved once and every terminal reads it, so shot counts carry sampling
    /// noise only and observables carry none.
    DensityMatrix,
    /// Stochastic Pauli propagation (SPP); serves marginal and observable
    /// queries only.
    StochasticPauli {
        num_samples: usize,
    },
    /// Deterministic sparse Pauli dynamics (SPD); serves marginal and
    /// observable queries only. [`SpdTruncation`] picks between the threshold
    /// policy and the budget one, which holds a fixed term count and so cannot
    /// die at the internal ceiling on growth the budget already caps.
    DeterministicPauli {
        truncation: SpdTruncation,
    },
    /// Heisenberg Pauli propagation through a noise model; serves expectation
    /// values and observable expectations only.
    ///
    /// Explicit-dispatch only: [`BackendKind::Auto`] never selects it. Noise
    /// enters as the channel's action on the Pauli basis, which shrinks
    /// coefficients while the circuit's rotations grow the term count, so the
    /// weighted sum stays small exactly when noise outpaces the branching
    /// rotations. Terms below `epsilon` are dropped once the sum exceeds
    /// `max_terms` (0 disables truncation and makes the run exact).
    ///
    /// A channel with no Pauli-basis form (custom Kraus, two-qubit Kraus,
    /// readout error) is rejected naming the density matrix.
    PauliPath {
        epsilon: f64,
        max_terms: usize,
    },
    /// Automatic backend selection with GPU acceleration opted in.
    ///
    /// Makes the same shape-based routing decisions as [`BackendKind::Auto`],
    /// but when the selected family (per sub-block, after subsystem
    /// decomposition) has a device capability row and the block clears the
    /// qubit-count crossover with VRAM to spare, that block runs on the
    /// supplied context. Every other choice, and every block that fails the
    /// crossover or VRAM check, runs on the identical CPU path Auto would
    /// take. Device paths resolve soft: an allocation that fails after the
    /// VRAM check degrades to the host, so a missing, unfit, or racing device
    /// stays on CPU rather than erroring.
    ///
    /// Applies to single runs, terminal shot and counts sampling, expectation
    /// values, temporal-Clifford tails, and non-Pauli noisy trajectories.
    #[cfg(feature = "gpu")]
    AutoGpu {
        context: Arc<GpuContext>,
    },
    /// Statevector backed by a CUDA GPU execution context.
    ///
    /// Circuits (or decomposed sub-blocks) with fewer than
    /// [`crate::gpu::min_qubits()`] qubits (tunable via
    /// `PRISM_GPU_MIN_QUBITS`, default [`crate::gpu::MIN_QUBITS_DEFAULT`])
    /// fall back to the host statevector path, since small states do not
    /// survive PCIe and launch-latency overhead.
    ///
    /// Compose with `simulate(...).backend(...).seed(...).run()` to get fusion
    /// plus independent-subsystem decomposition; each sub-block is evaluated
    /// against the crossover independently.
    #[cfg(feature = "gpu")]
    StatevectorGpu {
        context: Arc<GpuContext>,
    },
    /// Density matrix held in device memory on the supplied context.
    ///
    /// Explicit-dispatch only: neither [`BackendKind::Auto`] nor
    /// [`BackendKind::AutoGpu`] selects it. There is no crossover and no host
    /// fallback: every run allocates the `4^n` mixture on the device, after a
    /// budget check against the free VRAM that errors before allocating. An
    /// 11 GiB card holds 13 qubits (1 GiB at 13, 4 GiB at 14 plus scratch).
    /// Every channel, measurement, and readout sweep runs as a kernel; the
    /// noisy terminals answer from the exact mixture as
    /// [`BackendKind::DensityMatrix`] does.
    #[cfg(feature = "gpu")]
    DensityMatrixGpu {
        context: Arc<GpuContext>,
    },
    /// Stabilizer backend backed by a CUDA GPU tableau.
    ///
    /// Circuits (or decomposed sub-blocks) with fewer than
    /// [`crate::gpu::stabilizer_min_qubits()`] qubits (tunable via
    /// `PRISM_STABILIZER_GPU_MIN_QUBITS`, default
    /// [`crate::gpu::STABILIZER_MIN_QUBITS_DEFAULT`]) fall back to the CPU
    /// stabilizer path. Measurement and reset stay on device; probabilities and
    /// exports read back to the CPU algorithms.
    ///
    /// Compose with `simulate(...).backend(...).seed(...).run()` to pick up
    /// independent-subsystem decomposition; non-Clifford circuits are rejected
    /// at dispatch time with the same error shape as [`BackendKind::Stabilizer`].
    #[cfg(feature = "gpu")]
    StabilizerGpu {
        context: Arc<GpuContext>,
    },
    /// Exact state vector distributed across `2^p` ranks via a
    /// [`DistributedContext`]. The low `n - p` qubits are simulated locally with
    /// the standard SIMD kernels; the top `p` qubits select the rank.
    ///
    /// Results are independent of the rank count. With a single rank the path is
    /// identical to [`BackendKind::Statevector`].
    #[cfg(feature = "distributed")]
    StatevectorDistributed {
        context: Arc<DistributedContext>,
    },
}

impl BackendKind {
    /// Whether this kind uses automatic shape-based routing. True for
    /// [`BackendKind::Auto`] and its GPU-accelerated sibling
    /// [`BackendKind::AutoGpu`], which make identical routing decisions and
    /// differ only in whether a cleared statevector or stabilizer block runs on
    /// the device.
    #[inline]
    pub(crate) fn is_auto(&self) -> bool {
        match self {
            BackendKind::Auto => true,
            #[cfg(feature = "gpu")]
            BackendKind::AutoGpu { .. } => true,
            _ => false,
        }
    }

    /// False for the engines without a per-shot pure state: stabilizer rank,
    /// Pauli propagation, density matrix. The density matrix still serves a
    /// noise model, by evolving the mixture once and sampling the exact
    /// distribution, so its noisy terminals route around the trajectory engine
    /// rather than through it.
    pub fn supports_noisy_per_shot(&self) -> bool {
        !matches!(
            self,
            BackendKind::StabilizerRank
                | BackendKind::StochasticPauli { .. }
                | BackendKind::DeterministicPauli { .. }
                | BackendKind::PauliPath { .. }
        ) && !self.is_density_matrix()
    }

    /// True for the two kinds that hold the exact mixture, host or device.
    pub(crate) fn is_density_matrix(&self) -> bool {
        match self {
            BackendKind::DensityMatrix => true,
            #[cfg(feature = "gpu")]
            BackendKind::DensityMatrixGpu { .. } => true,
            _ => false,
        }
    }

    /// True for kinds that can run non-Pauli channels (damping, thermal
    /// relaxation, custom Kraus): the trajectory engine everywhere except the
    /// density matrix, which applies the channel to the mixture instead.
    pub fn supports_general_noise(&self) -> bool {
        match self {
            BackendKind::Auto
            | BackendKind::Statevector
            | BackendKind::Sparse
            | BackendKind::Mps { .. }
            | BackendKind::ProductState
            | BackendKind::Factored
            | BackendKind::TensorNetwork
            | BackendKind::TensorNetworkBounded { .. }
            | BackendKind::DensityMatrix => true,
            #[cfg(feature = "gpu")]
            BackendKind::AutoGpu { .. }
            | BackendKind::StatevectorGpu { .. }
            | BackendKind::DensityMatrixGpu { .. } => true,
            _ => false,
        }
    }

    pub(crate) fn is_stabilizer_family(&self) -> bool {
        matches!(
            self,
            BackendKind::Stabilizer | BackendKind::FactoredStabilizer
        ) || {
            #[cfg(feature = "gpu")]
            {
                matches!(self, BackendKind::StabilizerGpu { .. })
            }
            #[cfg(not(feature = "gpu"))]
            {
                false
            }
        }
    }

    pub(crate) fn general_noise_backend_names() -> &'static str {
        #[cfg(feature = "gpu")]
        {
            "Auto, Statevector, StatevectorGpu, Sparse, Mps, ProductState, Factored, TensorNetwork, DensityMatrix, or DensityMatrixGpu"
        }
        #[cfg(not(feature = "gpu"))]
        {
            "Auto, Statevector, Sparse, Mps, ProductState, Factored, TensorNetwork, or DensityMatrix"
        }
    }
}

pub(super) fn validate_explicit_backend(kind: &BackendKind, circuit: &Circuit) -> Result<()> {
    if kind.is_stabilizer_family() && !circuit.is_clifford_only() {
        return Err(PrismError::IncompatibleBackend {
            backend: "stabilizer".into(),
            reason: "circuit contains non-Clifford gates".into(),
        });
    }
    match kind {
        BackendKind::ProductState if circuit.has_entangling_gates() => {
            return Err(PrismError::IncompatibleBackend {
                backend: "productstate".into(),
                reason: "circuit contains entangling gates".into(),
            });
        }
        BackendKind::TensorNetworkBounded { tolerance }
            if !(tolerance.is_finite() && *tolerance > 0.0) =>
        {
            return Err(PrismError::InvalidParameter {
                message: format!(
                    "tensor-network tolerance must be a finite fraction above zero, got {tolerance}; \
                     TensorNetwork is the exact contraction"
                ),
            });
        }
        BackendKind::StabilizerRank if !circuit.has_t_gates() => {
            return Err(PrismError::IncompatibleBackend {
                backend: "stabilizer_rank".into(),
                reason: "circuit has no T gates; use Stabilizer instead".into(),
            });
        }
        BackendKind::DensityMatrix => {
            check_state_allocation(
                "density_matrix",
                circuit.num_qubits,
                max_density_matrix_qubits(),
                DM_QUBIT_CAP_ENV,
            )?;
        }
        _ => {}
    }
    Ok(())
}

/// Simulator family, independent of how it was selected (auto routing or an
/// explicit kind) and of where it executes (host or device).
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub(super) enum Family {
    ProductState,
    Stabilizer,
    Sparse,
    Mps,
    Factored,
    FactoredStabilizer,
    TensorNetwork,
    Statevector,
    DensityMatrix,
}

fn select_auto_backend_choice(circuit: &Circuit, has_partial_independence: bool) -> Family {
    if !circuit.has_entangling_gates() {
        Family::ProductState
    } else if circuit.is_clifford_only() {
        Family::Stabilizer
    } else if circuit.num_qubits > max_statevector_qubits() {
        if circuit.num_qubits <= MAX_SPARSE_INDEX_QUBITS
            && (circuit.is_sparse_friendly()
                || sparse_support_fits(circuit, AUTO_SPARSE_MAX_SUPPORT))
        {
            Family::Sparse
        } else {
            Family::Mps
        }
    } else if has_partial_independence {
        Family::Factored
    } else {
        Family::Statevector
    }
}

/// Run the circuit's gates on a sparse map and report whether it stays within
/// `cap` entries throughout. Declines resets, saves, and anything after a
/// measurement, whose branch could reach a different support.
fn sparse_support_fits(circuit: &Circuit, cap: usize) -> bool {
    if !circuit.has_terminal_measurements_only() {
        return false;
    }
    let mut probe = SparseBackend::new(0);
    if probe.init(circuit.num_qubits, 0).is_err() {
        return false;
    }
    circuit.instructions.iter().all(|inst| match inst {
        Instruction::Gate { .. } => probe.apply(inst).is_ok() && probe.entry_count() <= cap,
        Instruction::Measure { .. } | Instruction::Barrier { .. } => true,
        _ => false,
    })
}

pub(super) fn auto_selects_cpu_statevector(
    circuit: &Circuit,
    has_partial_independence: bool,
) -> bool {
    matches!(
        select_auto_backend_choice(circuit, has_partial_independence),
        Family::Statevector
    )
}

/// Execution target for a resolved family: host, or a device context with a
/// resolved failure mode. `soft` builds fall back to the host if device
/// allocation fails at `init`; hard builds surface the error.
#[derive(Clone)]
pub(super) enum Accel {
    Cpu,
    #[cfg(feature = "gpu")]
    Gpu {
        context: Arc<GpuContext>,
        soft: bool,
    },
}

/// Device eligibility for one family: the qubit crossover and the VRAM-fit
/// predicate, owned together. Families without device kernels are `None` rows
/// in [`gpu_capability`]; adding a device path for a family means adding its
/// row there, not touching dispatch.
#[cfg(feature = "gpu")]
struct GpuCapability {
    min_qubits: fn() -> usize,
    fits: fn(&Arc<GpuContext>, usize) -> bool,
}

/// Families with a `gpu_capability` row. Keep in sync with the match below;
/// a new device family needs an entry in both.
#[cfg(feature = "gpu")]
const GPU_CAPABLE_FAMILIES: [Family; 3] = [
    Family::Statevector,
    Family::Stabilizer,
    Family::DensityMatrix,
];

#[cfg(feature = "gpu")]
fn gpu_capability(family: Family) -> Option<GpuCapability> {
    match family {
        Family::Statevector => Some(GpuCapability {
            min_qubits: crate::gpu::min_qubits,
            fits: |ctx, n| ctx.fits_statevector_with_scratch(n).unwrap_or(false),
        }),
        Family::Stabilizer => Some(GpuCapability {
            min_qubits: crate::gpu::stabilizer_min_qubits,
            fits: |ctx, n| ctx.fits_tableau(n).unwrap_or(false),
        }),
        // Explicit-only: no crossover, and the `Auto` tree never selects the
        // family, so the soft fit gate is never consulted. The hard path budgets
        // the `4^n` buffer at `init`.
        Family::DensityMatrix => Some(GpuCapability {
            min_qubits: || 0,
            fits: |ctx, n| ctx.fits_statevector_with_scratch(2 * n).unwrap_or(false),
        }),
        Family::ProductState
        | Family::Sparse
        | Family::Mps
        | Family::Factored
        | Family::FactoredStabilizer
        | Family::TensorNetwork => None,
    }
}

/// The one decision point for CPU vs GPU execution of a family.
///
/// `AutoGpu` requires the family's capability row, the qubit crossover, and the
/// VRAM-fit gate, and resolves soft so an allocation race still degrades to the
/// host. An explicit GPU kind applies only the crossover for its own family and
/// resolves hard: the user asked for the device, so an unfit device errors
/// loudly. Every other kind, and every family without a capability row,
/// resolves to the host.
///
/// The verdict (including the VRAM query) is intended to be taken once per
/// user-level call and reused across shots; soft-mode init fallback covers
/// VRAM shrinking after the fact.
#[cfg(feature = "gpu")]
pub(super) fn accel_for(kind: &BackendKind, family: Family, num_qubits: usize) -> Accel {
    let Some((context, soft)) = gpu_request(kind, family) else {
        return Accel::Cpu;
    };
    let Some(cap) = gpu_capability(family) else {
        return Accel::Cpu;
    };
    if num_qubits < (cap.min_qubits)() {
        return Accel::Cpu;
    }
    if soft && !(cap.fits)(context, num_qubits) {
        return Accel::Cpu;
    }
    Accel::Gpu {
        context: context.clone(),
        soft,
    }
}

/// Which (kind, family) pairs request device execution, and whether the
/// request resolves soft. Shared by [`accel_for`] and [`may_resolve_to_gpu`]
/// so the kind-to-family matching lives in one place.
#[cfg(feature = "gpu")]
fn gpu_request(kind: &BackendKind, family: Family) -> Option<(&Arc<GpuContext>, bool)> {
    match (kind, family) {
        (BackendKind::AutoGpu { context }, _) => Some((context, true)),
        (BackendKind::StatevectorGpu { context }, Family::Statevector) => Some((context, false)),
        (BackendKind::StabilizerGpu { context }, Family::Stabilizer) => Some((context, false)),
        (BackendKind::DensityMatrixGpu { context }, Family::DensityMatrix) => {
            Some((context, false))
        }
        _ => None,
    }
}

/// Whether this kind could resolve any family to the device at the given
/// width. Over-approximates [`accel_for`]: it applies only the qubit
/// crossover and skips the VRAM-fit gate, so the verdict cannot flip between
/// a scheduling decision and the later per-block resolution. Multi-block
/// drivers consult it before spreading work across threads; GPU backends
/// share one CUDA stream per context and must not execute concurrently.
#[cfg(feature = "gpu")]
pub(super) fn may_resolve_to_gpu(kind: &BackendKind, num_qubits: usize) -> bool {
    GPU_CAPABLE_FAMILIES.into_iter().any(|family| {
        gpu_request(kind, family).is_some()
            && gpu_capability(family).is_some_and(|cap| num_qubits >= (cap.min_qubits)())
    })
}

#[cfg(not(feature = "gpu"))]
pub(super) fn accel_for(_kind: &BackendKind, _family: Family, _num_qubits: usize) -> Accel {
    Accel::Cpu
}

pub(super) fn build_statevector(accel: &Accel, seed: u64) -> StatevectorBackend {
    match accel {
        Accel::Cpu => StatevectorBackend::new(seed),
        #[cfg(feature = "gpu")]
        Accel::Gpu {
            context,
            soft: true,
        } => StatevectorBackend::new(seed).with_gpu_auto(context.clone()),
        #[cfg(feature = "gpu")]
        Accel::Gpu {
            context,
            soft: false,
        } => StatevectorBackend::new(seed).with_gpu(context.clone()),
    }
}

pub(super) fn build_density_matrix(accel: &Accel, seed: u64) -> DensityMatrixBackend {
    match accel {
        Accel::Cpu => DensityMatrixBackend::new(seed),
        #[cfg(feature = "gpu")]
        Accel::Gpu { context, .. } => DensityMatrixBackend::new(seed).with_gpu(context.clone()),
    }
}

/// True when no instruction can force destabilizer materialization: gates,
/// satisfied conditionals, and barriers only. Measurements, resets, and
/// guarded regions (whose bodies can hold either) all disqualify.
fn stabilizer_can_run_lazy(circuit: &Circuit) -> bool {
    circuit.instructions.iter().all(|inst| {
        matches!(
            inst,
            Instruction::Gate { .. }
                | Instruction::Conditional { .. }
                | Instruction::Barrier { .. }
        )
    })
}

fn build_stabilizer(accel: &Accel, lazy: bool, seed: u64) -> StabilizerBackend {
    // Lazy destabilizers: gates touch half the tableau. The plan sets `lazy`
    // only for circuits that never measure or reset, so the Gaussian
    // elimination that materialization costs is never paid; on a measuring
    // circuit the reconstruction can exceed the gate savings (the GHZ
    // measure-all shape loses outright, since its SGI gate cost is already
    // near zero). GPU init overrides the flag; the soft-fallback host path
    // keeps it.
    let stab = if lazy {
        StabilizerBackend::new_lazy(seed)
    } else {
        StabilizerBackend::new(seed)
    };
    match accel {
        Accel::Cpu => stab,
        #[cfg(feature = "gpu")]
        Accel::Gpu {
            context,
            soft: true,
        } => stab.with_gpu_auto(context.clone()),
        #[cfg(feature = "gpu")]
        Accel::Gpu {
            context,
            soft: false,
        } => stab.with_gpu(context.clone()),
    }
}

/// Resolved family plus acceleration for one user-level call. `build` is
/// cheap enough to call once per shot or per trajectory: a constructor, a
/// `Box::new`, and at most one `Arc` clone. All circuit analysis and every
/// driver VRAM query happen earlier, in [`resolve`].
#[derive(Clone)]
pub(super) enum BackendPlan {
    ProductState,
    Sparse,
    TensorNetwork {
        tolerance: Option<f64>,
    },
    Factored,
    FactoredStabilizer,
    Mps {
        max_bond_dim: usize,
    },
    Stabilizer {
        accel: Accel,
        lazy: bool,
    },
    Statevector {
        accel: Accel,
    },
    DensityMatrix {
        accel: Accel,
    },
    #[cfg(feature = "distributed")]
    Distributed(Arc<DistributedContext>),
}

/// True when two accelerator targets are the same device with the same failure
/// mode. A soft and a hard build of one context differ, because a device
/// allocation failure sends them to different places.
fn same_accel(a: &Accel, b: &Accel) -> bool {
    match (a, b) {
        (Accel::Cpu, Accel::Cpu) => true,
        #[cfg(feature = "gpu")]
        (
            Accel::Gpu {
                context: cx,
                soft: sx,
            },
            Accel::Gpu {
                context: cy,
                soft: sy,
            },
        ) => Arc::ptr_eq(cx, cy) && sx == sy,
        #[cfg(feature = "gpu")]
        _ => false,
    }
}

/// True when two plans build the same engine with the same settings, so one
/// backend can serve both.
///
/// Written out rather than derived because the distributed variant holds a
/// context that has no equality, and sharing a rank's backend across circuits
/// is not something this answers yes to.
pub(super) fn same_plan(a: &BackendPlan, b: &BackendPlan) -> bool {
    match (a, b) {
        (BackendPlan::ProductState, BackendPlan::ProductState)
        | (BackendPlan::Sparse, BackendPlan::Sparse)
        | (BackendPlan::Factored, BackendPlan::Factored)
        | (BackendPlan::FactoredStabilizer, BackendPlan::FactoredStabilizer) => true,
        (
            BackendPlan::TensorNetwork { tolerance: x },
            BackendPlan::TensorNetwork { tolerance: y },
        ) => x == y,
        (BackendPlan::Mps { max_bond_dim: x }, BackendPlan::Mps { max_bond_dim: y }) => x == y,
        (
            BackendPlan::Stabilizer {
                accel: ax,
                lazy: lx,
            },
            BackendPlan::Stabilizer {
                accel: ay,
                lazy: ly,
            },
        ) => same_accel(ax, ay) && lx == ly,
        (BackendPlan::Statevector { accel: x }, BackendPlan::Statevector { accel: y })
        | (BackendPlan::DensityMatrix { accel: x }, BackendPlan::DensityMatrix { accel: y }) => {
            same_accel(x, y)
        }
        _ => false,
    }
}

impl BackendPlan {
    /// Engine this plan builds, for a shot request that runs the plan zero
    /// times and so has no backend to read provenance off.
    pub(super) fn resolved(&self) -> ResolvedBackend {
        match self {
            BackendPlan::ProductState => ResolvedBackend::ProductState,
            BackendPlan::Sparse => ResolvedBackend::Sparse,
            BackendPlan::TensorNetwork { .. } => ResolvedBackend::TensorNetwork,
            BackendPlan::Factored => ResolvedBackend::Factored,
            BackendPlan::FactoredStabilizer => ResolvedBackend::FactoredStabilizer,
            BackendPlan::Mps { .. } => ResolvedBackend::Mps,
            BackendPlan::Stabilizer { .. } => ResolvedBackend::Stabilizer,
            BackendPlan::Statevector { .. } => ResolvedBackend::Statevector,
            BackendPlan::DensityMatrix { .. } => ResolvedBackend::DensityMatrix,
            #[cfg(feature = "distributed")]
            BackendPlan::Distributed(_) => ResolvedBackend::Distributed,
        }
    }

    pub(super) fn build(&self, seed: u64) -> Box<dyn Backend + Send> {
        match self {
            BackendPlan::ProductState => Box::new(ProductStateBackend::new(seed)),
            BackendPlan::Sparse => Box::new(SparseBackend::new(seed)),
            BackendPlan::TensorNetwork { tolerance } => Box::new(match *tolerance {
                Some(tolerance) => TensorNetworkBackend::with_tolerance(seed, tolerance),
                None => TensorNetworkBackend::new(seed),
            }),
            BackendPlan::Factored => Box::new(crate::backend::factored::FactoredBackend::new(seed)),
            BackendPlan::FactoredStabilizer => {
                Box::new(crate::backend::factored_stabilizer::FactoredStabilizerBackend::new(seed))
            }
            BackendPlan::Mps { max_bond_dim } => Box::new(MpsBackend::new(seed, *max_bond_dim)),
            BackendPlan::Stabilizer { accel, lazy } => {
                Box::new(build_stabilizer(accel, *lazy, seed))
            }
            BackendPlan::Statevector { accel } => Box::new(build_statevector(accel, seed)),
            BackendPlan::DensityMatrix { accel } => Box::new(build_density_matrix(accel, seed)),
            #[cfg(feature = "distributed")]
            BackendPlan::Distributed(context) => {
                Box::new(DistributedStatevectorBackend::new(context.clone(), seed))
            }
        }
    }

    pub(super) fn is_host_statevector(&self) -> bool {
        matches!(self, BackendPlan::Statevector { accel: Accel::Cpu })
    }

    pub(super) fn is_gpu(&self) -> bool {
        match self {
            BackendPlan::Stabilizer { accel, .. }
            | BackendPlan::Statevector { accel }
            | BackendPlan::DensityMatrix { accel } => !matches!(accel, Accel::Cpu),
            _ => false,
        }
    }

    #[cfg(test)]
    pub(super) fn family(&self) -> Family {
        match self {
            BackendPlan::ProductState => Family::ProductState,
            BackendPlan::Sparse => Family::Sparse,
            BackendPlan::TensorNetwork { .. } => Family::TensorNetwork,
            BackendPlan::Factored => Family::Factored,
            BackendPlan::FactoredStabilizer => Family::FactoredStabilizer,
            BackendPlan::Mps { .. } => Family::Mps,
            BackendPlan::Stabilizer { .. } => Family::Stabilizer,
            BackendPlan::Statevector { .. } => Family::Statevector,
            BackendPlan::DensityMatrix { .. } => Family::DensityMatrix,
            #[cfg(feature = "distributed")]
            BackendPlan::Distributed(_) => Family::Statevector,
        }
    }

    #[cfg(test)]
    pub(super) fn accel(&self) -> &Accel {
        match self {
            BackendPlan::Stabilizer { accel, .. }
            | BackendPlan::Statevector { accel }
            | BackendPlan::DensityMatrix { accel } => accel,
            _ => &Accel::Cpu,
        }
    }
}

/// Whole-call routing decision produced by [`resolve`]. Backend execution is
/// described by a buildable [`BackendPlan`]; the remaining variants are the
/// non-backend engines (stabilizer rank, Pauli propagation) that callers
/// handle outside the `Backend` world.
pub(super) enum ExecutionPlan {
    Backend(BackendPlan),
    StabilizerRank,
    StochasticPauli { num_samples: usize },
    DeterministicPauli { truncation: SpdTruncation },
    PauliPath,
}

/// Name the engine `kind` resolves to when that engine can discard state
/// weight or estimate by sampling, `None` when the route is exact.
///
/// Exactness is a property of the route and the circuit together, not of the
/// route alone: an engine that can truncate still answers exactly on a circuit
/// whose state it holds without loss, and only a predicate that proves that for
/// every circuit of the shape may say so here.
///
/// Resolved from the circuit rather than read off a finished run, so
/// [`Simulate::require_exact`] rejects before paying for the state it would
/// throw away. `sampling` marks a shot terminal with no noise model or start
/// state, where `Auto` takes the exact stabilizer-rank expansion for a
/// Clifford+T circuit before consulting the family tree.
///
/// [`Simulate::require_exact`]: crate::sim::Simulate::require_exact
pub(super) fn approximate_route_name(
    kind: &BackendKind,
    circuit: &Circuit,
    sampling: bool,
) -> Option<&'static str> {
    // Diagonal and permutation gates keep the all-zero start on a single basis
    // state, so an MPS holds it at bond 1 and has nothing to truncate. A start
    // state other than |0...0> never reaches an MPS route at all, since
    // `initial_state_plan` rejects every representation derived from that start.
    let mps_is_exact = circuit.is_sparse_friendly();
    match kind {
        BackendKind::Mps { .. } => {
            if mps_is_exact {
                return None;
            }
            return Some("Mps");
        }
        BackendKind::TensorNetworkBounded { .. } => return Some("TensorNetworkBounded"),
        BackendKind::StochasticPauli { .. } => return Some("StochasticPauli"),
        BackendKind::DeterministicPauli { truncation } => {
            return truncation.can_discard().then_some("DeterministicPauli");
        }
        BackendKind::PauliPath { epsilon, max_terms } => {
            let truncation = SpdTruncation::Threshold {
                epsilon: *epsilon,
                max_terms: *max_terms,
            };
            return truncation.can_discard().then_some("PauliPath");
        }
        _ => {}
    }
    if !kind.is_auto() {
        return None;
    }
    if sampling && crate::sim::auto_samples_by_stabilizer_rank(circuit) {
        return None;
    }
    let (_, has_partial_independence) = crate::sim::analyze_independence(circuit);
    match select_auto_backend_choice(circuit, has_partial_independence) {
        Family::Mps if !mps_is_exact => Some("Mps"),
        _ => None,
    }
}

pub(super) fn plan_for_family(
    kind: &BackendKind,
    family: Family,
    num_qubits: usize,
) -> BackendPlan {
    match family {
        Family::ProductState => BackendPlan::ProductState,
        Family::Sparse => BackendPlan::Sparse,
        Family::TensorNetwork => BackendPlan::TensorNetwork {
            tolerance: match kind {
                BackendKind::TensorNetworkBounded { tolerance } => Some(*tolerance),
                _ => None,
            },
        },
        Family::Factored => BackendPlan::Factored,
        Family::FactoredStabilizer => BackendPlan::FactoredStabilizer,
        Family::Mps => BackendPlan::Mps {
            max_bond_dim: AUTO_MPS_BOND_DIM,
        },
        Family::Stabilizer => BackendPlan::Stabilizer {
            accel: accel_for(kind, Family::Stabilizer, num_qubits),
            lazy: false,
        },
        Family::Statevector => BackendPlan::Statevector {
            accel: accel_for(kind, Family::Statevector, num_qubits),
        },
        Family::DensityMatrix => BackendPlan::DensityMatrix {
            accel: accel_for(kind, Family::DensityMatrix, num_qubits),
        },
    }
}

/// Resolve `kind` against `circuit` into an [`ExecutionPlan`]. Auto kinds run
/// the shape-based decision tree; explicit kinds map 1:1 onto their family. The
/// CPU-vs-GPU verdict, including the VRAM-fit query, is taken here through
/// [`accel_for`] and reused across shots; soft-mode init fallback covers VRAM
/// shrinking after resolution.
///
/// A shot loop resolves once and replays the plan, with one exception: a
/// decomposable circuit carrying a temporal Clifford prefix keeps the
/// per-shot route through the full pipeline, so this runs once per shot there.
/// The tree is pure and returns the same plan every time, so the repeat costs
/// its own instruction scans and changes nothing. Work added here is paid per
/// shot on that route, not once.
pub(super) fn resolve(
    kind: &BackendKind,
    circuit: &Circuit,
    has_partial_independence: bool,
) -> ExecutionPlan {
    let family = match kind {
        BackendKind::Auto => select_auto_backend_choice(circuit, has_partial_independence),
        #[cfg(feature = "gpu")]
        BackendKind::AutoGpu { .. } => {
            select_auto_backend_choice(circuit, has_partial_independence)
        }
        BackendKind::Statevector => Family::Statevector,
        BackendKind::Stabilizer => Family::Stabilizer,
        BackendKind::Sparse => Family::Sparse,
        BackendKind::Mps { max_bond_dim } => {
            return ExecutionPlan::Backend(BackendPlan::Mps {
                max_bond_dim: *max_bond_dim,
            });
        }
        BackendKind::ProductState => Family::ProductState,
        BackendKind::TensorNetwork | BackendKind::TensorNetworkBounded { .. } => {
            Family::TensorNetwork
        }
        BackendKind::Factored => Family::Factored,
        BackendKind::FactoredStabilizer => Family::FactoredStabilizer,
        BackendKind::DensityMatrix => Family::DensityMatrix,
        BackendKind::StabilizerRank => return ExecutionPlan::StabilizerRank,
        BackendKind::StochasticPauli { num_samples } => {
            return ExecutionPlan::StochasticPauli {
                num_samples: *num_samples,
            };
        }
        BackendKind::PauliPath { .. } => return ExecutionPlan::PauliPath,
        BackendKind::DeterministicPauli { truncation } => {
            return ExecutionPlan::DeterministicPauli {
                truncation: *truncation,
            };
        }
        #[cfg(feature = "gpu")]
        BackendKind::StatevectorGpu { .. } => Family::Statevector,
        #[cfg(feature = "gpu")]
        BackendKind::StabilizerGpu { .. } => Family::Stabilizer,
        #[cfg(feature = "gpu")]
        BackendKind::DensityMatrixGpu { .. } => Family::DensityMatrix,
        #[cfg(feature = "distributed")]
        BackendKind::StatevectorDistributed { context } => {
            return ExecutionPlan::Backend(BackendPlan::Distributed(context.clone()));
        }
    };
    let mut plan = plan_for_family(kind, family, circuit.num_qubits);
    if let BackendPlan::Stabilizer { lazy, .. } = &mut plan {
        *lazy = stabilizer_can_run_lazy(circuit);
    }
    ExecutionPlan::Backend(plan)
}

/// Backend plan for a run that starts from a caller-supplied state.
///
/// Routing is constrained here rather than consulted. Every shortcut the
/// shape-based tree can pick reads the circuit alone and is only valid from the
/// |0...0⟩ start: a Clifford circuit is a stabilizer state only when its input
/// is one, a product state and a subsystem split assume unentangled inputs, and
/// the Pauli engines propagate observables back to |0...0⟩. Selecting one of
/// them for an arbitrary start state is a wrong answer rather than an error, so
/// only the representations that can hold an arbitrary state are reachable: the
/// statevector, on the host, on a device, or sharded across ranks, and the
/// density matrix.
/// `Auto` lands on the statevector unconditionally: the caller already holds
/// `2^n` amplitudes, so the dense state is affordable by construction.
///
/// A start mixture narrows further, to the density matrix and its device
/// sibling: no other representation holds one, and `Auto` never selects the
/// density matrix, so it is declined by name like the rest.
pub(super) fn initial_state_plan(
    kind: &BackendKind,
    num_qubits: usize,
    start: StartState<'_>,
) -> Result<BackendPlan> {
    if let StartState::DensityMatrix(_) = start {
        return match kind {
            BackendKind::DensityMatrix => {
                Ok(plan_for_family(kind, Family::DensityMatrix, num_qubits))
            }
            #[cfg(feature = "gpu")]
            BackendKind::DensityMatrixGpu { .. } => {
                Ok(plan_for_family(kind, Family::DensityMatrix, num_qubits))
            }
            other => Err(PrismError::IncompatibleBackend {
                backend: format!("{other:?}"),
                reason: "a start density matrix runs on the density matrix, which holds a \
                         mixture; every other representation holds a pure state, so select \
                         `BackendKind::DensityMatrix` or its device sibling"
                    .into(),
            }),
        };
    }
    match kind {
        BackendKind::Auto | BackendKind::Statevector => {
            Ok(plan_for_family(kind, Family::Statevector, num_qubits))
        }
        #[cfg(feature = "gpu")]
        BackendKind::AutoGpu { .. } | BackendKind::StatevectorGpu { .. } => {
            Ok(plan_for_family(kind, Family::Statevector, num_qubits))
        }
        BackendKind::DensityMatrix => Ok(plan_for_family(kind, Family::DensityMatrix, num_qubits)),
        #[cfg(feature = "gpu")]
        BackendKind::DensityMatrixGpu { .. } => {
            Ok(plan_for_family(kind, Family::DensityMatrix, num_qubits))
        }
        #[cfg(feature = "distributed")]
        BackendKind::StatevectorDistributed { context } => {
            Ok(BackendPlan::Distributed(context.clone()))
        }
        other => Err(PrismError::IncompatibleBackend {
            backend: format!("{other:?}"),
            reason: "a start state other than |0...0> runs on the statevector or the density \
                     matrix; every other representation is derived from that start"
                .into(),
        }),
    }
}

pub(super) fn resolve_backend(
    kind: &BackendKind,
    circuit: &Circuit,
    has_partial_independence: bool,
) -> BackendPlan {
    match resolve(kind, circuit, has_partial_independence) {
        ExecutionPlan::Backend(plan) => plan,
        _ => unreachable!("non-backend dispatch should be handled by caller"),
    }
}

#[inline]
pub(super) fn min_clifford_prefix_gates(num_qubits: usize) -> usize {
    (num_qubits * 2).max(16)
}

pub(super) fn has_temporal_clifford_opportunity(kind: &BackendKind, circuit: &Circuit) -> bool {
    if !kind.is_auto() {
        return false;
    }
    if circuit.num_qubits > max_statevector_qubits() {
        return false;
    }
    // The split pays for itself only when the tail needs a dense state. A
    // Clifford-only circuit has no such tail, and exporting the tableau to run
    // measurements or a guarded region densely is strictly a loss.
    if circuit.is_clifford_only() {
        return false;
    }
    let min_gates = min_clifford_prefix_gates(circuit.num_qubits);
    let mut prefix_gates = 0;
    for inst in &circuit.instructions {
        match inst {
            Instruction::Gate { gate, .. } => {
                if !gate.is_clifford() {
                    break;
                }
                prefix_gates += 1;
            }
            Instruction::Measure { .. }
            | Instruction::Reset { .. }
            | Instruction::Conditional { .. }
            | Instruction::Region(_) => break,
            Instruction::Barrier { .. } | Instruction::Save { .. } => {}
        }
    }
    prefix_gates >= min_gates && prefix_gates < circuit.instructions.len()
}

/// Temporal-Clifford execution split into a seed-independent plan and a
/// per-seed run, so shot loops split and fuse the circuit once instead of
/// once per shot. The stabilizer prefix runs on the host tableau; the
/// crossover data in [`gpu_capability`] makes a device prefix unreachable
/// here (temporal-Clifford requires fitting the dense statevector, far below
/// the stabilizer crossover).
pub(super) struct TemporalCliffordPlan {
    prefix: Circuit,
    fused_tail: Circuit,
    tail_num_classical_bits: usize,
    tail_accel: Accel,
}

pub(super) fn plan_temporal_clifford(
    kind: &BackendKind,
    circuit: &Circuit,
) -> Option<TemporalCliffordPlan> {
    if !kind.is_auto() {
        return None;
    }
    if circuit.num_qubits > max_statevector_qubits() {
        return None;
    }
    if circuit.is_clifford_only() {
        return None;
    }
    let (prefix, tail) = circuit.clifford_prefix_split()?;
    if prefix.gate_count() < min_clifford_prefix_gates(circuit.num_qubits) {
        return None;
    }
    let tail_accel = accel_for(kind, Family::Statevector, circuit.num_qubits);
    let tail_num_classical_bits = tail.num_classical_bits;
    let fused_tail = match &tail_accel {
        Accel::Cpu => crate::circuit::fusion::fuse_circuit(&tail, true).into_owned(),
        #[cfg(feature = "gpu")]
        Accel::Gpu { .. } => {
            let expanded = crate::circuit::expand_qft_blocks(&tail);
            let expanded = crate::circuit::expand_pauli_rotations(&expanded).into_owned();
            crate::circuit::fusion::fuse_circuit(&expanded, true).into_owned()
        }
    };
    Some(TemporalCliffordPlan {
        prefix,
        fused_tail,
        tail_num_classical_bits,
        tail_accel,
    })
}

pub(super) fn run_temporal_clifford(
    plan: &TemporalCliffordPlan,
    seed: u64,
    want_probabilities: bool,
) -> Result<RunOutcome> {
    let mut stab = StabilizerBackend::new(seed);
    stab.init(plan.prefix.num_qubits, plan.prefix.num_classical_bits)?;
    stab.enable_lazy_destab();
    for inst in &plan.prefix.instructions {
        stab.apply(inst)?;
    }

    let state = stab.export_statevector()?;

    let mut sv = build_statevector(&plan.tail_accel, seed);
    sv.init_from_state(state, plan.tail_num_classical_bits)?;
    for inst in &plan.fused_tail.instructions {
        sv.apply(inst)?;
    }

    let probabilities = if want_probabilities {
        try_backend_probabilities(&sv)?
    } else {
        None
    };

    Ok(RunOutcome {
        classical_bits: sv.classical_results().to_vec(),
        probabilities,
        metadata: crate::sim::backend_metadata(&sv),
        saves: Vec::new(),
    })
}

#[cfg(all(test, feature = "gpu"))]
mod gpu_crossover_tests {
    use super::*;
    use crate::gates::Gate;

    fn stub_kind() -> BackendKind {
        BackendKind::StatevectorGpu {
            context: GpuContext::stub_for_tests(),
        }
    }

    fn run_query(kind: BackendKind, circuit: &Circuit, seed: u64) -> Result<RunOutcome> {
        crate::sim::simulate(circuit).backend(kind).seed(seed).run()
    }

    // The builder GPU shortcut must compose identically to constructing the
    // variant manually. Uses the stub context at a small circuit so crossover
    // fires and proves the composition is side-effect equivalent.
    #[test]
    fn builder_gpu_wraps_statevector_gpu_variant() {
        let ctx = GpuContext::stub_for_tests();
        let mut circuit = Circuit::new(4, 0);
        circuit.add_gate(Gate::H, &[0]);
        circuit.add_gate(Gate::Cx, &[0, 1]);

        let direct = crate::sim::simulate(&circuit)
            .gpu(ctx.clone())
            .seed(42)
            .run()
            .expect("builder GPU shortcut must honor crossover and route to CPU");
        let manual = crate::sim::simulate(&circuit)
            .backend(stub_kind())
            .seed(42)
            .run()
            .expect("manual variant reference");

        let dp = direct.probabilities.expect("direct probs").to_vec();
        let mp = manual.probabilities.expect("manual probs").to_vec();
        assert_eq!(dp, mp);
    }

    // A 4q circuit is far below the default 14q threshold. If the dispatch
    // layer were to build a GPU backend anyway, `GpuState::new` on the stub
    // context would return `BackendUnsupported`. Success proves the
    // crossover in `select_dispatch` is routing small circuits to the host
    // path.
    #[test]
    fn small_circuit_routes_to_cpu() {
        let mut circuit = Circuit::new(4, 0);
        circuit.add_gate(Gate::H, &[0]);
        circuit.add_gate(Gate::Cx, &[0, 1]);
        circuit.add_gate(Gate::H, &[2]);
        circuit.add_gate(Gate::Cx, &[2, 3]);

        let result = run_query(stub_kind(), &circuit, 42)
            .expect("stub context must not be touched for a 4q circuit");
        let probs = result
            .probabilities
            .expect("probabilities missing")
            .to_vec();

        let mut expected = [0.0_f64; 16];
        expected[0b0000] = 0.25;
        expected[0b0011] = 0.25;
        expected[0b1100] = 0.25;
        expected[0b1111] = 0.25;
        for (i, (p, e)) in probs.iter().zip(&expected).enumerate() {
            assert!((p - e).abs() < 1e-10, "p[{i}] = {p}, expected {e}");
        }
    }

    // `independent_bell_pairs(8)` spans 16 qubits but decomposes into 8
    // independent 2q blocks. With `BackendKind::StatevectorGpu`, each
    // sub-block is below the 14q threshold and must route to CPU. If
    // decomposition failed to fire, the 16q monolithic path would attempt
    // `GpuState::new` through the stub and return `BackendUnsupported`.
    // Success here proves decomposition survives across the GPU dispatch.
    #[test]
    fn decomposable_16q_circuit_runs_per_block_on_cpu() {
        let circuit = crate::circuits::independent_bell_pairs(8);
        assert_eq!(circuit.num_qubits, 16);

        let cpu = run_query(BackendKind::Statevector, &circuit, 42).expect("cpu baseline");
        let gpu = run_query(stub_kind(), &circuit, 42).expect("stub must stay out of the way");

        let cpu_p = cpu.probabilities.expect("cpu probs").to_vec();
        let gpu_p = gpu.probabilities.expect("gpu probs").to_vec();
        assert_eq!(cpu_p.len(), gpu_p.len());
        for (i, (c, g)) in cpu_p.iter().zip(gpu_p.iter()).enumerate() {
            assert!(
                (c - g).abs() < 1e-10,
                "prob[{i}] cpu={c}, gpu={g}, diff={}",
                (c - g).abs()
            );
        }
    }

    fn stabilizer_stub_kind() -> BackendKind {
        BackendKind::StabilizerGpu {
            context: GpuContext::stub_for_tests(),
        }
    }

    // A 4q Clifford circuit is far below the stabilizer GPU threshold, so the
    // stub context must never be touched. Produces the same measurement bits
    // as a plain CPU stabilizer run.
    #[test]
    fn stabilizer_gpu_small_circuit_routes_to_cpu() {
        let mut circuit = Circuit::new(4, 4);
        circuit.add_gate(Gate::H, &[0]);
        circuit.add_gate(Gate::Cx, &[0, 1]);
        circuit.add_gate(Gate::Cx, &[1, 2]);
        circuit.add_gate(Gate::Cx, &[2, 3]);
        circuit.add_measure(0, 0);
        circuit.add_measure(1, 1);
        circuit.add_measure(2, 2);
        circuit.add_measure(3, 3);

        let cpu_run = run_query(BackendKind::Stabilizer, &circuit, 42).expect("cpu baseline");
        let gpu_run = run_query(stabilizer_stub_kind(), &circuit, 42)
            .expect("stub must stay out of the way for small circuits");
        assert_eq!(cpu_run.classical_bits, gpu_run.classical_bits);
    }

    // Non-Clifford circuits are rejected at dispatch time with the same error
    // shape as `BackendKind::Stabilizer`.
    #[test]
    fn stabilizer_gpu_rejects_non_clifford_at_dispatch() {
        let mut circuit = Circuit::new(2, 0);
        circuit.add_gate(Gate::T, &[0]);
        let err = run_query(stabilizer_stub_kind(), &circuit, 42).unwrap_err();
        assert!(matches!(err, PrismError::IncompatibleBackend { .. }));
    }

    fn auto_gpu_stub_kind() -> BackendKind {
        BackendKind::AutoGpu {
            context: GpuContext::stub_for_tests(),
        }
    }

    fn assert_probs_match(a: &RunOutcome, b: &RunOutcome) {
        let ap = a.probabilities.as_ref().expect("probs a").to_vec();
        let bp = b.probabilities.as_ref().expect("probs b").to_vec();
        assert_eq!(ap.len(), bp.len());
        for (i, (x, y)) in ap.iter().zip(bp.iter()).enumerate() {
            assert!((x - y).abs() < 1e-10, "prob[{i}]: {x} vs {y}");
        }
    }

    // A small Clifford circuit selects the stabilizer choice, which sits below
    // the stabilizer GPU crossover, so `AutoGpu` builds a CPU stabilizer and
    // never touches the stub. Results match the plain `Auto` path.
    #[test]
    fn auto_gpu_small_clifford_routes_to_cpu() {
        let mut circuit = Circuit::new(4, 0);
        circuit.add_gate(Gate::H, &[0]);
        circuit.add_gate(Gate::Cx, &[0, 1]);
        circuit.add_gate(Gate::Cx, &[1, 2]);
        circuit.add_gate(Gate::Cx, &[2, 3]);

        let cpu = run_query(BackendKind::Auto, &circuit, 42).expect("cpu auto baseline");
        let gpu = run_query(auto_gpu_stub_kind(), &circuit, 42)
            .expect("stub must stay out of the way below the stabilizer crossover");
        assert_probs_match(&cpu, &gpu);
    }

    // `independent_bell_pairs(8)` spans 16 qubits but decomposes into 8
    // independent 2q blocks, each below the statevector crossover. Every block
    // stays on CPU under `AutoGpu`; if decomposition failed, the monolithic 16q
    // path would still hit the VRAM gate and fall back rather than error.
    #[test]
    fn auto_gpu_decomposable_16q_runs_per_block_on_cpu() {
        let circuit = crate::circuits::independent_bell_pairs(8);
        assert_eq!(circuit.num_qubits, 16);

        let cpu = run_query(BackendKind::Auto, &circuit, 42).expect("cpu auto baseline");
        let gpu = run_query(auto_gpu_stub_kind(), &circuit, 42).expect("stub stays out of the way");
        assert_probs_match(&cpu, &gpu);
    }

    // A 16q entangled non-Clifford circuit selects the statevector choice and
    // clears the 14q crossover, so `AutoGpu` reaches the GPU decision. Because
    // the stub cannot report VRAM, the fits check fails closed and the block
    // runs on CPU: same results as `Auto`, no error. The explicit
    // `StatevectorGpu` path has no VRAM gate, so the same circuit touches the
    // stub and surfaces `BackendUnsupported`, isolating the added fallback.
    #[test]
    fn auto_gpu_large_block_falls_back_to_cpu_on_stub() {
        let mut circuit = Circuit::new(16, 0);
        for q in 0..16 {
            circuit.add_gate(Gate::Rx(0.3), &[q]);
        }
        for q in 0..15 {
            circuit.add_gate(Gate::Cx, &[q, q + 1]);
        }

        let cpu = run_query(BackendKind::Auto, &circuit, 42).expect("cpu auto baseline");
        let gpu = run_query(auto_gpu_stub_kind(), &circuit, 42)
            .expect("stub VRAM query fails closed, so AutoGpu must fall back to CPU without error");
        assert_probs_match(&cpu, &gpu);

        let explicit = BackendKind::StatevectorGpu {
            context: GpuContext::stub_for_tests(),
        };
        assert!(matches!(
            run_query(explicit, &circuit, 42).unwrap_err(),
            PrismError::BackendUnsupported { .. }
        ));
    }
}

#[cfg(all(test, feature = "gpu"))]
mod accel_tests {
    use super::*;

    fn stub() -> Arc<GpuContext> {
        GpuContext::stub_for_tests()
    }

    fn is_cpu(accel: &Accel) -> bool {
        matches!(accel, Accel::Cpu)
    }

    #[test]
    fn auto_gpu_statevector_below_crossover_is_cpu() {
        let kind = BackendKind::AutoGpu { context: stub() };
        let n = crate::gpu::min_qubits() - 1;
        assert!(is_cpu(&accel_for(&kind, Family::Statevector, n)));
    }

    // Above the crossover the stub cannot report VRAM, so the fit gate fails
    // closed and the soft path resolves to the host.
    #[test]
    fn auto_gpu_statevector_fits_fails_closed_on_stub() {
        let kind = BackendKind::AutoGpu { context: stub() };
        let n = crate::gpu::min_qubits() + 2;
        assert!(is_cpu(&accel_for(&kind, Family::Statevector, n)));
    }

    #[test]
    fn auto_gpu_cpu_only_families_stay_cpu() {
        let kind = BackendKind::AutoGpu { context: stub() };
        for family in [
            Family::ProductState,
            Family::Sparse,
            Family::Mps,
            Family::Factored,
            Family::FactoredStabilizer,
            Family::TensorNetwork,
        ] {
            assert!(is_cpu(&accel_for(&kind, family, 1 << 10)), "{family:?}");
        }
    }

    #[test]
    fn explicit_statevector_gpu_is_hard_at_crossover_and_cpu_below() {
        let kind = BackendKind::StatevectorGpu { context: stub() };
        let n = crate::gpu::min_qubits();
        assert!(matches!(
            accel_for(&kind, Family::Statevector, n),
            Accel::Gpu { soft: false, .. }
        ));
        assert!(is_cpu(&accel_for(&kind, Family::Statevector, n - 1)));
    }

    #[test]
    fn explicit_stabilizer_gpu_is_hard_at_crossover_and_cpu_below() {
        let kind = BackendKind::StabilizerGpu { context: stub() };
        let n = crate::gpu::stabilizer_min_qubits();
        assert!(matches!(
            accel_for(&kind, Family::Stabilizer, n),
            Accel::Gpu { soft: false, .. }
        ));
        assert!(is_cpu(&accel_for(&kind, Family::Stabilizer, n - 1)));
    }

    // An explicit GPU kind accelerates only its own family; any other family
    // resolves to the host regardless of size.
    #[test]
    fn explicit_kind_other_family_is_cpu() {
        let sv = BackendKind::StatevectorGpu { context: stub() };
        assert!(is_cpu(&accel_for(&sv, Family::Stabilizer, 1 << 20)));
        let stab = BackendKind::StabilizerGpu { context: stub() };
        assert!(is_cpu(&accel_for(&stab, Family::Statevector, 1 << 20)));
    }

    #[test]
    fn cpu_kinds_are_cpu_everywhere() {
        for family in [Family::Statevector, Family::Stabilizer] {
            assert!(is_cpu(&accel_for(&BackendKind::Auto, family, 1 << 20)));
            assert!(is_cpu(&accel_for(
                &BackendKind::Statevector,
                family,
                1 << 20
            )));
        }
    }

    // `init_from_state` on a soft GPU backend degrades to the host when the
    // device upload fails, preserving the supplied amplitudes.
    #[test]
    fn init_from_state_soft_falls_back_to_host_on_stub() {
        use num_complex::Complex64;
        let mut sv = StatevectorBackend::new(42).with_gpu_auto(stub());
        let amp = Complex64::new(std::f64::consts::FRAC_1_SQRT_2, 0.0);
        let state = vec![amp, Complex64::new(0.0, 0.0), Complex64::new(0.0, 0.0), amp];
        sv.init_from_state(state.clone(), 0).unwrap();
        let exported = sv.export_statevector().unwrap();
        for (e, s) in exported.iter().zip(&state) {
            assert!((e - s).norm() < 1e-12);
        }
    }

    #[test]
    fn init_from_state_hard_errors_on_stub() {
        use num_complex::Complex64;
        let mut sv = StatevectorBackend::new(42).with_gpu(stub());
        let err = sv
            .init_from_state(vec![Complex64::new(1.0, 0.0), Complex64::new(0.0, 0.0)], 0)
            .unwrap_err();
        assert!(matches!(err, PrismError::BackendUnsupported { .. }));
    }
}

/// Table-driven coverage of [`resolve`]: every backend kind against every
/// circuit shape class, asserting the resolved family and execution target.
/// Extend it when a family gains a capability row.
#[cfg(test)]
mod dispatch_matrix_tests {
    use super::*;
    use crate::gates::Gate;

    fn product(n: usize) -> Circuit {
        let mut c = Circuit::new(n, 0);
        for q in 0..n {
            c.add_gate(Gate::Rx(0.3), &[q]);
        }
        c
    }

    fn clifford(n: usize) -> Circuit {
        let mut c = Circuit::new(n, 0);
        c.add_gate(Gate::H, &[0]);
        for q in 0..n - 1 {
            c.add_gate(Gate::Cx, &[q, q + 1]);
        }
        c
    }

    fn dense(n: usize) -> Circuit {
        let mut c = Circuit::new(n, 0);
        for q in 0..n {
            c.add_gate(Gate::Rx(0.3), &[q]);
        }
        for q in 0..n - 1 {
            c.add_gate(Gate::Cx, &[q, q + 1]);
        }
        c
    }

    /// One qubit past the statevector memory cap, or `None` when memory
    /// detection failed and the cap is disabled (no width is oversize then).
    fn oversize_qubits() -> Option<usize> {
        let cap = max_statevector_qubits();
        if cap >= usize::BITS as usize {
            eprintln!(
                "SKIP: statevector qubit cap disabled on this host; skipping oversize checks"
            );
            return None;
        }
        Some(cap + 1)
    }

    fn oversize_sparse(n: usize) -> Circuit {
        let mut c = Circuit::new(n, 0);
        for q in 0..n {
            c.add_gate(Gate::T, &[q]);
        }
        for q in 0..n - 1 {
            c.add_gate(Gate::Cx, &[q, q + 1]);
        }
        c
    }

    fn oversize_dense(n: usize) -> Circuit {
        let mut c = Circuit::new(n, 0);
        for q in 0..n {
            c.add_gate(Gate::Rx(0.3), &[q]);
        }
        for q in 0..n - 1 {
            c.add_gate(Gate::Cx, &[q, q + 1]);
        }
        c
    }

    fn resolved(kind: &BackendKind, circuit: &Circuit, hpi: bool) -> BackendPlan {
        match resolve(kind, circuit, hpi) {
            ExecutionPlan::Backend(plan) => plan,
            _ => panic!("expected a backend plan"),
        }
    }

    fn assert_cpu_family(kind: &BackendKind, circuit: &Circuit, hpi: bool, family: Family) {
        let plan = resolved(kind, circuit, hpi);
        assert_eq!(plan.family(), family, "kind {kind:?}");
        assert!(
            matches!(plan.accel(), Accel::Cpu),
            "kind {kind:?} family {family:?} must resolve to the host"
        );
    }

    fn auto_kinds() -> Vec<BackendKind> {
        #[cfg(feature = "gpu")]
        {
            vec![
                BackendKind::Auto,
                BackendKind::AutoGpu {
                    context: GpuContext::stub_for_tests(),
                },
            ]
        }
        #[cfg(not(feature = "gpu"))]
        {
            vec![BackendKind::Auto]
        }
    }

    // Auto and AutoGpu make identical family choices across every circuit
    // shape class; on the stub context every choice resolves to the host
    // (small circuits by crossover, large ones by the fail-closed VRAM gate).
    #[test]
    fn auto_family_matrix() {
        let oversize = oversize_qubits();
        for kind in auto_kinds() {
            assert_cpu_family(&kind, &product(6), false, Family::ProductState);
            assert_cpu_family(&kind, &clifford(6), false, Family::Stabilizer);
            assert_cpu_family(&kind, &clifford(16), false, Family::Stabilizer);
            assert_cpu_family(&kind, &dense(8), false, Family::Statevector);
            assert_cpu_family(&kind, &dense(16), false, Family::Statevector);
            assert_cpu_family(&kind, &dense(8), true, Family::Factored);
            if let Some(n) = oversize {
                assert_cpu_family(&kind, &oversize_sparse(n), false, Family::Sparse);
                assert_cpu_family(&kind, &oversize_dense(n), false, Family::Mps);
                let w_state = crate::circuits::w_state_circuit(n);
                assert_cpu_family(&kind, &w_state, false, Family::Sparse);
            }
        }
    }

    #[test]
    fn sparse_support_probe_bounds_every_gate_not_the_final_state() {
        // A W state ends on 64 basis states but holds 126 midway through its last
        // controlled rotation.
        let w_state = crate::circuits::w_state_circuit(64);
        assert!(sparse_support_fits(&w_state, 126));
        assert!(!sparse_support_fits(&w_state, 125));
        assert!(!sparse_support_fits(
            &oversize_dense(40),
            AUTO_SPARSE_MAX_SUPPORT
        ));

        let mut measured = Circuit::new(8, 1);
        measured.instructions = crate::circuits::w_state_circuit(8).instructions;
        measured.add_measure(0, 0);
        assert!(sparse_support_fits(&measured, AUTO_SPARSE_MAX_SUPPORT));
        measured.add_gate(Gate::H, &[1]);
        assert!(!sparse_support_fits(&measured, AUTO_SPARSE_MAX_SUPPORT));
    }

    #[test]
    fn auto_oversize_mps_uses_auto_bond_dim() {
        let Some(n) = oversize_qubits() else { return };
        for kind in auto_kinds() {
            let plan = resolved(&kind, &oversize_dense(n), false);
            assert!(matches!(
                plan,
                BackendPlan::Mps {
                    max_bond_dim: AUTO_MPS_BOND_DIM
                }
            ));
        }
    }

    // Explicit CPU kinds map 1:1 onto their family regardless of circuit
    // shape, always on the host.
    #[test]
    fn explicit_cpu_kind_matrix() {
        let circuit = dense(6);
        let cases = [
            (BackendKind::Statevector, Family::Statevector),
            (BackendKind::Stabilizer, Family::Stabilizer),
            (BackendKind::Sparse, Family::Sparse),
            (BackendKind::ProductState, Family::ProductState),
            (BackendKind::TensorNetwork, Family::TensorNetwork),
            (BackendKind::Factored, Family::Factored),
            (BackendKind::FactoredStabilizer, Family::FactoredStabilizer),
        ];
        for (kind, family) in cases {
            assert_cpu_family(&kind, &circuit, false, family);
        }

        let plan = resolved(&BackendKind::Mps { max_bond_dim: 77 }, &circuit, false);
        assert!(matches!(plan, BackendPlan::Mps { max_bond_dim: 77 }));
    }

    #[test]
    fn non_backend_kinds_resolve_to_their_engines() {
        let circuit = dense(6);
        assert!(matches!(
            resolve(&BackendKind::StabilizerRank, &circuit, false),
            ExecutionPlan::StabilizerRank
        ));
        assert!(matches!(
            resolve(
                &BackendKind::StochasticPauli { num_samples: 9 },
                &circuit,
                false
            ),
            ExecutionPlan::StochasticPauli { num_samples: 9 }
        ));
        assert!(matches!(
            resolve(
                &BackendKind::DeterministicPauli {
                    truncation: SpdTruncation::Threshold {
                        epsilon: 0.5,
                        max_terms: 3
                    }
                },
                &circuit,
                false
            ),
            ExecutionPlan::DeterministicPauli {
                truncation: SpdTruncation::Threshold {
                    epsilon: e,
                    max_terms: 3
                }
            } if e == 0.5
        ));
    }

    #[test]
    fn density_matrix_is_explicit_only() {
        // Auto never routes to the density-matrix family, whatever the shape.
        for circuit in [product(6), clifford(6), dense(6)] {
            for hpi in [false, true] {
                assert_ne!(
                    resolved(&BackendKind::Auto, &circuit, hpi).family(),
                    Family::DensityMatrix
                );
            }
        }
        // Explicit selection maps 1:1 onto the density-matrix plan, and the
        // backend it builds accepts fused payloads on every terminal.
        let plan = resolved(&BackendKind::DensityMatrix, &dense(6), false);
        assert_eq!(plan.family(), Family::DensityMatrix);
        assert!(plan.build(42).supports_fused_gates());
        assert!(matches!(plan.accel(), Accel::Cpu));
    }

    #[test]
    fn density_matrix_rejects_over_cap() {
        let cap = max_density_matrix_qubits();
        if cap >= usize::BITS as usize {
            eprintln!("SKIP: density-matrix cap disabled on this host");
            return;
        }
        let circuit = dense(cap + 1);
        let err = validate_explicit_backend(&BackendKind::DensityMatrix, &circuit).unwrap_err();
        assert!(matches!(err, PrismError::ResourceLimit(_)));
        // At or below the cap the same shape validates.
        assert!(validate_explicit_backend(&BackendKind::DensityMatrix, &dense(cap)).is_ok());
    }

    // Explicit GPU kinds resolve hard device execution at their crossover
    // (no VRAM gate) and host execution below it; the family choice never
    // changes with the target.
    #[cfg(feature = "gpu")]
    #[test]
    fn explicit_gpu_kind_matrix() {
        let sv_kind = BackendKind::StatevectorGpu {
            context: GpuContext::stub_for_tests(),
        };
        let large = dense(crate::gpu::min_qubits());
        let plan = resolved(&sv_kind, &large, false);
        assert_eq!(plan.family(), Family::Statevector);
        assert!(matches!(plan.accel(), Accel::Gpu { soft: false, .. }));
        assert_cpu_family(
            &sv_kind,
            &dense(crate::gpu::min_qubits() - 1),
            false,
            Family::Statevector,
        );

        let stab_kind = BackendKind::StabilizerGpu {
            context: GpuContext::stub_for_tests(),
        };
        let huge = Circuit::new(crate::gpu::stabilizer_min_qubits(), 0);
        let plan = resolved(&stab_kind, &huge, false);
        assert_eq!(plan.family(), Family::Stabilizer);
        assert!(matches!(plan.accel(), Accel::Gpu { soft: false, .. }));
        assert_cpu_family(&stab_kind, &clifford(8), false, Family::Stabilizer);
    }

    // Three CX on one pair form a block of rank at most 4, so the middle cut
    // reads bonds 2, 4, 4; a CX on a neighbouring pair in between ends the block
    // and the next CX doubles the bond from where it stands.
    #[test]
    fn mps_apply_cost_merges_a_pair_block_and_breaks_it_on_a_neighbour() {
        let overhead = MPS_GATE_CUT_OVERHEAD;
        let mut block = Circuit::new(8, 0);
        for _ in 0..3 {
            block.add_gate(Gate::Cx, &[3, 4]);
            block.add_gate(Gate::Ry(0.4), &[4]);
        }
        let cost = mps_apply_cost(&block, 256);
        assert_eq!(cost, 3.0 * overhead + 8.0 + 64.0 + 64.0);

        let mut broken = Circuit::new(8, 0);
        broken.add_gate(Gate::Cx, &[3, 4]);
        broken.add_gate(Gate::Cx, &[3, 4]);
        broken.add_gate(Gate::Cx, &[2, 3]);
        broken.add_gate(Gate::Cx, &[3, 4]);
        let cost = mps_apply_cost(&broken, 256);
        assert_eq!(cost, 4.0 * overhead + 8.0 + 64.0 + 8.0 + 512.0);
    }
}