prism-q 0.35.0

//! Stabilizer (Clifford) simulation backend.
//!
//! Uses the Gottesman-Knill theorem: Clifford circuits are simulated on a
//! stabilizer tableau of 2n Pauli strings in O(n) time per gate.
//!
//! # Memory layout
//!
//! Aaronson-Gottesman (2004) tableau with (2n+1) rows:
//! - Rows 0..n: destabilizer generators
//! - Rows n..2n: stabilizer generators
//! - Row 2n: scratch row (measurement computation)
//!
//! Each row stores n X-bits, n Z-bits (bit-packed into `Vec<u64>`), and one
//! phase bit (true = -1). Total memory: O(n²/8) bytes.
//!
//! # Gate support
//!
//! Clifford gates only: H, S, Sdg, SX, SXdg, X, Y, Z, Id, CX, CZ, SWAP.
//! Non-Clifford gates (T, Rx, Ry, Rz, Fused) return `BackendUnsupported`.
//! Pauli expectations are answered on the tableau (each string is +1, -1, or
//! 0), so marginals and observables run at any width.
//!
//! # When to prefer this backend
//!
//! - Clifford-only circuits (randomized benchmarking, error correction).
//! - Large qubit counts (1000+) past the statevector memory cap.
//! - Verification of Clifford subcircuits before layering T gates.
//!
//! # When NOT to use this backend
//!
//! - Any non-Clifford gate: automatic dispatch never routes here, and a
//!   hand-picked stabilizer backend fails at the first such gate. A small
//!   T count goes to the stabilizer-rank and Pauli-propagation engines.
//! - Amplitude or statevector output past about 20 qubits: the export is
//!   dense in 2^n even though the tableau is not.
//!
//! # Performance characteristics
//!
//! - Gate application: O(n) per gate (iterates 2n+1 rows, constant work per row)
//! - Measurement: O(n²) worst case (rowmul is O(n), applied to up to 2n rows)
//! - Memory: n=1000 → ~500 KB, n=10000 → ~50 MB
//! - Probability extraction: O(2^n * n), limited by dense output memory budget

use num_complex::Complex64;

#[cfg(feature = "gpu")]
use crate::backend::dense_probability_len;
use crate::backend::{Backend, overlap, reduced_density, schmidt};
use crate::circuit::Instruction;
#[cfg(feature = "gpu")]
use crate::error::PrismError;
use crate::error::Result;
use crate::gates::Gate;
use crate::sim::unified_pauli::PauliTerm;
use rand::SeedableRng;
use rand_chacha::ChaCha8Rng;

use std::borrow::Cow;
#[cfg(feature = "gpu")]
use std::sync::Arc;

pub(crate) mod dense;
pub(crate) mod diagnostics;
pub(crate) mod kernels;
#[cfg(test)]
mod tests;

use kernels::{MIN_WORDS_FOR_BATCH, rowmul_words};

#[cfg(feature = "gpu")]
use crate::gpu::kernels::stabilizer::CliffordBatchScratch;
#[cfg(feature = "gpu")]
use crate::gpu::{GpuContext, GpuTableau};

/// Clifford-only stabilizer simulation on an Aaronson-Gottesman tableau,
/// O(n²) space.
///
/// Cloning panics while a GPU tableau is attached: the device `CudaSlice`
/// cannot be duplicated, and a clone without it would corrupt later CPU-path
/// calls.
pub struct StabilizerBackend {
    pub(super) n: usize,
    pub(super) num_words: usize,
    pub(super) xz: Vec<u64>,
    pub(super) phase: Vec<bool>,
    pub(super) classical_bits: Vec<bool>,
    pub(super) rng: ChaCha8Rng,
    pub(super) qubit_active: Vec<Vec<u32>>,
    pub(super) total_weight: usize,
    pub(super) sgi_merge_buf: Vec<u32>,
    pub(super) sgi_new_a: Vec<u32>,
    pub(super) sgi_new_b: Vec<u32>,
    pub(super) sgi_max_active: usize,
    /// Set by every dense path that mutates the tableau without maintaining
    /// `qubit_active`. The counters `sgi_enabled` reads only describe the rows
    /// while the index is live, so a bulk entry rebuilds before consulting them.
    pub(super) sgi_stale: bool,
    pub(super) lazy_destab: bool,
    pub(super) gate_row_start: usize,
    #[cfg(feature = "gpu")]
    pub(super) gpu_context: Option<Arc<GpuContext>>,
    #[cfg(feature = "gpu")]
    pub(super) gpu_tableau: Option<GpuTableau>,
    /// Pending Clifford ops queued for the next batch launch. Flat layout of
    /// `[opcode, a, b, pad]` quads matching `CLIFOP_STRIDE` in
    /// `gpu::kernels::stabilizer`. Empty outside GPU mode.
    #[cfg(feature = "gpu")]
    pub(super) pending_gpu_ops: Vec<u32>,
    #[cfg(feature = "gpu")]
    pub(super) gpu_batch_scratch: CliffordBatchScratch,
    #[cfg(feature = "gpu")]
    pub(super) gpu_soft: bool,
}

