prism-q 0.35.0

//! Matrix Product State (MPS) simulation backend.
//!
//! Represents an n-qubit state as a chain of rank-3 tensors:
//!
//! ```text
//!   A[0] - A[1] - A[2] - ... - A[n-1]
//! ```
//!
//! Each tensor `A[k]` has shape (χ_left, 2, χ_right) where χ is the bond dimension.
//! Single-qubit gates are absorbed into the local tensor in O(χ²).
//! Two-qubit gates on adjacent sites use SVD truncation in O(χ³).
//! Non-adjacent two-qubit gates are routed via SWAP chains.
//!
//! # Memory layout
//!
//! - Each site tensor: contiguous `Vec<Complex64>` of length bond_left × 2 × bond_right.
//! - Element A[α, i, β] at index: `α * (2 * bond_right) + i * bond_right + β`.
//! - Bond dimension capped at `max_bond_dim`, which has no default: `MpsBackend::new`
//!   takes it, and automatic dispatch passes 256.
//! - SVD truncation applies both a relative tolerance (default 1e-12) and the
//!   bond cap.
//!
//! # Gate support
//!
//! Every single-qubit gate is absorbed into its site tensor. Adjacent
//! two-qubit gates (`Cx`, `Cz`, `Swap`, `Rzz`, `Cu`, `Fused2q`) go through the
//! SVD kernel, and non-adjacent pairs route via SWAP chains. `Mcu` expands to
//! a dense matrix over the gathered sites, `BatchPhase` sweeps its control
//! toward each target (a single phase reduces to the pair kernel), and the
//! other fused families (`Multi2q`, `BatchRzz`, `DiagonalBatch`, `MultiFused`)
//! decompose into the single-site and adjacent-pair kernels.
//!
//! # When to prefer this backend
//!
//! - 1D circuits with limited entanglement growth.
//! - Variational circuits with hardware-efficient ansatz.
//! - Large qubit counts (50+) with low entanglement.
//!
//! # When NOT to use this backend
//!
//! - Circuits with high entanglement across the chain (many Hadamards + random CX).
//! - Small qubit counts where statevector is faster and exact.

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

use crate::backend::{
    Backend, BasisSamples, NORM_CLAMP_MIN, dense_probability_len, dense_statevector_len, overlap,
    reserve_dense_output, schmidt, simd,
};
use crate::circuit::Instruction;
use crate::error::Result;
use crate::gates::Gate;
use crate::sim::unified_pauli::{PauliAxis, PauliTerm};

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

mod svd;

// Benchmark-only surface: benches link the crate externally and cannot see
// pub(crate), so the SVD entry points are pub; doc(hidden) keeps them out of
// the rendered API.
#[cfg(feature = "parallel")]
#[doc(hidden)]
pub use svd::svd_faer;
#[doc(hidden)]
pub use svd::{SvdResult, svd, svd_jacobi};

#[cfg(test)]
use svd::l2_norm;
use svd::{
    ThinQr, WeightSide, cut_rank, discards_resolvable_weight, fill_isometry_site_data,
    fill_scaled_u_site_data, fill_scaled_vt_data, fill_vt_site_data, scratch_slice,
    truncated_svd_rank, widest_block_cut,
};

const ZERO: Complex64 = Complex64::new(0.0, 0.0);
const ONE: Complex64 = Complex64::new(1.0, 0.0);
const DEFAULT_SVD_EPSILON: f64 = 1e-12;
/// Bond dimension at or above which a threshold cut takes the orthogonality
/// center whatever it drops. The environment of an ungauged chain amplifies
/// what a cut sheds by its own conditioning, which grows with the bonds:
/// brickwork chains at bonds to 128 read 1e-15 from the statevector at every
/// threshold to 1e-12, while at bond 256 a 16-qubit chain read 2.8e-7 from
/// values dropped at 8.6e-16 of the largest, under [`SVD_RESOLUTION`], and
/// another seed read 1.0e-8 with nothing dropped at all.
const GAUGE_RANK: usize = 16;
/// Relative size under which a value dropped by a cut on a chain under
/// [`GAUGE_RANK`] is taken as the factorization's own rounding.
const SVD_RESOLUTION: f64 = 4.0 * f64::EPSILON;
const MAX_SVD_SWEEPS: usize = 100;
/// Relative pair threshold of the one-sided Jacobi sweep: a pair rotates while
/// `|a_p^H a_q|` exceeds this times `|a_p| |a_q|`. Judged against the pair's
/// own norms rather than the whole matrix, so the columns of `U` come out
/// orthogonal on the small singular values too; an absolute residual left
/// them 0.56 apart on an 8x8 whose values span six decades.
const JACOBI_REL_TOL: f64 = 1e-15;
#[cfg(feature = "parallel")]
const MIN_BOND_FOR_PAR: usize = 32;
/// Shot count at or above which native sampling splits the shot loop across
/// the Rayon pool; per-shot substreams keep the output identical either way.
#[cfg(feature = "parallel")]
const MIN_SHOTS_FOR_PAR: usize = 32;
/// Leaves per Rayon task in the dense expansion, as a bit count. The split
/// depth is chosen to leave at least `2^6` amplitudes under each prefix, so
/// the gate is per-task work rather than total element count: each element
/// here costs a site contraction, not the few flops a statevector sweep pays.
#[cfg(feature = "parallel")]
const MIN_LEAVES_PER_TASK_BITS: usize = 6;
/// Cap on the split depth, bounding the serial prefix walk and the number of
/// scratch stacks allocated.
#[cfg(feature = "parallel")]
const MAX_SPLIT_BITS: usize = 6;

/// Raw output pointer for the dense expansion's Rayon split.
#[cfg(feature = "parallel")]
#[derive(Copy, Clone)]
struct DenseOutPtr<T>(*mut T);

#[cfg(feature = "parallel")]
// SAFETY: the expansion partitions basis indices by their top-site bits before
// entering Rayon, so each task writes a disjoint set of elements.
unsafe impl<T: Send> Send for DenseOutPtr<T> {}
#[cfg(feature = "parallel")]
// SAFETY: same partition as Send; sharing the wrapper grants no overlap.
unsafe impl<T: Send> Sync for DenseOutPtr<T> {}

#[cfg(feature = "parallel")]
impl<T> DenseOutPtr<T> {
    #[inline(always)]
    unsafe fn store(self, idx: usize, val: T) {
        // SAFETY: same contract as the enclosing unsafe fn.
        unsafe {
            *self.0.add(idx) = val;
        }
    }
}

#[derive(Clone)]
struct SiteTensor {
    bond_left: usize,
    bond_right: usize,
    data: Vec<Complex64>,
}

impl SiteTensor {
    fn new_zero_state() -> Self {
        Self {
            bond_left: 1,
            bond_right: 1,
            data: vec![ONE, ZERO],
        }
    }

    #[inline(always)]
    fn idx(&self, alpha: usize, i: usize, beta: usize) -> usize {
        alpha * (2 * self.bond_right) + i * self.bond_right + beta
    }
}

/// Controlled phase followed by SWAP, symmetric in the two qubits since the
/// phase acts on |11⟩ only.
fn cu_phase_swap_matrix(phase: Complex64) -> [[Complex64; 4]; 4] {
    let z = ZERO;
    let o = ONE;
    [[o, z, z, z], [z, z, o, z], [z, o, z, z], [z, z, z, phase]]
}

fn swap_gate_qubits(g: &[[Complex64; 4]; 4]) -> [[Complex64; 4]; 4] {
    let mut out = [[ZERO; 4]; 4];
    for i in 0..2 {
        for j in 0..2 {
            for k in 0..2 {
                for l in 0..2 {
                    out[j * 2 + i][l * 2 + k] = g[i * 2 + j][k * 2 + l];
                }
            }
        }
    }
    out
}

/// Build a 2^N × 2^N MCU gate matrix (row-major, flattened).
///
/// `qubit_order[i]` maps block-position `i` to the MCU role index:
/// roles 0..num_controls are controls, role num_controls is the target.
/// The matrix is identity except when all controls are |1⟩, where `mat` is
/// applied to the target qubit.
fn mcu_matrix(
    num_controls: usize,
    mat: &[[Complex64; 2]; 2],
    qubit_order: &[usize],
) -> Vec<Complex64> {
    let n = num_controls + 1;
    let dim = 1usize << n;
    let mut gate = vec![ZERO; dim * dim];

    // Block position p → bit (n-1-p) in state index (MSB = leftmost).
    let mut ctrl_mask: usize = 0;
    let mut tgt_bit: usize = 0;
    for (pos, &role) in qubit_order.iter().enumerate() {
        let bit = n - 1 - pos;
        if role < num_controls {
            ctrl_mask |= 1 << bit;
        } else {
            tgt_bit = bit;
        }
    }

    for s_in in 0..dim {
        if (s_in & ctrl_mask) == ctrl_mask {
            let t_in = (s_in >> tgt_bit) & 1;
            for (t_out, row) in mat.iter().enumerate() {
                let s_out = (s_in & !(1 << tgt_bit)) | (t_out << tgt_bit);
                gate[s_out * dim + s_in] = row[t_in];
            }
        } else {
            gate[s_in * dim + s_in] = ONE;
        }
    }
    gate
}

/// Isometry deviation [`MpsBackend::assert_gauge`] admits. A move refactorizes
/// the site it leaves, so the deviation is that factorization's own and does
/// not accumulate: 420 moves over the eight-site fixture in `mps_tests.rs`
/// take the worst site from 4.4e-16 to 1.1e-15.
#[cfg(test)]
const GAUGE_TOLERANCE: f64 = 1e-14;
/// Pauli axis for [`MpsBackend::pauli_expectation`]. Identity factors
/// are conveyed by omission from the factor list.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum MpsPauliAxis {
    X,
    Y,
    Z,
}

fn mps_pauli_factors(observable: &[PauliTerm]) -> Vec<(usize, MpsPauliAxis)> {
    observable
        .iter()
        .map(|term| {
            let axis = match term.axis {
                PauliAxis::X => MpsPauliAxis::X,
                PauliAxis::Y => MpsPauliAxis::Y,
                PauliAxis::Z => MpsPauliAxis::Z,
            };
            (term.qubit, axis)
        })
        .collect()
}

#[inline]
fn pauli_matrix_entry(axis: Option<MpsPauliAxis>, j: usize, k: usize) -> Complex64 {
    match (axis, j, k) {
        (None, 0, 0) | (None, 1, 1) => ONE,
        (None, _, _) => ZERO,
        (Some(MpsPauliAxis::X), 0, 1) | (Some(MpsPauliAxis::X), 1, 0) => ONE,
        (Some(MpsPauliAxis::X), _, _) => ZERO,
        (Some(MpsPauliAxis::Y), 0, 1) => Complex64::new(0.0, -1.0),
        (Some(MpsPauliAxis::Y), 1, 0) => Complex64::new(0.0, 1.0),
        (Some(MpsPauliAxis::Y), _, _) => ZERO,
        (Some(MpsPauliAxis::Z), 0, 0) => ONE,
        (Some(MpsPauliAxis::Z), 1, 1) => Complex64::new(-1.0, 0.0),
        (Some(MpsPauliAxis::Z), _, _) => ZERO,
    }
}

/// Matrix Product State backend with bounded bond dimension.
pub struct MpsBackend {
    num_qubits: usize,
    max_bond_dim: usize,
    svd_epsilon: f64,
    sites: Vec<SiteTensor>,
    logical_to_site: Vec<usize>,
    site_to_logical: Vec<usize>,
    classical_bits: Vec<bool>,
    rng: ChaCha8Rng,
    truncation_discarded: f64,
    /// Site the chain is canonical about, or `None` when no gauge is
    /// recorded. `Some` claims that every site left of it is a left isometry
    /// and every site right of it a right isometry, which is what makes a
    /// cut's discarded singular weight the state's 2-norm error. The two-site
    /// path establishes it where a cut is going to happen and maintains it
    /// from there; a write that cannot hold the claim clears it.
    center: Option<usize>,
    /// Widest bond any cut has written since [`Backend::init`], which decides
    /// with [`GAUGE_RANK`] whether a threshold cut takes the center and is the
    /// peak the result reports through [`Backend::bond_report`]. Read as a
    /// high-water mark so the decision costs a compare per gate rather than a
    /// scan of the chain.
    bond_high_water: usize,
    /// [`crate::backend::mps_workspace_cap_elements`] read once at
    /// construction, so the per-gate check is a field compare rather than an
    /// atomic load.
    workspace_cap: u128,
    /// Grow-only scratch for [`Self::apply_adjacent_two_qubit`], held across
    /// calls so the four contraction buffers cost no allocation per gate.
    /// Every element of each call's window is written by assignment before
    /// any read, so stale contents never leak between calls.
    scratch_theta: Vec<Complex64>,
    scratch_right_t: Vec<Complex64>,
    scratch_theta_prime: Vec<Complex64>,
    scratch_mat: Vec<Complex64>,
    /// Steps [`MpsBackend::move_center`] has taken since construction, so a
    /// test can price the walk a gate order asks for.
    #[cfg(test)]
    center_steps: usize,
    /// Relative squared weight the walk's factorizations have left behind in
    /// the columns they dropped, so a test can hold the walk to its own
    /// rounding.
    #[cfg(test)]
    qr_discarded: f64,
    /// Grow-only scratch for one step of [`MpsBackend::move_center`], which
    /// runs one to three times per two-qubit gate once the center is
    /// maintained: the factorization buffers and the neighbour the step
    /// absorbs `R` into.
    qr: ThinQr,
    scratch_absorb: Vec<Complex64>,
}

// Manual so clones start with empty scratch: the buffers are write-before-read
// and the stabilizer-rank path clones branch states per shot, so copying their
// contents would be pure memcpy waste.
impl Clone for MpsBackend {
    fn clone(&self) -> Self {
        Self {
            num_qubits: self.num_qubits,
            max_bond_dim: self.max_bond_dim,
            svd_epsilon: self.svd_epsilon,
            sites: self.sites.clone(),
            logical_to_site: self.logical_to_site.clone(),
            site_to_logical: self.site_to_logical.clone(),
            classical_bits: self.classical_bits.clone(),
            rng: self.rng.clone(),
            truncation_discarded: self.truncation_discarded,
            center: self.center,
            bond_high_water: self.bond_high_water,
            #[cfg(test)]
            center_steps: self.center_steps,
            #[cfg(test)]
            qr_discarded: self.qr_discarded,
            workspace_cap: self.workspace_cap,
            scratch_theta: Vec::new(),
            scratch_right_t: Vec::new(),
            scratch_theta_prime: Vec::new(),
            scratch_mat: Vec::new(),
            qr: ThinQr::default(),
            scratch_absorb: Vec::new(),
        }
    }
}