impl Clone for StabilizerBackend {
    fn clone(&self) -> Self {
        #[cfg(feature = "gpu")]
        if self.gpu_tableau.is_some() {
            // CPU host buffers are cleared while the tableau lives on device. A
            // silent clone would produce a backend with `n > 0`, empty `xz`, and
            // `gpu_tableau: None`, which would panic on the next CPU-path access.
            // No existing caller clones a GPU-attached backend; surface the
            // misuse loudly rather than corrupting state.
            panic!(
                "StabilizerBackend::clone is unsupported while a GPU tableau is attached; \
                 copy the tableau back to host first"
            );
        }
        Self {
            n: self.n,
            num_words: self.num_words,
            xz: self.xz.clone(),
            phase: self.phase.clone(),
            classical_bits: self.classical_bits.clone(),
            rng: self.rng.clone(),
            qubit_active: self.qubit_active.clone(),
            total_weight: self.total_weight,
            sgi_merge_buf: self.sgi_merge_buf.clone(),
            sgi_new_a: self.sgi_new_a.clone(),
            sgi_new_b: self.sgi_new_b.clone(),
            sgi_max_active: self.sgi_max_active,
            sgi_stale: self.sgi_stale,
            lazy_destab: self.lazy_destab,
            gate_row_start: self.gate_row_start,
            #[cfg(feature = "gpu")]
            gpu_context: self.gpu_context.clone(),
            #[cfg(feature = "gpu")]
            gpu_tableau: None,
            #[cfg(feature = "gpu")]
            pending_gpu_ops: Vec::new(),
            #[cfg(feature = "gpu")]
            gpu_batch_scratch: CliffordBatchScratch::default(),
            #[cfg(feature = "gpu")]
            gpu_soft: self.gpu_soft,
        }
    }
}

impl StabilizerBackend {
    pub fn new(seed: u64) -> Self {
        Self {
            n: 0,
            num_words: 0,
            xz: Vec::new(),
            phase: Vec::new(),
            classical_bits: Vec::new(),
            rng: ChaCha8Rng::seed_from_u64(seed),
            qubit_active: Vec::new(),
            total_weight: 0,
            sgi_merge_buf: Vec::new(),
            sgi_new_a: Vec::new(),
            sgi_new_b: Vec::new(),
            sgi_max_active: 0,
            sgi_stale: false,
            lazy_destab: false,
            gate_row_start: 0,
            #[cfg(feature = "gpu")]
            gpu_context: None,
            #[cfg(feature = "gpu")]
            gpu_tableau: None,
            #[cfg(feature = "gpu")]
            pending_gpu_ops: Vec::new(),
            #[cfg(feature = "gpu")]
            gpu_batch_scratch: CliffordBatchScratch::default(),
            #[cfg(feature = "gpu")]
            gpu_soft: false,
        }
    }

    /// Initialise the host-resident identity tableau. Shared by the CPU init
    /// path and the soft GPU fallback in [`Backend::init`].
    fn init_host(&mut self, n: usize, nw: usize, num_classical_bits: usize) -> Result<()> {
        self.n = n;
        self.num_words = nw;

        let total_rows = 2 * n + 1;
        let stride = 2 * nw;

        self.xz = vec![0u64; total_rows * stride];
        self.phase = vec![false; total_rows];

        for i in 0..n {
            let word = i / 64;
            let bit = i % 64;
            self.xz[i * stride + word] |= 1u64 << bit;
            self.xz[(i + n) * stride + nw + word] |= 1u64 << bit;
        }

        self.qubit_active = (0..n).map(|q| vec![q as u32, (n + q) as u32]).collect();
        self.total_weight = 2 * n;
        self.sgi_max_active = 2;
        self.sgi_stale = false;

        let want_lazy = self.lazy_destab;
        self.lazy_destab = false;
        self.gate_row_start = 0;
        self.finish_host_init(num_classical_bits);
        if want_lazy {
            self.enable_lazy_destab();
        }
        Ok(())
    }

    /// Start from a tableau exported by [`Self::export_tableau`] instead of
    /// |0...0⟩, so a prepared Clifford state can be handed to a second run.
    ///
    /// Layout, with `n = num_qubits` and `nw = ceil(n / 64)`: `words` holds
    /// `2 * n + 1` rows of `2 * nw` `u64` words, destabilizer rows `0..n`, then
    /// stabilizer rows `n..2n`, then the scratch row the export carries, whose
    /// contents are ignored. A row is its `nw` X words followed by its `nw` Z
    /// words, and qubit `q` is bit `q % 64` of word `q / 64` in each half.
    /// `phases[r]` is true when row `r` carries sign -1.
    ///
    /// The intended input is the export of an earlier run at the same width.
    /// Only the lengths and the diagonal pairing (destabilizer `i` anticommutes
    /// with stabilizer `i`) are checked; a row set that breaks the rest of the
    /// commutation structure is accepted and yields wrong measurement outcomes
    /// without an error. The random stream restarts from the seed given at
    /// construction.
    ///
    /// # Errors
    /// `InvalidParameter` when a length does not match `num_qubits` or a
    /// diagonal pair commutes.
    pub fn init_from_tableau(
        &mut self,
        num_qubits: usize,
        mut words: Vec<u64>,
        mut phases: Vec<bool>,
        num_classical_bits: usize,
    ) -> Result<()> {
        let n = num_qubits;
        let nw = n.div_ceil(64);
        kernels::rowops::check_imported_rows(n, &words, &phases)?;
        words[2 * n * 2 * nw..].fill(0);
        phases[2 * n] = false;

        #[cfg(feature = "gpu")]
        if let Some(ctx) = self.gpu_context.clone() {
            match GpuTableau::new(ctx, n) {
                Ok(mut tableau) => {
                    tableau.copy_from_host(&words, &phases)?;
                    self.attach_device_tableau(tableau, n, nw, num_classical_bits);
                    return Ok(());
                }
                Err(_) if self.gpu_soft => {
                    self.gpu_context = None;
                    self.gpu_tableau = None;
                }
                Err(e) => return Err(e),
            }
        }

        self.n = n;
        self.num_words = nw;
        self.xz = words;
        self.phase = phases;
        self.gate_row_start = if self.lazy_destab { n } else { 0 };
        self.qubit_active = vec![Vec::new(); n];
        self.rebuild_qubit_active();
        self.finish_host_init(num_classical_bits);
        Ok(())
    }

    /// Shared tail of the host init paths once the rows and the SGI index are
    /// in place: drops the GPU bookkeeping the new rows invalidate and sizes
    /// the classical register.
    fn finish_host_init(&mut self, num_classical_bits: usize) {
        #[cfg(feature = "gpu")]
        {
            self.pending_gpu_ops.clear();
            self.gpu_batch_scratch.clear();
        }
        crate::backend::init_classical_bits(&mut self.classical_bits, num_classical_bits);
    }