impl MpsBackend {
    pub fn new(seed: u64, max_bond_dim: usize) -> Self {
        Self {
            num_qubits: 0,
            max_bond_dim,
            svd_epsilon: DEFAULT_SVD_EPSILON,
            sites: Vec::new(),
            logical_to_site: Vec::new(),
            site_to_logical: Vec::new(),
            classical_bits: Vec::new(),
            rng: ChaCha8Rng::seed_from_u64(seed),
            truncation_discarded: 0.0,
            center: None,
            bond_high_water: 1,
            #[cfg(test)]
            center_steps: 0,
            #[cfg(test)]
            qr_discarded: 0.0,
            workspace_cap: crate::backend::mps_workspace_cap_elements(),
            scratch_theta: Vec::new(),
            scratch_right_t: Vec::new(),
            scratch_theta_prime: Vec::new(),
            scratch_mat: Vec::new(),
            qr: ThinQr::default(),
            scratch_absorb: Vec::new(),
        }
    }

    /// Create an MPS backend that keeps every non-zero singular value.
    pub(crate) fn new_exact(seed: u64) -> Self {
        let mut backend = Self::new(seed, usize::MAX);
        backend.svd_epsilon = 0.0;
        backend
    }

    /// Set the relative singular-value truncation threshold: each SVD drops
    /// values at or below `epsilon` times the largest, booked in
    /// [`Self::truncation_discarded`].
    ///
    /// The default 1e-12 keeps every numerically meaningful value. A larger
    /// threshold trades fidelity for bond dimension, and at a fixed cap can
    /// lower the total discarded weight by cutting low-weight Schmidt tails
    /// before the cap discards real weight.
    ///
    /// # Panics
    /// Panics unless `0 <= epsilon < 1`.
    pub fn set_svd_epsilon(&mut self, epsilon: f64) {
        assert!(
            (0.0..1.0).contains(&epsilon),
            "svd epsilon must lie in [0, 1)"
        );
        self.svd_epsilon = epsilon;
    }

    /// Zero the cumulative SVD-truncation weight. [`Backend::init`] does this,
    /// so a caller needs it only to scope the total to a sequence shorter than a
    /// run: the CAMPS T-gate path, which must reject silent truncation, brackets
    /// its sequence with this and [`Self::truncation_discarded`].
    ///
    /// [`Backend::init`]: crate::backend::Backend::init
    pub fn reset_truncation_tracking(&mut self) {
        self.truncation_discarded = 0.0;
    }

    /// Cumulative fraction of the state's squared weight discarded by SVD
    /// truncation since [`Backend::init`] or the last
    /// [`Self::reset_truncation_tracking`].
    ///
    /// The unit is squared: a cut that moves the state by `d` in the 2-norm
    /// books `d * d`, so a total of 1e-24 answers for a distance of 1e-12 and
    /// for an infidelity of 1e-24. Comparing it against a distance reads
    /// twelve orders too small.
    ///
    /// At the default threshold a meaningful value means the bond cap discarded
    /// real weight; after [`Self::set_svd_epsilon`] both sources contribute.
    ///
    /// A cut that can lose weight takes the orthogonality center onto its own
    /// sites first, so what it books is the error it made rather than a figure
    /// against a non-orthogonal environment: every cut that can reach the bond
    /// cap, every threshold cut once the chain's bonds have reached
    /// `GAUGE_RANK`, under that a two-site or block threshold cut that would
    /// drop a value above `SVD_RESOLUTION` of its largest, and every cut on a
    /// chain already carrying a center, which a measurement leaves behind on
    /// any chain whose cuts can lose weight at all. The center is parked by
    /// [`svd()`] rather than by the exact walk, so the environment is
    /// orthonormal to that factorization's isometry.
    ///
    /// The total sums one relative discard per SVD rather than measuring the
    /// final state, so a chain that truncates heavily can carry it past 1,
    /// where the fidelity it implies clamps to zero and certifies nothing.
    ///
    /// Truncation does not renormalize the chain. Every read rescales instead:
    /// expectation values, shot sampling, the probability terminal and
    /// [`Self::export_statevector`] all describe the normalized state.
    ///
    /// [`Backend::init`]: crate::backend::Backend::init
    pub fn truncation_discarded(&self) -> f64 {
        self.truncation_discarded
    }

    fn record_truncation(&mut self, singular_values: &[f64], chi_kept: usize) {
        if chi_kept >= singular_values.len() {
            return;
        }
        let total: f64 = singular_values.iter().map(|s| s * s).sum();
        if total <= 0.0 {
            return;
        }
        let discarded: f64 = singular_values[chi_kept..].iter().map(|s| s * s).sum();
        self.truncation_discarded += discarded / total;
    }

    /// Widest bond the chain holds now, 1 for a product state. The run's peak
    /// is `Backend::bond_report`.
    pub fn current_max_bond_dim(&self) -> usize {
        self.sites.iter().map(|s| s.bond_right).max().unwrap_or(1)
    }

    /// Restore the internal MPS site order to logical qubit order without
    /// changing the represented logical state.
    #[cfg(test)]
    pub(crate) fn canonicalize_logical_order(&mut self) -> Result<()> {
        let swap_mat = Gate::Swap.matrix_4x4();
        for target_site in 0..self.num_qubits {
            while self.site_to_logical[target_site] != target_site {
                let logical = target_site;
                let site = self.logical_to_site[logical];
                debug_assert!(site > target_site);
                self.apply_virtual_swap(site - 1, &swap_mat)?;
            }
        }
        Ok(())
    }

    /// `<self|other>` for two chains with the same site layout; an error when
    /// the qubit counts or layouts differ.
    ///
    /// Contracts left-to-right through the environment tensor in two passes
    /// per site, `O(n χ³)` against the naive four-index sum.
    pub fn inner_product(&self, other: &Self) -> Result<Complex64> {
        self.inner_product_with_scratch(other, &mut Vec::new(), &mut Vec::new())
    }

    /// [`Self::inner_product`] with caller-owned scratch. `tmp` and
    /// `next_env` are cleared and resized per site, so a caller evaluating
    /// many overlaps pays one small environment allocation per call instead
    /// of two per site.
    pub(crate) fn inner_product_with_scratch(
        &self,
        other: &Self,
        tmp: &mut Vec<Complex64>,
        next_env: &mut Vec<Complex64>,
    ) -> Result<Complex64> {
        if self.num_qubits != other.num_qubits {
            return Err(crate::error::PrismError::InvalidParameter {
                message: format!(
                    "MPS inner product requires equal qubit counts; got {} vs {}",
                    self.num_qubits, other.num_qubits,
                ),
            });
        }
        if self.logical_to_site != other.logical_to_site {
            return Err(crate::error::PrismError::InvalidParameter {
                message: "MPS inner product requires identical logical_to_site permutations"
                    .to_string(),
            });
        }
        if self.num_qubits == 0 {
            return Ok(ONE);
        }

        let mut env: Vec<Complex64> = vec![ONE; 1];
        let mut env_cols: usize = 1;

        for i in 0..self.num_qubits {
            let a = &self.sites[i];
            let b = &other.sites[i];
            let bl_a = a.bond_left;
            let br_a = a.bond_right;
            let bl_b = b.bond_left;
            let br_b = b.bond_right;

            // Step 1: tmp[β, j, α'] = Σ_α env[α, β] · conj(A[α, j, α'])
            tmp.clear();
            tmp.resize(bl_b * 2 * br_a, ZERO);
            for alpha in 0..bl_a {
                for j in 0..2 {
                    for alpha_p in 0..br_a {
                        let a_val = a.data[a.idx(alpha, j, alpha_p)].conj();
                        if a_val == ZERO {
                            continue;
                        }
                        for beta in 0..bl_b {
                            let e_ab = env[alpha * env_cols + beta];
                            tmp[beta * (2 * br_a) + j * br_a + alpha_p] += a_val * e_ab;
                        }
                    }
                }
            }

            // Step 2: next_env[α', β'] = Σ_{β, j} tmp[β, j, α'] · B[β, j, β']
            next_env.clear();
            next_env.resize(br_a * br_b, ZERO);
            for beta in 0..bl_b {
                for j in 0..2 {
                    for beta_p in 0..br_b {
                        let b_val = b.data[b.idx(beta, j, beta_p)];
                        if b_val == ZERO {
                            continue;
                        }
                        for alpha_p in 0..br_a {
                            let t = tmp[beta * (2 * br_a) + j * br_a + alpha_p];
                            next_env[alpha_p * br_b + beta_p] += t * b_val;
                        }
                    }
                }
            }

            std::mem::swap(&mut env, next_env);
            env_cols = br_b;
        }

        Ok(env[0])
    }

    /// `⟨ψ|P|ψ⟩` for the Pauli string with one factor per listed
    /// `(qubit, axis)`, identity elsewhere. Walks the contraction once
    /// (`O(N·χ³)`) and applies each site's Pauli matrix inline, with no MPS
    /// clone and no gate apply pass.
    /// Duplicate factors on the same qubit are rejected.
    pub fn pauli_expectation(&self, pauli_factors: &[(usize, MpsPauliAxis)]) -> Result<Complex64> {
        if self.num_qubits == 0 {
            return Ok(ONE);
        }

        let site_axis = self.pauli_site_axes(pauli_factors)?;

        let mut env: Vec<Complex64> = vec![ONE; 1];
        let mut env_cols: usize = 1;

        for (i, a) in self.sites.iter().enumerate().take(self.num_qubits) {
            env = self.contract_pauli_env(i, &env, env_cols, site_axis[i]);
            env_cols = a.bond_right;
        }

        Ok(env[0])
    }

    /// Site index of each requested Pauli factor, `None` where the string acts
    /// as identity. Rejects out-of-range qubits and duplicate factors.
    fn pauli_site_axes(
        &self,
        pauli_factors: &[(usize, MpsPauliAxis)],
    ) -> Result<Vec<Option<MpsPauliAxis>>> {
        let mut site_axis: Vec<Option<MpsPauliAxis>> = vec![None; self.num_qubits];
        for &(qubit, axis) in pauli_factors {
            if qubit >= self.num_qubits {
                return Err(crate::error::PrismError::InvalidParameter {
                    message: format!(
                        "MPS Pauli expectation: qubit {qubit} out of range (n={})",
                        self.num_qubits
                    ),
                });
            }
            let site = self.logical_to_site[qubit];
            if site_axis[site].is_some() {
                return Err(crate::error::PrismError::InvalidParameter {
                    message: format!("MPS Pauli expectation: duplicate factor on qubit {qubit}"),
                });
            }
            site_axis[site] = Some(axis);
        }
        Ok(site_axis)
    }

    /// Absorb one site into the `⟨ψ|P|ψ⟩` environment, applying `axis` as the
    /// site's Pauli matrix (identity when `None`).
    ///
    /// `env` is indexed `[bra bond, ket bond]` with `env_cols` columns; the
    /// result is indexed the same way with `bond_right` columns.
    fn contract_pauli_env(
        &self,
        site: usize,
        env: &[Complex64],
        env_cols: usize,
        axis: Option<MpsPauliAxis>,
    ) -> Vec<Complex64> {
        let a = &self.sites[site];
        let bl = a.bond_left;
        let br = a.bond_right;

        // Step 1: tmp[β, j, α'] = Σ_α env[α, β] · conj(A[α, j, α'])
        let mut tmp = vec![ZERO; bl * 2 * br];
        for alpha in 0..bl {
            for j in 0..2 {
                for alpha_p in 0..br {
                    let a_val = a.data[a.idx(alpha, j, alpha_p)].conj();
                    if a_val == ZERO {
                        continue;
                    }
                    for beta in 0..bl {
                        let e_ab = env[alpha * env_cols + beta];
                        tmp[beta * (2 * br) + j * br + alpha_p] += a_val * e_ab;
                    }
                }
            }
        }

        // Step 2: new_env[α', β'] = Σ_{β, j, k} tmp[β, j, α'] · M[j,k] · A[β, k, β']
        // with M = local Pauli matrix at this site (or identity if no factor).
        let mut new_env = vec![ZERO; br * br];
        for beta in 0..bl {
            for j in 0..2 {
                for k in 0..2 {
                    let m_jk = pauli_matrix_entry(axis, j, k);
                    if m_jk == ZERO {
                        continue;
                    }
                    for beta_p in 0..br {
                        let a_bk = a.data[a.idx(beta, k, beta_p)];
                        if a_bk == ZERO {
                            continue;
                        }
                        let ma_val = m_jk * a_bk;
                        for alpha_p in 0..br {
                            let t = tmp[beta * (2 * br) + j * br + alpha_p];
                            new_env[alpha_p * br + beta_p] += t * ma_val;
                        }
                    }
                }
            }
        }

        new_env
    }

    /// Identity environments left of every site, plus `⟨ψ|ψ⟩` as the last
    /// entry, in the `[bra bond, ket bond]` convention of
    /// [`MpsBackend::contract_pauli_env`].
    ///
    /// Entry `s` is the contraction of sites `0..s`, so it has `bond_left(s)`
    /// columns and entry `num_qubits` is the `1x1` norm.
    fn pauli_env_prefixes(&self) -> Vec<Vec<Complex64>> {
        let mut prefixes = Vec::with_capacity(self.num_qubits + 1);
        prefixes.push(vec![ONE; 1]);
        let mut env_cols = 1usize;
        for site in 0..self.num_qubits {
            let next = self.contract_pauli_env(site, prefixes.last().unwrap(), env_cols, None);
            env_cols = self.sites[site].bond_right;
            prefixes.push(next);
        }
        prefixes
    }

    /// Site currently hosting logical qubit `q`. SWAP routing changes this
    /// mapping: non-adjacent two-qubit gates run as chains of nearest-neighbour
    /// SWAPs and grow bond dimension in proportion to site distance. CAMPS
    /// anchor selection uses this to score candidate cascades.
    pub fn site_for_qubit(&self, q: usize) -> usize {
        self.logical_to_site[q]
    }

    /// Whether `P(q = |1⟩) = (1 − Re⟨Z_q⟩)/2` is below `tol`. For a qubit just
    /// disentangled by an OFD cascade, this is the qubit being in the pure `|0⟩`
    /// state on its site.
    pub fn is_qubit_in_zero_state(&self, q: usize, tol: f64) -> Result<bool> {
        let z = self.pauli_expectation(&[(q, MpsPauliAxis::Z)])?;
        let p_one = 0.5 * (1.0 - z.re);
        Ok(p_one < tol)
    }

    #[inline(always)]
    fn site_for_logical(&self, qubit: usize) -> usize {
        self.logical_to_site[qubit]
    }

    #[inline(always)]
    fn logical_for_site(&self, site: usize) -> usize {
        self.site_to_logical[site]
    }

    fn swap_layout_labels(&mut self, left_site: usize) {
        let left_logical = self.site_to_logical[left_site];
        let right_logical = self.site_to_logical[left_site + 1];
        self.site_to_logical.swap(left_site, left_site + 1);
        self.logical_to_site[left_logical] = left_site + 1;
        self.logical_to_site[right_logical] = left_site;
    }

    fn apply_virtual_swap(
        &mut self,
        left_site: usize,
        swap_mat: &[[Complex64; 4]; 4],
    ) -> Result<()> {
        self.apply_adjacent_two_qubit(swap_mat, left_site, true)?;
        self.swap_layout_labels(left_site);
        Ok(())
    }

    /// Apply a single-qubit unitary to `site` in place.
    ///
    /// A unitary leaves the gauge alone: `sum_i conj(u[i][j]) u[i][k]` is
    /// `delta_jk`, which collapses the site's isometry sum back to the one it
    /// had, so a recorded center survives untouched. An arbitrary matrix does
    /// not, which is why [`Backend::apply_1q_matrix`] drops the center first.
    #[inline(always)]
    fn apply_single_qubit_gate(&mut self, site: usize, u: &[[Complex64; 2]; 2]) {
        let t = &mut self.sites[site];
        let bl = t.bond_left;
        let br = t.bond_right;
        let prepared = simd::PreparedGate1q::new(u);
        for alpha in 0..bl {
            let base = alpha * (2 * br);
            let (lo, hi) = t.data[base..base + 2 * br].split_at_mut(br);
            prepared.apply_slice_pairs(lo, hi);
        }
    }

    /// Walk the orthogonality center to `to`, one site at a time.
    ///
    /// Each step refactorizes the site being left and absorbs the triangular
    /// factor into the neighbour it moves toward, so the walk is exact: it
    /// consults neither `svd_epsilon` nor the bond cap and books nothing
    /// through [`Self::truncation_discarded`]. Every site it crosses comes out
    /// an isometry whether or not it was one, and sites beyond `to` keep the
    /// gauge they had, so a chain canonical about [`Self::center`] comes out
    /// canonical about `to`. A chain that was not is repaired only over the
    /// span walked, which is why an unrecorded gauge goes through
    /// [`Self::establish_center`] instead. Bonds never grow: a step caps the
    /// bond it rewrites at the rank of the site it factorizes.
    ///
    /// # Panics
    /// Panics unless `to` indexes a site and a center is recorded.
    fn move_center(&mut self, to: usize) {
        assert!(
            to < self.sites.len(),
            "orthogonality center {to} is outside a chain of {} sites",
            self.sites.len()
        );
        let mut at = self.center.expect("the walk starts from a recorded center");
        #[cfg(test)]
        {
            self.center_steps += at.abs_diff(to);
        }
        while at < to {
            self.shift_center_right(at);
            at += 1;
        }
        while at > to {
            self.shift_center_left(at);
            at -= 1;
        }
        self.center = Some(to);
    }

    /// Record a center at `to` on a chain whose gauge is unknown.
    ///
    /// The first leg runs to the end away from `to`, so it crosses every site
    /// and the center it lands on is real whatever the sites held before. The
    /// second leg is then an ordinary move. Costs `n - 1` steps plus the
    /// distance from that end.
    fn establish_center(&mut self, to: usize) {
        let last = self.sites.len() - 1;
        let (from, via) = if 2 * to <= last { (last, 0) } else { (0, last) };
        self.center = Some(from);
        self.move_center(via);
        self.move_center(to);
    }

    /// Put the center on the pair starting at `left_site` and name the side
    /// the two-site update should leave it on.
    ///
    /// The update writes an isometry into one site of the pair and the
    /// singular weight into the other, and which is which is a write
    /// convention rather than extra work. Leaving the center on the far side
    /// in the direction of travel carries it through a monotone run of
    /// adjacent gates at no extra factorization.
    fn prepare_center(&mut self, left_site: usize) -> WeightSide {
        match self.center {
            Some(at) if at <= left_site => {
                self.move_center(left_site);
                WeightSide::Right
            }
            Some(_) => {
                self.move_center(left_site + 1);
                WeightSide::Left
            }
            None => {
                self.establish_center(left_site);
                WeightSide::Right
            }
        }
    }

    /// Whether a cut of `rank` singular values takes the center before it is
    /// made: it can reach the bond cap, or a threshold is set and the chain's
    /// bonds have reached [`GAUGE_RANK`]. Under that a threshold cut is judged
    /// from its spectrum by [`discards_resolvable_weight`].
    fn cut_can_lose(&self, rank: u128) -> bool {
        rank > self.max_bond_dim as u128
            || (self.svd_epsilon > 0.0
                && (rank >= GAUGE_RANK as u128 || self.bond_high_water >= GAUGE_RANK))
    }

    /// Whether any cut of this chain can lose weight, and so whether a gauge
    /// is worth recording at all. An unbounded cap with a threshold of zero
    /// keeps every singular value of every cut, which is what the exact
    /// constructor asks for, and a center there would put a walk on every
    /// later two-qubit gate for an error no cut makes.
    fn cuts_can_lose_weight(&self) -> bool {
        self.svd_epsilon > 0.0 || self.max_bond_dim < usize::MAX
    }

    /// Drop the recorded center unless the write lands on it.
    ///
    /// A write that is not an isometry on `site` breaks the gauge the record
    /// claims for every site but the center, and the next consumer then
    /// truncates against a non-orthogonal environment with nothing looking
    /// wrong. Rebuilding on demand is cheaper than repairing blind.
    fn invalidate_center_unless(&mut self, site: usize) {
        if self.center != Some(site) {
            self.center = None;
        }
    }

    /// One rightward step of [`Self::move_center`]: `site` keeps `Q` across
    /// its right bond and `R` multiplies into the site to its right.
    fn shift_center_right(&mut self, site: usize) {
        let bl = self.sites[site].bond_left;
        let bond = self.sites[site].bond_right;
        let rows = bl * 2;

        let data = &self.sites[site].data;
        let mat = scratch_slice(&mut self.scratch_mat, rows * bond);
        for alpha in 0..bl {
            for i in 0..2 {
                for gamma in 0..bond {
                    mat[gamma * rows + alpha * 2 + i] = data[alpha * (2 * bond) + i * bond + gamma];
                }
            }
        }

        self.qr.factorize(mat, rows, bond);
        let chi = self.qr.rank;
        #[cfg(test)]
        {
            self.qr_discarded += self.qr.discarded;
        }

        let br = self.sites[site + 1].bond_right;
        let next = &self.sites[site + 1].data;
        let r = &self.qr.r;
        let absorb = scratch_slice(&mut self.scratch_absorb, chi * 2 * br);
        let width = 2 * br;
        let row = |gamma: usize, out: &mut [Complex64]| {
            out.fill(ZERO);
            for delta in 0..bond {
                let weight = r[gamma * bond + delta];
                for (o, x) in out
                    .iter_mut()
                    .zip(&next[delta * width..delta * width + width])
                {
                    *o += weight * x;
                }
            }
        };

        #[cfg(feature = "parallel")]
        if chi >= MIN_BOND_FOR_PAR {
            absorb
                .par_chunks_mut(width)
                .enumerate()
                .for_each(|(gamma, out)| row(gamma, out));
        } else {
            absorb
                .chunks_mut(width)
                .enumerate()
                .for_each(|(gamma, out)| row(gamma, out));
        }

        #[cfg(not(feature = "parallel"))]
        absorb
            .chunks_mut(width)
            .enumerate()
            .for_each(|(gamma, out)| row(gamma, out));

        let right = &mut self.sites[site + 1];
        right.bond_left = chi;
        right.data.clear();
        right.data.extend_from_slice(absorb);

        fill_isometry_site_data(&mut self.sites[site].data, &self.qr.q, rows, bl, chi);
        self.sites[site].bond_right = chi;
    }

    /// One leftward step of [`Self::move_center`], the same factorization on
    /// the conjugate transpose: `site` keeps the rows of `Q†` across its left
    /// bond and `R†` multiplies into the site to its left.
    fn shift_center_left(&mut self, site: usize) {
        let bond = self.sites[site].bond_left;
        let br = self.sites[site].bond_right;
        let rows = 2 * br;

        let data = &self.sites[site].data;
        let mat = scratch_slice(&mut self.scratch_mat, rows * bond);
        for (out, x) in mat.iter_mut().zip(data.iter()) {
            *out = x.conj();
        }

        self.qr.factorize(mat, rows, bond);
        let chi = self.qr.rank;
        #[cfg(test)]
        {
            self.qr_discarded += self.qr.discarded;
        }

        let bl = self.sites[site - 1].bond_left;
        let prev = &self.sites[site - 1].data;
        let r = &self.qr.r;
        let absorb = scratch_slice(&mut self.scratch_absorb, bl * 2 * chi);
        // Summed over the old bond innermost, so the site being read and the
        // row of `R` both run stride 1 and each output element takes a single
        // store.
        let row = |alpha: usize, out: &mut [Complex64]| {
            for i in 0..2 {
                let base = alpha * (2 * bond) + i * bond;
                let amps = &prev[base..base + bond];
                for gamma in 0..chi {
                    let mut acc = ZERO;
                    for (amp, factor) in amps.iter().zip(&r[gamma * bond..gamma * bond + bond]) {
                        acc += amp * factor.conj();
                    }
                    out[i * chi + gamma] = acc;
                }
            }
        };

        #[cfg(feature = "parallel")]
        if bl >= MIN_BOND_FOR_PAR {
            absorb
                .par_chunks_mut(2 * chi)
                .enumerate()
                .for_each(|(alpha, out)| row(alpha, out));
        } else {
            absorb
                .chunks_mut(2 * chi)
                .enumerate()
                .for_each(|(alpha, out)| row(alpha, out));
        }

        #[cfg(not(feature = "parallel"))]
        absorb
            .chunks_mut(2 * chi)
            .enumerate()
            .for_each(|(alpha, out)| row(alpha, out));

        let left = &mut self.sites[site - 1];
        left.bond_right = chi;
        left.data.clear();
        left.data.extend_from_slice(absorb);

        let t = &mut self.sites[site];
        t.bond_left = chi;
        t.data.clear();
        t.data
            .extend(self.qr.q[..rows * chi].iter().map(|x| x.conj()));
    }

    /// Largest deviation from the orthogonality condition over the sites away
    /// from `center`: a site left of it owes `A† A = I` summed over its left
    /// bond and physical index, a site right of it owes `A A† = I` summed over
    /// its physical index and right bond.
    #[cfg(test)]
    fn gauge_deviation(&self, center: usize) -> f64 {
        let mut worst = 0.0f64;
        for (site, t) in self.sites.iter().enumerate() {
            if site < center {
                for g in 0..t.bond_right {
                    for h in 0..t.bond_right {
                        let mut sum = ZERO;
                        for alpha in 0..t.bond_left {
                            for i in 0..2 {
                                sum +=
                                    t.data[t.idx(alpha, i, g)].conj() * t.data[t.idx(alpha, i, h)];
                            }
                        }
                        let target = if g == h { ONE } else { ZERO };
                        worst = worst.max((sum - target).norm());
                    }
                }
            } else if site > center {
                for a in 0..t.bond_left {
                    for b in 0..t.bond_left {
                        let mut sum = ZERO;
                        for i in 0..2 {
                            for beta in 0..t.bond_right {
                                sum += t.data[t.idx(a, i, beta)] * t.data[t.idx(b, i, beta)].conj();
                            }
                        }
                        let target = if a == b { ONE } else { ZERO };
                        worst = worst.max((sum - target).norm());
                    }
                }
            }
        }
        worst
    }

    #[cfg(test)]
    fn assert_gauge(&self, center: usize) {
        let worst = self.gauge_deviation(center);
        assert!(
            worst <= GAUGE_TOLERANCE,
            "gauge deviation {worst:.3e} about center {center} exceeds {GAUGE_TOLERANCE:.0e}"
        );
    }

    /// Apply a two-qubit gate to the pair at `left_site` and `left_site + 1`.
    ///
    /// The cut this makes is exact in the 2-norm only against an orthonormal
    /// environment, so the pair takes the orthogonality center first wherever
    /// the cut can lose anything: before the factorization when
    /// [`Self::cut_can_lose`] says so from the rank and the chain's bonds, and
    /// after it when the spectrum shows a threshold cut dropping a resolvable
    /// value, in which case the pair is cut again in gauge. A pair under both
    /// marks stays off the walk entirely. Maintaining a center already
    /// recorded costs one to three steps, well under rebuilding it, so a chain
    /// that has one keeps it on every pair.
    fn apply_adjacent_two_qubit(
        &mut self,
        gate: &[[Complex64; 4]; 4],
        left_site: usize,
        left_is_first_qubit: bool,
    ) -> Result<()> {
        if !self.cut_adjacent_pair(gate, left_site, left_is_first_qubit)? {
            self.establish_center(left_site);
            self.cut_adjacent_pair(gate, left_site, left_is_first_qubit)?;
        }
        Ok(())
    }