    /// Commit to GPU mode on a live device tableau: the host rows are dropped
    /// and the SGI counters zeroed, since neither describes the tableau while
    /// it lives on the device.
    #[cfg(feature = "gpu")]
    fn attach_device_tableau(
        &mut self,
        tableau: GpuTableau,
        n: usize,
        nw: usize,
        num_classical_bits: usize,
    ) {
        self.n = n;
        self.num_words = nw;
        self.xz.clear();
        self.phase.clear();
        self.qubit_active = Vec::new();
        self.total_weight = 0;
        self.sgi_max_active = 0;
        self.lazy_destab = false;
        self.gate_row_start = 0;
        self.gpu_tableau = Some(tableau);
        self.pending_gpu_ops.clear();
        self.gpu_batch_scratch.clear();
        crate::backend::init_classical_bits(&mut self.classical_bits, num_classical_bits);
    }

    /// Opt into GPU acceleration using the given shared execution context.
    ///
    /// When set, [`Backend::init`] allocates a device-resident tableau instead of a
    /// host `Vec<u64>` and gate application routes through GPU kernels.
    #[cfg(feature = "gpu")]
    pub fn with_gpu(mut self, context: Arc<GpuContext>) -> Self {
        self.gpu_context = Some(context);
        self
    }

    /// Opt into GPU acceleration but fall back to the host tableau if the device
    /// allocation fails at [`Backend::init`], instead of erroring. Mirrors
    /// [`StatevectorBackend::with_gpu_auto`](crate::backend::statevector::StatevectorBackend::with_gpu_auto);
    /// used by the `Auto` GPU dispatch path, where an unfit or racing device must
    /// degrade to CPU rather than surface a failure.
    #[cfg(feature = "gpu")]
    pub fn with_gpu_auto(mut self, context: Arc<GpuContext>) -> Self {
        self.gpu_context = Some(context);
        self.gpu_soft = true;
        self
    }

    /// XOR `src_row` into `dst_row` and update `dst_row`'s phase per the
    /// Aaronson-Gottesman g-function. Dispatches to the GPU kernel when the
    /// tableau is device-resident, otherwise to the CPU `rowmul_words` SIMD
    /// helper.
    ///
    /// Not stable API: it lets `tests/golden_gpu.rs` drive the GPU rowmul
    /// kernel against the CPU reference, and may change in any release.
    #[doc(hidden)]
    pub fn rowmul_rows_for_testing(&mut self, src_row: usize, dst_row: usize) -> Result<()> {
        #[cfg(feature = "gpu")]
        if self.gpu_tableau.is_some() {
            self.flush_gpu_ops()?;
            let ctx = self
                .gpu_context
                .as_ref()
                .expect("gpu_tableau is_some but gpu_context is None")
                .clone();
            let Some(tableau) = self.gpu_tableau.as_mut() else {
                unreachable!("flush_gpu_ops does not drop gpu_tableau")
            };
            return crate::gpu::kernels::stabilizer::launch_rowmul_words(
                &ctx, tableau, src_row, dst_row,
            );
        }
        let nw = self.num_words;
        let stride = self.stride();
        // Source row must be copied (CPU `rowmul_words` takes `&[_]` for src).
        let src = self.xz[src_row * stride..(src_row + 1) * stride].to_vec();
        let sp = self.phase[src_row];
        let dp = self.phase[dst_row];
        let initial = if sp { 2u64 } else { 0 } + if dp { 2u64 } else { 0 };
        let dst = &mut self.xz[dst_row * stride..(dst_row + 1) * stride];
        let (dx, dz) = dst.split_at_mut(nw);
        let sum = rowmul_words(dx, &mut dz[..nw], &src[..nw], &src[nw..2 * nw], initial);
        self.phase[dst_row] = (sum & 3) >= 2;
        Ok(())
    }

    /// GPU measurement. Flush queued Cliffords, then run the
    /// Aaronson-Gottesman pivot search plus cascade entirely on-device.
    /// Only two small host roundtrips remain: the pivot sentinel (i32) and,
    /// on the deterministic branch, the outcome byte (u8).
    #[cfg(feature = "gpu")]
    fn apply_measure_gpu(&mut self, qubit: usize, classical_bit: usize) -> Result<()> {
        use crate::gpu::kernels::stabilizer as k;
        use rand::RngExt;
        self.flush_gpu_ops()?;
        let ctx = self
            .gpu_context
            .as_ref()
            .expect("apply_measure_gpu called without gpu_context")
            .clone();
        // Find the pivot first so the RNG is only advanced on the random
        // branch, matching the CPU `measure_random` draw pattern.
        let pivot = {
            let Some(tableau) = self.gpu_tableau.as_mut() else {
                unreachable!("apply_measure_gpu called without gpu_tableau")
            };
            k::launch_measure_find_pivot(&ctx, tableau, qubit)?
        };
        let outcome = if let Some(pivot_row) = pivot {
            let random_outcome: bool = self.rng.random();
            let Some(tableau) = self.gpu_tableau.as_mut() else {
                unreachable!("gpu_tableau was Some above the RNG draw")
            };
            k::launch_measure_cascade(&ctx, tableau, qubit, pivot_row)?;
            k::launch_measure_fixup(&ctx, tableau, qubit, pivot_row, random_outcome)?;
            random_outcome
        } else {
            let Some(tableau) = self.gpu_tableau.as_mut() else {
                unreachable!("gpu_tableau was Some in the pivot scope above")
            };
            k::launch_measure_deterministic(&ctx, tableau, qubit)?
        };
        self.classical_bits[classical_bit] = outcome;
        Ok(())
    }