    /// One attempt at [`Self::apply_adjacent_two_qubit`]. Returns `false`,
    /// having written nothing, when the chain has no center and the cut would
    /// drop a value the factorization resolves.
    fn cut_adjacent_pair(
        &mut self,
        gate: &[[Complex64; 4]; 4],
        left_site: usize,
        left_is_first_qubit: bool,
    ) -> Result<bool> {
        let right_site = left_site + 1;
        let bl = self.sites[left_site].bond_left;
        let bond_mid = self.sites[left_site].bond_right;
        let br = self.sites[right_site].bond_right;

        // Live contraction buffers: theta, theta_prime, and the SVD input plus
        // its internals at 4*bl*br each (about five at peak), and the right
        // transpose at 2*bond_mid*br.
        let workspace = 20u128 * bl as u128 * br as u128 + 2u128 * bond_mid as u128 * br as u128;
        if workspace > self.workspace_cap {
            return Err(crate::backend::workspace_allocation_error(
                "mps",
                "two-qubit gate application",
                workspace,
            ));
        }

        let maintained = self.center.is_some() || self.cut_can_lose(cut_rank(bl, bond_mid, br));
        let weight = if maintained {
            self.prepare_center(left_site)
        } else {
            WeightSide::Right
        };

        // A move rewrites the bonds it crosses at the rank it finds, so the
        // shapes are read again here.
        let bl = self.sites[left_site].bond_left;
        let bond_mid = self.sites[left_site].bond_right;
        let br = self.sites[right_site].bond_right;

        let g = if left_is_first_qubit {
            *gate
        } else {
            swap_gate_qubits(gate)
        };

        // 1. Contract: Θ[α, i, j, β] = Σ_γ A[k][α,i,γ] · A[k+1][γ,j,β]
        let left_data = &self.sites[left_site].data;
        let right_data = &self.sites[right_site].data;
        let chunk_size = 2 * 2 * br;
        let theta = scratch_slice(&mut self.scratch_theta, bl * chunk_size);

        // Transpose right: (bond_mid, 2, br) → (2, br, bond_mid) for sequential gamma access.
        let right_t = scratch_slice(&mut self.scratch_right_t, bond_mid * 2 * br);
        for gamma in 0..bond_mid {
            for j in 0..2usize {
                for beta in 0..br {
                    right_t[j * br * bond_mid + beta * bond_mid + gamma] =
                        right_data[gamma * (2 * br) + j * br + beta];
                }
            }
        }

        #[cfg(feature = "parallel")]
        if bl >= MIN_BOND_FOR_PAR {
            let right_ref = &right_t;
            theta
                .par_chunks_mut(chunk_size)
                .enumerate()
                .for_each(|(alpha, chunk)| {
                    for i in 0..2 {
                        for j in 0..2 {
                            for beta in 0..br {
                                let mut val = ZERO;
                                let rt_base = j * br * bond_mid + beta * bond_mid;
                                for gamma in 0..bond_mid {
                                    val += left_data[alpha * (2 * bond_mid) + i * bond_mid + gamma]
                                        * right_ref[rt_base + gamma];
                                }
                                chunk[i * (2 * br) + j * br + beta] = val;
                            }
                        }
                    }
                });
        } else {
            for alpha in 0..bl {
                for i in 0..2 {
                    for j in 0..2 {
                        for beta in 0..br {
                            let mut val = ZERO;
                            let rt_base = j * br * bond_mid + beta * bond_mid;
                            for gamma in 0..bond_mid {
                                val += left_data[alpha * (2 * bond_mid) + i * bond_mid + gamma]
                                    * right_t[rt_base + gamma];
                            }
                            theta[alpha * chunk_size + i * (2 * br) + j * br + beta] = val;
                        }
                    }
                }
            }
        }

        #[cfg(not(feature = "parallel"))]
        for alpha in 0..bl {
            for i in 0..2 {
                for j in 0..2 {
                    for beta in 0..br {
                        let mut val = ZERO;
                        let rt_base = j * br * bond_mid + beta * bond_mid;
                        for gamma in 0..bond_mid {
                            val += left_data[alpha * (2 * bond_mid) + i * bond_mid + gamma]
                                * right_t[rt_base + gamma];
                        }
                        theta[alpha * chunk_size + i * (2 * br) + j * br + beta] = val;
                    }
                }
            }
        }

        // 2. Apply gate: Θ'[α, i', j', β] = Σ_{i,j} G[i'2+j', i2+j] · Θ[α,i,j,β]
        let theta = &self.scratch_theta[..bl * chunk_size];
        let theta_prime = scratch_slice(&mut self.scratch_theta_prime, bl * chunk_size);

        #[cfg(feature = "parallel")]
        if bl >= MIN_BOND_FOR_PAR {
            let theta_ref = &theta;
            theta_prime
                .par_chunks_mut(chunk_size)
                .enumerate()
                .for_each(|(alpha, chunk)| {
                    for ip in 0..2 {
                        for jp in 0..2 {
                            for beta in 0..br {
                                let mut val = ZERO;
                                for i in 0..2 {
                                    for j in 0..2 {
                                        val += g[ip * 2 + jp][i * 2 + j]
                                            * theta_ref
                                                [alpha * chunk_size + i * (2 * br) + j * br + beta];
                                    }
                                }
                                chunk[ip * (2 * br) + jp * br + beta] = val;
                            }
                        }
                    }
                });
        } else {
            for alpha in 0..bl {
                for ip in 0..2 {
                    for jp in 0..2 {
                        for beta in 0..br {
                            let mut val = ZERO;
                            for i in 0..2 {
                                for j in 0..2 {
                                    val += g[ip * 2 + jp][i * 2 + j]
                                        * theta[alpha * chunk_size + i * (2 * br) + j * br + beta];
                                }
                            }
                            theta_prime[alpha * chunk_size + ip * (2 * br) + jp * br + beta] = val;
                        }
                    }
                }
            }
        }

        #[cfg(not(feature = "parallel"))]
        for alpha in 0..bl {
            for ip in 0..2 {
                for jp in 0..2 {
                    for beta in 0..br {
                        let mut val = ZERO;
                        for i in 0..2 {
                            for j in 0..2 {
                                val += g[ip * 2 + jp][i * 2 + j]
                                    * theta[alpha * chunk_size + i * (2 * br) + j * br + beta];
                            }
                        }
                        theta_prime[alpha * chunk_size + ip * (2 * br) + jp * br + beta] = val;
                    }
                }
            }
        }

        // 3. Reshape to matrix M: rows = bl*2, cols = 2*br
        //    M[(α*2+i), (j*br+β)] = Θ'[α, i, j, β]
        let rows = bl * 2;
        let cols = 2 * br;
        let theta_prime = &self.scratch_theta_prime[..bl * chunk_size];
        let mat = scratch_slice(&mut self.scratch_mat, rows * cols);
        for alpha in 0..bl {
            for i in 0..2 {
                for j in 0..2 {
                    for beta in 0..br {
                        let r = alpha * 2 + i;
                        let c = j * br + beta;
                        mat[c * rows + r] =
                            theta_prime[alpha * (2 * 2 * br) + i * (2 * br) + j * br + beta];
                    }
                }
            }
        }

        let svd_result = svd(mat, rows, cols);
        let chi_new = truncated_svd_rank(&svd_result.s, self.svd_epsilon, self.max_bond_dim);
        if !maintained && discards_resolvable_weight(&svd_result.s, chi_new) {
            return Ok(false);
        }
        self.record_truncation(&svd_result.s, chi_new);

        // A[k] has shape (bl, 2, chi_new) and A[k+1] shape (chi_new, 2, br).
        // U is column-major, U[col * rows + row]; V† is row-major,
        // Vt[row * cols + col]. Whichever side takes diag(S) carries the
        // weight, and the other comes out an isometry. Both sites are written
        // back through the buffers they already own, so a chain at a steady
        // bond allocates nothing per gate.
        let mut left_data = std::mem::take(&mut self.sites[left_site].data);
        let mut right_data = std::mem::take(&mut self.sites[right_site].data);
        match weight {
            WeightSide::Right => {
                fill_isometry_site_data(
                    &mut left_data,
                    &svd_result.u,
                    svd_result.u_rows,
                    bl,
                    chi_new,
                );
                fill_scaled_vt_data(&mut right_data, &svd_result, chi_new, 2, br);
            }
            WeightSide::Left => {
                fill_scaled_u_site_data(&mut left_data, &svd_result, bl, chi_new);
                fill_vt_site_data(&mut right_data, &svd_result, chi_new, 2, br);
            }
        }

        self.sites[left_site] = SiteTensor {
            bond_left: bl,
            bond_right: chi_new,
            data: left_data,
        };
        self.sites[right_site] = SiteTensor {
            bond_left: chi_new,
            bond_right: br,
            data: right_data,
        };
        self.bond_high_water = self.bond_high_water.max(chi_new);
        if maintained {
            self.center = Some(match weight {
                WeightSide::Right => right_site,
                WeightSide::Left => left_site,
            });
        }
        Ok(true)
    }

    fn apply_two_qubit_gate(
        &mut self,
        gate: &[[Complex64; 4]; 4],
        q0: usize,
        q1: usize,
    ) -> Result<()> {
        let p0 = self.site_for_logical(q0);
        let p1 = self.site_for_logical(q1);
        let k = p0.min(p1);
        let m = p0.max(p1);
        let left_is_first = p0 < p1;

        if m - k == 1 {
            self.apply_adjacent_two_qubit(gate, k, left_is_first)?;
        } else {
            let swap_mat = Gate::Swap.matrix_4x4();
            for s in (k + 1..m).rev() {
                self.apply_virtual_swap(s, &swap_mat)?;
            }
            self.apply_adjacent_two_qubit(gate, k, left_is_first)?;
        }
        Ok(())
    }

    /// Apply `instructions` in order, except that a run of two-qubit gates on
    /// adjacent, disjoint site pairs with increasing left sites enters from
    /// whichever end the recorded center is nearer. The gates of such a run
    /// commute, and the walk crosses the run once from either end, so the
    /// choice moves the entry cost and nothing else: on a brick layer that is
    /// one step against the width of the chain. A run is counted in gate
    /// entries by [`pair_entries`], so it spans instruction boundaries and
    /// reaches into `Multi2q` and `BatchRzz` lists, which is how fusion
    /// delivers a brick layer; an instruction without entries, a measurement
    /// or a barrier included, ends the run where it stands.
    fn apply_by_center(&mut self, instructions: &[Instruction]) -> Result<()> {
        let (mut i, mut k) = (0, 0);
        while i < instructions.len() {
            if pair_entries(&instructions[i]) == 0 {
                self.apply(&instructions[i])?;
                i += 1;
                continue;
            }
            let (mut ri, mut rk) = (i, k);
            let mut len = 0;
            let (mut first, mut last) = (0, 0);
            'scan: while ri < instructions.len() {
                let entries = pair_entries(&instructions[ri]);
                if entries == 0 {
                    break;
                }
                while rk < entries {
                    let (q0, q1) = pair_entry(&instructions[ri], rk);
                    let p0 = self.site_for_logical(q0);
                    let p1 = self.site_for_logical(q1);
                    let left = p0.min(p1);
                    if p0.abs_diff(p1) != 1 || (len > 0 && left < last + 2) {
                        break 'scan;
                    }
                    if len == 0 {
                        first = left;
                    }
                    last = left;
                    len += 1;
                    rk += 1;
                }
                ri += 1;
                rk = 0;
            }
            if len == 0 {
                self.apply_pair_entry(&instructions[i], k)?;
                k += 1;
                if k == pair_entries(&instructions[i]) {
                    i += 1;
                    k = 0;
                }
                continue;
            }
            if len > 1 && self.entry_steps(last) < self.entry_steps(first) {
                let (mut ii, mut kk) = (ri, rk);
                while (ii, kk) != (i, k) {
                    if kk == 0 {
                        ii -= 1;
                        kk = pair_entries(&instructions[ii]);
                    }
                    kk -= 1;
                    self.apply_pair_entry(&instructions[ii], kk)?;
                }
            } else {
                let (mut ii, mut kk) = (i, k);
                while (ii, kk) != (ri, rk) {
                    self.apply_pair_entry(&instructions[ii], kk)?;
                    kk += 1;
                    if kk == pair_entries(&instructions[ii]) {
                        ii += 1;
                        kk = 0;
                    }
                }
            }
            i = ri;
            k = rk;
        }
        Ok(())
    }

    /// Apply entry `k` of `instruction`, one of the [`pair_entries`] it owes.
    fn apply_pair_entry(&mut self, instruction: &Instruction, k: usize) -> Result<()> {
        match instruction {
            Instruction::Gate {
                gate: Gate::Multi2q(data),
                ..
            } => {
                let (q0, q1, mat) = &data.gates[k];
                self.apply_two_qubit_gate(mat, *q0, *q1)
            }
            Instruction::Gate {
                gate: Gate::BatchRzz(data),
                ..
            } => {
                let (q0, q1, theta) = data.edges[k];
                self.apply_two_qubit_gate(&Gate::Rzz(theta).matrix_4x4(), q0, q1)
            }
            Instruction::Gate {
                gate: Gate::Fused2q(mat),
                targets,
            } => self.apply_two_qubit_gate(mat, targets[0], targets[1]),
            Instruction::Gate { gate, targets } => {
                self.apply_two_qubit_gate(&gate.matrix_4x4(), targets[0], targets[1])
            }
            _ => unreachable!("an instruction without pair entries"),
        }
    }

    /// Steps [`Self::prepare_center`] walks to reach the pair at `left_site`.
    fn entry_steps(&self, left_site: usize) -> usize {
        match self.center {
            Some(at) if at <= left_site => left_site - at,
            Some(at) => at - left_site - 1,
            None => 0,
        }
    }

    /// Apply BatchPhase via bubble routing: sweep control toward targets,
    /// applying each CU-phase+SWAP as a single 2-site operation. O(k) tensor
    /// ops instead of O(k²) for k non-adjacent phases.
    fn apply_batch_phase_bubble(
        &mut self,
        control: usize,
        phases: &[(usize, Complex64)],
    ) -> Result<()> {
        if phases.is_empty() {
            return Ok(());
        }

        let mut right: Vec<(usize, Complex64)> = Vec::new();
        let mut left: Vec<(usize, Complex64)> = Vec::new();
        for &(target, phase) in phases {
            if target > control {
                right.push((target, phase));
            } else if target < control {
                left.push((target, phase));
            }
        }
        right.sort_by_key(|&(t, _)| t);
        left.sort_by_key(|&(t, _)| std::cmp::Reverse(t));

        let swap_mat = Gate::Swap.matrix_4x4();

        // Right sweep: bubble control rightward
        if !right.is_empty() {
            let mut cur_pos = control;
            let last_idx = right.len() - 1;
            for (idx, &(target, phase)) in right.iter().enumerate() {
                while cur_pos + 1 < target {
                    self.apply_adjacent_two_qubit(&swap_mat, cur_pos, true)?;
                    cur_pos += 1;
                }
                if idx < last_idx {
                    let combined = cu_phase_swap_matrix(phase);
                    self.apply_adjacent_two_qubit(&combined, cur_pos, true)?;
                    cur_pos += 1;
                } else {
                    let mat = [[ONE, ZERO], [ZERO, phase]];
                    let g = crate::gates::cu_matrix_4x4(&mat);
                    self.apply_adjacent_two_qubit(&g, cur_pos, true)?;
                }
            }
            while cur_pos > control {
                cur_pos -= 1;
                self.apply_adjacent_two_qubit(&swap_mat, cur_pos, true)?;
            }
        }

        // Left sweep: bubble control leftward
        if !left.is_empty() {
            let mut cur_pos = control;
            let last_idx = left.len() - 1;
            for (idx, &(target, phase)) in left.iter().enumerate() {
                while cur_pos > target + 1 {
                    cur_pos -= 1;
                    self.apply_adjacent_two_qubit(&swap_mat, cur_pos, true)?;
                }
                if idx < last_idx {
                    let combined = cu_phase_swap_matrix(phase);
                    self.apply_adjacent_two_qubit(&combined, cur_pos - 1, true)?;
                    cur_pos -= 1;
                } else {
                    let mat = [[ONE, ZERO], [ZERO, phase]];
                    let g = crate::gates::cu_matrix_4x4(&mat);
                    self.apply_adjacent_two_qubit(&g, cur_pos - 1, false)?;
                }
            }
            while cur_pos < control {
                self.apply_adjacent_two_qubit(&swap_mat, cur_pos, true)?;
                cur_pos += 1;
            }
        }
        Ok(())
    }

    /// Contract N adjacent site tensors into Θ with shape (χ_left, 2^N, χ_right).
    ///
    /// Returns (theta, bl, br) where theta is indexed as Θ[α, s, β] at
    /// `α * (dim * br) + s * br + β`, with dim = 2^n.
    fn contract_n_sites(&self, start: usize, n: usize) -> (Vec<Complex64>, usize, usize) {
        debug_assert!(n >= 2);
        let bl = self.sites[start].bond_left;
        let br = self.sites[start + n - 1].bond_right;

        let t0 = &self.sites[start];
        let mut theta = t0.data.clone();
        let mut cur_dim = 2usize;
        let mut cur_right = t0.bond_right;
        let mut new_theta = Vec::new();

        for k in 1..n {
            let t = &self.sites[start + k];
            let bond_mid = cur_right;
            let next_right = t.bond_right;
            let new_dim = cur_dim * 2;

            let nsite_chunk = new_dim * next_right;
            let needed = bl * nsite_chunk;
            new_theta.clear();
            new_theta.resize(needed, ZERO);
            let t_data = &t.data;
            let theta_ref = &theta;

            #[cfg(feature = "parallel")]
            if bl >= MIN_BOND_FOR_PAR {
                new_theta
                    .par_chunks_mut(nsite_chunk)
                    .enumerate()
                    .for_each(|(alpha, chunk)| {
                        for s in 0..cur_dim {
                            for j in 0..2 {
                                for beta in 0..next_right {
                                    let mut val = ZERO;
                                    for gamma in 0..bond_mid {
                                        val += theta_ref
                                            [alpha * (cur_dim * bond_mid) + s * bond_mid + gamma]
                                            * t_data
                                                [gamma * (2 * next_right) + j * next_right + beta];
                                    }
                                    let new_s = s * 2 + j;
                                    chunk[new_s * next_right + beta] = val;
                                }
                            }
                        }
                    });
            } else {
                for alpha in 0..bl {
                    for s in 0..cur_dim {
                        for j in 0..2 {
                            for beta in 0..next_right {
                                let mut val = ZERO;
                                for gamma in 0..bond_mid {
                                    val += theta_ref
                                        [alpha * (cur_dim * bond_mid) + s * bond_mid + gamma]
                                        * t_data[gamma * (2 * next_right) + j * next_right + beta];
                                }
                                let new_s = s * 2 + j;
                                new_theta[alpha * nsite_chunk + new_s * next_right + beta] = val;
                            }
                        }
                    }
                }
            }

            #[cfg(not(feature = "parallel"))]
            for alpha in 0..bl {
                for s in 0..cur_dim {
                    for j in 0..2 {
                        for beta in 0..next_right {
                            let mut val = ZERO;
                            for gamma in 0..bond_mid {
                                val += theta_ref
                                    [alpha * (cur_dim * bond_mid) + s * bond_mid + gamma]
                                    * t_data[gamma * (2 * next_right) + j * next_right + beta];
                            }
                            let new_s = s * 2 + j;
                            new_theta[alpha * nsite_chunk + new_s * next_right + beta] = val;
                        }
                    }
                }
            }
            std::mem::swap(&mut theta, &mut new_theta);
            cur_dim = new_dim;
            cur_right = next_right;
        }

        (theta, bl, br)
    }

    /// Apply a 2^N × 2^N gate matrix (row-major) to a contracted tensor in-place.
    fn apply_gate_to_theta(
        theta: &[Complex64],
        gate: &[Complex64],
        dim: usize,
        bl: usize,
        br: usize,
    ) -> Vec<Complex64> {
        let gate_chunk = dim * br;
        let mut out = vec![ZERO; bl * gate_chunk];

        #[cfg(feature = "parallel")]
        if bl >= MIN_BOND_FOR_PAR {
            out.par_chunks_mut(gate_chunk)
                .enumerate()
                .for_each(|(alpha, chunk)| {
                    for sp in 0..dim {
                        for beta in 0..br {
                            let mut val = ZERO;
                            for s in 0..dim {
                                val +=
                                    gate[sp * dim + s] * theta[alpha * gate_chunk + s * br + beta];
                            }
                            chunk[sp * br + beta] = val;
                        }
                    }
                });
            return out;
        }

        for alpha in 0..bl {
            for sp in 0..dim {
                for beta in 0..br {
                    let mut val = ZERO;
                    for s in 0..dim {
                        val += gate[sp * dim + s] * theta[alpha * gate_chunk + s * br + beta];
                    }
                    out[alpha * gate_chunk + sp * br + beta] = val;
                }
            }
        }
        out
    }

    /// Decompose a contracted tensor back into N site tensors via N-1 SVDs.
    ///
    /// Left-to-right sweep: at each step, reshape the remainder to
    /// (χ × 2) rows × (remaining_dim × χ_right) cols, SVD, truncate,
    /// extract left site from U, and continue with S·V† as the new remainder.
    fn decompose_n_sites(
        &mut self,
        theta: &[Complex64],
        start: usize,
        n: usize,
        bl: usize,
        br: usize,
        judge: bool,
    ) -> bool {
        let mut remaining = theta.to_vec();
        let mut cur_bl = bl;
        let mut remaining_dim = 1usize << n;
        let mut mat = Vec::new();
        let mut new_remaining = Vec::new();

        for k in 0..n - 1 {
            remaining_dim /= 2;
            let rows = cur_bl * 2;
            let cols = remaining_dim * br;

            // Reshape to column-major matrix for SVD: M[r, c] at mat[c * rows + r]
            // r = α * 2 + i, c = s_rest * br + β
            mat.clear();
            mat.resize(rows * cols, ZERO);
            #[cfg(feature = "parallel")]
            if cur_bl >= MIN_BOND_FOR_PAR * 2 {
                let rem = &remaining;
                let rd = remaining_dim;
                mat.par_iter_mut().enumerate().for_each(|(out_idx, elem)| {
                    let c = out_idx / rows;
                    let r = out_idx % rows;
                    let alpha = r / 2;
                    let i = r & 1;
                    let s_rest = c / br;
                    let beta = c % br;
                    let combined_s = i * rd + s_rest;
                    *elem = rem[alpha * (2 * rd * br) + combined_s * br + beta];
                });
            } else {
                for alpha in 0..cur_bl {
                    for i in 0..2 {
                        for s_rest in 0..remaining_dim {
                            for beta in 0..br {
                                let r = alpha * 2 + i;
                                let c = s_rest * br + beta;
                                let combined_s = i * remaining_dim + s_rest;
                                mat[c * rows + r] = remaining
                                    [alpha * (2 * remaining_dim * br) + combined_s * br + beta];
                            }
                        }
                    }
                }
            }

            #[cfg(not(feature = "parallel"))]
            for alpha in 0..cur_bl {
                for i in 0..2 {
                    for s_rest in 0..remaining_dim {
                        for beta in 0..br {
                            let r = alpha * 2 + i;
                            let c = s_rest * br + beta;
                            let combined_s = i * remaining_dim + s_rest;
                            mat[c * rows + r] = remaining
                                [alpha * (2 * remaining_dim * br) + combined_s * br + beta];
                        }
                    }
                }
            }

            let svd_result = svd(&mat, rows, cols);
            let chi_new = truncated_svd_rank(&svd_result.s, self.svd_epsilon, self.max_bond_dim);
            if judge && discards_resolvable_weight(&svd_result.s, chi_new) {
                return false;
            }
            self.record_truncation(&svd_result.s, chi_new);
            self.bond_high_water = self.bond_high_water.max(chi_new);

            // Extract left site from U: shape (cur_bl, 2, chi_new)
            let mut left_data = std::mem::take(&mut self.sites[start + k].data);
            fill_isometry_site_data(
                &mut left_data,
                &svd_result.u,
                svd_result.u_rows,
                cur_bl,
                chi_new,
            );

            self.sites[start + k] = SiteTensor {
                bond_left: cur_bl,
                bond_right: chi_new,
                data: left_data,
            };

            if k < n - 2 {
                // Build remainder from S·V†: shape (chi_new, remaining_dim, br)
                fill_scaled_vt_data(&mut new_remaining, &svd_result, chi_new, remaining_dim, br);
                std::mem::swap(&mut remaining, &mut new_remaining);
                cur_bl = chi_new;
            } else {
                // Last site: S·V† reshaped to (chi_new, 2, br)
                let mut right_data = std::mem::take(&mut self.sites[start + n - 1].data);
                fill_scaled_vt_data(&mut right_data, &svd_result, chi_new, 2, br);
                self.sites[start + n - 1] = SiteTensor {
                    bond_left: chi_new,
                    bond_right: br,
                    data: right_data,
                };
            }
        }
        true
    }

    /// Reject a dense `2^n x 2^n` gate matrix over the workspace budget before
    /// it is built.
    ///
    /// The N-site path holds two of them at peak, the assembled matrix and the
    /// role-reordered copy, and both are allocated before
    /// [`Self::apply_adjacent_n_qubit`] reaches its own theta check.
    fn check_n_qubit_gate_matrix(&self, n: usize) -> Result<()> {
        // Rejected outright rather than by comparison: `2^(2n + 1)` is not
        // representable here, and the cap is itself `u128::MAX` when memory
        // detection is unavailable, so a comparison would pass and leave the
        // caller to shift `1usize` by more than its width.
        if n >= 63 {
            return Err(crate::backend::workspace_allocation_error(
                "mps",
                "multi-qubit gate matrix",
                u128::MAX,
            ));
        }
        let elements = 2u128 << (2 * n);
        if elements > self.workspace_cap {
            return Err(crate::backend::workspace_allocation_error(
                "mps",
                "multi-qubit gate matrix",
                elements,
            ));
        }
        Ok(())
    }

    fn apply_adjacent_n_qubit(
        &mut self,
        gate: &[Complex64],
        dim: usize,
        start_site: usize,
    ) -> Result<()> {
        let n = dim.trailing_zeros() as usize; // dim = 2^n
        // Contraction, gate application, and decomposition each hold a
        // (bl, 2^n, br) tensor, with the decomposition remainder and SVD input
        // making about four live at peak.
        let bl = self.sites[start_site].bond_left;
        let br = self.sites[start_site + n - 1].bond_right;
        let workspace = 4u128 * dim as u128 * bl as u128 * br as u128;
        if workspace > self.workspace_cap {
            return Err(crate::backend::workspace_allocation_error(
                "mps",
                "multi-qubit gate application",
                workspace,
            ));
        }
        if !self.cut_block(gate, dim, start_site, n) {
            self.establish_center(start_site);
            let gauged = self.cut_block(gate, dim, start_site, n);
            debug_assert!(gauged, "a block cut in gauge judged its spectrum");
        }
        Ok(())
    }

    /// One attempt at [`Self::apply_adjacent_n_qubit`]. Returns `false`,
    /// having put the block's sites back, when the chain has no center and a
    /// cut of the decomposition drops a value the factorization resolves.
    ///
    /// The decomposition sweeps left to right and leaves the weight on the
    /// block's last site, so a center anywhere in the block ends there. One
    /// outside it would leave the sites between non-orthogonal, and the
    /// block's own cuts would be the ones paying for it. The two-site kernel
    /// judges its cut before writing anything, where the spectrum that shows
    /// it here is only reached once the sites before that cut are written, so
    /// the retry restores them from a copy. Only a threshold can drop a value
    /// on the ungauged arm, since a cut that could reach the cap gauges
    /// instead, so at a threshold of zero there is nothing to judge and no
    /// copy is taken.
    fn cut_block(&mut self, gate: &[Complex64], dim: usize, start_site: usize, n: usize) -> bool {
        let block = start_site..start_site + n;
        let bl = self.sites[start_site].bond_left;
        let br = self.sites[block.end - 1].bond_right;
        let maintained = self.center.is_some() || self.cut_can_lose(widest_block_cut(bl, br, n));
        if maintained {
            match self.center {
                Some(at) => self.move_center(at.clamp(start_site, block.end - 1)),
                None => self.establish_center(start_site),
            }
        }
        let judging = !maintained && self.svd_epsilon > 0.0;

        let (theta, bl, br) = self.contract_n_sites(start_site, n);
        let theta_prime = Self::apply_gate_to_theta(&theta, gate, dim, bl, br);
        if !judging {
            self.decompose_n_sites(&theta_prime, start_site, n, bl, br, false);
            if maintained {
                self.center = Some(block.end - 1);
            }
            return true;
        }

        let saved = self.sites[block.clone()].to_vec();
        let booked = self.truncation_discarded;
        let high_water = self.bond_high_water;
        if self.decompose_n_sites(&theta_prime, start_site, n, bl, br, true) {
            return true;
        }
        self.sites[block].clone_from_slice(&saved);
        self.truncation_discarded = booked;
        self.bond_high_water = high_water;
        false
    }

    /// Apply an N-qubit gate to arbitrary (possibly non-adjacent) qubits.
    ///
    /// SWAP-routes qubits into a contiguous block, applies the gate, then
    /// reverses the SWAPs.
    fn apply_n_qubit_gate(
        &mut self,
        gate: &[Complex64],
        dim: usize,
        qubits: &[usize],
    ) -> Result<()> {
        let n = qubits.len();
        debug_assert_eq!(dim, 1 << n);

        let mut indexed: Vec<(usize, usize)> = qubits
            .iter()
            .enumerate()
            .map(|(role, &pos)| (pos, role))
            .collect();
        indexed.sort_unstable_by_key(|&(pos, _)| pos);

        let sorted_positions: Vec<usize> = indexed.iter().map(|&(pos, _)| pos).collect();
        let start = sorted_positions[0];

        let contiguous = sorted_positions.last().unwrap() - start + 1 == n;

        if contiguous {
            let qubit_order: Vec<usize> = indexed.iter().map(|&(_, role)| role).collect();
            let reordered_gate = Self::reorder_n_gate(gate, dim, &qubit_order, n);
            self.apply_adjacent_n_qubit(&reordered_gate, dim, start)?;
        } else {
            let swap_mat = Gate::Swap.matrix_4x4();
            let mut current_positions = sorted_positions;
            let mut swap_log: Vec<usize> = Vec::new();

            for i in 1..n {
                let target_pos = start + i;
                while current_positions[i] > target_pos {
                    let s = current_positions[i] - 1;
                    self.apply_adjacent_two_qubit(&swap_mat, s, true)?;
                    swap_log.push(s);
                    // Update any tracked qubit that was displaced
                    for cp in &mut current_positions[(i + 1)..n] {
                        if *cp == s {
                            *cp = s + 1;
                        }
                    }
                    current_positions[i] -= 1;
                }
            }

            let qubit_order: Vec<usize> = indexed.iter().map(|&(_, role)| role).collect();
            let reordered_gate = Self::reorder_n_gate(gate, dim, &qubit_order, n);
            self.apply_adjacent_n_qubit(&reordered_gate, dim, start)?;

            for &s in swap_log.iter().rev() {
                self.apply_adjacent_two_qubit(&swap_mat, s, true)?;
            }
        }
        Ok(())
    }

    /// Reorder an N-qubit gate matrix to match the physical qubit ordering.
    ///
    /// `qubit_order[block_pos]` = the MCU role index for that block position.
    /// The original gate assumes role order 0,1,...,N-1 (controls then target).
    fn reorder_n_gate(
        gate: &[Complex64],
        dim: usize,
        qubit_order: &[usize],
        n: usize,
    ) -> Vec<Complex64> {
        let identity_order = qubit_order.iter().enumerate().all(|(i, &r)| r == i);
        if identity_order {
            return gate.to_vec();
        }

        let mut out = vec![ZERO; dim * dim];
        for s_in in 0..dim {
            for s_out in 0..dim {
                let mut logical_in = 0usize;
                let mut logical_out = 0usize;
                for (pos, &role) in qubit_order.iter().enumerate() {
                    let bit_pos = n - 1 - pos; // bit position for block pos
                    let role_bit = n - 1 - role; // bit position in role order
                    if s_in & (1 << bit_pos) != 0 {
                        logical_in |= 1 << role_bit;
                    }
                    if s_out & (1 << bit_pos) != 0 {
                        logical_out |= 1 << role_bit;
                    }
                }
                out[s_out * dim + s_in] = gate[logical_out * dim + logical_in];
            }
        }
        out
    }

    fn compute_left_env(&self, site: usize) -> Vec<Complex64> {
        let mut env = vec![ONE];
        let mut env_dim = 1usize;

        for s in 0..site {
            let t = &self.sites[s];
            let bl = t.bond_left;
            let br = t.bond_right;
            debug_assert_eq!(env_dim, bl);

            #[cfg(feature = "parallel")]
            if br >= MIN_BOND_FOR_PAR {
                let mut new_env = vec![ZERO; br * br];
                new_env.par_iter_mut().enumerate().for_each(|(idx, val)| {
                    let gamma = idx / br;
                    let gamma_p = idx % br;
                    let mut sum = ZERO;
                    for alpha in 0..bl {
                        for alpha_p in 0..bl {
                            let env_val = env[alpha * bl + alpha_p];
                            if env_val == ZERO {
                                continue;
                            }
                            for i in 0..2 {
                                sum += env_val
                                    * t.data[t.idx(alpha, i, gamma)]
                                    * t.data[t.idx(alpha_p, i, gamma_p)].conj();
                            }
                        }
                    }
                    *val = sum;
                });
                env = new_env;
                env_dim = br;
                continue;
            }

            let mut new_env = vec![ZERO; br * br];
            for gamma in 0..br {
                for gamma_p in 0..br {
                    let mut val = ZERO;
                    for alpha in 0..bl {
                        for alpha_p in 0..bl {
                            let env_val = env[alpha * bl + alpha_p];
                            if env_val == ZERO {
                                continue;
                            }
                            for i in 0..2 {
                                val += env_val
                                    * t.data[t.idx(alpha, i, gamma)]
                                    * t.data[t.idx(alpha_p, i, gamma_p)].conj();
                            }
                        }
                    }
                    new_env[gamma * br + gamma_p] = val;
                }
            }
            env = new_env;
            env_dim = br;
        }
        env
    }

    /// Absorb one site into a right environment, in the two passes
    /// [`Self::inner_product`] uses: `O(χ³)` instead of the four-index sum's
    /// `O(χ⁴)`. `tmp` carries the intermediate across a sweep so the inner
    /// buffer is allocated once rather than per site.
    ///
    /// [`Self::compute_right_env`] deliberately keeps the four-index form
    /// instead of calling this. It serves the single-site reduced density
    /// matrix, whose chains carry the zeros a projection leaves behind, and the
    /// `env_val == ZERO` skip there beat dense `O(χ³)` arithmetic by 14% on
    /// `mps/hotspots/measure_reset_32q_r3` (496 µs against 574 µs) while that
    /// row still went through the environments. The sweep below runs on
    /// unprojected chains where the environments are dense and the asymptotics
    /// decide instead. Measure both before unifying them.
    fn contract_right_env(
        &self,
        site: usize,
        env: &[Complex64],
        tmp: &mut Vec<Complex64>,
    ) -> Vec<Complex64> {
        let t = &self.sites[site];
        let bl = t.bond_left;
        let br = t.bond_right;

        // tmp[α, i, β'] = Σ_β A[α, i, β] · env[β, β']
        tmp.clear();
        tmp.resize(bl * 2 * br, ZERO);
        for alpha in 0..bl {
            for i in 0..2 {
                for beta in 0..br {
                    let a = t.data[t.idx(alpha, i, beta)];
                    if a == ZERO {
                        continue;
                    }
                    let dst = (alpha * 2 + i) * br;
                    let src = beta * br;
                    for beta_p in 0..br {
                        tmp[dst + beta_p] += a * env[src + beta_p];
                    }
                }
            }
        }

        // new_env[α, α'] = Σ_{i, β'} tmp[α, i, β'] · conj(A[α', i, β'])
        let mut new_env = vec![ZERO; bl * bl];
        for alpha_p in 0..bl {
            for i in 0..2 {
                for beta_p in 0..br {
                    let a = t.data[t.idx(alpha_p, i, beta_p)].conj();
                    if a == ZERO {
                        continue;
                    }
                    for alpha in 0..bl {
                        new_env[alpha * bl + alpha_p] += tmp[(alpha * 2 + i) * br + beta_p] * a;
                    }
                }
            }
        }
        new_env
    }

    fn compute_right_env(&self, site: usize) -> Vec<Complex64> {
        let mut env = vec![ONE];
        let mut env_dim = 1usize;

        for s in (site + 1..self.num_qubits).rev() {
            let t = &self.sites[s];
            let bl = t.bond_left;
            let br = t.bond_right;
            debug_assert_eq!(env_dim, br);

            #[cfg(feature = "parallel")]
            if bl >= MIN_BOND_FOR_PAR {
                let mut new_env = vec![ZERO; bl * bl];
                new_env.par_iter_mut().enumerate().for_each(|(idx, val)| {
                    let alpha = idx / bl;
                    let alpha_p = idx % bl;
                    let mut sum = ZERO;
                    for beta in 0..br {
                        for beta_p in 0..br {
                            let env_val = env[beta * br + beta_p];
                            if env_val == ZERO {
                                continue;
                            }
                            for i in 0..2 {
                                sum += t.data[t.idx(alpha, i, beta)]
                                    * t.data[t.idx(alpha_p, i, beta_p)].conj()
                                    * env_val;
                            }
                        }
                    }
                    *val = sum;
                });
                env = new_env;
                env_dim = bl;
                continue;
            }

            let mut new_env = vec![ZERO; bl * bl];
            for alpha in 0..bl {
                for alpha_p in 0..bl {
                    let mut val = ZERO;
                    for beta in 0..br {
                        for beta_p in 0..br {
                            let env_val = env[beta * br + beta_p];
                            if env_val == ZERO {
                                continue;
                            }
                            for i in 0..2 {
                                val += t.data[t.idx(alpha, i, beta)]
                                    * t.data[t.idx(alpha_p, i, beta_p)].conj()
                                    * env_val;
                            }
                        }
                    }
                    new_env[alpha * bl + alpha_p] = val;
                }
            }
            env = new_env;
            env_dim = bl;
        }
        env
    }

    /// Right environment of every site, in one backward sweep.
    ///
    /// `envs[s]` is what [`Self::compute_right_env`] returns for `s`, but the
    /// sweep shares each partial contraction with the site to its left, so all
    /// `n` cost one pass instead of `n`.
    fn right_environments(&self) -> Vec<Vec<Complex64>> {
        let n = self.num_qubits;
        let mut envs: Vec<Vec<Complex64>> = Vec::with_capacity(n);
        let mut tmp = Vec::new();
        envs.push(vec![ONE]);
        for site in (1..n).rev() {
            let next = self.contract_right_env(site, envs.last().unwrap(), &mut tmp);
            envs.push(next);
        }
        envs.reverse();
        envs
    }

    fn apply_reset(&mut self, qubit: usize) {
        let prob = self.born_weights(qubit);
        let outcome = usize::from(self.rng.random::<f64>() < prob[1].clamp(0.0, 1.0));
        self.collapse_site_to_zero(qubit, outcome, prob[outcome].clamp(0.0, 1.0));
    }

    fn apply_measure(&mut self, qubit: usize, classical_bit: usize) {
        let prob = self.born_weights(qubit);

        let measured = usize::from(self.rng.random::<f64>() < prob[1].clamp(0.0, 1.0));
        self.classical_bits[classical_bit] = measured == 1;
        self.project_site_outcome(qubit, measured, prob[measured].clamp(0.0, 1.0));
    }

    /// Unnormalized weight of each outcome at `site`, taken with the
    /// orthogonality center walked onto it.
    ///
    /// Every other site is then an isometry and both environments contract to
    /// identities, which leaves the weights as the site's own row sums: `O(χ²)`
    /// against the `O(n·χ⁴)` of the two environment sweeps a site off the
    /// center needs. The projection that follows writes `site` alone, the one
    /// site under no isometry claim, so the center survives it and the next
    /// measurement pays the distance between the two sites rather than a
    /// rebuild from the end of the chain. A chain
    /// [`Self::cuts_can_lose_weight`] rules out keeps no record past the
    /// projection: the walk is still the cheaper way to the weights, and the
    /// gauge would only tax the gates that follow.
    fn born_weights(&mut self, site: usize) -> [f64; 2] {
        match self.center {
            Some(_) => self.move_center(site),
            None => self.establish_center(site),
        }

        let t = &self.sites[site];
        let br = t.bond_right;
        let mut prob = [0.0f64; 2];
        for alpha in 0..t.bond_left {
            let row = &t.data[alpha * (2 * br)..(alpha + 1) * (2 * br)];
            for (outcome, weight) in prob.iter_mut().enumerate() {
                *weight += row[outcome * br..(outcome + 1) * br]
                    .iter()
                    .map(|x| x.norm_sqr())
                    .sum::<f64>();
            }
        }
        prob
    }

    /// Take the Born weight from the caller, which read it off the center in
    /// [`Self::born_weights`], the convention [`Self::collapse_site_to_zero`]
    /// already uses for reset. The write lands on the center, so the gauge
    /// comes through it intact.
    fn project_site_outcome(&mut self, site: usize, outcome: usize, prob: f64) -> f64 {
        debug_assert_eq!(self.center, Some(site), "a projection off the center");
        if !self.cuts_can_lose_weight() {
            self.center = None;
        }
        if prob <= NORM_CLAMP_MIN {
            return 0.0;
        }
        let inv_sqrt_prob = 1.0 / prob.sqrt();
        let scale = Complex64::new(inv_sqrt_prob, 0.0);
        let other = 1 - outcome;

        let t = &mut self.sites[site];
        let bl = t.bond_left;
        let br = t.bond_right;
        for alpha in 0..bl {
            for beta in 0..br {
                let idx_m = alpha * (2 * br) + outcome * br + beta;
                let idx_o = alpha * (2 * br) + other * br + beta;
                t.data[idx_m] *= scale;
                t.data[idx_o] = ZERO;
            }
        }
        prob
    }

    /// Collapse a site onto `outcome` and return it to the zero state.
    ///
    /// Writing the surviving component straight into the zero slot folds in the
    /// X that a trajectory reset would otherwise apply after projecting onto
    /// outcome 1, and takes the Born weight from the caller, which read it off
    /// the center in [`Self::born_weights`].
    fn collapse_site_to_zero(&mut self, site: usize, outcome: usize, prob: f64) {
        debug_assert_eq!(self.center, Some(site), "a reset off the center");
        if !self.cuts_can_lose_weight() {
            self.center = None;
        }
        if prob <= NORM_CLAMP_MIN {
            return;
        }
        let scale = Complex64::new(1.0 / prob.sqrt(), 0.0);
        let t = &mut self.sites[site];
        let br = t.bond_right;
        for alpha in 0..t.bond_left {
            let base = alpha * (2 * br);
            for beta in 0..br {
                let zero_idx = base + beta;
                t.data[zero_idx] = t.data[base + outcome * br + beta] * scale;
                t.data[zero_idx + br] = ZERO;
            }
        }
    }

    /// Project a logical qubit onto a requested computational-basis outcome.
    ///
    /// Returns the branch-local Born probability. A zero return means the
    /// requested branch should be discarded by the caller.
    pub(crate) fn project_z_outcome(&mut self, qubit: usize, outcome: bool) -> f64 {
        let site = self.site_for_logical(qubit);
        let outcome = usize::from(outcome);
        let prob = self.born_weights(site)[outcome].clamp(0.0, 1.0);
        self.project_site_outcome(site, outcome, prob)
    }

    fn reduced_density_site(&self, site: usize) -> [[Complex64; 2]; 2] {
        let l_env = self.compute_left_env(site);
        let r_env = self.compute_right_env(site);
        let t = &self.sites[site];
        let bl = t.bond_left;
        let br = t.bond_right;
        let mut rho = [[ZERO; 2]; 2];

        for (row, rho_row) in rho.iter_mut().enumerate() {
            for (col, rho_cell) in rho_row.iter_mut().enumerate() {
                let mut val = ZERO;
                for alpha in 0..bl {
                    for alpha_p in 0..bl {
                        let l_val = l_env[alpha * bl + alpha_p];
                        if l_val == ZERO {
                            continue;
                        }
                        for beta in 0..br {
                            for beta_p in 0..br {
                                let r_val = r_env[beta * br + beta_p];
                                if r_val == ZERO {
                                    continue;
                                }
                                val += l_val
                                    * t.data[t.idx(alpha, row, beta)]
                                    * t.data[t.idx(alpha_p, col, beta_p)].conj()
                                    * r_val;
                            }
                        }
                    }
                }
                *rho_cell = val;
            }
        }

        rho
    }

    /// Schmidt values across the bond between sites `bond` and `bond + 1`.
    ///
    /// With the center on either site of the bond, every other site is an
    /// isometry, so the singular values of the center site read as a matrix
    /// across that bond are the spectrum of the whole state. The site's data
    /// is already that matrix transposed in column-major order, and the
    /// transpose has the same singular values, so nothing is reshaped. The
    /// center walks to whichever site of the bond is nearer, or is
    /// established at the one nearer its end of the chain when none is
    /// recorded; the walk is exact and books nothing.
    fn schmidt_values_at_bond(&mut self, bond: usize) -> Vec<f64> {
        let last = self.sites.len() - 1;
        let site = match self.center {
            Some(at) if at > bond => bond + 1,
            Some(_) => bond,
            None if 2 * bond < last => bond,
            None => bond + 1,
        };
        match self.center {
            Some(_) => self.move_center(site),
            None => self.establish_center(site),
        }
        let t = &self.sites[site];
        let (rows, cols) = if site == bond {
            (t.bond_right, 2 * t.bond_left)
        } else {
            (2 * t.bond_right, t.bond_left)
        };
        schmidt::finish_schmidt_values(svd(&t.data, rows, cols).s)
    }

    /// Schmidt values from the eigenvalues of the reduced density matrix over
    /// the chain sites `sites`, ascending and at most half the chain. The
    /// center moves to the nearest site of the span so both environments are
    /// identities.
    fn schmidt_values_by_reduced_density(&mut self, sites: &[usize]) -> Result<Vec<f64>> {
        let what = "Schmidt values across a cut that is not contiguous in chain order";
        let k = sites.len();
        schmidt::check_schmidt_side(self.name(), what, k)?;
        let (first, last) = (sites[0], sites[k - 1]);
        let chi = self.sites[first..=last]
            .iter()
            .fold(1usize, |chi, t| chi.max(t.bond_left).max(t.bond_right));
        // The transfer object, its half-contracted form, and its successor at
        // the widest bond of the span, each over `4^k` open index pairs.
        let workspace = (7u128 * (chi as u128) * (chi as u128)) << (2 * k);
        if workspace > self.workspace_cap {
            return Err(crate::backend::workspace_allocation_error(
                self.name(),
                what,
                workspace,
            ));
        }
        let end = self.sites.len() - 1;
        match self.center {
            Some(at) => self.move_center(at.clamp(first, last)),
            None if first <= end - last => self.establish_center(first),
            None => self.establish_center(last),
        }
        let rho = self.reduced_density_sites(sites);
        Ok(schmidt::schmidt_values_from_density(&rho, 1 << k))
    }

    /// Reduced density matrix over the chain sites `sites` (ascending), row
    /// major `2^k x 2^k` with the leftmost site in the high bit of both
    /// indices. Requires the center inside `sites[0]..=sites[k - 1]`: the
    /// sites on either side of that span then contract to identities, and the
    /// sweep over the span carries a `(ket, bra, bond, bond)` transfer object
    /// that grows fourfold at each site it keeps open.
    fn reduced_density_sites(&self, sites: &[usize]) -> Vec<Complex64> {
        let (first, last) = (sites[0], sites[sites.len() - 1]);
        let bl = self.sites[first].bond_left;
        let mut env: Vec<Complex64> = (0..bl * bl)
            .map(|x| if x / bl == x % bl { ONE } else { ZERO })
            .collect();
        let mut dim = 1usize;
        let mut kept = sites.iter().peekable();
        for site in first..=last {
            let keep = kept.next_if_eq(&&site).is_some();
            let t = &self.sites[site];
            let (bl, br) = (t.bond_left, t.bond_right);
            let width = 2 * br;
            // half[o, a', i, b] = sum_a env[o, a, a'] A[a, i, b]
            let mut half = vec![ZERO; dim * dim * bl * width];
            for (o, rows) in half.chunks_exact_mut(bl * width).enumerate() {
                for (ap, out) in rows.chunks_exact_mut(width).enumerate() {
                    for a in 0..bl {
                        let e = env[(o * bl + a) * bl + ap];
                        for (x, amp) in out.iter_mut().zip(&t.data[a * width..(a + 1) * width]) {
                            *x += e * amp;
                        }
                    }
                }
            }
            // next[o', b, b'] = sum_{a', i, i'} half[o, a', i, b] conj(A[a', i', b'])
            // with `i == i'` traced, or `i` and `i'` appended to `o` as the
            // low bit of the ket and bra when the site is kept.
            let next_dim = if keep { 2 * dim } else { dim };
            let mut next = vec![ZERO; next_dim * next_dim * br * br];
            for (o, rows) in half.chunks_exact(bl * width).enumerate() {
                let (ket, bra) = (o / dim, o % dim);
                for (ap, row) in rows.chunks_exact(width).enumerate() {
                    for i in 0..2 {
                        let ket_row = &row[i * br..(i + 1) * br];
                        for ip in 0..2 {
                            if !keep && ip != i {
                                continue;
                            }
                            let o_next = if keep {
                                (2 * ket + i) * next_dim + 2 * bra + ip
                            } else {
                                o
                            };
                            let bra_row = &t.data[ap * width + ip * br..ap * width + (ip + 1) * br];
                            let dst = &mut next[o_next * br * br..(o_next + 1) * br * br];
                            for (b, &x) in ket_row.iter().enumerate() {
                                for (y, d) in bra_row.iter().zip(&mut dst[b * br..(b + 1) * br]) {
                                    *d += x * y.conj();
                                }
                            }
                        }
                    }
                }
            }
            env = next;
            dim = next_dim;
        }
        let br = self.sites[last].bond_right;
        env.chunks_exact(br * br)
            .map(|block| (0..br).map(|b| block[b * br + b]).sum())
            .collect()
    }

    /// Conditional outcome weights at `site`, given the boundary vector `left`
    /// that carries the bits already fixed to its left.
    ///
    /// Writes `w[i * bond_right + β] = Σ_α left[α] · A[α, i, β]` into `w` and
    /// returns the unnormalized weight of each outcome. The conditioned left
    /// environment is the rank-one `left[α]·conj(left[α'])`, which collapses
    /// the four-index contraction a general environment owes to `O(χ²)`.
    fn site_conditional_weights(
        &self,
        site: usize,
        left: &[Complex64],
        right: &[Complex64],
        w: &mut [Complex64],
    ) -> [f64; 2] {
        let t = &self.sites[site];
        let bl = t.bond_left;
        let br = t.bond_right;

        let w = &mut w[..2 * br];
        w.fill(ZERO);
        for (alpha, &l) in left[..bl].iter().enumerate() {
            if l == ZERO {
                continue;
            }
            let row = &t.data[alpha * (2 * br)..(alpha + 1) * (2 * br)];
            for (dst, &src) in w.iter_mut().zip(row) {
                *dst += l * src;
            }
        }

        let mut prob = [0.0f64; 2];
        for (i, prob_out) in prob.iter_mut().enumerate() {
            let wi = &w[i * br..(i + 1) * br];
            let mut acc = ZERO;
            for (beta, &w_beta) in wi.iter().enumerate() {
                if w_beta == ZERO {
                    continue;
                }
                let row = &right[beta * br..(beta + 1) * br];
                let inner: Complex64 = row.iter().zip(wi).map(|(&r, &v)| r * v.conj()).sum();
                acc += w_beta * inner;
            }
            *prob_out = acc.re;
        }
        prob
    }

    /// Contract site `site` against the right-side vector `src`, selecting
    /// physical index `bit`, into `dst` of length `bond_left`.
    #[inline(always)]
    fn contract_site_right(
        &self,
        site: usize,
        bit: usize,
        src: &[Complex64],
        dst: &mut [Complex64],
    ) {
        let t = &self.sites[site];
        debug_assert_eq!(src.len(), t.bond_right);
        debug_assert_eq!(dst.len(), t.bond_left);
        for (alpha, slot) in dst.iter_mut().enumerate() {
            let mut acc = ZERO;
            for (beta, &v) in src.iter().enumerate() {
                acc += t.data[t.idx(alpha, bit, beta)] * v;
            }
            *slot = acc;
        }
    }

    /// Depth-first expansion of the subtree below `depth`, writing one
    /// amplitude per leaf through `store`.
    ///
    /// Sites are consumed right to left, so `scratch[d]` holds the vector on
    /// the left bond of site `n - d` and the deepest level decides site 0.
    /// With an unpermuted layout that makes consecutive leaves differ in the
    /// lowest basis bit, so leaf writes run in address order.
    fn expand_subtree(
        &self,
        depth: usize,
        basis: usize,
        scratch: &mut [Vec<Complex64>],
        store: &mut impl FnMut(usize, Complex64),
    ) {
        let n = self.num_qubits;
        if depth == n {
            store(basis, scratch[n][0]);
            return;
        }
        let site = n - 1 - depth;
        let logical = self.site_to_logical[site];
        for bit in 0..2 {
            let (lo, hi) = scratch.split_at_mut(depth + 1);
            self.contract_site_right(site, bit, &lo[depth], &mut hi[0]);
            self.expand_subtree(depth + 1, basis | (bit << logical), scratch, store);
        }
    }

    /// Scratch stack for [`MpsBackend::expand_subtree`], one buffer per depth.
    fn expansion_scratch(&self) -> Vec<Vec<Complex64>> {
        let n = self.num_qubits;
        let mut scratch = Vec::with_capacity(n + 1);
        scratch.push(vec![ONE]);
        for depth in 1..=n {
            scratch.push(vec![ZERO; self.sites[n - depth].bond_left]);
        }
        scratch
    }

    /// Every `(basis prefix, left-bond vector)` pair at `depth`, in DFS order,
    /// grown from `root`.
    #[cfg(feature = "parallel")]
    fn expansion_prefixes(&self, depth: usize, root: Complex64) -> Vec<(usize, Vec<Complex64>)> {
        let mut scratch = self.expansion_scratch();
        scratch[0][0] = root;
        let mut out = Vec::with_capacity(1 << depth);
        self.collect_prefixes(0, depth, 0, &mut scratch, &mut out);
        out
    }

    #[cfg(feature = "parallel")]
    fn collect_prefixes(
        &self,
        depth: usize,
        stop: usize,
        basis: usize,
        scratch: &mut [Vec<Complex64>],
        out: &mut Vec<(usize, Vec<Complex64>)>,
    ) {
        if depth == stop {
            out.push((basis, scratch[depth].clone()));
            return;
        }
        let site = self.num_qubits - 1 - depth;
        let logical = self.site_to_logical[site];
        for bit in 0..2 {
            let (lo, hi) = scratch.split_at_mut(depth + 1);
            self.contract_site_right(site, bit, &lo[depth], &mut hi[0]);
            self.collect_prefixes(depth + 1, stop, basis | (bit << logical), scratch, out);
        }
    }

    /// Fill `out` with `convert(amplitude)` per basis index, sharing every
    /// prefix contraction across the `2^(n-k)` basis states below it.
    ///
    /// [`MpsBackend::chain_amplitude`] recomputes that prefix per amplitude and
    /// allocates once per site, so a dense fill costs `2^n * sum_k bl_k*br_k`
    /// against `sum_k 2^(k+1)*bl_k*br_k` here. The parallel split fixes the top
    /// `m` sites, whose basis-bit sets are disjoint across prefixes.
    ///
    /// The root vector carries `1/sqrt(⟨ψ|ψ⟩)`, so a truncated chain expands to
    /// the normalized state at no cost per leaf.
    fn expand_dense<T: Copy + Send>(&self, out: &mut [T], convert: impl Fn(Complex64) -> T + Sync) {
        let root = self.normalized_root();
        #[cfg(feature = "parallel")]
        {
            let split = self
                .num_qubits
                .saturating_sub(MIN_LEAVES_PER_TASK_BITS)
                .min(MAX_SPLIT_BITS);
            if split > 0 {
                let ptr = DenseOutPtr(out.as_mut_ptr());
                self.expansion_prefixes(split, root)
                    .into_par_iter()
                    .for_each(|(basis, vector)| {
                        let mut scratch = self.expansion_scratch();
                        scratch[split].copy_from_slice(&vector);
                        self.expand_subtree(split, basis, &mut scratch, &mut |idx, amp| {
                            // SAFETY: this task owns every basis index whose top
                            // `split` site bits match its prefix, and those sets
                            // partition the output, so no two tasks reach the
                            // same element.
                            unsafe { ptr.store(idx, convert(amp)) }
                        });
                    });
                return;
            }
        }

        let mut scratch = self.expansion_scratch();
        scratch[0][0] = root;
        self.expand_subtree(0, 0, &mut scratch, &mut |idx, amp| out[idx] = convert(amp));
    }

    /// A zero-norm chain keeps the unit root, the same guard
    /// [`Backend::pauli_expectations`] applies.
    fn normalized_root(&self) -> Complex64 {
        let norm = self.pauli_env_prefixes()[self.num_qubits][0].re;
        if norm > 0.0 {
            Complex64::new(1.0 / norm.sqrt(), 0.0)
        } else {
            ONE
        }
    }

    fn dispatch_gate(&mut self, gate: &Gate, targets: &[usize]) -> Result<()> {
        match gate {
            Gate::Rzz(_) => {
                let g = gate.matrix_4x4();
                self.apply_two_qubit_gate(&g, targets[0], targets[1])?;
            }
            Gate::Cx => {
                let g = gate.matrix_4x4();
                self.apply_two_qubit_gate(&g, targets[0], targets[1])?;
            }
            Gate::Cz => {
                let g = gate.matrix_4x4();
                self.apply_two_qubit_gate(&g, targets[0], targets[1])?;
            }
            Gate::Swap => {
                let g = Gate::Swap.matrix_4x4();
                self.apply_two_qubit_gate(&g, targets[0], targets[1])?;
            }
            Gate::Cu(_) => {
                let g = gate.matrix_4x4();
                self.apply_two_qubit_gate(&g, targets[0], targets[1])?;
            }
            Gate::Mcu(data) => {
                let num_ctrl = data.num_controls as usize;
                let n = num_ctrl + 1;
                self.check_n_qubit_gate_matrix(n)?;
                let all_qubits: Vec<usize> = targets
                    .iter()
                    .map(|&qubit| self.site_for_logical(qubit))
                    .collect();
                let dim = 1usize << n;
                let role_order: Vec<usize> = (0..n).collect();
                let gate_mat = mcu_matrix(num_ctrl, &data.mat, &role_order);
                self.apply_n_qubit_gate(&gate_mat, dim, &all_qubits)?;
            }
            Gate::Unitary(data) => {
                let n = data.num_qubits();
                self.check_n_qubit_gate_matrix(n)?;
                let sites: Vec<usize> = targets
                    .iter()
                    .map(|&qubit| self.site_for_logical(qubit))
                    .collect();
                self.apply_n_qubit_gate(data.matrix(), 1 << n, &sites)?;
            }
            Gate::BatchPhase(data) => {
                if let [(target, phase)] = data.phases[..] {
                    let g = crate::gates::cu_matrix_4x4(&[[ONE, ZERO], [ZERO, phase]]);
                    self.apply_two_qubit_gate(&g, targets[0], target)?;
                } else {
                    let control = self.site_for_logical(targets[0]);
                    let phases: Vec<(usize, Complex64)> = data
                        .phases
                        .iter()
                        .map(|&(qubit, phase)| (self.site_for_logical(qubit), phase))
                        .collect();
                    self.apply_batch_phase_bubble(control, &phases)?;
                }
            }
            Gate::BatchRzz(data) => {
                for &(q0, q1, theta) in &data.edges {
                    let g = Gate::Rzz(theta).matrix_4x4();
                    self.apply_two_qubit_gate(&g, q0, q1)?;
                }
            }
            Gate::DiagonalBatch(data) => {
                for entry in &data.entries {
                    if let Some((q, mat)) = entry.as_1q_matrix() {
                        self.apply_single_qubit_gate(self.site_for_logical(q), &mat);
                    } else if let Some((q0, q1, mat)) = entry.as_2q_matrix() {
                        self.apply_two_qubit_gate(&mat, q0, q1)?;
                    }
                }
            }
            Gate::MultiFused(data) => {
                for &(target, ref mat) in &data.gates {
                    self.apply_single_qubit_gate(self.site_for_logical(target), mat);
                }
            }
            Gate::Fused2q(mat) => {
                self.apply_two_qubit_gate(mat, targets[0], targets[1])?;
            }
            Gate::Multi2q(data) => {
                for &(q0, q1, ref mat) in &data.gates {
                    self.apply_two_qubit_gate(mat, q0, q1)?;
                }
            }
            single_qubit => {
                let mat = single_qubit.matrix_2x2();
                self.apply_single_qubit_gate(self.site_for_logical(targets[0]), &mat);
            }
        }
        Ok(())
    }
}