    /// GPU reset. Mirrors the CPU strategy: measure on-device into a scratch
    /// classical slot, then queue an X on the Clifford batch when the
    /// outcome is 1 so the next flush flips the qubit.
    #[cfg(feature = "gpu")]
    fn apply_reset_gpu(&mut self, qubit: usize) -> Result<()> {
        use crate::gpu::kernels::stabilizer::op;
        let prev_len = self.classical_bits.len();
        self.classical_bits.push(false);
        let scratch = prev_len;
        self.apply_measure_gpu(qubit, scratch)?;
        let outcome = self.classical_bits[scratch];
        self.classical_bits.truncate(prev_len);
        if outcome {
            self.queue_1q_gpu(op::X, qubit);
        }
        Ok(())
    }

    /// Copy the device tableau back to host and replay any queued Clifford
    /// ops onto the copied buffers.
    ///
    /// Exists so `&self` read paths (probabilities, export_tableau,
    /// export_statevector) can observe post-flush state without mutating the
    /// device tableau or the queue. The queue is left intact so the next
    /// `&mut self` entry point flushes normally. If the queue is empty the
    /// copied buffers are returned unchanged.
    #[cfg(feature = "gpu")]
    fn copy_device_tableau_with_pending(&self) -> Result<(Vec<u64>, Vec<bool>)> {
        use crate::gpu::kernels::stabilizer::{CLIFOP_STRIDE, op};
        let tableau = self
            .gpu_tableau
            .as_ref()
            .expect("copy_device_tableau_with_pending called without gpu_tableau");
        let (xz, phase) = tableau.copy_to_host()?;
        if self.pending_gpu_ops.is_empty() {
            return Ok((xz, phase));
        }
        let mut cpu = StabilizerBackend::new(0);
        cpu.n = self.n;
        cpu.num_words = self.num_words;
        cpu.xz = xz;
        cpu.phase = phase;
        for chunk in self.pending_gpu_ops.as_chunks::<CLIFOP_STRIDE>().0 {
            let opcode = chunk[0];
            let a = chunk[1] as usize;
            let b = chunk[2] as usize;
            match opcode {
                op::H => cpu.dispatch_gate(&Gate::H, &[a])?,
                op::S => cpu.dispatch_gate(&Gate::S, &[a])?,
                op::SDG => cpu.dispatch_gate(&Gate::Sdg, &[a])?,
                op::X => cpu.dispatch_gate(&Gate::X, &[a])?,
                op::Y => cpu.dispatch_gate(&Gate::Y, &[a])?,
                op::Z => cpu.dispatch_gate(&Gate::Z, &[a])?,
                op::SX => cpu.dispatch_gate(&Gate::SX, &[a])?,
                op::SXDG => cpu.dispatch_gate(&Gate::SXdg, &[a])?,
                op::CX => cpu.dispatch_gate(&Gate::Cx, &[a, b])?,
                op::CZ => cpu.dispatch_gate(&Gate::Cz, &[a, b])?,
                op::SWAP => cpu.dispatch_gate(&Gate::Swap, &[a, b])?,
                _ => {
                    return Err(PrismError::BackendUnsupported {
                        backend: self.name().to_string(),
                        operation: format!("unknown queued Clifford opcode {opcode}"),
                    });
                }
            }
        }
        Ok((cpu.xz, cpu.phase))
    }

    fn dense_probabilities(&self, xz: &[u64], phase: &[bool]) -> Result<Vec<f64>> {
        dense::dense_probabilities(
            xz,
            phase,
            self.n,
            self.num_words,
            self.stride(),
            self.name(),
        )
    }

    /// GPU probabilities. Check the dense output width first, then copy the
    /// tableau back, replaying queued Clifford ops that have not been flushed to
    /// the device so the result matches a fully-flushed state.
    #[cfg(feature = "gpu")]
    fn probabilities_gpu(&self) -> Result<Vec<f64>> {
        dense_probability_len(self.name(), self.n)?;
        let (xz, phase) = self.copy_device_tableau_with_pending()?;
        self.dense_probabilities(&xz, &phase)
    }

    #[cfg(feature = "gpu")]
    fn queue_1q_gpu(&mut self, opcode: u32, target: usize) {
        self.pending_gpu_ops
            .extend_from_slice(&[opcode, target as u32, 0, 0]);
    }

    #[cfg(feature = "gpu")]
    fn queue_2q_gpu(&mut self, opcode: u32, a: usize, b: usize) {
        self.pending_gpu_ops
            .extend_from_slice(&[opcode, a as u32, b as u32, 0]);
    }

    /// Drain and apply the queued Clifford op list in a single kernel launch.
    ///
    /// Must be called before any code path that reads the device tableau or
    /// hands control back to the CPU algorithms (copy-back, probabilities,
    /// measurement, reset, the testing-only rowmul helper, statevector
    /// export). Cheap no-op when the queue is empty.
    #[cfg(feature = "gpu")]
    pub(super) fn flush_gpu_ops(&mut self) -> Result<()> {
        if self.pending_gpu_ops.is_empty() {
            return Ok(());
        }
        let ctx = self
            .gpu_context
            .as_ref()
            .expect("flush_gpu_ops called without gpu_context")
            .clone();
        let tableau = self
            .gpu_tableau
            .as_mut()
            .expect("flush_gpu_ops called without gpu_tableau");
        crate::gpu::kernels::stabilizer::launch_clifford_batch(
            &ctx,
            tableau,
            &self.pending_gpu_ops,
            &mut self.gpu_batch_scratch,
        )?;
        self.pending_gpu_ops.clear();
        Ok(())
    }