/// Two-qubit gate entries `instruction` contributes to a run of
/// [`MpsBackend::apply_by_center`]: one for a gate the adjacent kernel takes
/// whole, one per entry of a `Multi2q` or `BatchRzz`, none for anything else.
fn pair_entries(instruction: &Instruction) -> usize {
    match instruction {
        Instruction::Gate { gate, .. } => match gate {
            Gate::Rzz(_) | Gate::Cx | Gate::Cz | Gate::Swap | Gate::Cu(_) | Gate::Fused2q(_) => 1,
            Gate::Multi2q(data) => data.gates.len(),
            Gate::BatchRzz(data) => data.edges.len(),
            _ => 0,
        },
        _ => 0,
    }
}

/// Logical qubits of entry `k` of `instruction`, one of the
/// [`pair_entries`] it owes.
fn pair_entry(instruction: &Instruction, k: usize) -> (usize, usize) {
    match instruction {
        Instruction::Gate {
            gate: Gate::Multi2q(data),
            ..
        } => {
            let (q0, q1, _) = &data.gates[k];
            (*q0, *q1)
        }
        Instruction::Gate {
            gate: Gate::BatchRzz(data),
            ..
        } => {
            let (q0, q1, _) = data.edges[k];
            (q0, q1)
        }
        Instruction::Gate { targets, .. } => (targets[0], targets[1]),
        _ => unreachable!("an instruction without pair entries"),
    }
}

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

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

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

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

    /// A chain whose bonds are all 1 certifies itself: it is a product state,
    /// and with nothing discarded on the way to it there was no approximation
    /// to make. Any wider chain reports `Approximate` whether or not this run
    /// truncated, since the bond cap can discard weight. The claim comes from
    /// the chain rather than from the cut history, which records what a cut
    /// wrote and so passes over a bond a center move widened.
    fn exactness(&self) -> crate::sim::Exactness {
        if self.current_max_bond_dim() == 1 && self.truncation_discarded <= 0.0 {
            return crate::sim::Exactness::Exact;
        }
        crate::sim::Exactness::Approximate {
            fidelity_lower_bound: Some((1.0 - self.truncation_discarded).max(0.0)),
        }
    }

    fn bond_report(&self) -> Option<crate::sim::BondReport> {
        Some(crate::sim::BondReport {
            peak: self.bond_high_water,
            cap: self.max_bond_dim,
        })
    }

    fn init(&mut self, num_qubits: usize, num_classical_bits: usize) -> Result<()> {
        self.num_qubits = num_qubits;
        self.truncation_discarded = 0.0;
        self.center = None;
        self.bond_high_water = 1;
        crate::backend::init_classical_bits(&mut self.classical_bits, num_classical_bits);
        self.sites = (0..num_qubits)
            .map(|_| SiteTensor::new_zero_state())
            .collect();
        self.logical_to_site = (0..num_qubits).collect();
        self.site_to_logical = (0..num_qubits).collect();
        Ok(())
    }

    fn apply(&mut self, instruction: &Instruction) -> Result<()> {
        match instruction {
            Instruction::Gate { gate, targets } => {
                let list = matches!(gate, Gate::Multi2q(_) | Gate::BatchRzz(_));
                if list && pair_entries(instruction) > 0 {
                    self.apply_by_center(std::slice::from_ref(instruction))?;
                } else {
                    self.dispatch_gate(gate, targets)?;
                }
            }
            Instruction::Measure {
                qubit,
                classical_bit,
            } => {
                self.apply_measure(self.site_for_logical(*qubit), *classical_bit);
            }
            Instruction::Reset { qubit } => {
                self.apply_reset(self.site_for_logical(*qubit));
            }
            Instruction::Barrier { .. } => {}
            Instruction::Save { label, .. } => {
                return Err(crate::backend::save_not_applied("Mps", 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 apply_instructions(&mut self, instructions: &[Instruction]) -> Result<()> {
        self.apply_by_center(instructions)
    }

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

    fn apply_1q_matrix(&mut self, qubit: usize, matrix: &[[Complex64; 2]; 2]) -> Result<()> {
        let site = self.site_for_logical(qubit);
        self.invalidate_center_unless(site);
        self.apply_single_qubit_gate(site, matrix);
        Ok(())
    }

    fn reduced_density_matrix_1q(&self, qubit: usize) -> Result<[[Complex64; 2]; 2]> {
        Ok(self.reduced_density_site(self.site_for_logical(qubit)))
    }

    /// Contracted out of the chain rather than traced out of a vector: with the
    /// center inside the span, the sites on either side of it collapse to
    /// identities and one sweep carries a transfer object that grows fourfold
    /// per kept site. The sweep indexes by chain position, so a SWAP-routed
    /// layout is permuted back into the requested qubit order at the end.
    fn reduced_density_matrix(&mut self, subsystem: &[usize]) -> Result<Vec<Complex64>> {
        schmidt::validate_qubit_set(subsystem, self.num_qubits)?;
        let k = subsystem.len();
        let dim = crate::backend::reduced_density::reduced_density_side(self.name(), k)?;

        let mut order: Vec<usize> = (0..k).collect();
        order.sort_by_key(|&i| self.site_for_logical(subsystem[i]));
        let sites: Vec<usize> = order
            .iter()
            .map(|&i| self.site_for_logical(subsystem[i]))
            .collect();

        let (first, last) = (sites[0], sites[k - 1]);
        let chi = self.sites[first..=last]
            .iter()
            .fold(1usize, |chi, t| chi.max(t.bond_left).max(t.bond_right));
        let workspace = (7u128 * (chi as u128) * (chi as u128)) << (2 * k);
        if workspace > self.workspace_cap {
            return Err(crate::backend::workspace_allocation_error(
                self.name(),
                "a reduced density matrix over the span these qubits occupy",
                workspace,
            ));
        }

        let end = self.sites.len() - 1;
        match self.center {
            Some(at) => self.move_center(at.clamp(first, last)),
            None if first <= end - last => self.establish_center(first),
            None => self.establish_center(last),
        }
        let contracted = self.reduced_density_sites(&sites);

        // Bit `k - 1 - j` of a contracted index is the site at position `j`,
        // which hosts `subsystem[order[j]]` and so becomes bit `order[j]`.
        let permuted: Vec<usize> = (0..dim)
            .map(|t| {
                (0..k)
                    .filter(|&j| (t >> (k - 1 - j)) & 1 == 1)
                    .fold(0, |index, j| index | (1 << order[j]))
            })
            .collect();
        let mut rho = vec![ZERO; dim * dim];
        for (t, row) in contracted.chunks_exact(dim).enumerate() {
            for (tp, &amp) in row.iter().enumerate() {
                rho[permuted[t] * dim + permuted[tp]] = amp;
            }
        }
        crate::backend::reduced_density::normalize_trace(&mut rho, dim);
        Ok(rho)
    }

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

    /// The one-SVD route when the subsystem is a run of sites at either end
    /// of the chain; otherwise the eigenvalues of the reduced density matrix
    /// over the smaller side of the cut, under the dense cap. The chain order
    /// is never re-sorted for it.
    fn schmidt_values(&mut self, subsystem: &[usize]) -> Result<Vec<f64>> {
        schmidt::validate_subsystem(subsystem, self.num_qubits)?;
        let n = self.num_qubits;
        let mut sites: Vec<usize> = subsystem
            .iter()
            .map(|&q| self.site_for_logical(q))
            .collect();
        sites.sort_unstable();
        let k = sites.len();
        if sites[k - 1] - sites[0] + 1 == k {
            if sites[0] == 0 {
                return Ok(self.schmidt_values_at_bond(k - 1));
            }
            if sites[k - 1] == n - 1 {
                return Ok(self.schmidt_values_at_bond(sites[0] - 1));
            }
        }
        if 2 * k > n {
            let mut inside = vec![false; n];
            for &site in &sites {
                inside[site] = true;
            }
            sites = (0..n).filter(|&site| !inside[site]).collect();
        }
        self.schmidt_values_by_reduced_density(&sites)
    }

    /// [`MpsBackend::inner_product`] when both chains carry the same site
    /// layout, `O(n chi^3)` and no width cap. The contraction pairs sites by
    /// chain position, so a chain whose virtual swaps left a different order
    /// takes the dense route instead.
    fn overlap_sq(&self, other: &dyn Backend) -> Result<f64> {
        if let Some(chain) = other
            .as_any()
            .and_then(|any| any.downcast_ref::<MpsBackend>())
            && chain.num_qubits == self.num_qubits
            && chain.logical_to_site == self.logical_to_site
        {
            let inner = self.inner_product(chain)?;
            let left = self.inner_product(self)?.re;
            let right = chain.inner_product(chain)?.re;
            return Ok(overlap::normalized(inner.norm_sqr(), left, right));
        }
        overlap::export_overlap_sq(
            self.name(),
            self.num_qubits(),
            || self.export_statevector(),
            other,
        )
    }

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

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

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

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

    /// Sequential conditional sampling over the chain: `O(n·χ²)` per shot after
    /// one `O(n·χ³)` sweep for the right environments.
    ///
    /// Each shot draws from its own ChaCha8 substream keyed on `(seed, shot)`
    /// (streams 2 and up: 0 is the backend's own draw stream, 1 the
    /// trajectory noise stream), so the sampled words are identical at any
    /// thread count and the shot loop parallelizes at or above
    /// `MIN_SHOTS_FOR_PAR`.
    ///
    /// Sites are visited left to right and each bit is recorded against the
    /// logical qubit currently hosted there, so a layout permuted by SWAP
    /// routing needs no canonicalization pass.
    fn sample_basis_states(&mut self, num_shots: usize, seed: u64) -> Result<BasisSamples> {
        let n = self.num_qubits;
        let mut samples = BasisSamples::new(num_shots, n);
        if n == 0 {
            return Ok(samples);
        }

        let right = self.right_environments();
        let max_bond = self
            .sites
            .iter()
            .map(|site| site.bond_left.max(site.bond_right))
            .max()
            .unwrap_or(1);

        let words_per_shot = samples.words_per_shot();
        let sample_shot =
            |shot: usize, shot_words: &mut [u64], left: &mut [Complex64], w: &mut [Complex64]| {
                let mut rng = ChaCha8Rng::seed_from_u64(seed);
                rng.set_stream(shot as u64 + 2);
                left[0] = ONE;
                left[1..].fill(ZERO);
                for (site, right_env) in right.iter().enumerate() {
                    let br = self.sites[site].bond_right;
                    let prob = self.site_conditional_weights(site, left, right_env, w);

                    let total = prob[0] + prob[1];
                    let prob_one = if total > 0.0 {
                        (prob[1] / total).clamp(0.0, 1.0)
                    } else {
                        0.0
                    };
                    let outcome = usize::from(rng.random::<f64>() < prob_one);
                    if outcome == 1 {
                        let qubit = self.logical_for_site(site);
                        shot_words[qubit / 64] |= 1u64 << (qubit % 64);
                    }

                    let scale = 1.0 / prob[outcome].max(NORM_CLAMP_MIN).sqrt();
                    left[..br].copy_from_slice(&w[outcome * br..(outcome + 1) * br]);
                    for value in &mut left[..br] {
                        *value *= scale;
                    }
                }
            };

        #[cfg(feature = "parallel")]
        if num_shots >= MIN_SHOTS_FOR_PAR {
            samples
                .words_mut()
                .par_chunks_mut(words_per_shot)
                .enumerate()
                .for_each_init(
                    || (vec![ZERO; max_bond], vec![ZERO; 2 * max_bond]),
                    |(left, w), (shot, shot_words)| sample_shot(shot, shot_words, left, w),
                );
            return Ok(samples);
        }

        let mut left = vec![ZERO; max_bond];
        let mut w = vec![ZERO; 2 * max_bond];
        for shot in 0..num_shots {
            let start = shot * words_per_shot;
            sample_shot(
                shot,
                &mut samples.words_mut()[start..start + words_per_shot],
                &mut left,
                &mut w,
            );
        }
        Ok(samples)
    }

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

    /// Identity environments are swept once and shared, so an observable
    /// contracts only the sites between its first and last non-identity
    /// factor: a single-qubit observable costs one site instead of a full
    /// chain walk. Divided by `⟨ψ|ψ⟩` because a truncated chain is not
    /// exactly normalized.
    fn pauli_expectations(&self, observables: &[Vec<PauliTerm>]) -> Result<Vec<f64>> {
        if self.num_qubits == 0 {
            return Ok(vec![1.0; observables.len()]);
        }
        if observables.len() < 2 {
            // One observable costs the same two sweeps either way, and keeping
            // the plain contraction here leaves a route that shares no code
            // with the windowed one for a test to compare against.
            let norm = self.pauli_expectation(&[])?.re;
            return observables
                .iter()
                .map(|observable| {
                    let value = self.pauli_expectation(&mps_pauli_factors(observable))?;
                    Ok(if norm == 0.0 { 0.0 } else { value.re / norm })
                })
                .collect();
        }

        let prefixes = self.pauli_env_prefixes();
        let rights = self.right_environments();
        let norm = prefixes[self.num_qubits][0].re;

        observables
            .iter()
            .map(|observable| {
                let site_axis = self.pauli_site_axes(&mps_pauli_factors(observable))?;
                if norm == 0.0 {
                    return Ok(0.0);
                }
                let Some(first) = site_axis.iter().position(Option::is_some) else {
                    return Ok(1.0);
                };
                let last = site_axis.iter().rposition(Option::is_some).unwrap();

                let mut env = self.contract_pauli_env(
                    first,
                    &prefixes[first],
                    self.sites[first].bond_left,
                    site_axis[first],
                );
                let mut cols = self.sites[first].bond_right;
                for (site, axis) in site_axis.iter().enumerate().take(last + 1).skip(first + 1) {
                    env = self.contract_pauli_env(site, &env, cols, *axis);
                    cols = self.sites[site].bond_right;
                }

                // The identity cap right of `last` carries the ket bond first,
                // so closing the two chains is a trace of the two matrices.
                let right = &rights[last];
                let mut value = ZERO;
                for bra in 0..cols {
                    for ket in 0..cols {
                        value += env[bra * cols + ket] * right[ket * cols + bra];
                    }
                }
                Ok(value.re / norm)
            })
            .collect()
    }

    fn export_statevector(&self) -> Result<Vec<Complex64>> {
        let dim = dense_statevector_len(self.name(), "statevector export", self.num_qubits)?;
        let mut state = Vec::new();
        reserve_dense_output(&mut state, dim, self.name(), "statevector export")?;
        state.resize(dim, ZERO);
        self.expand_dense(&mut state, |amp| amp);
        Ok(state)
    }
}

#[cfg(test)]
#[path = "mps_tests.rs"]
mod tests;