    /// Dispatch a Clifford gate for batched GPU execution. Only callable when
    /// `gpu_tableau` is Some and a `gpu_context` is attached; the public
    /// `apply` guard enforces both invariants before this is called.
    ///
    /// Each call appends onto `pending_gpu_ops`; no kernel is launched until
    /// the next `flush_gpu_ops`. Non-Clifford gates fail loudly and leave the
    /// queue intact so diagnostic state is preserved.
    #[cfg(feature = "gpu")]
    fn dispatch_gate_gpu(&mut self, gate: &Gate, targets: &[usize]) -> Result<()> {
        use crate::gpu::kernels::stabilizer::op;
        match gate {
            Gate::Id => Ok(()),
            Gate::H => {
                self.queue_1q_gpu(op::H, targets[0]);
                Ok(())
            }
            Gate::S => {
                self.queue_1q_gpu(op::S, targets[0]);
                Ok(())
            }
            Gate::Sdg => {
                self.queue_1q_gpu(op::SDG, targets[0]);
                Ok(())
            }
            Gate::X => {
                self.queue_1q_gpu(op::X, targets[0]);
                Ok(())
            }
            Gate::Y => {
                self.queue_1q_gpu(op::Y, targets[0]);
                Ok(())
            }
            Gate::Z => {
                self.queue_1q_gpu(op::Z, targets[0]);
                Ok(())
            }
            Gate::SX => {
                self.queue_1q_gpu(op::SX, targets[0]);
                Ok(())
            }
            Gate::SXdg => {
                self.queue_1q_gpu(op::SXDG, targets[0]);
                Ok(())
            }
            Gate::Cx => {
                self.queue_2q_gpu(op::CX, targets[0], targets[1]);
                Ok(())
            }
            Gate::Cz => {
                self.queue_2q_gpu(op::CZ, targets[0], targets[1]);
                Ok(())
            }
            Gate::Swap => {
                self.queue_2q_gpu(op::SWAP, targets[0], targets[1]);
                Ok(())
            }
            Gate::T
            | Gate::Tdg
            | Gate::Rx(_)
            | Gate::Ry(_)
            | Gate::Rz(_)
            | Gate::P(_)
            | Gate::Rzz(_)
            | Gate::Cu(_)
            | Gate::Mcu(_)
            | Gate::Fused(_)
            | Gate::BatchPhase(_)
            | Gate::BatchRzz(_)
            | Gate::DiagonalBatch(_)
            | Gate::MultiFused(_)
            | Gate::Fused2q(_)
            | Gate::Multi2q(_)
            | Gate::QftBlock { .. }
            | Gate::PauliRot(_)
            | Gate::MultiPauliRot(_)
            | Gate::Unitary(_) => Err(PrismError::BackendUnsupported {
                backend: self.name().to_string(),
                operation: format!(
                    "non-Clifford gate `{}` (stabilizer backend supports Clifford gates only)",
                    gate.name()
                ),
            }),
        }
    }

    /// Create a backend that starts in lazy destabilizer mode (see
    /// [`Self::enable_lazy_destab`]).
    pub fn new_lazy(seed: u64) -> Self {
        let mut s = Self::new(seed);
        s.lazy_destab = true;
        s
    }

    /// Defer destabilizer rows: gates update only the n stabilizer rows until
    /// a measurement or reset materializes the destabilizers via Gaussian
    /// elimination. Halves gate cost for long Clifford prefixes with few
    /// measurements. Both the per-instruction [`Backend::apply`] path and the
    /// bulk apply paths preserve laziness across gates.
    ///
    /// Call this on a freshly initialized tableau: the sparse generator index
    /// it installs takes the identity generators' supports as read rather than
    /// scanning the rows.
    pub fn enable_lazy_destab(&mut self) {
        if self.lazy_destab || self.n == 0 {
            return;
        }
        self.lazy_destab = true;
        self.gate_row_start = self.n;
        let n = self.n;
        self.qubit_active = (0..n).map(|q| vec![(n + q) as u32]).collect();
        self.total_weight = n;
        self.sgi_max_active = 1;
        self.sgi_stale = false;
    }

    pub(super) fn ensure_destabilizers(&mut self) {
        if !self.lazy_destab {
            return;
        }
        let sgi_live = self.sgi_enabled();
        kernels::rowops::materialize_destabilizers(
            &mut self.xz,
            &mut self.phase,
            self.n,
            self.num_words,
        );
        self.lazy_destab = false;
        self.gate_row_start = 0;
        if sgi_live {
            // Materialization replaced the stabilizer half with its reduced
            // basis, so the per-row supports the index tracked no longer
            // describe the rows. Rebuild eagerly here: the SGI loop reads the
            // index on the very next instruction and does not consult
            // `sgi_stale`.
            self.rebuild_qubit_active();
        } else {
            self.sgi_stale = true;
        }
    }

    pub fn raw_tableau(&self) -> (&[u64], &[bool]) {
        (&self.xz, &self.phase)
    }

    /// The tableau as a read-only query sees it. Rows `n..2n` are current
    /// whatever the destabilizers are doing, so unlike
    /// [`StabilizerBackend::rows_with_destabilizers`] this borrows the host
    /// rows instead of materializing a copy of the other half.
    fn stabilizer_rows(&self) -> Result<kernels::rowops::TableauRows<'_>> {
        #[cfg(feature = "gpu")]
        if self.gpu_tableau.is_some() {
            let (xz, phase) = self.copy_device_tableau_with_pending()?;
            return Ok((Cow::Owned(xz), Cow::Owned(phase)));
        }
        let (xz, phase) = self.raw_tableau();
        Ok((Cow::Borrowed(xz), Cow::Borrowed(phase)))
    }

    /// Ebits across the cut at `subsystem`, the rank its diagnostics share.
    fn cut_rank(&self, subsystem: &[usize]) -> Result<usize> {
        let (xz, _) = self.stabilizer_rows()?;
        Ok(diagnostics::subsystem_rank(
            &xz,
            self.n,
            self.num_words,
            subsystem,
        ))
    }

    /// Whether the host rows are the live tableau, which they are not while a
    /// device copy holds it.
    #[cfg(feature = "gpu")]
    fn host_rows_live(&self) -> bool {
        self.gpu_tableau.is_none()
    }

    #[cfg(not(feature = "gpu"))]
    fn host_rows_live(&self) -> bool {
        true
    }

    /// Tableau rows with current destabilizers, for a read-only query: the
    /// host rows as they stand, a materialized copy under lazy destabilizers,
    /// or the device tableau brought back with its queued gates applied.
    fn rows_with_destabilizers(&self) -> Result<kernels::rowops::TableauRows<'_>> {
        #[cfg(feature = "gpu")]
        if self.gpu_tableau.is_some() {
            let (xz, phase) = self.copy_device_tableau_with_pending()?;
            return Ok((Cow::Owned(xz), Cow::Owned(phase)));
        }
        if self.lazy_destab {
            let mut xz = self.xz.clone();
            let mut phase = self.phase.clone();
            kernels::rowops::materialize_destabilizers(&mut xz, &mut phase, self.n, self.num_words);
            return Ok((Cow::Owned(xz), Cow::Owned(phase)));
        }
        Ok((Cow::Borrowed(&self.xz), Cow::Borrowed(&self.phase)))
    }

    /// Apply gate and satisfied-conditional instructions, skipping
    /// measurements, resets, and barriers. Routes through the SGI or
    /// word-batch bulk paths at the same thresholds as
    /// [`Backend::apply_instructions`].
    ///
    /// # Errors
    /// A guarded region can hold a measurement or reset, which this path drops
    /// by contract, so one is rejected rather than silently skipped.
    pub fn apply_gates_only(&mut self, instructions: &[Instruction]) -> Result<()> {
        if instructions
            .iter()
            .any(|inst| matches!(inst, Instruction::Region(_)))
        {
            return Err(crate::error::PrismError::IncompatibleBackend {
                backend: "Stabilizer".to_string(),
                reason: "the gates-only path skips measurement and reset, so it cannot carry a \
                         guarded region"
                    .to_string(),
            });
        }
        let nw = self.num_words;
        if nw < MIN_WORDS_FOR_BATCH {
            for instruction in instructions {
                match instruction {
                    Instruction::Gate { gate, targets } => self.dispatch_gate(gate, targets)?,
                    Instruction::Conditional {
                        condition,
                        gate,
                        targets,
                    } if condition.evaluate(&self.classical_bits) => {
                        self.dispatch_gate(gate, targets)?;
                    }
                    _ => {}
                }
            }
            return Ok(());
        }

        if self.sgi_ready() {
            return self.apply_gates_only_sgi(instructions);
        }

        self.apply_gates_only_word_batch(instructions)
    }

    pub(super) fn rowmul(&mut self, h: usize, i: usize) {
        kernels::rowops::rowmul(&mut self.xz, &mut self.phase, self.num_words, h, i);
    }

    pub(super) fn copy_row(&mut self, dst: usize, src: usize) {
        kernels::rowops::copy_row(&mut self.xz, &mut self.phase, self.num_words, dst, src);
    }

    pub(super) fn zero_row(&mut self, r: usize) {
        kernels::rowops::zero_row(&mut self.xz, &mut self.phase, self.num_words, r);
    }

    pub(super) fn apply_measure(&mut self, qubit: usize, classical_bit: usize) {
        self.ensure_destabilizers();
        self.apply_measure_with_info(qubit, classical_bit);
    }

    pub(super) fn apply_reset(&mut self, qubit: usize) -> Result<()> {
        self.ensure_destabilizers();
        let prev_len = self.classical_bits.len();
        self.classical_bits.push(false);
        let scratch = prev_len;
        self.apply_measure_with_info(qubit, scratch);
        let outcome = self.classical_bits[scratch];
        self.classical_bits.truncate(prev_len);
        if outcome {
            self.dispatch_gate(&Gate::X, &[qubit])?;
        }
        Ok(())
    }

    pub(crate) fn apply_measure_with_info(
        &mut self,
        qubit: usize,
        classical_bit: usize,
    ) -> (bool, Vec<usize>) {
        self.sgi_stale = true;
        let n = self.n;
        let word = qubit / 64;
        let bit_mask = 1u64 << (qubit % 64);
        let stride = self.stride();

        let mut p: Option<usize> = None;
        for i in n..2 * n {
            if self.xz[i * stride + word] & bit_mask != 0 {
                p = Some(i);
                break;
            }
        }

        if let Some(p_row) = p {
            let p_base = p_row * stride;
            let mut support = Vec::new();
            for q in 0..n {
                if self.xz[p_base + q / 64] & (1u64 << (q % 64)) != 0 {
                    support.push(q);
                }
            }
            self.measure_random(p_row, word, bit_mask, classical_bit);
            (true, support)
        } else {
            let scratch = 2 * n;
            self.zero_row(scratch);

            for i in 0..n {
                if self.xz[i * stride + word] & bit_mask != 0 {
                    self.rowmul(scratch, i + n);
                }
            }

            let outcome = self.phase[scratch];
            self.classical_bits[classical_bit] = outcome;
            (false, Vec::new())
        }
    }

    /// Measure a run of distinct qubits into distinct classical bits using one
    /// tableau index pass instead of two strided full-tableau scans per
    /// measurement.
    ///
    /// Bit-identical to applying the measurements sequentially: a random
    /// outcome invalidates the index, which is rebuilt before continuing.
    /// Returns, per measurement, whether the outcome was random, the X-support
    /// of the pivot row for random outcomes, and the outcome bit.
    pub(crate) fn batch_measure_ref_info(
        &mut self,
        measurements: &[(usize, usize)],
    ) -> (Vec<bool>, Vec<Vec<usize>>, Vec<bool>) {
        self.sgi_stale = true;
        self.ensure_destabilizers();
        let num_meas = measurements.len();
        let n = self.n;
        let nw = self.num_words;
        let stride = self.stride();

        let mut is_random = vec![false; num_meas];
        let mut random_x_support: Vec<Vec<usize>> = vec![Vec::new(); num_meas];
        let mut outcomes = vec![false; num_meas];

        let mut qubit_to_meas: Vec<usize> = vec![usize::MAX; n];
        for (mi, &(qubit, _)) in measurements.iter().enumerate() {
            qubit_to_meas[qubit] = mi;
        }

        let mut first_destab = vec![usize::MAX; num_meas];

        // Only the stabilizer half can set `first_destab`, so the destabilizer
        // rows carry no information this index reads.
        let build_index = |first_destab: &mut [usize],
                           xz: &[u64],
                           qubit_to_meas: &[usize],
                           n: usize,
                           nw: usize,
                           stride: usize,
                           num_meas: usize| {
            first_destab.fill(usize::MAX);
            for r in n..2 * n {
                let r_base = r * stride;
                for w in 0..nw {
                    let x_word = xz[r_base + w];
                    if x_word == 0 {
                        continue;
                    }
                    let mut bits = x_word;
                    while bits != 0 {
                        let b = bits.trailing_zeros() as usize;
                        let q = w * 64 + b;
                        if q < n {
                            let mi = qubit_to_meas[q];
                            if mi < num_meas && first_destab[mi] == usize::MAX {
                                first_destab[mi] = r;
                            }
                        }
                        bits &= bits - 1;
                    }
                }
            }
        };

        build_index(
            &mut first_destab,
            &self.xz,
            &qubit_to_meas,
            n,
            nw,
            stride,
            num_meas,
        );

        for mi in 0..num_meas {
            let (qubit, classical_bit) = measurements[mi];
            if first_destab[mi] != usize::MAX {
                let (_, support) = self.apply_measure_with_info(qubit, classical_bit);
                is_random[mi] = true;
                outcomes[mi] = self.classical_bits[classical_bit];
                random_x_support[mi] = support;
                build_index(
                    &mut first_destab,
                    &self.xz,
                    &qubit_to_meas,
                    n,
                    nw,
                    stride,
                    num_meas,
                );
            }
        }

        let mut all_diagonal = true;
        'diag_check: for i in 0..n {
            let base = (i + n) * stride;
            for w in 0..nw {
                if self.xz[base + w] != 0 {
                    all_diagonal = false;
                    break 'diag_check;
                }
            }
        }

        if all_diagonal {
            for i in 0..n {
                let phase_i = self.phase[i + n];
                if !phase_i {
                    continue;
                }
                let base = i * stride;
                for w in 0..nw {
                    let x_word = self.xz[base + w];
                    if x_word == 0 {
                        continue;
                    }
                    let mut bits = x_word;
                    while bits != 0 {
                        let b = bits.trailing_zeros() as usize;
                        let q = w * 64 + b;
                        if q < n {
                            let mi = qubit_to_meas[q];
                            if mi < num_meas && !is_random[mi] {
                                outcomes[mi] ^= true;
                            }
                        }
                        bits &= bits - 1;
                    }
                }
            }
            for mi in 0..num_meas {
                if !is_random[mi] {
                    self.classical_bits[measurements[mi].1] = outcomes[mi];
                }
            }
        } else {
            for mi in 0..num_meas {
                if is_random[mi] {
                    continue;
                }
                let (qubit, classical_bit) = measurements[mi];
                let word = qubit / 64;
                let bit_mask = 1u64 << (qubit % 64);
                let scratch = 2 * n;
                self.zero_row(scratch);
                for i in 0..n {
                    if self.xz[i * stride + word] & bit_mask != 0 {
                        self.rowmul(scratch, i + n);
                    }
                }
                outcomes[mi] = self.phase[scratch];
                self.classical_bits[classical_bit] = outcomes[mi];
            }
        }
        (is_random, random_x_support, outcomes)
    }

    /// Export a host-readable snapshot of the tableau: bit-packed `xz` and
    /// per-row `phase`, in the layout [`Self::init_from_tableau`] documents
    /// and accepts. Destabilizers are current in the copy: a lazy backend
    /// materializes them, and a device tableau is copied back with queued ops
    /// replayed.
    pub fn export_tableau(&self) -> Result<(Vec<u64>, Vec<bool>)> {
        let (xz, phase) = self.rows_with_destabilizers()?;
        Ok((xz.into_owned(), phase.into_owned()))
    }

    /// Constructs the 2^n amplitude vector by projecting a support basis state
    /// through each stabilizer generator: |ψ⟩ = ∏_i (I + g_i)/2 |seed⟩, normalized.
    ///
    /// Each projection applies Pauli string g_i to the dense vector in O(2^n),
    /// giving O(n x 2^n) total.
    pub fn export_statevector(&self) -> Result<Vec<Complex64>> {
        #[cfg(feature = "gpu")]
        if self.gpu_tableau.is_some() {
            // Host copy-back path: device tableau plus any queued Clifford ops.
            let (xz, phase) = self.copy_device_tableau_with_pending()?;
            return self.dense_statevector(&xz, &phase);
        }
        self.dense_statevector(&self.xz, &self.phase)
    }

    fn dense_statevector(&self, xz: &[u64], phase: &[bool]) -> Result<Vec<Complex64>> {
        dense::dense_statevector(
            xz,
            phase,
            self.n,
            self.num_words,
            self.stride(),
            self.name(),
            "statevector export",
        )
    }
}

impl crate::backend::sealed::Sealed for StabilizerBackend {}

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

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

    fn init(&mut self, num_qubits: usize, num_classical_bits: usize) -> Result<()> {
        let n = num_qubits;
        let nw = n.div_ceil(64);

        #[cfg(feature = "gpu")]
        if let Some(ctx) = self.gpu_context.clone() {
            // Allocate the device tableau first so an allocation failure returns
            // cleanly without touching any existing state. The transition to GPU
            // mode commits only once the tableau is live.
            let tableau = match GpuTableau::new(ctx, n) {
                Ok(tableau) => tableau,
                Err(_) if self.gpu_soft => {
                    self.gpu_context = None;
                    self.gpu_tableau = None;
                    return self.init_host(n, nw, num_classical_bits);
                }
                Err(e) => return Err(e),
            };
            self.attach_device_tableau(tableau, n, nw, num_classical_bits);
            return Ok(());
        }

        self.init_host(n, nw, num_classical_bits)
    }

    fn apply(&mut self, instruction: &Instruction) -> Result<()> {
        #[cfg(feature = "gpu")]
        if self.gpu_tableau.is_some() {
            match instruction {
                Instruction::Barrier { .. } => return Ok(()),
                Instruction::Save { label, .. } => {
                    return Err(crate::backend::save_not_applied("Stabilizer", label));
                }
                Instruction::Conditional { condition, .. }
                    if !condition.evaluate(&self.classical_bits) =>
                {
                    return Ok(());
                }
                Instruction::Gate { gate, targets } => {
                    return self.dispatch_gate_gpu(gate, targets);
                }
                Instruction::Conditional { gate, targets, .. } => {
                    return self.dispatch_gate_gpu(gate, targets);
                }
                Instruction::Measure {
                    qubit,
                    classical_bit,
                } => {
                    return self.apply_measure_gpu(*qubit, *classical_bit);
                }
                Instruction::Reset { qubit } => {
                    return self.apply_reset_gpu(*qubit);
                }
                Instruction::Region(region) => {
                    return self.apply_region(region);
                }
            }
        }
        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::Save { label, .. } => {
                return Err(crate::backend::save_not_applied("Stabilizer", label));
            }
            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<()> {
        #[cfg(feature = "gpu")]
        if self.gpu_tableau.is_some() {
            return self.apply_reset_gpu(qubit);
        }
        self.apply_reset(qubit)
    }

    fn apply_instructions(&mut self, instructions: &[Instruction]) -> Result<()> {
        #[cfg(feature = "gpu")]
        if self.gpu_tableau.is_some() {
            for instruction in instructions {
                self.apply(instruction)?;
            }
            // Flush here so the device tableau matches the queued ops before
            // the next read (probabilities, export, classical bit access).
            // Measurement and reset flush mid-sequence via `gpu_sync_to_host`.
            self.flush_gpu_ops()?;
            return Ok(());
        }
        let nw = self.num_words;
        if nw < MIN_WORDS_FOR_BATCH {
            for instruction in instructions {
                self.apply(instruction)?;
            }
            return Ok(());
        }

        if self.sgi_ready() {
            return self.apply_instructions_sgi(instructions);
        }

        self.apply_instructions_word_batch(instructions)
    }

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

    fn probabilities(&self) -> Result<Vec<f64>> {
        #[cfg(feature = "gpu")]
        if self.gpu_tableau.is_some() {
            return self.probabilities_gpu();
        }
        self.dense_probabilities(&self.xz, &self.phase)
    }

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

    fn supports_fused_gates(&self) -> bool {
        false
    }

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

    /// Zero when the string anticommutes with a generator, otherwise the sign
    /// that puts it in the stabilizer group. One pass over the rows per string,
    /// with both scratch rows allocated once for the whole request.
    fn pauli_expectations(&self, observables: &[Vec<PauliTerm>]) -> Result<Vec<f64>> {
        let (xz, phase) = self.rows_with_destabilizers()?;
        let nw = self.num_words;
        let mut scratch = vec![0u64; 4 * nw];
        let (term, acc) = scratch.split_at_mut(2 * nw);
        observables
            .iter()
            .map(|observable| {
                kernels::rowops::pack_pauli_string(observable, self.n, nw, term)?;
                Ok(kernels::rowops::pauli_expectation(
                    &xz, &phase, self.n, nw, term, acc,
                ))
            })
            .collect()
    }

    /// `2^r` copies of `2^(-r/2)` for the rank `r` of the cut: a stabilizer
    /// state's spectrum is flat, so the rank is the whole of it. The list is
    /// priced against the dense export cap, past which only
    /// [`Backend::entanglement_entropy`] answers.
    fn schmidt_values(&mut self, subsystem: &[usize]) -> Result<Vec<f64>> {
        schmidt::validate_subsystem(subsystem, self.n)?;
        diagnostics::flat_spectrum(self.name(), self.cut_rank(subsystem)?)
    }

    /// `(rank of the generators restricted to the cut) - |A|` in units of
    /// `ln 2`, at any width: the rank is read off one elimination and the
    /// spectrum it stands for is never built.
    fn entanglement_entropy(&mut self, subsystem: &[usize]) -> Result<f64> {
        schmidt::validate_subsystem(subsystem, self.n)?;
        Ok(self.cut_rank(subsystem)? as f64 * std::f64::consts::LN_2)
    }

    /// `2^-k` times the sum over the stabilizer subgroup supported inside
    /// `subsystem`, which is the projector a stabilizer marginal is. The
    /// group is read off one elimination over the complement's columns, so
    /// the `2^n` export is never taken.
    fn reduced_density_matrix(&mut self, subsystem: &[usize]) -> Result<Vec<Complex64>> {
        schmidt::validate_qubit_set(subsystem, self.n)?;
        let dim = reduced_density::reduced_density_side(self.name(), subsystem.len())?;
        let (xz, phase) = self.stabilizer_rows()?;
        Ok(diagnostics::subsystem_density(
            &xz,
            &phase,
            self.n,
            self.num_words,
            subsystem,
            dim,
        ))
    }

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

    /// `2^(n - r)` for the rank `r` of the two stabilizer groups merged, and
    /// `0` on a sign conflict, at any width. A device-resident tableau and
    /// any other representation take the dense route.
    fn overlap_sq(&self, other: &dyn Backend) -> Result<f64> {
        if let Some(tableau) = other
            .as_any()
            .and_then(|any| any.downcast_ref::<StabilizerBackend>())
            && tableau.n == self.n
            && self.host_rows_live()
            && tableau.host_rows_live()
        {
            return Ok(crate::sim::stabilizer_rank::stabilizer_overlap_sq(
                self, tableau, self.n,
            ));
        }
        overlap::export_overlap_sq(
            self.name(),
            self.num_qubits(),
            || self.export_statevector(),
            other,
        )
    }

    fn export_statevector(&self) -> Result<Vec<Complex64>> {
        // Delegate to the inherent method; it handles both CPU and GPU paths.
        StabilizerBackend::export_statevector(self)
    }
}