prism-q 0.35.0

//! Gate application kernels for the statevector backend.
//!
//! Each kernel has a sequential implementation and, when the `parallel` feature
//! is enabled, a Rayon-parallelized variant dispatched above
//! `PARALLEL_THRESHOLD_QUBITS`.

use num_complex::Complex64;
use rand::RngExt;
use smallvec::SmallVec;
use std::sync::atomic::{AtomicU64, Ordering};
use std::sync::{Arc, OnceLock, RwLock};

use super::StatevectorBackend;
use super::insert_zero_bit;
#[cfg(feature = "parallel")]
use super::{MIN_PAR_ELEMS, PARALLEL_THRESHOLD_QUBITS, SendPtr};
use crate::backend::simd;
use crate::backend::{MCU_QUBIT_BUF, is_phase_one, measurement_inv_norm, sorted_mcu_qubits};
use crate::circuit::{QftTextbookStep, qft_textbook_steps};
use crate::gates::{
    BatchPhaseData, BatchRzzData, DiagEntry, DiagonalBatchData, Gate, MAX_UNITARY_QUBITS,
    MULTI_2Q_HIGH_BUDGET, diag_entries_phase, multi_2q_high_budget_for, multi_2q_low_bits,
    multi_2q_tile_bits_for, pauli_rot_masks,
};
use crate::sim::unified_pauli::PauliAxis;
#[cfg(feature = "parallel")]
use rayon::prelude::*;

use crate::backend::cache;
#[cfg(feature = "parallel")]
use crate::backend::{MIN_PAR_ITERS, chunk_min_len};

/// Largest qubit target whose full period (2^(t+1) elements) fits within `tile_size`.
#[inline(always)]
const fn max_target_for_tile(tile_size: usize) -> usize {
    tile_size.ilog2() as usize - 1
}

/// Amplitudes per L2 tile: the tile budget, 2^14 on the x86 cores measured.
#[inline(always)]
fn multi_gate_l2_tile() -> usize {
    cache::tile_budget_bytes() / size_of::<Complex64>()
}

/// Highest target the L2 tier takes in a sequential pass.
#[inline(always)]
pub(crate) fn multi_gate_max_l2_target() -> usize {
    max_target_for_tile(multi_gate_l2_tile())
}

/// Amplitudes per L2 tile for a parallel pass over `state_len` amplitudes: the budget,
/// cut to leave the pool its tile count (see [`cache::parallel_tile_len`]).
#[cfg(feature = "parallel")]
#[inline(always)]
fn multi_gate_l2_tile_par(state_len: usize) -> usize {
    cache::parallel_tile_len(state_len, multi_gate_l2_tile(), MIN_PAR_ELEMS)
}

/// The in-place tile a `Multi2q` tiered pass over a `num_qubits` state runs with: the
/// subcube tile, so a batch that sits below it lands in one tier and keeps its order.
#[inline(always)]
fn multi_2q_small_tile(num_qubits: usize) -> usize {
    1 << multi_2q_tile_bits_for(num_qubits)
}

const MULTI_GATE_L3_TILE: usize = 131_072;
const MULTI_GATE_MAX_L3_TARGET: usize = max_target_for_tile(MULTI_GATE_L3_TILE);

/// Length of the longest prefix of `gates` that
/// [`StatevectorBackend::apply_multi_2q`] applies in list order on a
/// `num_qubits` state: every gate sits in the in-place tile, so the tiered pass
/// runs one tier, or the prefix fits a [`subcube_plan`]. Past both, the tiered
/// pass reorders gates across tiers, so callers that need list order beyond
/// one batch (the density-matrix registers) cut the list into runs of this
/// length. At least 1 for a non-empty list, since one gate is always in order.
pub(crate) fn multi_2q_ordered_prefix(
    gates: &[(usize, usize, [[Complex64; 4]; 4])],
    num_qubits: usize,
) -> usize {
    let tile_bits = multi_2q_tile_bits_for(num_qubits);
    let budget = multi_2q_high_budget_for(num_qubits);
    let low_bits = multi_2q_low_bits();
    let mut high: SmallVec<[usize; MULTI_2Q_HIGH_BUDGET]> = SmallVec::new();
    let mut in_tile = true;
    let mut over_budget = false;
    for (i, &(q0, q1, _)) in gates.iter().enumerate() {
        in_tile &= q0.max(q1) < tile_bits;
        for q in [q0, q1] {
            if q >= low_bits && !high.contains(&q) {
                if high.len() == budget {
                    over_budget = true;
                } else {
                    high.push(q);
                }
            }
        }
        if over_budget && !in_tile {
            return i.max(1);
        }
    }
    gates.len()
}

/// Tile geometry for a `Multi2q` batch that reaches past the lowest tile bits:
/// a tile is `2^high.len()` contiguous runs of `2^low` amplitudes, one run per
/// setting of the `high` qubits, and `low + high.len() == tile_bits`.
/// Gathered into one buffer, run `c` occupies `[c << low, (c + 1) << low)`, so
/// a low qubit keeps its bit and the `j`th high qubit becomes bit `low + j`.
struct SubcubePlan {
    tile_bits: usize,
    low: usize,
    high: SmallVec<[usize; MULTI_2Q_HIGH_BUDGET]>,
    /// Qubit index to bit position inside the gathered tile.
    map: Vec<usize>,
}

impl SubcubePlan {
    /// Index bits outside the tile, ascending: the ones an outer index walks.
    fn rest(&self, num_qubits: usize) -> Vec<usize> {
        (self.low..num_qubits)
            .filter(|q| !self.high.contains(q))
            .collect()
    }
}

/// The subcube tile for `gates`, or `None` when the batch is served better by
/// the in-place tiles (every qubit below the tile bits, or the state is one
/// tile) or cannot be tiled at all (more than [`multi_2q_high_budget`]
/// distinct qubits at or above [`multi_2q_low_bits`]).
fn subcube_plan(
    gates: &[(usize, usize, [[Complex64; 4]; 4])],
    num_qubits: usize,
) -> Option<SubcubePlan> {
    let tile_bits = multi_2q_tile_bits_for(num_qubits);
    if num_qubits <= tile_bits {
        return None;
    }
    let max_q = gates.iter().map(|&(q0, q1, _)| q0.max(q1)).max()?;
    if max_q < tile_bits {
        return None;
    }
    let budget = multi_2q_high_budget_for(num_qubits);
    let low_bits = multi_2q_low_bits();
    let mut high: SmallVec<[usize; MULTI_2Q_HIGH_BUDGET]> = SmallVec::new();
    for &(q0, q1, _) in gates {
        for q in [q0, q1] {
            if q >= low_bits && !high.contains(&q) {
                if high.len() == budget {
                    return None;
                }
                high.push(q);
            }
        }
    }
    // A high qubit that lands inside the contiguous run rides there instead,
    // which lengthens the run and can free another; iterate to a fixed point.
    let mut low = tile_bits - high.len();
    loop {
        let before = high.len();
        high.retain(|&mut q| q >= low);
        if high.len() == before {
            break;
        }
        low = tile_bits - high.len();
    }
    high.sort_unstable();
    let mut map = vec![usize::MAX; num_qubits];
    for (q, slot) in map.iter_mut().enumerate().take(low) {
        *slot = q;
    }
    for (j, &h) in high.iter().enumerate() {
        map[h] = low + j;
    }
    Some(SubcubePlan {
        tile_bits,
        low,
        high,
        map,
    })
}

/// [`subcube_plan`] over `qubits`, the qubits a batch needs inside the tile (the X and
/// Y qubits of a `MultiPauliRot` batch, the targets of a `MultiFused` one), with the
/// same `None` cases.
fn qubit_subcube_plan(
    qubits: impl Iterator<Item = usize> + Clone,
    num_qubits: usize,
) -> Option<SubcubePlan> {
    let tile_bits = multi_2q_tile_bits_for(num_qubits);
    if num_qubits <= tile_bits {
        return None;
    }
    let max_q = qubits.clone().max()?;
    if max_q < tile_bits {
        return None;
    }
    let budget = multi_2q_high_budget_for(num_qubits);
    let low_bits = multi_2q_low_bits();
    let mut high: SmallVec<[usize; MULTI_2Q_HIGH_BUDGET]> = SmallVec::new();
    for q in qubits {
        if q >= low_bits && !high.contains(&q) {
            if high.len() == budget {
                return None;
            }
            high.push(q);
        }
    }
    // A high qubit that lands inside the contiguous run rides there instead,
    // which lengthens the run and can free another; iterate to a fixed point.
    let mut low = tile_bits - high.len();
    loop {
        let before = high.len();
        high.retain(|&mut q| q >= low);
        if high.len() == before {
            break;
        }
        low = tile_bits - high.len();
    }
    high.sort_unstable();
    let mut map = vec![usize::MAX; num_qubits];
    for (q, slot) in map.iter_mut().enumerate().take(low) {
        *slot = q;
    }
    for (j, &h) in high.iter().enumerate() {
        map[h] = low + j;
    }
    Some(SubcubePlan {
        tile_bits,
        low,
        high,
        map,
    })
}

fn prepare_2q(
    gates: &[(usize, usize, [[Complex64; 4]; 4])],
) -> Vec<(usize, usize, simd::PreparedGate2q)> {
    gates
        .iter()
        .map(|&(q0, q1, ref mat)| (q0, q1, simd::PreparedGate2q::new(mat)))
        .collect()
}

/// Prepare the gates whose high qubit fits the L2 tile (up to `max_l2_target`) or the
/// L3 tile, by tier. Gates above the L3 tile are left out.
#[inline(always)]
fn split_multi_2q_tiers(
    gates: &[(usize, usize, [[Complex64; 4]; 4])],
    max_l2_target: usize,
    small_gates: &mut SmallVec<[(usize, usize, simd::PreparedGate2q); 2]>,
    medium_gates: &mut SmallVec<[(usize, usize, simd::PreparedGate2q); 2]>,
) {
    for &(q0, q1, ref mat) in gates {
        let max_q = q0.max(q1);
        if max_q <= max_l2_target {
            small_gates.push((q0, q1, simd::PreparedGate2q::new(mat)));
        } else if max_q <= MULTI_GATE_MAX_L3_TARGET {
            medium_gates.push((q0, q1, simd::PreparedGate2q::new(mat)));
        }
    }
}

thread_local! {
    /// The gathered tile, one per thread and kept between batches: a fresh
    /// allocation per subcube paid more in page faults than the sweeps it
    /// saved.
    static SUBCUBE_TILE: std::cell::RefCell<Vec<Complex64>> =
        const { std::cell::RefCell::new(Vec::new()) };
}

fn with_subcube_tile<R>(tile_bits: usize, f: impl FnOnce(&mut [Complex64]) -> R) -> R {
    SUBCUBE_TILE.with(|tile| {
        let mut tile = tile.borrow_mut();
        if tile.len() != 1 << tile_bits {
            tile.resize(1 << tile_bits, Complex64::new(0.0, 0.0));
        }
        f(&mut tile)
    })
}

/// Gather the subcube `outer` selects into this thread's tile, apply `gates`
/// there in order through the contiguous tiled kernel, and scatter it back.
/// Applying in place across the runs instead was measured slower at every
/// width: the strided group walk gives up the paired AVX2 kernel.
///
/// # Safety
/// `state` must point at `2^num_qubits` amplitudes laid out for `plan` and
/// `rest`, and no other access to this subcube may overlap the call.
unsafe fn apply_subcube(
    state: *mut Complex64,
    plan: &SubcubePlan,
    rest: &[usize],
    outer: usize,
    gates: &[(usize, usize, simd::PreparedGate2q)],
) {
    let run = 1usize << plan.low;
    let mut base = 0usize;
    for (j, &p) in rest.iter().enumerate() {
        base |= ((outer >> j) & 1) << p;
    }
    let run_offset = |c: usize| {
        let mut off = base;
        for (j, &h) in plan.high.iter().enumerate() {
            off |= ((c >> j) & 1) << h;
        }
        off
    };
    let runs = 1usize << plan.high.len();
    with_subcube_tile(plan.tile_bits, |tile| {
        for c in 0..runs {
            // SAFETY: same contract as the enclosing unsafe fn; the run is
            // inside the subcube and the tile holds `runs << low` elements.
            unsafe {
                std::ptr::copy_nonoverlapping(
                    state.add(run_offset(c)),
                    tile.as_mut_ptr().add(c << plan.low),
                    run,
                );
            }
        }
        // A zero subcube stays zero under any unitary, so it needs no gates and no
        // write back. Early in a circuit most subcubes are zero; a dense one exits
        // this scan at its first amplitude.
        if tile[..runs << plan.low]
            .iter()
            .all(|a| a.re == 0.0 && a.im == 0.0)
        {
            return;
        }
        for &(q0, q1, ref prepared) in gates {
            prepared.apply_tiled(tile, plan.tile_bits, plan.map[q0], plan.map[q1]);
        }
        for c in 0..runs {
            // SAFETY: same contract as the enclosing unsafe fn.
            unsafe {
                std::ptr::copy_nonoverlapping(
                    tile.as_ptr().add(c << plan.low),
                    state.add(run_offset(c)),
                    run,
                );
            }
        }
    });
}

/// [`apply_subcube`] for a `MultiPauliRot` batch: gather the subcube, apply
/// `rotations` there in order, and scatter it back. Kept apart from the `Multi2q`
/// gather: one generic body for both read 5% slower on `statevector/qv/16`.
///
/// # Safety
/// `state` must point at `2^num_qubits` amplitudes laid out for `plan` and
/// `rest`, and no other access to this subcube may overlap the call.
unsafe fn apply_pauli_rot_subcube(
    state: *mut Complex64,
    plan: &SubcubePlan,
    rest: &[usize],
    outer: usize,
    rotations: &[TiledPauliRot],
) {
    let run = 1usize << plan.low;
    let mut base = 0usize;
    for (j, &p) in rest.iter().enumerate() {
        base |= ((outer >> j) & 1) << p;
    }
    let run_offset = |c: usize| {
        let mut off = base;
        for (j, &h) in plan.high.iter().enumerate() {
            off |= ((c >> j) & 1) << h;
        }
        off
    };
    let runs = 1usize << plan.high.len();
    with_subcube_tile(plan.tile_bits, |tile| {
        for c in 0..runs {
            // SAFETY: same contract as the enclosing unsafe fn; the run is
            // inside the subcube and the tile holds `runs << low` elements.
            unsafe {
                std::ptr::copy_nonoverlapping(
                    state.add(run_offset(c)),
                    tile.as_mut_ptr().add(c << plan.low),
                    run,
                );
            }
        }
        // A zero subcube stays zero under any unitary, so it needs no gates and no
        // write back. Early in a circuit most subcubes are zero; a dense one exits
        // this scan at its first amplitude.
        if tile[..runs << plan.low]
            .iter()
            .all(|a| a.re == 0.0 && a.im == 0.0)
        {
            return;
        }
        for rotation in rotations {
            rotation.apply(tile, base);
        }
        for c in 0..runs {
            // SAFETY: same contract as the enclosing unsafe fn.
            unsafe {
                std::ptr::copy_nonoverlapping(
                    tile.as_ptr().add(c << plan.low),
                    state.add(run_offset(c)),
                    run,
                );
            }
        }
    });
}

/// [`apply_subcube`] for a `MultiFused` batch, each gate keyed by its bit in the
/// gathered tile. Kept apart from the `Multi2q` gather for the reason
/// [`apply_pauli_rot_subcube`] gives.
///
/// # Safety
/// `state` must point at `2^num_qubits` amplitudes laid out for `plan` and
/// `rest`, and no other access to this subcube may overlap the call.
unsafe fn apply_multi_1q_subcube(
    state: *mut Complex64,
    plan: &SubcubePlan,
    rest: &[usize],
    outer: usize,
    gates: &[(usize, simd::PreparedGate1q)],
) {
    let run = 1usize << plan.low;
    let mut base = 0usize;
    for (j, &p) in rest.iter().enumerate() {
        base |= ((outer >> j) & 1) << p;
    }
    let run_offset = |c: usize| {
        let mut off = base;
        for (j, &h) in plan.high.iter().enumerate() {
            off |= ((c >> j) & 1) << h;
        }
        off
    };
    let runs = 1usize << plan.high.len();
    with_subcube_tile(plan.tile_bits, |tile| {
        for c in 0..runs {
            // SAFETY: same contract as the enclosing unsafe fn; the run is
            // inside the subcube and the tile holds `runs << low` elements.
            unsafe {
                std::ptr::copy_nonoverlapping(
                    state.add(run_offset(c)),
                    tile.as_mut_ptr().add(c << plan.low),
                    run,
                );
            }
        }
        // A zero subcube stays zero under any unitary.
        if tile[..runs << plan.low]
            .iter()
            .all(|a| a.re == 0.0 && a.im == 0.0)
        {
            return;
        }
        for &(bit, ref prepared) in gates {
            prepared.apply_tiled(tile, bit);
        }
        for c in 0..runs {
            // SAFETY: same contract as the enclosing unsafe fn.
            unsafe {
                std::ptr::copy_nonoverlapping(
                    tile.as_ptr().add(c << plan.low),
                    state.add(run_offset(c)),
                    run,
                );
            }
        }
    });
}

/// Run `gates`, keyed by tile bit, over every subcube of `state` that `plan` selects:
/// one pass over the state.
fn apply_multi_1q_subcubes(
    state: &mut [Complex64],
    plan: &SubcubePlan,
    gates: &[(usize, simd::PreparedGate1q)],
) {
    let num_qubits = state.len().trailing_zeros() as usize;
    let rest = plan.rest(num_qubits);
    let subcubes = 1usize << (num_qubits - plan.tile_bits);
    #[cfg(feature = "parallel")]
    if num_qubits >= PARALLEL_THRESHOLD_QUBITS {
        let ptr = SendPtr(state.as_mut_ptr());
        (0..subcubes).into_par_iter().for_each(|outer| {
            // SAFETY: `ptr` holds 2^num_qubits amplitudes and distinct `outer`
            // values address disjoint subcubes, so no two tasks touch the same
            // amplitude.
            unsafe { apply_multi_1q_subcube(ptr.as_complex_ptr(), plan, &rest, outer, gates) };
        });
        return;
    }
    let ptr = state.as_mut_ptr();
    for outer in 0..subcubes {
        // SAFETY: `ptr` holds 2^num_qubits amplitudes and each `outer` addresses
        // its own subcube, applied one at a time here.
        unsafe { apply_multi_1q_subcube(ptr, plan, &rest, outer, gates) };
    }
}

/// Targets folded into one shared high-target traversal. `lanes * tile_len` is
/// fixed at the L2 budget, so each lane's run shrinks as `2^-k` while the
/// stream count grows; past three the tiles are too short to prefetch.
#[cfg(feature = "parallel")]
const MAX_SHARED_1Q_TARGETS: usize = 3;

#[inline(always)]
fn adjacent_2q_indices(offset: usize, stride: usize, q0_is_lo: bool) -> [usize; 4] {
    if q0_is_lo {
        [
            offset,
            offset + (stride << 1),
            offset + stride,
            offset + (stride * 3),
        ]
    } else {
        [
            offset,
            offset + stride,
            offset + (stride << 1),
            offset + (stride * 3),
        ]
    }
}

#[inline(always)]
fn swap_slices_kernel(a: &mut [Complex64], b: &mut [Complex64]) {
    debug_assert_eq!(a.len(), b.len());
    if a.len() < 4 {
        for (x, y) in a.iter_mut().zip(b.iter_mut()) {
            std::mem::swap(x, y);
        }
    } else {
        simd::swap_slices(a, b);
    }
}

#[inline(always)]
fn negate_slice_kernel(slice: &mut [Complex64]) {
    if slice.len() < 4 {
        for amp in slice {
            *amp = -*amp;
        }
    } else {
        simd::negate_slice(slice);
    }
}

/// Fold the `|1>` branch of a reset onto the `|0>` branch over one aligned
/// `(lo, hi)` pair, then clear `hi`. Shared by the sequential reset and the
/// block and sub-tile arms of [`StatevectorBackend::fold_reset_par`].
#[inline(always)]
fn fold_reset_pair(lo: &mut [Complex64], hi: &mut [Complex64], outcome: bool) {
    if outcome {
        lo.copy_from_slice(hi);
    }
    simd::zero_slice(hi);
}

/// Diagonal `(d0, d1)` on `target` when the state holds fewer than four
/// `2^(target+1)` blocks, where tiling by the block leaves at most three Rayon
/// tasks. Scales each half of every block in parallel sub-tiles instead.
#[cfg(feature = "parallel")]
#[inline(always)]
fn apply_diagonal_high_target_par(
    state: &mut [Complex64],
    target: usize,
    d0: Complex64,
    d1: Complex64,
    skip_lo: bool,
) {
    let half = 1usize << target;
    for block in state.chunks_mut(half << 1) {
        let (lo, hi) = block.split_at_mut(half);
        lo.par_chunks_mut(MIN_PAR_ELEMS)
            .zip(hi.par_chunks_mut(MIN_PAR_ELEMS))
            .for_each(|(lo_t, hi_t)| {
                if !skip_lo {
                    simd::scale_complex_slice(lo_t, d0);
                }
                simd::scale_complex_slice(hi_t, d1);
            });
    }
}

pub(crate) const BATCH_PHASE_GROUP_SIZE: usize = 10;
pub(crate) const BATCH_PHASE_TABLE_SIZE: usize = 1024;
pub(crate) const MAX_BATCH_PHASE_GROUPS: usize = 4;

pub(crate) const BATCH_RZZ_GROUP_SIZE: usize = 8;
pub(crate) const BATCH_RZZ_TABLE_SIZE: usize = 256;
pub(crate) const MAX_BATCH_RZZ_GROUPS: usize = 4;
#[cfg(target_arch = "x86_64")]
const BATCH_RZZ_BMI2_MAX_UNIQUE: usize = 10;
#[cfg(target_arch = "x86_64")]
const BATCH_RZZ_BMI2_TABLE_SIZE: usize = 1024;
#[cfg(target_arch = "x86_64")]
const BATCH_RZZ_BMI2_MAX_GROUPS: usize = 8;

const _: () = assert!(
    MAX_BATCH_PHASE_GROUPS * BATCH_PHASE_GROUP_SIZE == BatchPhaseData::MAX_PHASES
        && MAX_BATCH_RZZ_GROUPS * BATCH_RZZ_GROUP_SIZE == BatchRzzData::MAX_EDGES,
    "fusion splits these runs against the gate-level cap, which must match the table shape"
);

pub(crate) const DIAG_BATCH_MAX_QUBITS_PER_GROUP: usize = DiagonalBatchData::MAX_QUBITS_PER_GROUP;
pub(crate) const DIAG_BATCH_TABLE_SIZE: usize = 1 << DIAG_BATCH_MAX_QUBITS_PER_GROUP;
pub(crate) const MAX_DIAG_BATCH_GROUPS: usize = DiagonalBatchData::MAX_GROUPS;

type QftTwiddleTable = Arc<[Complex64]>;

struct CachedQftTwiddles {
    table: QftTwiddleTable,
    last_used: AtomicU64,
}

// Hits take the read lock only: the LRU stamp and the tick are atomics so a
// lookup never needs exclusive access. Misses build the table outside any lock
// and take the write lock to insert and evict. Relaxed suffices on the atomics
// because the write guard excludes every reader: the guard's release/acquire
// pair orders each stamp before the eviction scan that reads it.
struct QftTwiddleCache {
    entries: Vec<Option<CachedQftTwiddles>>,
    bytes: usize,
    limit_bytes: usize,
    tick: AtomicU64,
    /// Holds the conjugate tables the inverse transform reads.
    inverse: bool,
}

impl QftTwiddleCache {
    fn new(limit_bytes: usize) -> Self {
        Self {
            entries: Vec::new(),
            bytes: 0,
            limit_bytes,
            tick: AtomicU64::new(0),
            inverse: false,
        }
    }

    #[inline]
    fn next_tick(&self) -> u64 {
        self.tick.fetch_add(1, Ordering::Relaxed).wrapping_add(1)
    }
}

const QFT_TWIDDLE_CACHE_DEFAULT_LIMIT_MB: usize = 256;

// Soft LRU cap for QFT twiddles. One oversized table may occupy the cache
// alone, which keeps large repeated QFT runs from rebuilding twiddles every
// sample. `PRISM_QFT_TWIDDLE_CACHE_LIMIT_MB=0` disables caching.
#[inline]
fn qft_twiddle_cache_limit_bytes() -> usize {
    static CACHED: OnceLock<usize> = OnceLock::new();
    *CACHED.get_or_init(|| {
        crate::env_knobs::usize_knob(
            "PRISM_QFT_TWIDDLE_CACHE_LIMIT_MB",
            QFT_TWIDDLE_CACHE_DEFAULT_LIMIT_MB,
            0,
        )
        .saturating_mul(1024 * 1024)
    })
}

#[inline(always)]
fn qft_twiddle_table_bytes(table_len: usize) -> usize {
    table_len.saturating_mul(std::mem::size_of::<Complex64>())
}

/// DIF twiddles `exp(2 pi i k / 2^n) / sqrt(2)`, conjugated for `inverse`.
fn qft_twiddles_scaled(n: usize, inverse: bool) -> QftTwiddleTable {
    static FORWARD: OnceLock<RwLock<QftTwiddleCache>> = OnceLock::new();
    static INVERSE: OnceLock<RwLock<QftTwiddleCache>> = OnceLock::new();
    let cache = if inverse {
        INVERSE.get_or_init(|| {
            RwLock::new(QftTwiddleCache {
                inverse: true,
                ..QftTwiddleCache::new(qft_twiddle_cache_limit_bytes())
            })
        })
    } else {
        FORWARD.get_or_init(|| RwLock::new(QftTwiddleCache::new(qft_twiddle_cache_limit_bytes())))
    };
    cached_qft_twiddles(cache, n)
}

fn cached_qft_twiddles(cache: &RwLock<QftTwiddleCache>, n: usize) -> QftTwiddleTable {
    let (cache_limit, inverse) = {
        let guard = cache.read().unwrap();
        if let Some(Some(entry)) = guard.entries.get(n) {
            entry
                .last_used
                .fetch_max(guard.next_tick(), Ordering::Relaxed);
            return Arc::clone(&entry.table);
        }
        (guard.limit_bytes, guard.inverse)
    };

    let inv_sqrt2 = std::f64::consts::FRAC_1_SQRT_2;
    let qft_size = 1usize << n;
    let mut twiddles = Vec::with_capacity(qft_size / 2);
    for k in 0..qft_size / 2 {
        let angle = std::f64::consts::TAU * k as f64 / qft_size as f64;
        let w = Complex64::from_polar(inv_sqrt2, angle);
        twiddles.push(if inverse { w.conj() } else { w });
    }
    let twiddles: Arc<[Complex64]> = twiddles.into();
    let table_bytes = qft_twiddle_table_bytes(twiddles.len());
    if cache_limit == 0 {
        return twiddles;
    }

    let mut guard = cache.write().unwrap();
    let tick = guard.next_tick();
    if guard.entries.len() <= n {
        guard.entries.resize_with(n + 1, || None);
    }
    if let Some(existing) = &guard.entries[n] {
        existing.last_used.fetch_max(tick, Ordering::Relaxed);
        return Arc::clone(&existing.table);
    }

    while guard.bytes.saturating_add(table_bytes) > cache_limit {
        let Some((idx, _)) = guard
            .entries
            .iter()
            .enumerate()
            .filter_map(|(idx, entry)| {
                entry
                    .as_ref()
                    .map(|entry| (idx, entry.last_used.load(Ordering::Relaxed)))
            })
            .min_by_key(|&(_, last_used)| last_used)
        else {
            break;
        };
        if let Some(evicted) = guard.entries[idx].take() {
            guard.bytes = guard
                .bytes
                .saturating_sub(qft_twiddle_table_bytes(evicted.table.len()));
        }
    }

    guard.entries[n] = Some(CachedQftTwiddles {
        table: Arc::clone(&twiddles),
        last_used: AtomicU64::new(tick),
    });
    guard.bytes = guard.bytes.saturating_add(table_bytes);
    twiddles
}

/// Fewest low qubits [`StatevectorBackend::apply_qft_block`] runs together in one tile.
const QFT_MIN_TILE_BITS: usize = 13;

/// Most low qubits per tile: twice the tile budget, 512 KiB on the x86 cores measured.
#[inline(always)]
fn qft_max_tile_bits() -> usize {
    (2 * cache::tile_budget_bytes() / size_of::<Complex64>()).ilog2() as usize
}

/// Tile width for an `n`-qubit QFT over a `state_bits`-qubit state: as wide as leaves
/// the pool its tile count ([`cache::min_parallel_tiles`]), and above the minimum one
/// bit narrower where that leaves an even count of stages above the tile, which then
/// all run as radix-4 pairs instead of paying a radix-2 pass.
fn qft_tile_bits(n: usize, state_bits: usize) -> usize {
    qft_tile_bits_for(n, qft_max_tile_bits(), qft_min_tiles_log2(n, state_bits))
}

/// Log2 of the tiles each `2^n` chunk of a `state_bits`-qubit state owes the pool.
fn qft_min_tiles_log2(n: usize, state_bits: usize) -> usize {
    (cache::min_parallel_tiles().ilog2() as usize).saturating_sub(state_bits - n)
}

/// Tile width before the radix-4 pairing adjustment of [`qft_tile_bits_for`], which
/// the gathered high-stage passes do not need.
fn qft_widest_tile_bits_for(n: usize, max_tile_bits: usize, min_tiles_log2: usize) -> usize {
    max_tile_bits
        .min(n.saturating_sub(min_tiles_log2))
        .max(QFT_MIN_TILE_BITS.min(n))
}

fn qft_tile_bits_for(n: usize, max_tile_bits: usize, min_tiles_log2: usize) -> usize {
    let widest = qft_widest_tile_bits_for(n, max_tile_bits, min_tiles_log2);
    if widest > QFT_MIN_TILE_BITS && (n - widest) % 2 == 1 {
        widest - 1
    } else {
        widest
    }
}

/// Smallest QFT the gathered high-stage passes run on, where a whole-register transform
/// on an 8-thread host first leaves [`QFT_GATHER_MIN_STAGES`] stages above its tile.
/// Provisional until swept.
const QFT_GATHER_MIN_QUBITS: usize = 16;

/// Fewest stages above the tile that take the gathered passes. From three, one gathered
/// pass replaces at least two radix-2 or radix-4 passes; at two or fewer both take one.
const QFT_GATHER_MIN_STAGES: usize = 3;

/// Fewest amplitudes per gathered row: four cache lines, so a strided row read is a
/// burst rather than a lone line.
const QFT_GATHER_MIN_ROW: usize = 16;

/// Amplitudes one gathered group holds: half the tile budget, so the buffers of two SMT
/// siblings share one L2.
#[inline(always)]
fn qft_gather_len() -> usize {
    cache::tile_budget_bytes() / 2 / size_of::<Complex64>()
}

/// Most stages one gathered pass runs: 9 on a 256 KiB tile budget.
#[inline(always)]
fn qft_gather_max_bits() -> usize {
    (qft_gather_len() / QFT_GATHER_MIN_ROW).ilog2() as usize
}

/// Unit phases `exp(2 pi i x / 2^total_bits)` for `x < 2^total_bits`, conjugated for
/// `inverse`, as a product of three tables. The low table spans a gathered group, so a
/// phase below that length is one read.
struct QftGatherPhases {
    low: Box<[Complex64]>,
    mid: Box<[Complex64]>,
    high: Box<[Complex64]>,
    low_bits: usize,
    mid_bits: usize,
}

impl QftGatherPhases {
    fn new(total_bits: usize, inverse: bool) -> Self {
        let low_bits = (qft_gather_len().ilog2() as usize).min(total_bits);
        let rest = total_bits - low_bits;
        let mid_bits = rest.div_ceil(2);
        let scale = std::f64::consts::TAU / (1u64 << total_bits) as f64;
        let table = |bits: usize, shift: usize| -> Box<[Complex64]> {
            (0..1usize << bits)
                .map(|x| {
                    let w = Complex64::from_polar(1.0, (x << shift) as f64 * scale);
                    if inverse { w.conj() } else { w }
                })
                .collect()
        };
        Self {
            low: table(low_bits, 0),
            mid: table(mid_bits, low_bits),
            high: table(rest - mid_bits, low_bits + mid_bits),
            low_bits,
            mid_bits,
        }
    }

    #[inline(always)]
    fn phase(&self, x: usize) -> Complex64 {
        let mid = (x >> self.low_bits) & ((1 << self.mid_bits) - 1);
        self.low[x & (self.low.len() - 1)]
            * self.mid[mid]
            * self.high[x >> (self.low_bits + self.mid_bits)]
    }
}

fn qft_gather_phases(total_bits: usize, inverse: bool) -> &'static QftGatherPhases {
    const WIDTHS: usize = usize::BITS as usize;
    static FORWARD: [OnceLock<QftGatherPhases>; WIDTHS] = [const { OnceLock::new() }; WIDTHS];
    static INVERSE: [OnceLock<QftGatherPhases>; WIDTHS] = [const { OnceLock::new() }; WIDTHS];
    let slots = if inverse { &INVERSE } else { &FORWARD };
    slots[total_bits].get_or_init(|| QftGatherPhases::new(total_bits, inverse))
}

/// Per-stride QFT twiddles for the in-tile stages: entries `[s, 2s)` hold
/// `exp(2 pi i k / 2s) / sqrt(2)` for each power of two `s`, conjugated for `inverse`.
/// They are bit-identical to the whole-state table's entries `k * total / 2s`, which a
/// small stride would read one page apart.
fn qft_compact_twiddles(inverse: bool) -> &'static [Complex64] {
    static FORWARD: OnceLock<Box<[Complex64]>> = OnceLock::new();
    static INVERSE: OnceLock<Box<[Complex64]>> = OnceLock::new();
    let build = || {
        let inv_sqrt2 = std::f64::consts::FRAC_1_SQRT_2;
        let mut table = vec![Complex64::new(0.0, 0.0); 1 << qft_max_tile_bits()];
        let mut s = 1;
        while s < table.len() {
            for k in 0..s {
                let angle = std::f64::consts::TAU * k as f64 / (2 * s) as f64;
                let w = Complex64::from_polar(inv_sqrt2, angle);
                table[s + k] = if inverse { w.conj() } else { w };
            }
            s <<= 1;
        }
        table.into_boxed_slice()
    };
    if inverse {
        INVERSE.get_or_init(build)
    } else {
        FORWARD.get_or_init(build)
    }
}

#[inline(always)]
fn reverse_low_bits(x: usize, bits: usize) -> usize {
    if bits == 0 {
        return 0;
    }
    x.reverse_bits() >> (usize::BITS as usize - bits)
}

/// Bits per side of the square tiles the bit reversal transposes: 32 rows of 32
/// amplitudes, so a tile pair takes 32 KiB.
const BITREV_TILE_BITS: usize = 5;

/// Move tile `b` of a `bits`-bit state to tile `reverse(b)` and back, viewing an index
/// as `(a, b, c)` with `a` the high and `c` the low [`BITREV_TILE_BITS`] bits.
/// Amplitude `(a, b, c)` moves to `(rev(c), rev(b), rev(a))`, so every row read or
/// written is a contiguous run. Out of line, so its 32 KiB of tiles live only while a
/// task runs instead of in every frame of Rayon's recursive split.
///
/// # Safety
///
/// `state` must be valid for `1 << bits` amplitudes, `bits >= 2 * BITREV_TILE_BITS`,
/// and no other thread may access tile `b` or its reverse during the call.
#[inline(never)]
unsafe fn bit_reverse_tile_pair(state: *mut Complex64, bits: usize, b: usize) {
    const Q: usize = BITREV_TILE_BITS;
    const SIDE: usize = 1 << Q;
    let rb = reverse_low_bits(b, bits - 2 * Q);
    let row_stride = 1usize << (bits - Q);
    let rev: [usize; SIDE] = std::array::from_fn(|x| reverse_low_bits(x, Q));
    let mut tile_b = [Complex64::new(0.0, 0.0); SIDE * SIDE];
    let mut tile_rb = [Complex64::new(0.0, 0.0); SIDE * SIDE];
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        for a in 0..SIDE {
            let row = a * row_stride;
            std::ptr::copy_nonoverlapping(
                state.add(row + (b << Q)),
                tile_b.as_mut_ptr().add(a * SIDE),
                SIDE,
            );
            if rb != b {
                std::ptr::copy_nonoverlapping(
                    state.add(row + (rb << Q)),
                    tile_rb.as_mut_ptr().add(a * SIDE),
                    SIDE,
                );
            }
        }
        for a in 0..SIDE {
            let row = a * row_stride;
            let out_rb = state.add(row + (rb << Q));
            for c in 0..SIDE {
                *out_rb.add(c) = tile_b[rev[c] * SIDE + rev[a]];
            }
            if rb != b {
                let out_b = state.add(row + (b << Q));
                for c in 0..SIDE {
                    *out_b.add(c) = tile_rb[rev[c] * SIDE + rev[a]];
                }
            }
        }
    }
}

#[inline]
fn apply_bit_reverse_permutation(state: &mut [Complex64], bits: usize) {
    let len = state.len();
    debug_assert_eq!(len, 1usize << bits);
    if len <= 2 {
        return;
    }
    if bits < 2 * BITREV_TILE_BITS {
        for i in 0..len {
            let j = reverse_low_bits(i, bits);
            if j > i {
                state.swap(i, j);
            }
        }
        return;
    }

    let middle = bits - 2 * BITREV_TILE_BITS;
    let owns_pair = move |&b: &usize| reverse_low_bits(b, middle) >= b;

    #[cfg(feature = "parallel")]
    if bits >= PARALLEL_THRESHOLD_QUBITS {
        let ptr = SendPtr(state.as_mut_ptr());
        (0..1usize << middle)
            .into_par_iter()
            .filter(owns_pair)
            .for_each(move |b| {
                // SAFETY: reversal of the middle bits is an involution, so the pairs
                // `{b, reverse(b)}` with `b <= reverse(b)` partition the tiles and each
                // task owns its pair. Tiles are disjoint sets of rows of the state.
                unsafe { bit_reverse_tile_pair(ptr.as_complex_ptr(), bits, b) }
            });
        return;
    }

    let ptr = state.as_mut_ptr();
    for b in (0..1usize << middle).filter(owns_pair) {
        // SAFETY: `ptr` spans the whole state, borrowed mutably for this loop alone.
        unsafe { bit_reverse_tile_pair(ptr, bits, b) }
    }
}

/// Reverse the low `bits` index bits within every `2^bits` chunk of `state`.
fn apply_bit_reverse_chunks(state: &mut [Complex64], bits: usize) {
    let chunk = 1usize << bits;
    if state.len() == chunk {
        apply_bit_reverse_permutation(state, bits);
        return;
    }
    #[cfg(feature = "parallel")]
    if state.len() >= 1usize << PARALLEL_THRESHOLD_QUBITS {
        state
            .par_chunks_mut(chunk.max(MIN_PAR_ELEMS))
            .for_each(|task| {
                for c in task.chunks_mut(chunk) {
                    apply_bit_reverse_permutation(c, bits);
                }
            });
        return;
    }
    for c in state.chunks_mut(chunk) {
        apply_bit_reverse_permutation(c, bits);
    }
}

#[derive(Clone, Copy)]
pub(crate) struct BatchRzzGroup {
    pub(crate) table: [Complex64; BATCH_RZZ_TABLE_SIZE],
    pub(crate) q0s: [usize; BATCH_RZZ_GROUP_SIZE],
    pub(crate) q1s: [usize; BATCH_RZZ_GROUP_SIZE],
    pub(crate) len: usize,
}

pub(crate) fn build_batch_rzz_tables(
    edges: &[(usize, usize, f64)],
    groups: &mut [BatchRzzGroup; MAX_BATCH_RZZ_GROUPS],
) -> usize {
    let num_groups = edges.len().div_ceil(BATCH_RZZ_GROUP_SIZE);
    assert!(
        num_groups <= MAX_BATCH_RZZ_GROUPS,
        "BatchRzz carries {} edges, past the {} the group tables hold",
        edges.len(),
        BatchRzzData::MAX_EDGES
    );

    for (g, group) in groups.iter_mut().enumerate().take(num_groups) {
        let start = g * BATCH_RZZ_GROUP_SIZE;
        let end = (start + BATCH_RZZ_GROUP_SIZE).min(edges.len());
        let group_len = end - start;
        group.len = group_len;

        for (k, &(q0, q1, _)) in edges[start..end].iter().enumerate() {
            group.q0s[k] = q0;
            group.q1s[k] = q1;
        }

        for bits in 0..BATCH_RZZ_TABLE_SIZE {
            let mut angle = 0.0f64;
            for k in 0..group_len {
                let theta = edges[start + k].2;
                let parity = (bits >> k) & 1;
                angle += if parity == 0 {
                    -theta / 2.0
                } else {
                    theta / 2.0
                };
            }
            group.table[bits] = Complex64::from_polar(1.0, angle);
        }
    }
    num_groups
}

#[inline(always)]
fn extract_rzz_bits(i: usize, group: &BatchRzzGroup) -> usize {
    let mut bits = 0usize;
    for k in 0..group.len {
        bits |= (((i >> group.q0s[k]) ^ (i >> group.q1s[k])) & 1) << k;
    }
    bits
}

#[cfg(target_arch = "x86_64")]
#[derive(Clone, Copy)]
struct BatchRzzBmi2Group {
    table: [Complex64; BATCH_RZZ_BMI2_TABLE_SIZE],
    pext_mask: u64,
}

/// Pack `edges` into lookup groups of at most [`BATCH_RZZ_BMI2_MAX_UNIQUE`] qubits
/// each, placing every edge in the group it widens least. A group takes any number of
/// edges among its qubits, so a ring plus a random pairing fits in a few groups where
/// runs of consecutive edges would not. `None` when the edges need more than
/// [`BATCH_RZZ_BMI2_MAX_GROUPS`] groups.
#[cfg(target_arch = "x86_64")]
fn build_batch_rzz_bmi2_tables(
    edges: &[(usize, usize, f64)],
    groups: &mut [BatchRzzBmi2Group; BATCH_RZZ_BMI2_MAX_GROUPS],
) -> Option<usize> {
    debug_assert!(edges.len() <= u64::BITS as usize);
    let mut masks = [0u64; BATCH_RZZ_BMI2_MAX_GROUPS];
    let mut members = [0u64; BATCH_RZZ_BMI2_MAX_GROUPS];
    let mut num_groups = 0;
    for (e, &(q0, q1, _)) in edges.iter().enumerate() {
        let pair = (1u64 << q0) | (1u64 << q1);
        let best = (0..num_groups)
            .filter(|&g| (masks[g] | pair).count_ones() as usize <= BATCH_RZZ_BMI2_MAX_UNIQUE)
            .min_by_key(|&g| (masks[g] | pair).count_ones() - masks[g].count_ones());
        let g = match best {
            Some(g) => g,
            None if num_groups < BATCH_RZZ_BMI2_MAX_GROUPS => {
                num_groups += 1;
                num_groups - 1
            }
            None => return None,
        };
        masks[g] |= pair;
        members[g] |= 1u64 << e;
    }

    for (g, group) in groups.iter_mut().enumerate().take(num_groups) {
        let mask = masks[g];
        group.pext_mask = mask;
        let unique = mask.count_ones() as usize;
        // Flipping table bit `k` flips the parity of every edge on it, which
        // multiplies the phase by `e^{i theta}` or its conjugate.
        let mut incident = [[(0u8, ONE); BatchRzzData::MAX_EDGES]; BATCH_RZZ_BMI2_MAX_UNIQUE];
        let mut degree = [0usize; BATCH_RZZ_BMI2_MAX_UNIQUE];
        let mut base = ONE;
        let mut bits = members[g];
        while bits != 0 {
            let (q0, q1, theta) = edges[bits.trailing_zeros() as usize];
            let p0 = (mask & ((1u64 << q0) - 1)).count_ones() as usize;
            let p1 = (mask & ((1u64 << q1) - 1)).count_ones() as usize;
            let flip = Complex64::from_polar(1.0, theta);
            incident[p0][degree[p0]] = (p1 as u8, flip);
            degree[p0] += 1;
            incident[p1][degree[p1]] = (p0 as u8, flip);
            degree[p1] += 1;
            base *= Complex64::from_polar(1.0, -theta / 2.0);
            bits &= bits - 1;
        }

        group.table[0] = base;
        for c in 1..1usize << unique {
            let k = c.trailing_zeros() as usize;
            let prev = c & (c - 1);
            let mut phase = group.table[prev];
            for &(other, flip) in &incident[k][..degree[k]] {
                phase *= if (prev >> other) & 1 == 0 {
                    flip
                } else {
                    flip.conj()
                };
            }
            group.table[c] = phase;
        }
    }
    Some(num_groups)
}

#[cfg(target_arch = "x86_64")]
#[target_feature(enable = "bmi2,fma")]
unsafe fn batch_rzz_tile_bmi2(
    tile: &mut [Complex64],
    base_idx: usize,
    groups: &[BatchRzzBmi2Group; BATCH_RZZ_BMI2_MAX_GROUPS],
    num_groups: usize,
) {
    use std::arch::x86_64::_pext_u64;
    let one = Complex64::new(1.0, 0.0);
    for (j, amp) in tile.iter_mut().enumerate() {
        let i = base_idx + j;
        let mut combined = one;
        for group in groups.iter().take(num_groups) {
            let bits = _pext_u64(i as u64, group.pext_mask) as usize;
            combined *= group.table[bits];
        }
        *amp *= combined;
    }
}

#[derive(Clone, Copy)]
pub(crate) struct BatchPhaseGroup {
    pub(crate) table: [Complex64; BATCH_PHASE_TABLE_SIZE],
    pub(crate) shifts: [usize; BATCH_PHASE_GROUP_SIZE],
    pub(crate) len: usize,
    pub(crate) pext_mask: u64,
}

const ONE: Complex64 = Complex64::new(1.0, 0.0);

impl BatchRzzGroup {
    const EMPTY: Self = Self {
        table: [ONE; BATCH_RZZ_TABLE_SIZE],
        q0s: [0; BATCH_RZZ_GROUP_SIZE],
        q1s: [0; BATCH_RZZ_GROUP_SIZE],
        len: 0,
    };
}

#[cfg(target_arch = "x86_64")]
impl BatchRzzBmi2Group {
    const EMPTY: Self = Self {
        table: [ONE; BATCH_RZZ_BMI2_TABLE_SIZE],
        pext_mask: 0,
    };
}

impl BatchPhaseGroup {
    const EMPTY: Self = Self {
        table: [ONE; BATCH_PHASE_TABLE_SIZE],
        shifts: [0; BATCH_PHASE_GROUP_SIZE],
        len: 0,
        pext_mask: 0,
    };
}

/// Lookup tables the batched diagonal kernels fill per gate.
///
/// Heap-resident and reused across gates, so a gate rewrites only the table prefix its
/// groups index instead of initialising a 64 KB array on the stack. Each table allocates
/// on the first gate of its kind.
#[derive(Default)]
pub(crate) struct BatchTableScratch {
    rzz: Option<Box<[BatchRzzGroup; MAX_BATCH_RZZ_GROUPS]>>,
    #[cfg(target_arch = "x86_64")]
    rzz_bmi2: Option<Box<[BatchRzzBmi2Group; BATCH_RZZ_BMI2_MAX_GROUPS]>>,
    phase: Option<Box<[BatchPhaseGroup; MAX_BATCH_PHASE_GROUPS]>>,
    diag: Option<Box<DiagBatchTables>>,
}

impl BatchTableScratch {
    fn rzz(&mut self) -> &mut [BatchRzzGroup; MAX_BATCH_RZZ_GROUPS] {
        self.rzz
            .get_or_insert_with(|| Box::new([BatchRzzGroup::EMPTY; MAX_BATCH_RZZ_GROUPS]))
    }

    #[cfg(target_arch = "x86_64")]
    fn rzz_bmi2(&mut self) -> &mut [BatchRzzBmi2Group; BATCH_RZZ_BMI2_MAX_GROUPS] {
        self.rzz_bmi2
            .get_or_insert_with(|| Box::new([BatchRzzBmi2Group::EMPTY; BATCH_RZZ_BMI2_MAX_GROUPS]))
    }

    fn phase(&mut self) -> &mut [BatchPhaseGroup; MAX_BATCH_PHASE_GROUPS] {
        self.phase
            .get_or_insert_with(|| Box::new([BatchPhaseGroup::EMPTY; MAX_BATCH_PHASE_GROUPS]))
    }

    fn diag(&mut self) -> &mut DiagBatchTables {
        self.diag
            .get_or_insert_with(|| Box::new(DiagBatchTables::empty()))
    }
}

/// Build lookup tables for batched controlled-phase application.
///
/// Partitions `phases` into groups of up to 10. For each group, builds a
/// 1024-entry table mapping every combination of target-bit patterns to the
/// combined phase. Returns the number of groups filled.
pub(crate) fn build_batch_phase_tables(
    phases: &[(usize, Complex64)],
    groups: &mut [BatchPhaseGroup; MAX_BATCH_PHASE_GROUPS],
) -> usize {
    let one = Complex64::new(1.0, 0.0);
    let num_groups = phases.len().div_ceil(BATCH_PHASE_GROUP_SIZE);
    assert!(
        num_groups <= MAX_BATCH_PHASE_GROUPS,
        "BatchPhase carries {} entries, past the {} the group tables hold",
        phases.len(),
        BatchPhaseData::MAX_PHASES
    );

    for (g, group) in groups.iter_mut().enumerate().take(num_groups) {
        let start = g * BATCH_PHASE_GROUP_SIZE;
        let end = (start + BATCH_PHASE_GROUP_SIZE).min(phases.len());
        let group_len = end - start;
        group.len = group_len;

        let mut sorted: [(usize, Complex64); BATCH_PHASE_GROUP_SIZE] =
            [(0, one); BATCH_PHASE_GROUP_SIZE];
        sorted[..group_len].copy_from_slice(&phases[start..end]);
        sorted[..group_len].sort_unstable_by_key(|&(q, _)| q);

        let mut mask = 0u64;
        for (j, &(q, _)) in sorted.iter().enumerate().take(group_len) {
            group.shifts[j] = q;
            mask |= 1u64 << q;
        }
        group.pext_mask = mask;

        group.table[0] = one;
        for (j, &(_, phase)) in sorted.iter().enumerate().take(group_len) {
            let prev_entries = 1usize << j;
            for k in 0..prev_entries {
                group.table[k | prev_entries] = group.table[k] * phase;
            }
        }
    }
    num_groups
}

/// Extract bits at specified positions from `idx` and pack them into a
/// contiguous low-order pattern for table lookup.
#[inline(always)]
fn extract_bits(idx: usize, shifts: &[usize; BATCH_PHASE_GROUP_SIZE], len: usize) -> usize {
    let mut bits = 0usize;
    for (j, &shift) in shifts.iter().enumerate().take(len) {
        bits |= ((idx >> shift) & 1) << j;
    }
    bits
}

/// Pre-built LUT data for `apply_diagonal_batch`.
///
/// Returned by [`build_diagonal_batch_tables`]. The statevector backend consumes this
/// directly; the GPU backend uploads the tables + per-group metadata to the device and
/// launches a batched kernel.
pub(crate) struct DiagBatchTables {
    pub(crate) tables: [[Complex64; DIAG_BATCH_TABLE_SIZE]; MAX_DIAG_BATCH_GROUPS],
    pub(crate) unique_qubits: SmallVec<[usize; 32]>,
    pub(crate) group_sizes: [usize; MAX_DIAG_BATCH_GROUPS],
    /// Per-group PEXT masks used by the x86_64 BMI2 CPU path only. The GPU path
    /// rebuilds the equivalent `group_shifts` array from `unique_qubits + group_sizes`
    /// because PTX has no PEXT; non-x86_64 CPU builds don't use this either. Kept
    /// populated regardless (the fill is trivial) so the struct has one layout.
    #[cfg_attr(not(target_arch = "x86_64"), allow(dead_code))]
    pub(crate) group_pext_masks: [u64; MAX_DIAG_BATCH_GROUPS],
    pub(crate) num_groups: usize,
}

impl DiagBatchTables {
    fn empty() -> Self {
        Self {
            tables: [[ONE; DIAG_BATCH_TABLE_SIZE]; MAX_DIAG_BATCH_GROUPS],
            unique_qubits: SmallVec::new(),
            group_sizes: [0; MAX_DIAG_BATCH_GROUPS],
            group_pext_masks: [0; MAX_DIAG_BATCH_GROUPS],
            num_groups: 0,
        }
    }
}

/// Root of `x` in the qubit-affinity forest, with path halving.
fn find_root(parent: &mut [u8; 64], mut x: u8) -> u8 {
    while parent[x as usize] != x {
        parent[x as usize] = parent[parent[x as usize] as usize];
        x = parent[x as usize];
    }
    x
}

/// Group `entries` by qubit affinity and build per-group phase LUTs.
///
/// A two-qubit entry folds its phase into a single group's table, so both of its qubits
/// have to share a group. Qubits joined by such entries therefore form components, and the
/// components are packed into groups as units, largest first.
///
/// Returns `None` when a component is wider than `DIAG_BATCH_MAX_QUBITS_PER_GROUP` or the
/// packing needs more than `MAX_DIAG_BATCH_GROUPS` groups. Both cases need the per-element
/// fallback path. Returns `Some` with `num_groups == 0` when `entries` is empty (no-op).
#[cfg(any(test, feature = "gpu"))]
pub(crate) fn build_diagonal_batch_tables(entries: &[DiagEntry]) -> Option<DiagBatchTables> {
    let mut out = DiagBatchTables::empty();
    fill_diagonal_batch_tables(entries, &mut out).then_some(out)
}

/// [`build_diagonal_batch_tables`] into a reused `out`, touching only the table prefix
/// each group indexes. Returns false when the entries are not groupable; `out` is then
/// stale and the caller takes the per-element fallback.
fn fill_diagonal_batch_tables(entries: &[DiagEntry], out: &mut DiagBatchTables) -> bool {
    let mut unique_qubits = SmallVec::<[usize; 32]>::new();
    let add_qubit = |q: usize, uq: &mut SmallVec<[usize; 32]>| {
        if !uq.contains(&q) {
            uq.push(q);
        }
    };
    for e in entries {
        match e {
            DiagEntry::Phase1q { qubit, .. } => add_qubit(*qubit, &mut unique_qubits),
            DiagEntry::Phase2q { q0, q1, .. } | DiagEntry::Parity2q { q0, q1, .. } => {
                add_qubit(*q0, &mut unique_qubits);
                add_qubit(*q1, &mut unique_qubits);
            }
        }
    }
    unique_qubits.sort_unstable();
    let num_unique = unique_qubits.len();

    if num_unique == 0 {
        out.unique_qubits.clear();
        out.group_sizes = [0; MAX_DIAG_BATCH_GROUPS];
        out.group_pext_masks = [0; MAX_DIAG_BATCH_GROUPS];
        out.num_groups = 0;
        return true;
    }

    let mut dense = [0u8; 64];
    for (idx, &q) in unique_qubits.iter().enumerate() {
        dense[q] = idx as u8;
    }
    let mut parent: [u8; 64] = std::array::from_fn(|i| i as u8);
    for e in entries {
        if let DiagEntry::Phase2q { q0, q1, .. } | DiagEntry::Parity2q { q0, q1, .. } = e {
            let ra = find_root(&mut parent, dense[*q0]);
            let rb = find_root(&mut parent, dense[*q1]);
            if ra != rb {
                parent[ra as usize] = rb;
            }
        }
    }

    let mut comp_of = [0u8; 64];
    let mut comp_len: SmallVec<[usize; 64]> = SmallVec::new();
    let mut slot_of = [u8::MAX; 64];
    for &q in &unique_qubits {
        let root = find_root(&mut parent, dense[q]) as usize;
        if slot_of[root] == u8::MAX {
            slot_of[root] = comp_len.len() as u8;
            comp_len.push(0);
        }
        comp_of[q] = slot_of[root];
        comp_len[slot_of[root] as usize] += 1;
    }
    let mut order: SmallVec<[u8; 64]> = (0..comp_len.len() as u8).collect();
    order.sort_unstable_by_key(|&c| std::cmp::Reverse(comp_len[c as usize]));

    let mut group_of_comp = [0u8; 64];
    let mut group_sizes = [0usize; MAX_DIAG_BATCH_GROUPS];
    for &c in &order {
        let len = comp_len[c as usize];
        let Some(g) = group_sizes
            .iter()
            .position(|&used| used + len <= DIAG_BATCH_MAX_QUBITS_PER_GROUP)
        else {
            return false;
        };
        group_of_comp[c as usize] = g as u8;
        group_sizes[g] += len;
    }
    let num_groups = group_sizes.iter().filter(|&&size| size > 0).count();

    let mut qubit_group = [0u8; 64];
    let mut qubit_pos = [0u8; 64];
    let mut group_pext_masks = [0u64; MAX_DIAG_BATCH_GROUPS];
    let sorted = std::mem::take(&mut unique_qubits);
    for (g, mask) in group_pext_masks.iter_mut().enumerate().take(num_groups) {
        let members = sorted
            .iter()
            .filter(|&&q| group_of_comp[comp_of[q] as usize] as usize == g);
        for (pos, &q) in members.enumerate() {
            qubit_group[q] = g as u8;
            qubit_pos[q] = pos as u8;
            *mask |= 1u64 << q;
            unique_qubits.push(q);
        }
    }

    let tables = &mut out.tables;
    for (table, &k) in tables.iter_mut().zip(&group_sizes).take(num_groups) {
        table[..1 << k].fill(ONE);
    }

    for e in entries {
        match e {
            DiagEntry::Phase1q { qubit, d0, d1 } => {
                let g = qubit_group[*qubit] as usize;
                let p = qubit_pos[*qubit] as usize;
                let k = group_sizes[g];
                for (bits, entry) in tables[g][..1 << k].iter_mut().enumerate() {
                    if (bits >> p) & 1 == 1 {
                        *entry *= d1;
                    } else {
                        *entry *= d0;
                    }
                }
            }
            DiagEntry::Phase2q { q0, q1, phase } => {
                let g = qubit_group[*q0] as usize;
                let p0 = qubit_pos[*q0] as usize;
                let p1 = qubit_pos[*q1] as usize;
                debug_assert_eq!(g, qubit_group[*q1] as usize);
                let k = group_sizes[g];
                for (bits, entry) in tables[g][..1 << k].iter_mut().enumerate() {
                    if ((bits >> p0) & 1 == 1) && ((bits >> p1) & 1 == 1) {
                        *entry *= phase;
                    }
                }
            }
            DiagEntry::Parity2q { q0, q1, same, diff } => {
                let g = qubit_group[*q0] as usize;
                let p0 = qubit_pos[*q0] as usize;
                let p1 = qubit_pos[*q1] as usize;
                debug_assert_eq!(g, qubit_group[*q1] as usize);
                let k = group_sizes[g];
                for (bits, entry) in tables[g][..1 << k].iter_mut().enumerate() {
                    let parity = ((bits >> p0) ^ (bits >> p1)) & 1;
                    *entry *= if parity == 0 { *same } else { *diff };
                }
            }
        }
    }

    out.unique_qubits = unique_qubits;
    out.group_sizes = group_sizes;
    out.group_pext_masks = group_pext_masks;
    out.num_groups = num_groups;
    true
}

#[cfg(target_arch = "x86_64")]
#[target_feature(enable = "bmi2,fma")]
unsafe fn apply_batch_phase_bmi2(
    state: &mut [Complex64],
    ctrl_bit: usize,
    groups: &[BatchPhaseGroup; MAX_BATCH_PHASE_GROUPS],
    num_groups: usize,
) {
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        use std::arch::x86_64::_pext_u64;
        let ctrl_mask = 1usize << ctrl_bit;
        let half = state.len() >> 1;
        for k in 0..half {
            let i = insert_zero_bit(k, ctrl_bit) | ctrl_mask;
            let mut combined = Complex64::new(1.0, 0.0);
            for group in groups.iter().take(num_groups) {
                let bits = _pext_u64(i as u64, group.pext_mask) as usize;
                combined *= group.table[bits];
            }
            // SAFETY: insert_zero_bit(k, ctrl_bit) | ctrl_mask < state.len() for k < half
            *state.get_unchecked_mut(i) *= combined;
        }
    }
}

#[cfg(all(target_arch = "x86_64", feature = "parallel"))]
#[target_feature(enable = "bmi2,fma")]
unsafe fn batch_phase_tile_bmi2(
    tile: &mut [Complex64],
    base_idx: usize,
    ctrl_bit: usize,
    groups: &[BatchPhaseGroup; MAX_BATCH_PHASE_GROUPS],
    num_groups: usize,
) {
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        use std::arch::x86_64::_pext_u64;
        let ctrl_mask = 1usize << ctrl_bit;
        let tile_bits = tile.len().trailing_zeros() as usize;

        if ctrl_bit >= tile_bits {
            if (base_idx & ctrl_mask) == 0 {
                return;
            }
            for (j, amp) in tile.iter_mut().enumerate() {
                let i = base_idx + j;
                let mut combined = Complex64::new(1.0, 0.0);
                for group in groups.iter().take(num_groups) {
                    let bits = _pext_u64(i as u64, group.pext_mask) as usize;
                    combined *= group.table[bits];
                }
                *amp *= combined;
            }
        } else {
            let half = tile.len() >> 1;
            for k in 0..half {
                let j = insert_zero_bit(k, ctrl_bit) | ctrl_mask;
                let i = base_idx + j;
                let mut combined = Complex64::new(1.0, 0.0);
                for group in groups.iter().take(num_groups) {
                    let bits = _pext_u64(i as u64, group.pext_mask) as usize;
                    combined *= group.table[bits];
                }
                // SAFETY: j < tile.len() guaranteed by insert_zero_bit mapping
                *tile.get_unchecked_mut(j) *= combined;
            }
        }
    }
}

/// Run one FFT stage in parallel.
///
/// Large group counts split by group. High-stride stages split inside each
/// group so all cores still get work.
#[cfg(feature = "parallel")]
#[inline]
fn fft_stage_par(
    state: &mut [Complex64],
    stride: usize,
    block_size: usize,
    twiddle_step: usize,
    twiddles_scaled: &[Complex64],
    inv_sqrt2: f64,
) {
    let total = state.len();
    let num_groups = total / block_size;

    let apply_pairs = |lo: &mut [Complex64], hi: &mut [Complex64], k_base: usize| {
        debug_assert_eq!(lo.len(), hi.len());
        for j in 0..lo.len() {
            let k = k_base + j;
            let w = twiddles_scaled[k * twiddle_step];
            let (out_lo, out_hi) = radix2_butterfly_values(lo[j], hi[j], w, inv_sqrt2);
            lo[j] = out_lo;
            hi[j] = out_hi;
        }
    };

    if num_groups >= 4 && block_size >= MIN_PAR_ELEMS {
        // Many large groups: split state by group.
        state.par_chunks_mut(block_size).for_each(|group| {
            let (lo, hi) = group.split_at_mut(stride);
            apply_pairs(lo, hi, 0);
        });
    } else if block_size < MIN_PAR_ELEMS {
        // Small groups: bundle work so each Rayon job has enough elements.
        let task_size = MIN_PAR_ELEMS;
        state.par_chunks_mut(task_size).for_each(|task_chunk| {
            let mut g = 0;
            while g + block_size <= task_chunk.len() {
                let group = &mut task_chunk[g..g + block_size];
                let (lo, hi) = group.split_at_mut(stride);
                apply_pairs(lo, hi, 0);
                g += block_size;
            }
        });
    } else {
        // Few groups: split each large group into zipped chunks.
        const MIN_PAR_PAIRS: usize = MIN_PAR_ELEMS / 2;
        for group_start in (0..total).step_by(block_size) {
            let group = &mut state[group_start..group_start + block_size];
            let (lo, hi) = group.split_at_mut(stride);
            let sub = MIN_PAR_PAIRS.min(stride).max(1);
            lo.par_chunks_mut(sub)
                .zip(hi.par_chunks_mut(sub))
                .enumerate()
                .for_each(|(chunk_idx, (lo_c, hi_c))| {
                    let k_base = chunk_idx * sub;
                    apply_pairs(lo_c, hi_c, k_base);
                });
        }
    }
}

/// Scalar radix-2 DIF butterfly. `w` includes the `1/sqrt(2)` factor.
#[inline(always)]
fn radix2_butterfly_values(
    a: Complex64,
    b: Complex64,
    w: Complex64,
    inv_sqrt2: f64,
) -> (Complex64, Complex64) {
    let lo = Complex64::new((a.re + b.re) * inv_sqrt2, (a.im + b.im) * inv_sqrt2);
    let hi = (a - b) * w;
    (lo, hi)
}

/// Apply one full radix-2 DIF stage in place to `state` at the given stride.
/// Butterfly `k` takes `twiddles_scaled[k * twiddle_step]`: step `total / (2 * stride)`
/// on the whole-state table, 1 on a per-stride run.
#[inline]
fn run_radix2_stage_seq(
    state: &mut [Complex64],
    stride: usize,
    twiddles_scaled: &[Complex64],
    twiddle_step: usize,
    inv_sqrt2: f64,
) {
    let block_size = stride << 1;
    let total = state.len();
    let mut group_start = 0;
    while group_start < total {
        for k in 0..stride {
            let i0 = group_start + k;
            let i1 = i0 + stride;
            let w = twiddles_scaled[k * twiddle_step];
            let (lo, hi) = radix2_butterfly_values(state[i0], state[i1], w, inv_sqrt2);
            state[i0] = lo;
            state[i1] = hi;
        }
        group_start += block_size;
    }
}

/// Scalar radix-4 butterfly math. Twiddles include the `1/sqrt(2)` factor.
#[inline(always)]
#[allow(clippy::too_many_arguments)]
fn radix4_butterfly_values(
    a: Complex64,
    b: Complex64,
    c: Complex64,
    d: Complex64,
    w_2s_k: Complex64,
    w_2s_kps: Complex64,
    w_s_k: Complex64,
    inv_sqrt2: f64,
) -> (Complex64, Complex64, Complex64, Complex64) {
    let t_ac_sum = Complex64::new((a.re + c.re) * inv_sqrt2, (a.im + c.im) * inv_sqrt2);
    let t_ac_diff = (a - c) * w_2s_k;
    let t_bd_sum = Complex64::new((b.re + d.re) * inv_sqrt2, (b.im + d.im) * inv_sqrt2);
    let t_bd_diff = (b - d) * w_2s_kps;

    let oa = Complex64::new(
        (t_ac_sum.re + t_bd_sum.re) * inv_sqrt2,
        (t_ac_sum.im + t_bd_sum.im) * inv_sqrt2,
    );
    let ob = (t_ac_sum - t_bd_sum) * w_s_k;
    let oc = Complex64::new(
        (t_ac_diff.re + t_bd_diff.re) * inv_sqrt2,
        (t_ac_diff.im + t_bd_diff.im) * inv_sqrt2,
    );
    let od = (t_ac_diff - t_bd_diff) * w_s_k;
    (oa, ob, oc, od)
}

/// Apply one fused radix-4 butterfly at state index `i`.
#[inline(always)]
#[cfg_attr(target_arch = "aarch64", allow(dead_code))]
fn radix4_butterfly_scalar(
    state: &mut [Complex64],
    i: usize,
    inner_stride: usize,
    w_2s_k: Complex64,
    w_2s_kps: Complex64,
    w_s_k: Complex64,
    inv_sqrt2: f64,
) {
    let s = inner_stride;
    let (oa, ob, oc, od) = radix4_butterfly_values(
        state[i],
        state[i + s],
        state[i + 2 * s],
        state[i + 3 * s],
        w_2s_k,
        w_2s_kps,
        w_s_k,
        inv_sqrt2,
    );
    state[i] = oa;
    state[i + s] = ob;
    state[i + 2 * s] = oc;
    state[i + 3 * s] = od;
}

/// Apply one radix-4 butterfly across four disjoint quarter slices.
#[inline(always)]
#[cfg_attr(
    any(target_arch = "aarch64", not(feature = "parallel")),
    allow(dead_code)
)]
#[allow(clippy::too_many_arguments)]
fn radix4_butterfly_quartet_scalar(
    qa: &mut [Complex64],
    qb: &mut [Complex64],
    qc: &mut [Complex64],
    qd: &mut [Complex64],
    j: usize,
    w_2s_k: Complex64,
    w_2s_kps: Complex64,
    w_s_k: Complex64,
    inv_sqrt2: f64,
) {
    let (oa, ob, oc, od) = radix4_butterfly_values(
        qa[j], qb[j], qc[j], qd[j], w_2s_k, w_2s_kps, w_s_k, inv_sqrt2,
    );
    qa[j] = oa;
    qb[j] = ob;
    qc[j] = oc;
    qd[j] = od;
}

/// FMA variant for one radix-4 butterfly.
#[cfg(target_arch = "x86_64")]
#[target_feature(enable = "fma")]
unsafe fn radix4_butterfly_fma(
    base_ptr: *mut f64,
    inner_stride: usize,
    w_2s_k: Complex64,
    w_2s_kps: Complex64,
    w_s_k: Complex64,
) {
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        use std::arch::x86_64::{
            _mm_add_pd, _mm_loadu_pd, _mm_mul_pd, _mm_set1_pd, _mm_storeu_pd, _mm_sub_pd,
        };
        let s = inner_stride;
        let inv_sqrt2_v = _mm_set1_pd(std::f64::consts::FRAC_1_SQRT_2);

        let a_ptr = base_ptr;
        let b_ptr = base_ptr.add(s * 2);
        let c_ptr = base_ptr.add(2 * s * 2);
        let d_ptr = base_ptr.add(3 * s * 2);

        let a = _mm_loadu_pd(a_ptr);
        let b = _mm_loadu_pd(b_ptr);
        let c = _mm_loadu_pd(c_ptr);
        let d = _mm_loadu_pd(d_ptr);

        let w_2s_k_rr = _mm_set1_pd(w_2s_k.re);
        let w_2s_k_ii = _mm_set1_pd(w_2s_k.im);
        let w_2s_kps_rr = _mm_set1_pd(w_2s_kps.re);
        let w_2s_kps_ii = _mm_set1_pd(w_2s_kps.im);
        let w_s_k_rr = _mm_set1_pd(w_s_k.re);
        let w_s_k_ii = _mm_set1_pd(w_s_k.im);

        let t_ac_sum = _mm_mul_pd(_mm_add_pd(a, c), inv_sqrt2_v);
        let t_ac_diff = simd::complex_mul_fma(w_2s_k_rr, w_2s_k_ii, _mm_sub_pd(a, c));
        let t_bd_sum = _mm_mul_pd(_mm_add_pd(b, d), inv_sqrt2_v);
        let t_bd_diff = simd::complex_mul_fma(w_2s_kps_rr, w_2s_kps_ii, _mm_sub_pd(b, d));

        _mm_storeu_pd(
            a_ptr,
            _mm_mul_pd(_mm_add_pd(t_ac_sum, t_bd_sum), inv_sqrt2_v),
        );
        _mm_storeu_pd(
            b_ptr,
            simd::complex_mul_fma(w_s_k_rr, w_s_k_ii, _mm_sub_pd(t_ac_sum, t_bd_sum)),
        );
        _mm_storeu_pd(
            c_ptr,
            _mm_mul_pd(_mm_add_pd(t_ac_diff, t_bd_diff), inv_sqrt2_v),
        );
        _mm_storeu_pd(
            d_ptr,
            simd::complex_mul_fma(w_s_k_rr, w_s_k_ii, _mm_sub_pd(t_ac_diff, t_bd_diff)),
        );
    }
}

/// AVX2-256 SIMD variant: processes TWO radix-4 butterflies per call (at
/// consecutive `k` and `k+1`, four complex amplitudes each, eight elements
/// total). Each load/store is one 256-bit op covering two adjacent complex
/// values (which are 32 bytes contiguous in memory). Twiddles are loaded
/// scalar-by-scalar since the inner twiddle step is 2.
#[cfg(target_arch = "x86_64")]
#[target_feature(enable = "avx2,fma")]
#[allow(clippy::too_many_arguments)]
unsafe fn radix4_butterfly_pair_avx2fma(
    base_ptr: *mut f64,
    inner_stride: usize,
    w_2s_k0: Complex64,
    w_2s_k1: Complex64,
    w_2s_kps0: Complex64,
    w_2s_kps1: Complex64,
    w_s_k0: Complex64,
    w_s_k1: Complex64,
) {
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        use std::arch::x86_64::{
            _mm256_add_pd, _mm256_loadu_pd, _mm256_mul_pd, _mm256_set_pd, _mm256_set1_pd,
            _mm256_storeu_pd, _mm256_sub_pd,
        };
        let s = inner_stride;
        let inv_sqrt2_v = _mm256_set1_pd(std::f64::consts::FRAC_1_SQRT_2);

        let a_ptr = base_ptr;
        let b_ptr = base_ptr.add(s * 2);
        let c_ptr = base_ptr.add(2 * s * 2);
        let d_ptr = base_ptr.add(3 * s * 2);

        let a = _mm256_loadu_pd(a_ptr);
        let b = _mm256_loadu_pd(b_ptr);
        let c = _mm256_loadu_pd(c_ptr);
        let d = _mm256_loadu_pd(d_ptr);

        // c_rr layout: [w0.re, w0.re, w1.re, w1.re], broadcasts each twiddle's
        // real part to the two lanes corresponding to that butterfly's complex.
        let w_2s_k_rr = _mm256_set_pd(w_2s_k1.re, w_2s_k1.re, w_2s_k0.re, w_2s_k0.re);
        let w_2s_k_ii = _mm256_set_pd(w_2s_k1.im, w_2s_k1.im, w_2s_k0.im, w_2s_k0.im);
        let w_2s_kps_rr = _mm256_set_pd(w_2s_kps1.re, w_2s_kps1.re, w_2s_kps0.re, w_2s_kps0.re);
        let w_2s_kps_ii = _mm256_set_pd(w_2s_kps1.im, w_2s_kps1.im, w_2s_kps0.im, w_2s_kps0.im);
        let w_s_k_rr = _mm256_set_pd(w_s_k1.re, w_s_k1.re, w_s_k0.re, w_s_k0.re);
        let w_s_k_ii = _mm256_set_pd(w_s_k1.im, w_s_k1.im, w_s_k0.im, w_s_k0.im);

        let t_ac_sum = _mm256_mul_pd(_mm256_add_pd(a, c), inv_sqrt2_v);
        let t_ac_diff = simd::complex_mul_avx2fma(w_2s_k_rr, w_2s_k_ii, _mm256_sub_pd(a, c));
        let t_bd_sum = _mm256_mul_pd(_mm256_add_pd(b, d), inv_sqrt2_v);
        let t_bd_diff = simd::complex_mul_avx2fma(w_2s_kps_rr, w_2s_kps_ii, _mm256_sub_pd(b, d));

        _mm256_storeu_pd(
            a_ptr,
            _mm256_mul_pd(_mm256_add_pd(t_ac_sum, t_bd_sum), inv_sqrt2_v),
        );
        _mm256_storeu_pd(
            b_ptr,
            simd::complex_mul_avx2fma(w_s_k_rr, w_s_k_ii, _mm256_sub_pd(t_ac_sum, t_bd_sum)),
        );
        _mm256_storeu_pd(
            c_ptr,
            _mm256_mul_pd(_mm256_add_pd(t_ac_diff, t_bd_diff), inv_sqrt2_v),
        );
        _mm256_storeu_pd(
            d_ptr,
            simd::complex_mul_avx2fma(w_s_k_rr, w_s_k_ii, _mm256_sub_pd(t_ac_diff, t_bd_diff)),
        );
    }
}

#[cfg(target_arch = "x86_64")]
#[target_feature(enable = "avx2,fma")]
#[cfg_attr(not(feature = "parallel"), allow(dead_code))]
#[allow(clippy::too_many_arguments)]
unsafe fn radix4_butterfly_pair_slices_avx2fma(
    a_ptr: *mut f64,
    b_ptr: *mut f64,
    c_ptr: *mut f64,
    d_ptr: *mut f64,
    w_2s_k0: Complex64,
    w_2s_k1: Complex64,
    w_2s_kps0: Complex64,
    w_2s_kps1: Complex64,
    w_s_k0: Complex64,
    w_s_k1: Complex64,
) {
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        use std::arch::x86_64::{
            _mm256_add_pd, _mm256_loadu_pd, _mm256_mul_pd, _mm256_set_pd, _mm256_set1_pd,
            _mm256_storeu_pd, _mm256_sub_pd,
        };
        let inv_sqrt2_v = _mm256_set1_pd(std::f64::consts::FRAC_1_SQRT_2);

        let a = _mm256_loadu_pd(a_ptr);
        let b = _mm256_loadu_pd(b_ptr);
        let c = _mm256_loadu_pd(c_ptr);
        let d = _mm256_loadu_pd(d_ptr);

        let w_2s_k_rr = _mm256_set_pd(w_2s_k1.re, w_2s_k1.re, w_2s_k0.re, w_2s_k0.re);
        let w_2s_k_ii = _mm256_set_pd(w_2s_k1.im, w_2s_k1.im, w_2s_k0.im, w_2s_k0.im);
        let w_2s_kps_rr = _mm256_set_pd(w_2s_kps1.re, w_2s_kps1.re, w_2s_kps0.re, w_2s_kps0.re);
        let w_2s_kps_ii = _mm256_set_pd(w_2s_kps1.im, w_2s_kps1.im, w_2s_kps0.im, w_2s_kps0.im);
        let w_s_k_rr = _mm256_set_pd(w_s_k1.re, w_s_k1.re, w_s_k0.re, w_s_k0.re);
        let w_s_k_ii = _mm256_set_pd(w_s_k1.im, w_s_k1.im, w_s_k0.im, w_s_k0.im);

        let t_ac_sum = _mm256_mul_pd(_mm256_add_pd(a, c), inv_sqrt2_v);
        let t_ac_diff = simd::complex_mul_avx2fma(w_2s_k_rr, w_2s_k_ii, _mm256_sub_pd(a, c));
        let t_bd_sum = _mm256_mul_pd(_mm256_add_pd(b, d), inv_sqrt2_v);
        let t_bd_diff = simd::complex_mul_avx2fma(w_2s_kps_rr, w_2s_kps_ii, _mm256_sub_pd(b, d));

        _mm256_storeu_pd(
            a_ptr,
            _mm256_mul_pd(_mm256_add_pd(t_ac_sum, t_bd_sum), inv_sqrt2_v),
        );
        _mm256_storeu_pd(
            b_ptr,
            simd::complex_mul_avx2fma(w_s_k_rr, w_s_k_ii, _mm256_sub_pd(t_ac_sum, t_bd_sum)),
        );
        _mm256_storeu_pd(
            c_ptr,
            _mm256_mul_pd(_mm256_add_pd(t_ac_diff, t_bd_diff), inv_sqrt2_v),
        );
        _mm256_storeu_pd(
            d_ptr,
            simd::complex_mul_avx2fma(w_s_k_rr, w_s_k_ii, _mm256_sub_pd(t_ac_diff, t_bd_diff)),
        );
    }
}

#[cfg(target_arch = "aarch64")]
#[inline(always)]
unsafe fn neon_twiddle_parts(
    w: Complex64,
) -> (
    std::arch::aarch64::float64x2_t,
    std::arch::aarch64::float64x2_t,
) {
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        use std::arch::aarch64::{vcombine_f64, vdup_n_f64, vdupq_n_f64};
        (
            vdupq_n_f64(w.re),
            vcombine_f64(vdup_n_f64(-w.im), vdup_n_f64(w.im)),
        )
    }
}

#[cfg(target_arch = "aarch64")]
#[target_feature(enable = "neon")]
#[allow(clippy::too_many_arguments)]
unsafe fn radix4_butterfly_neon(
    base_ptr: *mut f64,
    inner_stride: usize,
    w_2s_k: Complex64,
    w_2s_kps: Complex64,
    w_s_k: Complex64,
) {
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        use std::arch::aarch64::{
            vaddq_f64, vdupq_n_f64, vld1q_f64, vmulq_f64, vst1q_f64, vsubq_f64,
        };
        let s = inner_stride;
        let inv_sqrt2_v = vdupq_n_f64(std::f64::consts::FRAC_1_SQRT_2);

        let a_ptr = base_ptr;
        let b_ptr = base_ptr.add(s * 2);
        let c_ptr = base_ptr.add(2 * s * 2);
        let d_ptr = base_ptr.add(3 * s * 2);

        let a = vld1q_f64(a_ptr);
        let b = vld1q_f64(b_ptr);
        let c = vld1q_f64(c_ptr);
        let d = vld1q_f64(d_ptr);

        let (w_2s_k_rr, w_2s_k_ii_as) = neon_twiddle_parts(w_2s_k);
        let (w_2s_kps_rr, w_2s_kps_ii_as) = neon_twiddle_parts(w_2s_kps);
        let (w_s_k_rr, w_s_k_ii_as) = neon_twiddle_parts(w_s_k);

        let t_ac_sum = vmulq_f64(vaddq_f64(a, c), inv_sqrt2_v);
        let t_ac_diff = simd::complex_mul_neon(w_2s_k_rr, w_2s_k_ii_as, vsubq_f64(a, c));
        let t_bd_sum = vmulq_f64(vaddq_f64(b, d), inv_sqrt2_v);
        let t_bd_diff = simd::complex_mul_neon(w_2s_kps_rr, w_2s_kps_ii_as, vsubq_f64(b, d));

        vst1q_f64(a_ptr, vmulq_f64(vaddq_f64(t_ac_sum, t_bd_sum), inv_sqrt2_v));
        vst1q_f64(
            b_ptr,
            simd::complex_mul_neon(w_s_k_rr, w_s_k_ii_as, vsubq_f64(t_ac_sum, t_bd_sum)),
        );
        vst1q_f64(
            c_ptr,
            vmulq_f64(vaddq_f64(t_ac_diff, t_bd_diff), inv_sqrt2_v),
        );
        vst1q_f64(
            d_ptr,
            simd::complex_mul_neon(w_s_k_rr, w_s_k_ii_as, vsubq_f64(t_ac_diff, t_bd_diff)),
        );
    }
}

#[cfg(target_arch = "aarch64")]
#[target_feature(enable = "neon")]
#[cfg_attr(not(feature = "parallel"), allow(dead_code))]
#[allow(clippy::too_many_arguments)]
unsafe fn radix4_butterfly_slices_neon(
    a_ptr: *mut f64,
    b_ptr: *mut f64,
    c_ptr: *mut f64,
    d_ptr: *mut f64,
    w_2s_k: Complex64,
    w_2s_kps: Complex64,
    w_s_k: Complex64,
) {
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        use std::arch::aarch64::{
            vaddq_f64, vdupq_n_f64, vld1q_f64, vmulq_f64, vst1q_f64, vsubq_f64,
        };
        let inv_sqrt2_v = vdupq_n_f64(std::f64::consts::FRAC_1_SQRT_2);

        let a = vld1q_f64(a_ptr);
        let b = vld1q_f64(b_ptr);
        let c = vld1q_f64(c_ptr);
        let d = vld1q_f64(d_ptr);

        let (w_2s_k_rr, w_2s_k_ii_as) = neon_twiddle_parts(w_2s_k);
        let (w_2s_kps_rr, w_2s_kps_ii_as) = neon_twiddle_parts(w_2s_kps);
        let (w_s_k_rr, w_s_k_ii_as) = neon_twiddle_parts(w_s_k);

        let t_ac_sum = vmulq_f64(vaddq_f64(a, c), inv_sqrt2_v);
        let t_ac_diff = simd::complex_mul_neon(w_2s_k_rr, w_2s_k_ii_as, vsubq_f64(a, c));
        let t_bd_sum = vmulq_f64(vaddq_f64(b, d), inv_sqrt2_v);
        let t_bd_diff = simd::complex_mul_neon(w_2s_kps_rr, w_2s_kps_ii_as, vsubq_f64(b, d));

        vst1q_f64(a_ptr, vmulq_f64(vaddq_f64(t_ac_sum, t_bd_sum), inv_sqrt2_v));
        vst1q_f64(
            b_ptr,
            simd::complex_mul_neon(w_s_k_rr, w_s_k_ii_as, vsubq_f64(t_ac_sum, t_bd_sum)),
        );
        vst1q_f64(
            c_ptr,
            vmulq_f64(vaddq_f64(t_ac_diff, t_bd_diff), inv_sqrt2_v),
        );
        vst1q_f64(
            d_ptr,
            simd::complex_mul_neon(w_s_k_rr, w_s_k_ii_as, vsubq_f64(t_ac_diff, t_bd_diff)),
        );
    }
}

#[inline]
fn apply_radix4_groups(
    group: &mut [Complex64],
    s: usize,
    k_base: usize,
    step_outer: usize,
    step_inner: usize,
    twiddles_scaled: &[Complex64],
    inv_sqrt2: f64,
) {
    debug_assert_eq!(group.len() % (s * 4), 0);
    let groups = group.len() / (s * 4);
    #[cfg(target_arch = "aarch64")]
    let _ = inv_sqrt2;

    #[cfg(target_arch = "x86_64")]
    {
        if simd::has_avx2_fma() && s >= 2 {
            // SAFETY: has_avx2_fma() is runtime checked. Adjacent pairs stay
            // inside the same radix-4 group because s >= 2. The safe
            // alternative is the scalar radix-4 loop below; this SIMD path is
            // kept because the measured QFT speedup depends on avoiding scalar
            // Complex64 arithmetic in the fused stage-pair hot loop.
            unsafe {
                let base_ptr = group.as_mut_ptr() as *mut f64;
                for g in 0..groups {
                    let group_off = g * s * 4;
                    let mut k_offset = 0usize;
                    while k_offset + 2 <= s {
                        let k = k_base + k_offset;
                        let w_2s_k0 = twiddles_scaled[k * step_outer];
                        let w_2s_k1 = twiddles_scaled[(k + 1) * step_outer];
                        let w_2s_kps0 = twiddles_scaled[(k + s) * step_outer];
                        let w_2s_kps1 = twiddles_scaled[(k + s + 1) * step_outer];
                        let w_s_k0 = twiddles_scaled[k * step_inner];
                        let w_s_k1 = twiddles_scaled[(k + 1) * step_inner];
                        let pair_base = base_ptr.add((group_off + k_offset) * 2);
                        radix4_butterfly_pair_avx2fma(
                            pair_base, s, w_2s_k0, w_2s_k1, w_2s_kps0, w_2s_kps1, w_s_k0, w_s_k1,
                        );
                        k_offset += 2;
                    }
                    if k_offset < s {
                        let k = k_base + k_offset;
                        let w_2s_k = twiddles_scaled[k * step_outer];
                        let w_2s_kps = twiddles_scaled[(k + s) * step_outer];
                        let w_s_k = twiddles_scaled[k * step_inner];
                        let pair_base = base_ptr.add((group_off + k_offset) * 2);
                        radix4_butterfly_fma(pair_base, s, w_2s_k, w_2s_kps, w_s_k);
                    }
                }
            }
            return;
        }
        if simd::has_fma() {
            // SAFETY: has_fma() is runtime checked. The safe alternative is
            // the scalar radix-4 loop below; this SIMD path is kept because the
            // measured QFT speedup depends on avoiding scalar Complex64
            // arithmetic in the fused stage-pair hot loop.
            unsafe {
                let base_ptr = group.as_mut_ptr() as *mut f64;
                for g in 0..groups {
                    let group_off = g * s * 4;
                    for k_offset in 0..s {
                        let k = k_base + k_offset;
                        let w_2s_k = twiddles_scaled[k * step_outer];
                        let w_2s_kps = twiddles_scaled[(k + s) * step_outer];
                        let w_s_k = twiddles_scaled[k * step_inner];
                        let pair_base = base_ptr.add((group_off + k_offset) * 2);
                        radix4_butterfly_fma(pair_base, s, w_2s_k, w_2s_kps, w_s_k);
                    }
                }
            }
            return;
        }
    }

    #[cfg(target_arch = "aarch64")]
    {
        // SAFETY: NEON is baseline on aarch64. The loop maps each radix-4
        // butterfly to four in-bounds amplitudes inside `group`; the safe
        // alternative is the scalar fallback below. This SIMD counterpart is
        // required so Apple Silicon does not pay scalar Complex64 overhead in
        // the same measured QFT hot loop as x86.
        unsafe {
            let base_ptr = group.as_mut_ptr() as *mut f64;
            for g in 0..groups {
                let group_off = g * s * 4;
                for k_offset in 0..s {
                    let k = k_base + k_offset;
                    let w_2s_k = twiddles_scaled[k * step_outer];
                    let w_2s_kps = twiddles_scaled[(k + s) * step_outer];
                    let w_s_k = twiddles_scaled[k * step_inner];
                    let pair_base = base_ptr.add((group_off + k_offset) * 2);
                    radix4_butterfly_neon(pair_base, s, w_2s_k, w_2s_kps, w_s_k);
                }
            }
        }
    }

    #[cfg(not(target_arch = "aarch64"))]
    {
        for g in 0..groups {
            let base = g * s * 4;
            for k_offset in 0..s {
                let k = k_base + k_offset;
                let w_2s_k = twiddles_scaled[k * step_outer];
                let w_2s_kps = twiddles_scaled[(k + s) * step_outer];
                let w_s_k = twiddles_scaled[k * step_inner];
                radix4_butterfly_scalar(
                    group,
                    base + k_offset,
                    s,
                    w_2s_k,
                    w_2s_kps,
                    w_s_k,
                    inv_sqrt2,
                );
            }
        }
    }
}

/// Run two DIF stages at inner stride `s` within `slice`, reading [`qft_compact_twiddles`].
#[inline]
fn fft_stage_pair_in_slice(
    slice: &mut [Complex64],
    s: usize,
    compact: &[Complex64],
    inv_sqrt2: f64,
) {
    apply_radix4_groups(slice, s, 0, 1, 2, &compact[s << 1..s << 2], inv_sqrt2);
}

/// Offsets per quarter that one task of [`fft_stage_pair_chunked`] covers.
#[cfg(feature = "parallel")]
const QFT_PAIR_CHUNK: usize = 1024;
/// Radix-4 groups that share one gathered twiddle chunk within a task.
#[cfg(feature = "parallel")]
const QFT_PAIR_GROUPS_PER_TASK: usize = 16;

/// Run a radix-4 stage pair as tasks of one twiddle chunk over a batch of groups.
///
/// Below the top pair a stride's twiddles sit `step_outer` and `step_inner` entries
/// apart in the whole-state table, up to a page each. A task gathers its chunk once
/// and reuses it across the batch, and splitting by chunk keeps every thread busy when
/// a stage has fewer groups than threads. `s` must be a multiple of the chunk.
#[cfg(feature = "parallel")]
fn fft_stage_pair_chunked(
    state: &mut [Complex64],
    s: usize,
    step_outer: usize,
    step_inner: usize,
    twiddles_scaled: &[Complex64],
    inv_sqrt2: f64,
) {
    const CHUNK: usize = QFT_PAIR_CHUNK;
    debug_assert_eq!(s % CHUNK, 0);
    let chunks = s / CHUNK;
    let groups = state.len() / (s << 2);
    let batch = QFT_PAIR_GROUPS_PER_TASK.min(groups);
    let batches = groups / batch;
    let ptr = SendPtr(state.as_mut_ptr());
    (0..chunks * batches).into_par_iter().for_each(move |task| {
        let k0 = task / batches * CHUNK;
        let first_group = task % batches * batch;
        // SAFETY: a task owns offsets `[k0, k0 + CHUNK)` of the four quarters of each
        // group in its batch. Tasks differ in chunk or batch, so no two share a run.
        unsafe {
            fft_stage_pair_chunk(
                ptr.as_complex_ptr(),
                s,
                k0,
                first_group..first_group + batch,
                step_outer,
                step_inner,
                twiddles_scaled,
                inv_sqrt2,
            );
        }
    });
}

/// Gather the twiddles for offsets `[k0, k0 + QFT_PAIR_CHUNK)` and apply them to
/// `groups`. Out of line, so its 48 KiB of twiddles live only while a task runs
/// instead of in every frame of Rayon's recursive split.
///
/// # Safety
///
/// `state` must be valid for every group in `groups`, and no other thread may access
/// those offsets of any quarter of these groups during the call.
#[cfg(feature = "parallel")]
#[inline(never)]
#[allow(clippy::too_many_arguments)]
unsafe fn fft_stage_pair_chunk(
    state: *mut Complex64,
    s: usize,
    k0: usize,
    groups: std::ops::Range<usize>,
    step_outer: usize,
    step_inner: usize,
    twiddles_scaled: &[Complex64],
    inv_sqrt2: f64,
) {
    const CHUNK: usize = QFT_PAIR_CHUNK;
    let w_outer: [Complex64; CHUNK] =
        std::array::from_fn(|j| twiddles_scaled[(k0 + j) * step_outer]);
    let w_outer_s: [Complex64; CHUNK] =
        std::array::from_fn(|j| twiddles_scaled[(k0 + j + s) * step_outer]);
    let w_inner: [Complex64; CHUNK] =
        std::array::from_fn(|j| twiddles_scaled[(k0 + j) * step_inner]);
    for g in groups {
        let base = g * (s << 2) + k0;
        // SAFETY: same contract as the enclosing unsafe fn.
        let (qa, qb, qc, qd) = unsafe {
            (
                std::slice::from_raw_parts_mut(state.add(base), CHUNK),
                std::slice::from_raw_parts_mut(state.add(base + s), CHUNK),
                std::slice::from_raw_parts_mut(state.add(base + 2 * s), CHUNK),
                std::slice::from_raw_parts_mut(state.add(base + 3 * s), CHUNK),
            )
        };
        radix4_runs(qa, qb, qc, qd, &w_outer, &w_outer_s, &w_inner, inv_sqrt2);
    }
}

/// Apply radix-4 butterflies at matching offsets of four quarter slices, with offset
/// `j` taking index `j` of each twiddle slice.
#[cfg(feature = "parallel")]
#[inline(always)]
#[allow(clippy::too_many_arguments)]
fn radix4_runs(
    qa: &mut [Complex64],
    qb: &mut [Complex64],
    qc: &mut [Complex64],
    qd: &mut [Complex64],
    w_outer: &[Complex64],
    w_outer_s: &[Complex64],
    w_inner: &[Complex64],
    inv_sqrt2: f64,
) {
    let n = qa.len();
    #[cfg(target_arch = "aarch64")]
    let _ = inv_sqrt2;

    #[cfg(target_arch = "x86_64")]
    if simd::has_avx2_fma() && n >= 2 {
        let mut j = 0usize;
        // SAFETY: AVX2 detected. Each call reads and writes two adjacent amplitudes
        // at offset `j` of four disjoint slices of length `n`, with `j + 2 <= n`.
        unsafe {
            let (pa, pb, pc, pd) = (
                qa.as_mut_ptr() as *mut f64,
                qb.as_mut_ptr() as *mut f64,
                qc.as_mut_ptr() as *mut f64,
                qd.as_mut_ptr() as *mut f64,
            );
            while j + 2 <= n {
                radix4_butterfly_pair_slices_avx2fma(
                    pa.add(j * 2),
                    pb.add(j * 2),
                    pc.add(j * 2),
                    pd.add(j * 2),
                    w_outer[j],
                    w_outer[j + 1],
                    w_outer_s[j],
                    w_outer_s[j + 1],
                    w_inner[j],
                    w_inner[j + 1],
                );
                j += 2;
            }
        }
        if j < n {
            radix4_butterfly_quartet_scalar(
                qa,
                qb,
                qc,
                qd,
                j,
                w_outer[j],
                w_outer_s[j],
                w_inner[j],
                inv_sqrt2,
            );
        }
        return;
    }

    #[cfg(target_arch = "aarch64")]
    // SAFETY: NEON is baseline on aarch64. Each call reads and writes the amplitude at
    // offset `j < n` of four disjoint slices of length `n`.
    unsafe {
        let (pa, pb, pc, pd) = (
            qa.as_mut_ptr() as *mut f64,
            qb.as_mut_ptr() as *mut f64,
            qc.as_mut_ptr() as *mut f64,
            qd.as_mut_ptr() as *mut f64,
        );
        for j in 0..n {
            radix4_butterfly_slices_neon(
                pa.add(j * 2),
                pb.add(j * 2),
                pc.add(j * 2),
                pd.add(j * 2),
                w_outer[j],
                w_outer_s[j],
                w_inner[j],
            );
        }
    }

    #[cfg(not(target_arch = "aarch64"))]
    for j in 0..n {
        radix4_butterfly_quartet_scalar(
            qa,
            qb,
            qc,
            qd,
            j,
            w_outer[j],
            w_outer_s[j],
            w_inner[j],
            inv_sqrt2,
        );
    }
}

/// Run two DIF FFT stages as one radix-4 pass over every `2 * twiddles_scaled.len()`
/// amplitude chunk of `state`.
///
/// `inner_stride = 2^(outer_stage - 1)`. Falls back to a sequential loop when
/// the `parallel` feature is off.
#[inline]
fn fft_stage_pair_par(
    state: &mut [Complex64],
    inner_stride: usize,
    twiddles_scaled: &[Complex64],
    inv_sqrt2: f64,
) {
    let s = inner_stride;
    let fft_len = 2 * twiddles_scaled.len();
    let step_outer = fft_len / (s << 2);
    let step_inner = fft_len / (s << 1);

    #[cfg(feature = "parallel")]
    fft_stage_pair_chunked(state, s, step_outer, step_inner, twiddles_scaled, inv_sqrt2);

    #[cfg(not(feature = "parallel"))]
    apply_radix4_groups(
        state,
        s,
        0,
        step_outer,
        step_inner,
        twiddles_scaled,
        inv_sqrt2,
    );
}

/// Run the DIF stages above `tile_bits` of every `2^n` chunk as gathered passes of at
/// most [`qft_gather_max_bits`] stages each, split evenly.
fn qft_gather_high_stages(state: &mut [Complex64], n: usize, tile_bits: usize, inverse: bool) {
    let passes = (n - tile_bits).div_ceil(qft_gather_max_bits());
    let mut top = n;
    for left in (1..=passes).rev() {
        let bits = (top - tile_bits).div_ceil(left);
        qft_gather_pass(state, QftGatherPass::new(top - bits, bits, inverse));
        top -= bits;
    }
}

fn qft_gather_pass(state: &mut [Complex64], pass: QftGatherPass) {
    let groups = state.len() >> (pass.bits + pass.row_len.ilog2() as usize);
    let len = pass.row_len << pass.bits;

    #[cfg(feature = "parallel")]
    if state.len() >= 1usize << PARALLEL_THRESHOLD_QUBITS {
        let ptr = SendPtr(state.as_mut_ptr());
        (0..groups).into_par_iter().for_each_init(
            || vec![Complex64::new(0.0, 0.0); len],
            move |buf, g| {
                // SAFETY: a group owns `row_len` offsets of each of its rows, and the
                // groups tile the state, so no two tasks share an amplitude.
                unsafe { pass.run_group(ptr.as_complex_ptr(), g, buf) }
            },
        );
        return;
    }

    let mut buf = vec![Complex64::new(0.0, 0.0); len];
    let ptr = state.as_mut_ptr();
    for g in 0..groups {
        // SAFETY: `ptr` spans the whole state, borrowed mutably for this loop alone.
        unsafe { pass.run_group(ptr, g, &mut buf) }
    }
}

/// One gathered pass of the QFT's DIF network: stages `lo..lo + bits` of every
/// `2^(lo + bits)` chunk, a group of `2^bits` rows at a time.
///
/// A group copies, for `row_len` consecutive offsets `L` below `lo`, the amplitudes at
/// every value of the pass bits into a row each. The stage at local stride `2^u`
/// twiddles position `k 2^lo + L` by `exp(2 pi i (k 2^lo + L) / 2^(lo + u + 1))`, the
/// local DIF twiddle for `k` times a factor in `L` alone. Both inputs of every later
/// butterfly carry the same factor, so the rows run a plain `2^bits`-point DIF on
/// [`qft_compact_twiddles`] and the factors collect into one phase
/// `exp(2 pi i L rev(r) / 2^(lo + bits))` on output row `r`, applied on the way back.
#[derive(Clone, Copy)]
struct QftGatherPass {
    lo: usize,
    bits: usize,
    row_len: usize,
    compact: &'static [Complex64],
    phases: &'static QftGatherPhases,
}

impl QftGatherPass {
    fn new(lo: usize, bits: usize, inverse: bool) -> Self {
        Self {
            lo,
            bits,
            row_len: (qft_gather_len() >> bits).min(1 << lo),
            compact: qft_compact_twiddles(inverse),
            phases: qft_gather_phases(lo + bits, inverse),
        }
    }

    /// Gather group `g` into `buf`, run its stages, and scatter it back phased.
    ///
    /// # Safety
    ///
    /// `state` must be valid for every amplitude of group `g`, and no other thread may
    /// access them during the call.
    #[inline]
    unsafe fn run_group(&self, state: *mut Complex64, g: usize, buf: &mut [Complex64]) {
        let blocks_log2 = self.lo - self.row_len.ilog2() as usize;
        let col0 = (g & ((1 << blocks_log2) - 1)) * self.row_len;
        let base = ((g >> blocks_log2) << (self.lo + self.bits)) + col0;
        for (r, row) in buf.chunks_exact_mut(self.row_len).enumerate() {
            // SAFETY: same contract as the enclosing unsafe fn.
            unsafe {
                std::ptr::copy_nonoverlapping(
                    state.add(base + (r << self.lo)),
                    row.as_mut_ptr(),
                    self.row_len,
                );
            }
        }
        qft_dif_rows(buf, self.row_len, self.bits, self.compact);
        for (r, row) in buf.chunks_exact(self.row_len).enumerate() {
            let rev = reverse_low_bits(r, self.bits);
            // SAFETY: same contract as the enclosing unsafe fn.
            let out = unsafe {
                std::slice::from_raw_parts_mut(state.add(base + (r << self.lo)), self.row_len)
            };
            qft_scatter_row(
                out,
                row,
                self.phases.phase(col0 * rev),
                &self.phases.low,
                rev,
            );
        }
    }
}

/// Run a `2^bits`-point DIF down the rows of `buf`, every column of `row_len` an
/// independent transform, in the stage order of the in-tile stages.
fn qft_dif_rows(buf: &mut [Complex64], row_len: usize, bits: usize, compact: &[Complex64]) {
    let mut top = bits;
    if bits % 2 == 1 {
        top -= 1;
        let s = 1usize << top;
        for block in buf.chunks_exact_mut(2 * s * row_len) {
            let (a, b) = block.split_at_mut(s * row_len);
            for (k, (ra, rb)) in a
                .chunks_exact_mut(row_len)
                .zip(b.chunks_exact_mut(row_len))
                .enumerate()
            {
                radix2_rows(ra, rb, compact[s + k]);
            }
        }
    }
    while top > 0 {
        let s = 1usize << (top - 2);
        let quarter = s * row_len;
        for block in buf.chunks_exact_mut(4 * quarter) {
            let (a, rest) = block.split_at_mut(quarter);
            let (b, rest) = rest.split_at_mut(quarter);
            let (c, d) = rest.split_at_mut(quarter);
            let rows = a
                .chunks_exact_mut(row_len)
                .zip(b.chunks_exact_mut(row_len))
                .zip(c.chunks_exact_mut(row_len))
                .zip(d.chunks_exact_mut(row_len));
            for (k, (((ra, rb), rc), rd)) in rows.enumerate() {
                radix4_rows(
                    ra,
                    rb,
                    rc,
                    rd,
                    compact[2 * s + k],
                    compact[3 * s + k],
                    compact[s + k],
                );
            }
        }
        top -= 2;
    }
}

/// Radix-2 DIF butterflies between two rows of one even length sharing twiddle `w`.
#[inline(always)]
fn radix2_rows(lo: &mut [Complex64], hi: &mut [Complex64], w: Complex64) {
    debug_assert_eq!(lo.len(), hi.len());
    #[cfg(target_arch = "x86_64")]
    if simd::has_avx2_fma() {
        // SAFETY: AVX2 detected. The rows share one even length.
        unsafe { radix2_rows_avx2fma(lo, hi, w) };
        return;
    }
    let inv_sqrt2 = std::f64::consts::FRAC_1_SQRT_2;
    for (a, b) in lo.iter_mut().zip(hi.iter_mut()) {
        (*a, *b) = radix2_butterfly_values(*a, *b, w, inv_sqrt2);
    }
}

#[cfg(target_arch = "x86_64")]
#[target_feature(enable = "avx2,fma")]
unsafe fn radix2_rows_avx2fma(lo: &mut [Complex64], hi: &mut [Complex64], w: Complex64) {
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        use std::arch::x86_64::{
            _mm256_add_pd, _mm256_loadu_pd, _mm256_mul_pd, _mm256_set1_pd, _mm256_storeu_pd,
            _mm256_sub_pd,
        };
        let inv_sqrt2_v = _mm256_set1_pd(std::f64::consts::FRAC_1_SQRT_2);
        let (w_rr, w_ii) = (_mm256_set1_pd(w.re), _mm256_set1_pd(w.im));
        let (pl, ph) = (lo.as_mut_ptr() as *mut f64, hi.as_mut_ptr() as *mut f64);
        for off in (0..2 * lo.len()).step_by(4) {
            let a = _mm256_loadu_pd(pl.add(off));
            let b = _mm256_loadu_pd(ph.add(off));
            _mm256_storeu_pd(pl.add(off), _mm256_mul_pd(_mm256_add_pd(a, b), inv_sqrt2_v));
            _mm256_storeu_pd(
                ph.add(off),
                simd::complex_mul_avx2fma(w_rr, w_ii, _mm256_sub_pd(a, b)),
            );
        }
    }
}

/// Radix-4 DIF butterflies across four rows of one even length sharing one twiddle
/// triple, in the argument order of [`radix4_butterfly_values`].
#[inline(always)]
fn radix4_rows(
    qa: &mut [Complex64],
    qb: &mut [Complex64],
    qc: &mut [Complex64],
    qd: &mut [Complex64],
    w_outer: Complex64,
    w_outer_s: Complex64,
    w_inner: Complex64,
) {
    let n = qa.len();
    debug_assert!(qb.len() == n && qc.len() == n && qd.len() == n);

    #[cfg(target_arch = "x86_64")]
    if simd::has_avx2_fma() {
        // SAFETY: AVX2 detected. The four rows share one even length.
        unsafe { radix4_rows_avx2fma(qa, qb, qc, qd, w_outer, w_outer_s, w_inner) };
        return;
    }

    #[cfg(target_arch = "aarch64")]
    // SAFETY: NEON is baseline on aarch64. Each call reads and writes the amplitude at
    // offset `j < n` of four disjoint rows of length `n`.
    unsafe {
        let (pa, pb, pc, pd) = (
            qa.as_mut_ptr() as *mut f64,
            qb.as_mut_ptr() as *mut f64,
            qc.as_mut_ptr() as *mut f64,
            qd.as_mut_ptr() as *mut f64,
        );
        for j in 0..n {
            radix4_butterfly_slices_neon(
                pa.add(j * 2),
                pb.add(j * 2),
                pc.add(j * 2),
                pd.add(j * 2),
                w_outer,
                w_outer_s,
                w_inner,
            );
        }
    }

    #[cfg(not(target_arch = "aarch64"))]
    for j in 0..n {
        radix4_butterfly_quartet_scalar(
            qa,
            qb,
            qc,
            qd,
            j,
            w_outer,
            w_outer_s,
            w_inner,
            std::f64::consts::FRAC_1_SQRT_2,
        );
    }
}

#[cfg(target_arch = "x86_64")]
#[target_feature(enable = "avx2,fma")]
unsafe fn radix4_rows_avx2fma(
    qa: &mut [Complex64],
    qb: &mut [Complex64],
    qc: &mut [Complex64],
    qd: &mut [Complex64],
    w_outer: Complex64,
    w_outer_s: Complex64,
    w_inner: Complex64,
) {
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        use std::arch::x86_64::{
            _mm256_add_pd, _mm256_loadu_pd, _mm256_mul_pd, _mm256_set1_pd, _mm256_storeu_pd,
            _mm256_sub_pd,
        };
        let inv_sqrt2_v = _mm256_set1_pd(std::f64::consts::FRAC_1_SQRT_2);
        let (wo_rr, wo_ii) = (_mm256_set1_pd(w_outer.re), _mm256_set1_pd(w_outer.im));
        let (ws_rr, ws_ii) = (_mm256_set1_pd(w_outer_s.re), _mm256_set1_pd(w_outer_s.im));
        let (wi_rr, wi_ii) = (_mm256_set1_pd(w_inner.re), _mm256_set1_pd(w_inner.im));
        let (pa, pb, pc, pd) = (
            qa.as_mut_ptr() as *mut f64,
            qb.as_mut_ptr() as *mut f64,
            qc.as_mut_ptr() as *mut f64,
            qd.as_mut_ptr() as *mut f64,
        );
        for off in (0..2 * qa.len()).step_by(4) {
            let a = _mm256_loadu_pd(pa.add(off));
            let b = _mm256_loadu_pd(pb.add(off));
            let c = _mm256_loadu_pd(pc.add(off));
            let d = _mm256_loadu_pd(pd.add(off));

            let t_ac_sum = _mm256_mul_pd(_mm256_add_pd(a, c), inv_sqrt2_v);
            let t_ac_diff = simd::complex_mul_avx2fma(wo_rr, wo_ii, _mm256_sub_pd(a, c));
            let t_bd_sum = _mm256_mul_pd(_mm256_add_pd(b, d), inv_sqrt2_v);
            let t_bd_diff = simd::complex_mul_avx2fma(ws_rr, ws_ii, _mm256_sub_pd(b, d));

            _mm256_storeu_pd(
                pa.add(off),
                _mm256_mul_pd(_mm256_add_pd(t_ac_sum, t_bd_sum), inv_sqrt2_v),
            );
            _mm256_storeu_pd(
                pb.add(off),
                simd::complex_mul_avx2fma(wi_rr, wi_ii, _mm256_sub_pd(t_ac_sum, t_bd_sum)),
            );
            _mm256_storeu_pd(
                pc.add(off),
                _mm256_mul_pd(_mm256_add_pd(t_ac_diff, t_bd_diff), inv_sqrt2_v),
            );
            _mm256_storeu_pd(
                pd.add(off),
                simd::complex_mul_avx2fma(wi_rr, wi_ii, _mm256_sub_pd(t_ac_diff, t_bd_diff)),
            );
        }
    }
}

/// Write `row` into `out` with amplitude `j` times `anchor * low[j * step]`.
#[inline(always)]
fn qft_scatter_row(
    out: &mut [Complex64],
    row: &[Complex64],
    anchor: Complex64,
    low: &[Complex64],
    step: usize,
) {
    debug_assert_eq!(out.len(), row.len());
    #[cfg(target_arch = "x86_64")]
    if simd::has_avx2_fma() {
        // SAFETY: AVX2 detected. `out` and `row` share one even length.
        unsafe { qft_scatter_row_avx2fma(out, row, anchor, low, step) };
        return;
    }
    for (j, (o, &a)) in out.iter_mut().zip(row).enumerate() {
        *o = a * (anchor * low[j * step]);
    }
}

#[cfg(target_arch = "x86_64")]
#[target_feature(enable = "avx2,fma")]
unsafe fn qft_scatter_row_avx2fma(
    out: &mut [Complex64],
    row: &[Complex64],
    anchor: Complex64,
    low: &[Complex64],
    step: usize,
) {
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        use std::arch::x86_64::{
            _mm256_loadu_pd, _mm256_movedup_pd, _mm256_permute_pd, _mm256_set_pd, _mm256_set1_pd,
            _mm256_storeu_pd,
        };
        let (a_rr, a_ii) = (_mm256_set1_pd(anchor.re), _mm256_set1_pd(anchor.im));
        let (src, dst) = (row.as_ptr() as *const f64, out.as_mut_ptr() as *mut f64);
        let mut x = 0;
        for off in (0..2 * row.len()).step_by(4) {
            let (p0, p1) = (low[x], low[x + step]);
            let w =
                simd::complex_mul_avx2fma(a_rr, a_ii, _mm256_set_pd(p1.im, p1.re, p0.im, p0.re));
            let v = _mm256_loadu_pd(src.add(off));
            _mm256_storeu_pd(
                dst.add(off),
                simd::complex_mul_avx2fma(_mm256_movedup_pd(w), _mm256_permute_pd(w, 0b1111), v),
            );
            x += 2 * step;
        }
    }
}

/// Shortest run of contiguous pairs handed to the SIMD pair kernel. Below it the
/// AVX2 loop runs at most one 256-bit iteration per call, so the per-run dispatch
/// outweighs it and the scalar loop takes the block.
const MIN_PAULI_ROT_SIMD_RUN: usize = 4;

/// The pair arithmetic of `exp(-i θ P / 2)` for [`pauli_rot_pair_halves`].
///
/// `m_lo` and `m_hi` are the parity-even cross coefficients toward the pivot-clear
/// and pivot-set side; `even` and `odd` hold the pair matrix
/// `[[c, ±m_lo], [±m_hi, c]]` for each parity with its SIMD tier chosen once.
struct PauliRotPair {
    c: f64,
    m_lo: Complex64,
    m_hi: Complex64,
    even: simd::PreparedGate1q,
    odd: simd::PreparedGate1q,
}

impl PauliRotPair {
    /// The pair arithmetic for a rotation by `theta` whose string has `num_y` Y letters.
    ///
    /// `m_hi` points toward the pivot-set side of a pair; the opposite direction
    /// differs by the parity of `xmask & zmask`, which is the Y count.
    fn new(theta: f64, num_y: u32) -> Self {
        let (s, c) = (theta / 2.0).sin_cos();
        let m_hi = match num_y % 4 {
            0 => Complex64::new(0.0, -s),
            1 => Complex64::new(s, 0.0),
            2 => Complex64::new(0.0, s),
            _ => Complex64::new(-s, 0.0),
        };
        let m_lo = if num_y % 2 == 1 { -m_hi } else { m_hi };
        let diag = Complex64::new(c, 0.0);
        Self {
            c,
            m_lo,
            m_hi,
            even: simd::PreparedGate1q::new(&[[diag, m_lo], [m_hi, diag]]),
            odd: simd::PreparedGate1q::new(&[[diag, -m_lo], [-m_hi, diag]]),
        }
    }
}

/// Mix the pairs `(lo[i], hi[i ^ xlow])` of one Pauli-rotation pivot block.
///
/// `base` is the global index of `lo[0]` and must be a multiple of `lo.len()`,
/// so `(base + i) & zmask` reads the pair's parity sign directly; `flip` (0 or 1)
/// is added to that parity. Pairs whose index bits below the lowest X or Z letter
/// vary form contiguous runs on both sides with one parity, and runs of at least
/// [`MIN_PAULI_ROT_SIMD_RUN`] go through the SIMD pair kernel.
#[inline(always)]
fn pauli_rot_pair_halves(
    lo: &mut [Complex64],
    hi: &mut [Complex64],
    base: usize,
    xlow: usize,
    zmask: usize,
    pair: &PauliRotPair,
    flip: u32,
) {
    let run = 1usize << (((xlow | zmask) & (lo.len() - 1)) | lo.len()).trailing_zeros();
    if run >= MIN_PAULI_ROT_SIMD_RUN {
        for start in (0..lo.len()).step_by(run) {
            let k = start ^ xlow;
            let prepared = if ((base + start) & zmask).count_ones() & 1 == flip {
                &pair.even
            } else {
                &pair.odd
            };
            prepared.apply(&mut lo[start..start + run], &mut hi[k..k + run]);
        }
        return;
    }

    let (c, m_lo, m_hi) = (pair.c, pair.m_lo, pair.m_hi);
    for (i, amp_lo) in lo.iter_mut().enumerate() {
        let k = i ^ xlow;
        let a = *amp_lo;
        let b = hi[k];
        if ((base + i) & zmask).count_ones() & 1 == flip {
            *amp_lo = a * c + b * m_lo;
            hi[k] = b * c + a * m_hi;
        } else {
            *amp_lo = a * c - b * m_lo;
            hi[k] = b * c - a * m_hi;
        }
    }
}

/// How one `MultiPauliRot` rotation acts inside a tile.
#[allow(clippy::large_enum_variant)]
enum TiledAction {
    /// A diagonal string: the phase for an even and for an odd Z parity.
    Parity([Complex64; 2]),
    Pairs(PauliRotPair),
}

/// One rotation of a `MultiPauliRot` batch in tile coordinates.
///
/// `xmask` and `zmask` cover tile bits. `zrest` keeps the Z letters on qubits
/// outside the tile: their parity is fixed across a tile and, when odd, negates
/// both cross terms or swaps the two phases.
struct TiledPauliRot {
    xmask: usize,
    zmask: usize,
    zrest: usize,
    action: TiledAction,
}

impl TiledPauliRot {
    /// Move a rotation's global masks into the tile `tile_bit` describes: the tile
    /// bit a qubit lands on, or `None` for a qubit outside the tile, which must
    /// carry no X or Y letter.
    fn new(
        xmask: usize,
        zmask: usize,
        theta: f64,
        tile_bit: impl Fn(usize) -> Option<usize>,
    ) -> Self {
        let (mut tile_x, mut tile_z, mut zrest) = (0usize, 0usize, 0usize);
        let mut letters = xmask | zmask;
        while letters != 0 {
            let q = letters.trailing_zeros() as usize;
            letters &= letters - 1;
            match tile_bit(q) {
                Some(bit) => {
                    tile_x |= (xmask >> q & 1) << bit;
                    tile_z |= (zmask >> q & 1) << bit;
                }
                None => {
                    debug_assert_eq!(xmask >> q & 1, 0, "X or Y letter outside the tile");
                    zrest |= 1 << q;
                }
            }
        }
        let action = if xmask == 0 {
            let (s, c) = (theta / 2.0).sin_cos();
            TiledAction::Parity([Complex64::new(c, -s), Complex64::new(c, s)])
        } else {
            TiledAction::Pairs(PauliRotPair::new(theta, (xmask & zmask).count_ones()))
        };
        Self {
            xmask: tile_x,
            zmask: tile_z,
            zrest,
            action,
        }
    }

    /// Apply to `tile`, whose first amplitude has global index `base` once the
    /// tile bits are cleared.
    #[inline(always)]
    fn apply(&self, tile: &mut [Complex64], base: usize) {
        let flip = (base & self.zrest).count_ones() & 1;
        match &self.action {
            TiledAction::Parity([even, odd]) => {
                let phases = if flip == 1 {
                    [*odd, *even]
                } else {
                    [*even, *odd]
                };
                for (i, amp) in tile.iter_mut().enumerate() {
                    *amp *= phases[((i & self.zmask).count_ones() & 1) as usize];
                }
            }
            TiledAction::Pairs(pair) => {
                let pivot = usize::BITS as usize - 1 - self.xmask.leading_zeros() as usize;
                let half = 1usize << pivot;
                let block = half << 1;
                let xlow = self.xmask ^ half;
                for (k, chunk) in tile.chunks_mut(block).enumerate() {
                    let (lo, hi) = chunk.split_at_mut(half);
                    pauli_rot_pair_halves(lo, hi, k * block, xlow, self.zmask, pair, flip);
                }
            }
        }
    }
}

/// Apply one gate over the whole state, tiling or splitting by block size.
#[cfg(feature = "parallel")]
#[inline(always)]
pub(crate) fn apply_single_gate_par(
    state: &mut [Complex64],
    target: usize,
    mat: &[[Complex64; 2]; 2],
) {
    let half = 1usize << target;
    let block_size = half << 1;
    let prepared = simd::PreparedGate1q::new(mat);

    const MIN_TILE: usize = 8192;

    let tile_size = MIN_TILE.max(block_size);
    let num_tiles = state.len() / tile_size;

    if block_size <= MIN_TILE && num_tiles >= 4 {
        state.par_chunks_mut(MIN_TILE).for_each(|tile| {
            prepared.apply_full_sequential(tile, target);
        });
    } else if num_tiles >= 4 {
        state.par_chunks_mut(block_size).for_each(|chunk| {
            let (lo, hi) = chunk.split_at_mut(half);
            prepared.apply(lo, hi);
        });
    } else {
        let sub_tile = MIN_TILE.min(half);
        for block in state.chunks_mut(block_size) {
            let (lo, hi) = block.split_at_mut(half);
            lo.par_chunks_mut(sub_tile)
                .zip(hi.par_chunks_mut(sub_tile))
                .for_each(|(lo_t, hi_t)| {
                    prepared.apply(lo_t, hi_t);
                });
        }
    }
}

/// Apply gates on `k` distinct targets in one traversal instead of `k`.
///
/// Gates on disjoint qubits commute, so the stages run in target order.
/// For targets above the L2 tier only: below it the tiled pass keeps a tile
/// cache resident across gates, which this shape cannot. The L3 tier's targets
/// join it where their tiles would not stay resident either.
#[cfg(feature = "parallel")]
fn apply_multi_1q_shared(state: &mut [Complex64], gates: &[(usize, [[Complex64; 2]; 2])]) {
    let mut sorted: SmallVec<[(usize, [[Complex64; 2]; 2]); 4]> = gates.into();
    sorted.sort_unstable_by_key(|&(target, _)| target);
    let k = sorted.len();
    let lanes = 1usize << k;
    let low_target = sorted[0].0;

    let prepared: SmallVec<[simd::PreparedGate1q; 4]> = sorted
        .iter()
        .map(|(_, mat)| simd::PreparedGate1q::new(mat))
        .collect();

    let mut offsets: SmallVec<[usize; 16]> = SmallVec::with_capacity(lanes);
    for j in 0..lanes {
        let mut offset = 0usize;
        for (axis, &(target, _)) in sorted.iter().enumerate() {
            if (j >> axis) & 1 == 1 {
                offset |= 1usize << target;
            }
        }
        offsets.push(offset);
    }

    let run = 1usize << low_target;
    let tile_len = (multi_gate_l2_tile() / lanes).max(1);
    let tiles_per_run = run / tile_len;
    let blocks = (state.len() >> k) >> low_target;

    let base_of = |block: usize| {
        let mut base = block << low_target;
        for &(target, _) in sorted.iter() {
            base = insert_zero_bit(base, target);
        }
        base
    };

    let start_of = |unit: usize| base_of(unit / tiles_per_run) + (unit % tiles_per_run) * tile_len;
    let units = blocks * tiles_per_run;

    let ptr = SendPtr(state.as_mut_ptr());
    (0..units).into_par_iter().for_each(|unit| {
        let start = start_of(unit);
        for (axis, prep) in prepared.iter().enumerate() {
            let stride = 1usize << axis;
            for j in 0..lanes {
                if j & stride != 0 {
                    continue;
                }
                // SAFETY: the target weights are super-increasing, so
                // the target bits and tile index partition the buffer
                // into disjoint `tile_len` runs; no two lanes or units
                // alias.
                unsafe {
                    let base = ptr.as_complex_ptr();
                    let lo = std::slice::from_raw_parts_mut(base.add(start + offsets[j]), tile_len);
                    let hi = std::slice::from_raw_parts_mut(
                        base.add(start + offsets[j | stride]),
                        tile_len,
                    );
                    prep.apply(lo, hi);
                }
            }
        }
    });
}

/// Gather offsets of a dense unitary's `2^k` amplitudes from a base index:
/// matrix index `m` reads bit `k - 1 - i` for `targets[i]`, so `targets[0]`
/// is the most significant matrix-index bit.
#[inline(always)]
fn unitary_offsets(targets: &[usize]) -> [usize; 1 << MAX_UNITARY_QUBITS] {
    let k = targets.len();
    let mut offsets = [0usize; 1 << MAX_UNITARY_QUBITS];
    for (m, offset) in offsets.iter_mut().enumerate().take(1 << k) {
        for (i, &target) in targets.iter().enumerate() {
            *offset |= ((m >> (k - 1 - i)) & 1) << target;
        }
    }
    offsets
}

/// Targets sorted ascending, the order `insert_zero_bit` needs to place a base
/// index around every target bit.
#[inline(always)]
fn sorted_unitary_targets(targets: &[usize]) -> [usize; MAX_UNITARY_QUBITS] {
    let mut sorted = [0usize; MAX_UNITARY_QUBITS];
    sorted[..targets.len()].copy_from_slice(targets);
    sorted[..targets.len()].sort_unstable();
    sorted
}

/// `state[base + offsets[r]] = sum_c mat[r * dim + c] * state[base + offsets[c]]`
/// over the `dim` amplitudes one base index selects, `amps` the gather buffer.
///
/// # Safety
/// `base + off` must index the buffer behind `ptr` for every entry of
/// `offsets`, and nothing else may touch those amplitudes meanwhile.
#[inline(always)]
unsafe fn apply_unitary_at(
    ptr: *mut Complex64,
    base: usize,
    offsets: &[usize],
    mat: &[Complex64],
    amps: &mut [Complex64],
) {
    let dim = offsets.len();
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        for (slot, &off) in amps.iter_mut().zip(offsets) {
            *slot = *ptr.add(base + off);
        }
        for (row, &off) in offsets.iter().enumerate() {
            let coeffs = mat.get_unchecked(row * dim..(row + 1) * dim);
            let mut acc = Complex64::new(0.0, 0.0);
            for (c, a) in coeffs.iter().zip(amps.iter()) {
                acc += c * a;
            }
            *ptr.add(base + off) = acc;
        }
    }
}

/// Apply a dense `2^k x 2^k` matrix over `targets` in one gather-scatter pass,
/// `targets[0]` the most significant matrix-index bit. Shared with the
/// factored backend, whose blocks are bare statevector slices.
#[inline(always)]
pub(crate) fn apply_unitary_seq(state: &mut [Complex64], targets: &[usize], mat: &[Complex64]) {
    let k = targets.len();
    let dim = 1usize << k;
    let sorted = sorted_unitary_targets(targets);
    let offsets = unitary_offsets(targets);
    let num_bases = state.len() >> k;
    let ptr = state.as_mut_ptr();
    let mut amps = [Complex64::new(0.0, 0.0); 1 << MAX_UNITARY_QUBITS];
    for i in 0..num_bases {
        let mut base = i;
        for &q in &sorted[..k] {
            base = insert_zero_bit(base, q);
        }
        // SAFETY: insert_zero_bit over the ascending targets maps each base
        // index onto an in-bounds index with every target bit clear, so the
        // `dim` offsets select in-bounds amplitudes disjoint across bases.
        unsafe {
            apply_unitary_at(ptr, base, &offsets[..dim], mat, &mut amps[..dim]);
        }
    }
}

/// [`apply_unitary_seq`] with the base indices split across Rayon tasks.
#[cfg(feature = "parallel")]
#[inline(always)]
pub(crate) fn par_apply_unitary(state: &mut [Complex64], targets: &[usize], mat: &[Complex64]) {
    let k = targets.len();
    let dim = 1usize << k;
    let sorted = sorted_unitary_targets(targets);
    let offsets = unitary_offsets(targets);
    let num_bases = state.len() >> k;
    let ptr = SendPtr(state.as_mut_ptr());
    (0..num_bases)
        .into_par_iter()
        .with_min_len(MIN_PAR_ITERS)
        .for_each(move |i| {
            let mut amps = [Complex64::new(0.0, 0.0); 1 << MAX_UNITARY_QUBITS];
            let mut base = i;
            for &q in &sorted[..k] {
                base = insert_zero_bit(base, q);
            }
            // SAFETY: insert_zero_bit over the ascending targets maps each
            // base index onto an in-bounds index with every target bit clear,
            // so the `dim` offsets select in-bounds amplitudes disjoint across
            // bases, and no two Rayon tasks share an amplitude.
            unsafe {
                apply_unitary_at(
                    ptr.as_complex_ptr(),
                    base,
                    &offsets[..dim],
                    mat,
                    &mut amps[..dim],
                );
            }
        });
}

type TargetedGate1q = (usize, [[Complex64; 2]; 2]);

/// Sort `gates` into the L2 tile, L3 tile and untiled tiers by target. Both upper
/// tiers run in ascending target order, so folding the L3 tier into the untiled one
/// leaves every amplitude bitwise unchanged.
#[inline(always)]
fn split_multi_1q_tiers(
    gates: &[TargetedGate1q],
    max_l2_target: usize,
    medium_tiled: bool,
    small_gates: &mut SmallVec<[(usize, simd::PreparedGate1q); 16]>,
    medium_gates: &mut SmallVec<[(usize, simd::PreparedGate1q); 4]>,
    large_gates: &mut SmallVec<[TargetedGate1q; 4]>,
) {
    for &(target, mat) in gates {
        if target <= max_l2_target {
            small_gates.push((target, simd::PreparedGate1q::new(&mat)));
        } else if target <= MULTI_GATE_MAX_L3_TARGET && medium_tiled {
            medium_gates.push((target, simd::PreparedGate1q::new(&mat)));
        } else {
            large_gates.push((target, mat));
        }
    }
    medium_gates.sort_unstable_by_key(|&(target, _)| target);
    large_gates.sort_unstable_by_key(|&(target, _)| target);
}

/// Whether the L3 tier beats folding its targets into the shared passes: every
/// worker needs its own tile, and the concurrent tiles must fit the last-level cache.
/// An unreported cache keeps the tier.
#[cfg(feature = "parallel")]
fn keep_l3_tier(state_len: usize) -> bool {
    let Some(cache) = crate::backend::cache::last_level_cache() else {
        return true;
    };
    let tiles = state_len / MULTI_GATE_L3_TILE;
    let workers = rayon::current_num_threads();
    let tile_bytes = MULTI_GATE_L3_TILE * std::mem::size_of::<Complex64>();
    tiles >= workers && crate::backend::cache::tiles_fit(cache, tile_bytes, workers)
}

/// Sort diagonal `gates` into the L2 tile, L3 tile and untiled tiers by target, as
/// `(target, d0, d1)` with the tiled tiers also carrying `is_phase_one(d0)`.
#[inline(always)]
fn split_multi_1q_diagonal_tiers(
    gates: &[TargetedGate1q],
    max_l2_target: usize,
    small_gates: &mut SmallVec<[(usize, Complex64, Complex64, bool); 16]>,
    medium_gates: &mut SmallVec<[(usize, Complex64, Complex64, bool); 4]>,
    large_gates: &mut SmallVec<[(usize, Complex64, Complex64); 4]>,
) {
    for &(target, mat) in gates {
        let d0 = mat[0][0];
        let d1 = mat[1][1];
        let skip_lo = is_phase_one(d0);
        if target <= max_l2_target {
            small_gates.push((target, d0, d1, skip_lo));
        } else if target <= MULTI_GATE_MAX_L3_TARGET {
            medium_gates.push((target, d0, d1, skip_lo));
        } else {
            large_gates.push((target, d0, d1));
        }
    }
}

/// Passes over a `num_qubits` state that [`StatevectorBackend::apply_multi_1q`] makes
/// for a non-diagonal `gates` batch, on the route [`multi_1q_route`] picks. Callers
/// weigh it against the passes of another route.
pub(crate) fn multi_1q_sweeps(gates: &[(usize, [[Complex64; 2]; 2])], num_qubits: usize) -> usize {
    #[cfg(feature = "parallel")]
    if num_qubits >= PARALLEL_THRESHOLD_QUBITS {
        let state_len = 1usize << num_qubits;
        return multi_1q_route(
            gates,
            num_qubits,
            max_target_for_tile(multi_gate_l2_tile_par(state_len)),
            keep_l3_tier(state_len),
            MAX_SHARED_1Q_TARGETS,
        )
        .1;
    }
    multi_1q_route(gates, num_qubits, multi_gate_max_l2_target(), true, 1).1
}

/// The passes [`split_multi_1q_tiers`] leads to with these tier bounds, when the gates
/// above the tiles share traversals `shared` at a time.
fn tier_sweeps(
    gates: &[(usize, [[Complex64; 2]; 2])],
    max_l2_target: usize,
    medium_tiled: bool,
    shared: usize,
) -> usize {
    let (mut small, mut medium, mut large) = (false, false, 0usize);
    for &(target, _) in gates {
        if target <= max_l2_target {
            small = true;
        } else if target <= MULTI_GATE_MAX_L3_TARGET && medium_tiled {
            medium = true;
        } else {
            large += 1;
        }
    }
    small as usize + medium as usize + large.div_ceil(shared)
}

/// How a `MultiFused` batch runs. Every route applies the targets inside the
/// in-place tile first, in list order, then the rest by ascending target, so the
/// routes agree bitwise.
enum Multi1qRoute {
    /// The tile tiers, then shared traversals of the targets above them.
    Tiered,
    /// One subcube gather over every target.
    Whole(SubcubePlan),
    /// The in-place tile pass, then the targets above it cut into runs of the
    /// `Multi2q` high budget, each `(len, plan)`: one subcube gather per run with a
    /// plan, shared traversals for a run whose gather would save no pass.
    Split(SmallVec<[(usize, Option<SubcubePlan>); 2]>),
}

/// The route for `gates` on a `num_qubits` state and its pass count. A gather route
/// is taken only when it makes strictly fewer passes than the tiers, whose targets
/// above `max_l2_target` share traversals `shared` at a time.
fn multi_1q_route(
    gates: &[TargetedGate1q],
    num_qubits: usize,
    max_l2_target: usize,
    medium_tiled: bool,
    shared: usize,
) -> (Multi1qRoute, usize) {
    let tiered = tier_sweeps(gates, max_l2_target, medium_tiled, shared);
    if tiered <= 1 {
        return (Multi1qRoute::Tiered, tiered);
    }
    if let Some(plan) = qubit_subcube_plan(gates.iter().map(|&(target, _)| target), num_qubits) {
        return (Multi1qRoute::Whole(plan), 1);
    }
    let mut above: SmallVec<[usize; 16]> = gates
        .iter()
        .map(|&(target, _)| target)
        .filter(|&target| target > max_l2_target)
        .collect();
    above.sort_unstable();
    let mut passes = usize::from(above.len() < gates.len());
    let mut runs = SmallVec::new();
    for run in above.chunks(multi_2q_high_budget_for(num_qubits)) {
        let plan = if run.len() > shared {
            qubit_subcube_plan(run.iter().copied(), num_qubits)
        } else {
            None
        };
        passes += if plan.is_some() {
            1
        } else {
            run.len().div_ceil(shared)
        };
        runs.push((run.len(), plan));
    }
    if passes < tiered {
        (Multi1qRoute::Split(runs), passes)
    } else {
        (Multi1qRoute::Tiered, tiered)
    }
}

/// The gates of a [`Multi1qRoute::Whole`] gather keyed by their bit in `plan`'s tile:
/// `small` in its list order, then `above` by ascending target.
fn whole_subcube_gates(
    plan: &SubcubePlan,
    small: SmallVec<[(usize, simd::PreparedGate1q); 16]>,
    above: &[TargetedGate1q],
) -> SmallVec<[(usize, simd::PreparedGate1q); 16]> {
    small
        .into_iter()
        .map(|(target, prepared)| (plan.map[target], prepared))
        .chain(
            above
                .iter()
                .map(|&(target, ref mat)| (plan.map[target], simd::PreparedGate1q::new(mat))),
        )
        .collect()
}

/// Apply `above`, the gates past the in-place tile by ascending target, in the runs
/// of a [`Multi1qRoute::Split`].
fn apply_multi_1q_runs(
    state: &mut [Complex64],
    above: &[TargetedGate1q],
    runs: &[(usize, Option<SubcubePlan>)],
) {
    let mut rest = above;
    for &(len, ref plan) in runs {
        let (run, tail) = rest.split_at(len);
        rest = tail;
        if let Some(plan) = plan {
            let tiled: SmallVec<[(usize, simd::PreparedGate1q); 16]> = run
                .iter()
                .map(|&(target, ref mat)| (plan.map[target], simd::PreparedGate1q::new(mat)))
                .collect();
            apply_multi_1q_subcubes(state, plan, &tiled);
            continue;
        }
        #[cfg(feature = "parallel")]
        if state.len() >= 1 << PARALLEL_THRESHOLD_QUBITS {
            for chunk in run.chunks(MAX_SHARED_1Q_TARGETS) {
                if chunk.len() > 1 {
                    apply_multi_1q_shared(state, chunk);
                } else {
                    apply_single_gate_par(state, chunk[0].0, &chunk[0].1);
                }
            }
            continue;
        }
        for &(target, ref mat) in run {
            simd::PreparedGate1q::new(mat).apply_full_sequential(state, target);
        }
    }
}

/// Apply a `MultiFused` batch on the route [`multi_1q_route`] picks. Shared with the
/// factored backend, whose blocks are bare statevector slices.
#[cfg(feature = "parallel")]
#[inline(always)]
pub(crate) fn apply_multi_1q_par(state: &mut [Complex64], gates: &[(usize, [[Complex64; 2]; 2])]) {
    if gates.is_empty() {
        return;
    }
    if gates.len() == 1 {
        apply_single_gate_par(state, gates[0].0, &gates[0].1);
        return;
    }
    let medium_tiled = keep_l3_tier(state.len());
    let (route, _) = multi_1q_route(
        gates,
        state.len().trailing_zeros() as usize,
        max_target_for_tile(multi_gate_l2_tile_par(state.len())),
        medium_tiled,
        MAX_SHARED_1Q_TARGETS,
    );
    apply_multi_1q_par_routed(state, gates, medium_tiled, &route);
}

/// [`apply_multi_1q_par`] on `route`, with the L3 tier on or folded into the shared
/// passes when the route is [`Multi1qRoute::Tiered`].
#[cfg(feature = "parallel")]
#[inline(always)]
fn apply_multi_1q_par_routed(
    state: &mut [Complex64],
    gates: &[(usize, [[Complex64; 2]; 2])],
    medium_tiled: bool,
    route: &Multi1qRoute,
) {
    let tile_size = multi_gate_l2_tile_par(state.len());
    let mut small_gates = SmallVec::new();
    let mut medium_gates = SmallVec::new();
    let mut large_gates = SmallVec::new();
    split_multi_1q_tiers(
        gates,
        max_target_for_tile(tile_size),
        medium_tiled && matches!(route, Multi1qRoute::Tiered),
        &mut small_gates,
        &mut medium_gates,
        &mut large_gates,
    );

    if let Multi1qRoute::Whole(plan) = route {
        let tiled = whole_subcube_gates(plan, small_gates, &large_gates);
        apply_multi_1q_subcubes(state, plan, &tiled);
        return;
    }

    if !small_gates.is_empty() {
        state
            .par_chunks_mut(tile_size)
            .with_min_len(chunk_min_len(tile_size))
            .for_each(|tile| {
                for &(target, ref prepared) in &small_gates {
                    prepared.apply_tiled(tile, target);
                }
            });
    }

    if !medium_gates.is_empty() {
        let outer_block = 1usize << (MULTI_GATE_MAX_L3_TARGET + 1);
        let tile_size = MULTI_GATE_L3_TILE.max(outer_block);
        state
            .par_chunks_mut(tile_size)
            .with_min_len(chunk_min_len(tile_size))
            .for_each(|tile| {
                for &(target, ref prepared) in &medium_gates {
                    prepared.apply_tiled(tile, target);
                }
            });
    }

    if let Multi1qRoute::Split(runs) = route {
        apply_multi_1q_runs(state, &large_gates, runs);
        return;
    }

    for chunk in large_gates.chunks(MAX_SHARED_1Q_TARGETS) {
        if chunk.len() > 1 {
            apply_multi_1q_shared(state, chunk);
        } else {
            apply_single_gate_par(state, chunk[0].0, &chunk[0].1);
        }
    }
}

impl StatevectorBackend {
    #[inline(always)]
    pub(super) fn apply_single_gate(&mut self, target: usize, mat: [[Complex64; 2]; 2]) {
        #[cfg(feature = "parallel")]
        if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
            apply_single_gate_par(&mut self.state, target, &mat);
            return;
        }

        let prepared = simd::PreparedGate1q::new(&mat);
        prepared.apply_full_sequential(&mut self.state, target);
    }

    #[inline(always)]
    pub(super) fn apply_diagonal_gate(&mut self, target: usize, d0: Complex64, d1: Complex64) {
        #[cfg(feature = "parallel")]
        if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
            self.apply_diagonal_gate_par(target, d0, d1);
            return;
        }

        let skip_lo = is_phase_one(d0);
        simd::apply_diagonal_sequential(&mut self.state, target, d0, d1, skip_lo);
    }

    #[cfg(feature = "parallel")]
    #[inline(always)]
    fn apply_diagonal_gate_par(&mut self, target: usize, d0: Complex64, d1: Complex64) {
        let skip_lo = is_phase_one(d0);

        let half = 1usize << target;
        let block_size = half << 1;
        let num_blocks = self.state.len() / block_size;

        if num_blocks < 4 {
            apply_diagonal_high_target_par(&mut self.state, target, d0, d1, skip_lo);
            return;
        }

        const MIN_TILE: usize = 8192;
        let tile_size = MIN_TILE.max(block_size);
        self.state.par_chunks_mut(tile_size).for_each(|tile| {
            simd::apply_diagonal_sequential(tile, target, d0, d1, skip_lo);
        });
    }

    #[inline(always)]
    pub(super) fn apply_cx(&mut self, control: usize, target: usize) {
        #[cfg(feature = "parallel")]
        if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
            self.apply_cx_par(control, target);
            return;
        }

        if control > target {
            let ctrl_half = 1usize << control;
            let block_size = ctrl_half << 1;
            let tgt_half = 1usize << target;
            let tgt_block = tgt_half << 1;

            for chunk in self.state.chunks_mut(block_size) {
                let (_, hi) = chunk.split_at_mut(ctrl_half);
                for sub in hi.chunks_mut(tgt_block) {
                    let (sub_lo, sub_hi) = sub.split_at_mut(tgt_half);
                    swap_slices_kernel(sub_lo, sub_hi);
                }
            }
        } else {
            let ctrl_half = 1usize << control;
            let ctrl_block = ctrl_half << 1;
            let tgt_half = 1usize << target;
            let block_size = tgt_half << 1;

            for chunk in self.state.chunks_mut(block_size) {
                let (lo, hi) = chunk.split_at_mut(tgt_half);
                for (lo_sub, hi_sub) in lo.chunks_mut(ctrl_block).zip(hi.chunks_mut(ctrl_block)) {
                    let (_, lo_active) = lo_sub.split_at_mut(ctrl_half);
                    let (_, hi_active) = hi_sub.split_at_mut(ctrl_half);
                    swap_slices_kernel(lo_active, hi_active);
                }
            }
        }
    }

    #[cfg(feature = "parallel")]
    #[inline(always)]
    fn apply_cx_par(&mut self, control: usize, target: usize) {
        if control > target {
            let ctrl_half = 1usize << control;
            let block_size = ctrl_half << 1;
            let tgt_half = 1usize << target;
            let tgt_block = tgt_half << 1;
            let num_blocks = self.state.len() / block_size;

            if num_blocks >= 4 {
                self.state
                    .par_chunks_mut(block_size)
                    .with_min_len(chunk_min_len(block_size))
                    .for_each(|chunk| {
                        let (_, hi) = chunk.split_at_mut(ctrl_half);
                        for sub in hi.chunks_mut(tgt_block) {
                            let (sub_lo, sub_hi) = sub.split_at_mut(tgt_half);
                            swap_slices_kernel(sub_lo, sub_hi);
                        }
                    });
            } else {
                let inner_tile = MIN_PAR_ELEMS.max(tgt_block);
                for chunk in self.state.chunks_mut(block_size) {
                    let (_, hi) = chunk.split_at_mut(ctrl_half);
                    hi.par_chunks_mut(inner_tile).for_each(|tile| {
                        for sub in tile.chunks_mut(tgt_block) {
                            let (sub_lo, sub_hi) = sub.split_at_mut(tgt_half);
                            swap_slices_kernel(sub_lo, sub_hi);
                        }
                    });
                }
            }
        } else {
            let ctrl_half = 1usize << control;
            let ctrl_block = ctrl_half << 1;
            let tgt_half = 1usize << target;
            let block_size = tgt_half << 1;
            let num_blocks = self.state.len() / block_size;

            if num_blocks >= 4 {
                self.state
                    .par_chunks_mut(block_size)
                    .with_min_len(chunk_min_len(block_size))
                    .for_each(|chunk| {
                        let (lo, hi) = chunk.split_at_mut(tgt_half);
                        for (lo_sub, hi_sub) in
                            lo.chunks_mut(ctrl_block).zip(hi.chunks_mut(ctrl_block))
                        {
                            let (_, lo_active) = lo_sub.split_at_mut(ctrl_half);
                            let (_, hi_active) = hi_sub.split_at_mut(ctrl_half);
                            swap_slices_kernel(lo_active, hi_active);
                        }
                    });
            } else {
                let inner_tile = MIN_PAR_ELEMS.max(ctrl_block);
                for chunk in self.state.chunks_mut(block_size) {
                    let (lo, hi) = chunk.split_at_mut(tgt_half);
                    lo.par_chunks_mut(inner_tile)
                        .zip(hi.par_chunks_mut(inner_tile))
                        .for_each(|(lo_tile, hi_tile)| {
                            for (lo_sub, hi_sub) in lo_tile
                                .chunks_mut(ctrl_block)
                                .zip(hi_tile.chunks_mut(ctrl_block))
                            {
                                let (_, lo_active) = lo_sub.split_at_mut(ctrl_half);
                                let (_, hi_active) = hi_sub.split_at_mut(ctrl_half);
                                swap_slices_kernel(lo_active, hi_active);
                            }
                        });
                }
            }
        }
    }

    #[inline(always)]
    pub(super) fn apply_cz(&mut self, q0: usize, q1: usize) {
        #[cfg(feature = "parallel")]
        if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
            self.apply_cz_par(q0, q1);
            return;
        }

        let (lo_q, hi_q) = if q0 < q1 { (q0, q1) } else { (q1, q0) };
        let lo_half = 1usize << lo_q;
        let lo_block = lo_half << 1;
        let hi_half = 1usize << hi_q;
        let block_size = hi_half << 1;

        for chunk in self.state.chunks_mut(block_size) {
            let (_, hi_group) = chunk.split_at_mut(hi_half);
            for sub in hi_group.chunks_mut(lo_block) {
                let (_, sub_hi) = sub.split_at_mut(lo_half);
                negate_slice_kernel(sub_hi);
            }
        }
    }

    #[cfg(feature = "parallel")]
    #[inline(always)]
    fn apply_cz_par(&mut self, q0: usize, q1: usize) {
        let (lo_q, hi_q) = if q0 < q1 { (q0, q1) } else { (q1, q0) };
        let lo_half = 1usize << lo_q;
        let lo_block = lo_half << 1;
        let hi_half = 1usize << hi_q;
        let block_size = hi_half << 1;
        let num_blocks = self.state.len() / block_size;

        if num_blocks >= 4 {
            self.state
                .par_chunks_mut(block_size)
                .with_min_len(chunk_min_len(block_size))
                .for_each(|chunk| {
                    let (_, hi_group) = chunk.split_at_mut(hi_half);
                    for sub in hi_group.chunks_mut(lo_block) {
                        let (_, sub_hi) = sub.split_at_mut(lo_half);
                        negate_slice_kernel(sub_hi);
                    }
                });
        } else {
            let inner_tile = MIN_PAR_ELEMS.max(lo_block);
            for chunk in self.state.chunks_mut(block_size) {
                let (_, hi_group) = chunk.split_at_mut(hi_half);
                hi_group.par_chunks_mut(inner_tile).for_each(|tile| {
                    for sub in tile.chunks_mut(lo_block) {
                        let (_, sub_hi) = sub.split_at_mut(lo_half);
                        negate_slice_kernel(sub_hi);
                    }
                });
            }
        }
    }

    #[inline(always)]
    pub(super) fn apply_rzz(&mut self, q0: usize, q1: usize, theta: f64) {
        #[cfg(feature = "parallel")]
        if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
            self.apply_rzz_par(q0, q1, theta);
            return;
        }

        let phase_same = Complex64::from_polar(1.0, -theta / 2.0);
        let phase_diff = Complex64::from_polar(1.0, theta / 2.0);
        let phases = [phase_same, phase_diff];

        for (i, amp) in self.state.iter_mut().enumerate() {
            let parity = ((i >> q0) ^ (i >> q1)) & 1;
            *amp *= phases[parity];
        }
    }

    #[cfg(feature = "parallel")]
    #[inline(always)]
    fn apply_rzz_par(&mut self, q0: usize, q1: usize, theta: f64) {
        let phase_same = Complex64::from_polar(1.0, -theta / 2.0);
        let phase_diff = Complex64::from_polar(1.0, theta / 2.0);
        let phases = [phase_same, phase_diff];

        self.state
            .par_chunks_mut(MIN_PAR_ELEMS)
            .enumerate()
            .for_each(|(chunk_idx, chunk)| {
                let base = chunk_idx * MIN_PAR_ELEMS;
                for (j, amp) in chunk.iter_mut().enumerate() {
                    let i = base + j;
                    let parity = ((i >> q0) ^ (i >> q1)) & 1;
                    *amp *= phases[parity];
                }
            });
    }

    /// Apply `exp(-i θ P / 2)` in one pass, `P` given as one letter per target.
    ///
    /// `P|j> = i^{num_y} (-1)^{popcount(j & zmask)} |j ^ xmask>`, so the gate
    /// mixes each pair `(j, j ^ xmask)` with `cos(θ/2)` on the diagonal and a
    /// sign-resolved `-i sin(θ/2) i^{num_y}` cross term. A Z-only string has
    /// `xmask == 0` and reduces to parity phases on the diagonal.
    #[inline(always)]
    pub(super) fn apply_pauli_rot(&mut self, targets: &[usize], theta: f64, axes: &[PauliAxis]) {
        let (xmask, zmask, num_y) = pauli_rot_masks(targets, axes);
        self.apply_pauli_rot_masks(xmask, zmask, num_y, theta);
    }

    #[inline(always)]
    fn apply_pauli_rot_masks(&mut self, xmask: usize, zmask: usize, num_y: u32, theta: f64) {
        if xmask == 0 {
            let (s, c) = (theta / 2.0).sin_cos();
            self.apply_parity_phase(zmask, Complex64::new(c, -s), Complex64::new(c, s));
            return;
        }

        let pair = PauliRotPair::new(theta, num_y);
        let pivot = usize::BITS as usize - 1 - xmask.leading_zeros() as usize;
        let half = 1usize << pivot;
        let block = half << 1;
        let xlow = xmask ^ half;

        #[cfg(feature = "parallel")]
        if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
            self.apply_pauli_rot_pairs_par(pivot, xlow, zmask, &pair);
            return;
        }

        for (chunk_idx, chunk) in self.state.chunks_mut(block).enumerate() {
            let (lo, hi) = chunk.split_at_mut(half);
            pauli_rot_pair_halves(lo, hi, chunk_idx * block, xlow, zmask, &pair, 0);
        }
    }

    #[cfg(feature = "parallel")]
    #[inline(always)]
    fn apply_pauli_rot_pairs_par(
        &mut self,
        pivot: usize,
        xlow: usize,
        zmask: usize,
        pair: &PauliRotPair,
    ) {
        let half = 1usize << pivot;
        let block = half << 1;
        let num_blocks = self.state.len() / block;

        if num_blocks >= 4 {
            self.state
                .par_chunks_mut(block)
                .with_min_len(chunk_min_len(block))
                .enumerate()
                .for_each(|(chunk_idx, chunk)| {
                    let (lo, hi) = chunk.split_at_mut(half);
                    pauli_rot_pair_halves(lo, hi, chunk_idx * block, xlow, zmask, pair, 0);
                });
        } else {
            // Tiles are a power of two above `xlow`, so `i ^ xlow` stays
            // inside its tile and each (lo, hi) tile pair is disjoint.
            let tile = MIN_PAR_ELEMS.max((xlow + 1).next_power_of_two());
            for (chunk_idx, chunk) in self.state.chunks_mut(block).enumerate() {
                let chunk_base = chunk_idx * block;
                let (lo, hi) = chunk.split_at_mut(half);
                lo.par_chunks_mut(tile)
                    .zip(hi.par_chunks_mut(tile))
                    .enumerate()
                    .for_each(|(t, (lo_t, hi_t))| {
                        pauli_rot_pair_halves(
                            lo_t,
                            hi_t,
                            chunk_base + t * tile,
                            xlow,
                            zmask,
                            pair,
                            0,
                        );
                    });
            }
        }
    }

    /// Multiply each amplitude by `even` or `odd` from the parity of its
    /// index bits under `zmask`.
    #[inline(always)]
    fn apply_parity_phase(&mut self, zmask: usize, even: Complex64, odd: Complex64) {
        #[cfg(feature = "parallel")]
        if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
            self.apply_parity_phase_par(zmask, even, odd);
            return;
        }

        let phases = [even, odd];
        for (i, amp) in self.state.iter_mut().enumerate() {
            *amp *= phases[((i & zmask).count_ones() & 1) as usize];
        }
    }

    #[cfg(feature = "parallel")]
    #[inline(always)]
    fn apply_parity_phase_par(&mut self, zmask: usize, even: Complex64, odd: Complex64) {
        let phases = [even, odd];
        self.state
            .par_chunks_mut(MIN_PAR_ELEMS)
            .enumerate()
            .for_each(|(chunk_idx, chunk)| {
                let base = chunk_idx * MIN_PAR_ELEMS;
                for (j, amp) in chunk.iter_mut().enumerate() {
                    *amp *= phases[(((base + j) & zmask).count_ones() & 1) as usize];
                }
            });
    }

    /// Apply `(xmask, zmask, theta)` Pauli rotations in order in one cache-tiled pass.
    ///
    /// Strings whose X and Y letters all sit below the tile bits run on the state in
    /// place; others gather the subcube holding their X and Y qubits, as
    /// [`apply_multi_2q`](Self::apply_multi_2q) does. Z letters never need the tile.
    /// A batch no tile holds takes one pass per rotation.
    pub(super) fn apply_multi_pauli_rot(&mut self, rotations: &[(usize, usize, f64)]) {
        let xall = rotations.iter().fold(0usize, |acc, &(x, _, _)| acc | x);
        let xqubits = (0..self.num_qubits).filter(|&q| xall >> q & 1 == 1);
        if let Some(plan) = qubit_subcube_plan(xqubits, self.num_qubits) {
            let tiled: Vec<TiledPauliRot> = rotations
                .iter()
                .map(|&(x, z, theta)| {
                    TiledPauliRot::new(x, z, theta, |q| {
                        Some(plan.map[q]).filter(|&bit| bit != usize::MAX)
                    })
                })
                .collect();
            let rest = plan.rest(self.num_qubits);
            let subcubes = 1usize << (self.num_qubits - plan.tile_bits);
            #[cfg(feature = "parallel")]
            if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
                let ptr = SendPtr(self.state.as_mut_ptr());
                (0..subcubes).into_par_iter().for_each(|outer| {
                    // SAFETY: `ptr` holds 2^num_qubits amplitudes and distinct
                    // `outer` values address disjoint subcubes, so no two tasks
                    // touch the same amplitude.
                    unsafe {
                        apply_pauli_rot_subcube(ptr.as_complex_ptr(), &plan, &rest, outer, &tiled)
                    };
                });
                return;
            }
            let state = self.state.as_mut_ptr();
            for outer in 0..subcubes {
                // SAFETY: `state` holds 2^num_qubits amplitudes and each `outer`
                // addresses its own subcube, applied one at a time here.
                unsafe { apply_pauli_rot_subcube(state, &plan, &rest, outer, &tiled) };
            }
            return;
        }

        let tile_bits = multi_2q_tile_bits_for(self.num_qubits).min(self.num_qubits);
        if xall >> tile_bits != 0 {
            for &(x, z, theta) in rotations {
                self.apply_pauli_rot_masks(x, z, (x & z).count_ones(), theta);
            }
            return;
        }
        let tiled: Vec<TiledPauliRot> = rotations
            .iter()
            .map(|&(x, z, theta)| TiledPauliRot::new(x, z, theta, |q| (q < tile_bits).then_some(q)))
            .collect();
        let tile_len = 1usize << tile_bits;

        #[cfg(feature = "parallel")]
        if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
            self.state
                .par_chunks_mut(tile_len)
                .with_min_len(chunk_min_len(tile_len))
                .enumerate()
                .for_each(|(k, tile)| {
                    for rotation in &tiled {
                        rotation.apply(tile, k * tile_len);
                    }
                });
            return;
        }

        for (k, tile) in self.state.chunks_mut(tile_len).enumerate() {
            for rotation in &tiled {
                rotation.apply(tile, k * tile_len);
            }
        }
    }

    #[inline(always)]
    pub(super) fn apply_batch_rzz(&mut self, edges: &[(usize, usize, f64)]) {
        #[cfg(feature = "parallel")]
        if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
            self.apply_batch_rzz_par(edges);
            return;
        }

        #[cfg(target_arch = "x86_64")]
        if simd::has_bmi2() && simd::has_fma() {
            let bmi2_groups = self.batch_tables.rzz_bmi2();
            if let Some(num_groups) = build_batch_rzz_bmi2_tables(edges, bmi2_groups) {
                // SAFETY: has_bmi2() && has_fma() verified above
                unsafe { batch_rzz_tile_bmi2(&mut self.state, 0, bmi2_groups, num_groups) };
                return;
            }
        }

        let groups = self.batch_tables.rzz();
        let num_groups = build_batch_rzz_tables(edges, groups);

        for (i, amp) in self.state.iter_mut().enumerate() {
            let mut combined = ONE;
            for group in groups.iter().take(num_groups) {
                let bits = extract_rzz_bits(i, group);
                combined *= group.table[bits];
            }
            *amp *= combined;
        }
    }

    #[cfg(feature = "parallel")]
    #[inline(always)]
    fn apply_batch_rzz_par(&mut self, edges: &[(usize, usize, f64)]) {
        #[cfg(target_arch = "x86_64")]
        if simd::has_bmi2() && simd::has_fma() {
            let bmi2_groups = self.batch_tables.rzz_bmi2();
            if let Some(num_groups) = build_batch_rzz_bmi2_tables(edges, bmi2_groups) {
                let bmi2_groups = &*bmi2_groups;
                self.state
                    .par_chunks_mut(MIN_PAR_ELEMS)
                    .enumerate()
                    .for_each(|(chunk_idx, chunk)| {
                        let base = chunk_idx * MIN_PAR_ELEMS;
                        // SAFETY: has_bmi2() && has_fma() verified above
                        unsafe {
                            batch_rzz_tile_bmi2(chunk, base, bmi2_groups, num_groups);
                        }
                    });
                return;
            }
        }

        let groups = self.batch_tables.rzz();
        let num_groups = build_batch_rzz_tables(edges, groups);
        let groups = &*groups;

        self.state
            .par_chunks_mut(MIN_PAR_ELEMS)
            .enumerate()
            .for_each(|(chunk_idx, chunk)| {
                let base = chunk_idx * MIN_PAR_ELEMS;
                for (j, amp) in chunk.iter_mut().enumerate() {
                    let i = base + j;
                    let mut combined = Complex64::new(1.0, 0.0);
                    for group in groups.iter().take(num_groups) {
                        let bits = extract_rzz_bits(i, group);
                        combined *= group.table[bits];
                    }
                    *amp *= combined;
                }
            });
    }

    #[inline(always)]
    pub(super) fn apply_cu(&mut self, control: usize, target: usize, mat: [[Complex64; 2]; 2]) {
        #[cfg(feature = "parallel")]
        if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
            self.apply_cu_par(control, target, mat);
            return;
        }

        let prepared = simd::PreparedGate1q::new(&mat);

        if control > target {
            let ctrl_half = 1usize << control;
            let block_size = ctrl_half << 1;
            for chunk in self.state.chunks_mut(block_size) {
                let (_, hi) = chunk.split_at_mut(ctrl_half);
                prepared.apply_full_sequential(hi, target);
            }
        } else {
            let ctrl_mask = 1usize << control;
            let tgt_mask = 1usize << target;
            let num_iters = 1usize << (self.num_qubits - 2);
            let base_ptr = self.state.as_mut_ptr() as *mut f64;
            for i in 0..num_iters {
                let base = insert_zero_bit(insert_zero_bit(i, control), target);
                let idx0 = base | ctrl_mask;
                let idx1 = idx0 | tgt_mask;
                // SAFETY: indices from insert_zero_bit bijection are in-bounds and disjoint.
                unsafe {
                    prepared.apply_pair_ptr(base_ptr.add(idx0 * 2), base_ptr.add(idx1 * 2));
                }
            }
        }
    }

    #[cfg(feature = "parallel")]
    #[inline(always)]
    fn apply_cu_par(&mut self, control: usize, target: usize, mat: [[Complex64; 2]; 2]) {
        if control > target {
            let ctrl_half = 1usize << control;
            let block_size = ctrl_half << 1;
            let prepared = simd::PreparedGate1q::new(&mat);
            let num_blocks = self.state.len() / block_size;

            if num_blocks >= 4 {
                self.state
                    .par_chunks_mut(block_size)
                    .with_min_len(chunk_min_len(block_size))
                    .for_each(|chunk| {
                        let (_, hi) = chunk.split_at_mut(ctrl_half);
                        prepared.apply_full_sequential(hi, target);
                    });
            } else {
                let tgt_half = 1usize << target;
                let tgt_block = tgt_half << 1;
                let inner_tile = MIN_PAR_ELEMS.max(tgt_block);
                for chunk in self.state.chunks_mut(block_size) {
                    let (_, hi) = chunk.split_at_mut(ctrl_half);
                    hi.par_chunks_mut(inner_tile).for_each(|tile| {
                        prepared.apply_full_sequential(tile, target);
                    });
                }
            }
        } else {
            let ctrl_mask = 1usize << control;
            let tgt_mask = 1usize << target;
            let lo = control;
            let hi = target;
            let num_iters = 1usize << (self.num_qubits - 2);
            let ptr = SendPtr(self.state.as_mut_ptr());
            let prepared = simd::PreparedGate1q::new(&mat);

            // SAFETY: insert_zero_bit bijection gives disjoint index pairs per iteration.
            (0..num_iters)
                .into_par_iter()
                .with_min_len(MIN_PAR_ITERS)
                .for_each(move |i| {
                    let base = insert_zero_bit(insert_zero_bit(i, lo), hi);
                    let idx0 = base | ctrl_mask;
                    let idx1 = idx0 | tgt_mask;
                    // SAFETY: idx0 and idx1 are in bounds and unique for this iteration.
                    // The outer iterator maps every pair once, so Rayon tasks do not alias.
                    unsafe {
                        let fp = ptr.as_f64_ptr();
                        prepared.apply_pair_ptr(fp.add(idx0 * 2), fp.add(idx1 * 2));
                    }
                });
        }
    }

    #[inline(always)]
    pub(super) fn apply_unitary(&mut self, targets: &[usize], mat: &[Complex64]) {
        #[cfg(feature = "parallel")]
        if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
            par_apply_unitary(&mut self.state, targets, mat);
            return;
        }
        apply_unitary_seq(&mut self.state, targets, mat);
    }

    #[inline(always)]
    pub(super) fn apply_mcu(
        &mut self,
        controls: &[usize],
        target: usize,
        mat: [[Complex64; 2]; 2],
    ) {
        #[cfg(feature = "parallel")]
        if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
            self.apply_mcu_par(controls, target, mat);
            return;
        }

        let ctrl_mask: usize = controls.iter().map(|&q| 1usize << q).fold(0, |a, b| a | b);
        let tgt_mask = 1usize << target;
        let mut sorted_buf = [0usize; MCU_QUBIT_BUF];
        let num_special = sorted_mcu_qubits(controls, target, &mut sorted_buf);
        let sorted = &sorted_buf[..num_special];

        let num_iters = 1usize << (self.num_qubits - num_special);
        let prepared = simd::PreparedGate1q::new(&mat);
        let base_ptr = self.state.as_mut_ptr() as *mut f64;

        for i in 0..num_iters {
            let mut base = i;
            for &q in sorted {
                base = insert_zero_bit(base, q);
            }
            let idx0 = base | ctrl_mask;
            let idx1 = idx0 | tgt_mask;
            // SAFETY: indices from insert_zero_bit bijection are in-bounds and disjoint.
            unsafe {
                prepared.apply_pair_ptr(base_ptr.add(idx0 * 2), base_ptr.add(idx1 * 2));
            }
        }
    }

    #[cfg(feature = "parallel")]
    #[inline(always)]
    fn apply_mcu_par(&mut self, controls: &[usize], target: usize, mat: [[Complex64; 2]; 2]) {
        let ctrl_mask: usize = controls.iter().map(|&q| 1usize << q).fold(0, |a, b| a | b);
        let tgt_mask = 1usize << target;
        let mut sorted_buf = [0usize; MCU_QUBIT_BUF];
        let num_special = sorted_mcu_qubits(controls, target, &mut sorted_buf);
        let sorted = &sorted_buf[..num_special];

        let num_iters = 1usize << (self.num_qubits - num_special);
        let ptr = SendPtr(self.state.as_mut_ptr());
        let prepared = simd::PreparedGate1q::new(&mat);

        // SAFETY: insert_zero_bit bijection gives disjoint index pairs per iteration.
        (0..num_iters)
            .into_par_iter()
            .with_min_len(MIN_PAR_ITERS)
            .for_each(move |i| {
                let mut base = i;
                for &q in sorted {
                    base = insert_zero_bit(base, q);
                }
                let idx0 = base | ctrl_mask;
                let idx1 = idx0 | tgt_mask;
                // SAFETY: idx0 and idx1 are in bounds and unique for this iteration.
                // The outer iterator maps every pair once, so Rayon tasks do not alias.
                unsafe {
                    let fp = ptr.as_f64_ptr();
                    prepared.apply_pair_ptr(fp.add(idx0 * 2), fp.add(idx1 * 2));
                }
            });
    }

    #[inline(always)]
    pub(super) fn apply_cu_phase(&mut self, control: usize, target: usize, phase: Complex64) {
        #[cfg(feature = "parallel")]
        if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
            self.apply_cu_phase_par(control, target, phase);
            return;
        }

        let (lo, hi) = if control < target {
            (control, target)
        } else {
            (target, control)
        };
        let lo_half = 1usize << lo;
        let lo_block = lo_half << 1;
        let hi_half = 1usize << hi;
        let block_size = hi_half << 1;

        for chunk in self.state.chunks_mut(block_size) {
            let hi_group = &mut chunk[hi_half..];
            for sub in hi_group.chunks_mut(lo_block) {
                simd::scale_complex_slice(&mut sub[lo_half..], phase);
            }
        }
    }

    #[cfg(feature = "parallel")]
    #[inline(always)]
    fn apply_cu_phase_par(&mut self, control: usize, target: usize, phase: Complex64) {
        let (lo, hi) = if control < target {
            (control, target)
        } else {
            (target, control)
        };
        let lo_half = 1usize << lo;
        let lo_block = lo_half << 1;
        let hi_half = 1usize << hi;
        let block_size = hi_half << 1;
        let num_blocks = self.state.len() / block_size;

        if num_blocks >= 4 {
            self.state
                .par_chunks_mut(block_size)
                .with_min_len(chunk_min_len(block_size))
                .for_each(|chunk| {
                    let hi_group = &mut chunk[hi_half..];
                    for sub in hi_group.chunks_mut(lo_block) {
                        simd::scale_complex_slice(&mut sub[lo_half..], phase);
                    }
                });
        } else {
            let inner_tile = MIN_PAR_ELEMS.max(lo_block);
            for chunk in self.state.chunks_mut(block_size) {
                let hi_group = &mut chunk[hi_half..];
                hi_group.par_chunks_mut(inner_tile).for_each(|tile| {
                    for sub in tile.chunks_mut(lo_block) {
                        simd::scale_complex_slice(&mut sub[lo_half..], phase);
                    }
                });
            }
        }
    }

    #[inline(always)]
    pub(super) fn apply_mcu_phase(&mut self, controls: &[usize], target: usize, phase: Complex64) {
        #[cfg(feature = "parallel")]
        if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
            self.apply_mcu_phase_par(controls, target, phase);
            return;
        }

        let all_mask: usize = controls
            .iter()
            .chain(std::iter::once(&target))
            .map(|&q| 1usize << q)
            .fold(0, |a, b| a | b);
        let mut sorted_buf = [0usize; MCU_QUBIT_BUF];
        let num_special = sorted_mcu_qubits(controls, target, &mut sorted_buf);
        let sorted = &sorted_buf[..num_special];

        let num_iters = 1usize << (self.num_qubits - num_special);
        for i in 0..num_iters {
            let mut base = i;
            for &q in sorted {
                base = insert_zero_bit(base, q);
            }
            self.state[base | all_mask] *= phase;
        }
    }

    #[cfg(feature = "parallel")]
    #[inline(always)]
    fn apply_mcu_phase_par(&mut self, controls: &[usize], target: usize, phase: Complex64) {
        let all_mask: usize = controls
            .iter()
            .chain(std::iter::once(&target))
            .map(|&q| 1usize << q)
            .fold(0, |a, b| a | b);
        let mut sorted_buf = [0usize; MCU_QUBIT_BUF];
        let num_special = sorted_mcu_qubits(controls, target, &mut sorted_buf);
        let sorted = &sorted_buf[..num_special];

        let num_iters = 1usize << (self.num_qubits - num_special);
        let ptr = SendPtr(self.state.as_mut_ptr());

        // SAFETY: insert_zero_bit bijection gives disjoint indices per iteration.
        (0..num_iters)
            .into_par_iter()
            .with_min_len(MIN_PAR_ITERS)
            .for_each(move |i| {
                let mut base = i;
                for &q in sorted {
                    base = insert_zero_bit(base, q);
                }
                let idx = base | all_mask;
                // SAFETY: idx is in bounds and unique for this iteration. The
                // insert_zero_bit mapping excludes all special qubits before masks are set.
                unsafe {
                    let val = ptr.load(idx);
                    ptr.store(idx, val * phase);
                }
            });
    }

    #[inline(always)]
    pub(super) fn apply_swap(&mut self, q0: usize, q1: usize) {
        #[cfg(feature = "parallel")]
        if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
            self.apply_swap_par(q0, q1);
            return;
        }

        let (lo, hi) = if q0 < q1 { (q0, q1) } else { (q1, q0) };
        let lo_half = 1usize << lo;
        let lo_block = lo_half << 1;
        let hi_half = 1usize << hi;
        let block_size = hi_half << 1;

        for chunk in self.state.chunks_mut(block_size) {
            let (lo_group, hi_group) = chunk.split_at_mut(hi_half);
            for (lo_sub, hi_sub) in lo_group
                .chunks_mut(lo_block)
                .zip(hi_group.chunks_mut(lo_block))
            {
                let (_, lo_sub_hi) = lo_sub.split_at_mut(lo_half);
                let (hi_sub_lo, _) = hi_sub.split_at_mut(lo_half);
                simd::swap_slices(lo_sub_hi, hi_sub_lo);
            }
        }
    }

    #[cfg(feature = "parallel")]
    #[inline(always)]
    fn apply_swap_par(&mut self, q0: usize, q1: usize) {
        let (lo, hi) = if q0 < q1 { (q0, q1) } else { (q1, q0) };
        let lo_half = 1usize << lo;
        let lo_block = lo_half << 1;
        let hi_half = 1usize << hi;
        let block_size = hi_half << 1;
        let num_blocks = self.state.len() / block_size;

        if num_blocks >= 4 {
            self.state
                .par_chunks_mut(block_size)
                .with_min_len(chunk_min_len(block_size))
                .for_each(|chunk| {
                    let (lo_group, hi_group) = chunk.split_at_mut(hi_half);
                    let lo_subs = lo_group.chunks_mut(lo_block);
                    let hi_subs = hi_group.chunks_mut(lo_block);
                    for (lo_sub, hi_sub) in lo_subs.zip(hi_subs) {
                        let (_, lo_sub_hi) = lo_sub.split_at_mut(lo_half);
                        let (hi_sub_lo, _) = hi_sub.split_at_mut(lo_half);
                        simd::swap_slices(lo_sub_hi, hi_sub_lo);
                    }
                });
        } else {
            let inner_tile = MIN_PAR_ELEMS.max(lo_block);
            for chunk in self.state.chunks_mut(block_size) {
                let (lo_group, hi_group) = chunk.split_at_mut(hi_half);
                lo_group
                    .par_chunks_mut(inner_tile)
                    .zip(hi_group.par_chunks_mut(inner_tile))
                    .for_each(|(lo_tile, hi_tile)| {
                        for (lo_sub, hi_sub) in lo_tile
                            .chunks_mut(lo_block)
                            .zip(hi_tile.chunks_mut(lo_block))
                        {
                            let (_, lo_sub_hi) = lo_sub.split_at_mut(lo_half);
                            let (hi_sub_lo, _) = hi_sub.split_at_mut(lo_half);
                            simd::swap_slices(lo_sub_hi, hi_sub_lo);
                        }
                    });
            }
        }
    }

    #[inline(always)]
    pub(crate) fn apply_fused_2q(&mut self, q0: usize, q1: usize, mat: &[[Complex64; 4]; 4]) {
        if q0.max(q1) - q0.min(q1) == 1 {
            #[cfg(feature = "parallel")]
            if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
                self.apply_fused_2q_adjacent_par(q0, q1, mat);
                return;
            }

            self.apply_fused_2q_adjacent(q0, q1, mat);
            return;
        }

        #[cfg(feature = "parallel")]
        if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
            self.apply_fused_2q_par(q0, q1, mat);
            return;
        }

        let prepared = simd::PreparedGate2q::new(mat);
        prepared.apply_full(&mut self.state, self.num_qubits, q0, q1);
    }

    #[inline(always)]
    fn apply_fused_2q_adjacent(&mut self, q0: usize, q1: usize, mat: &[[Complex64; 4]; 4]) {
        let lo = q0.min(q1);
        let stride = 1usize << lo;
        let block_size = stride << 2;
        let q0_is_lo = q0 == lo;
        let prepared = simd::PreparedGate2q::new(mat);

        for chunk in self.state.chunks_mut(block_size) {
            let ptr = chunk.as_mut_ptr() as *mut f64;
            for offset in 0..stride {
                let i = adjacent_2q_indices(offset, stride, q0_is_lo);
                // SAFETY: each offset maps to one disjoint 4-amplitude group in
                // this chunk, and all indices are below block_size.
                unsafe {
                    prepared.apply_group_ptr(ptr, i);
                }
            }
        }
    }

    #[cfg(feature = "parallel")]
    #[inline(always)]
    fn apply_fused_2q_adjacent_group(
        ptr: SendPtr,
        group: usize,
        lo: usize,
        stride: usize,
        block_size: usize,
        q0_is_lo: bool,
        prepared: &simd::PreparedGate2q,
    ) {
        let block = group >> lo;
        let offset = group & (stride - 1);
        let base = block * block_size;
        let local = adjacent_2q_indices(offset, stride, q0_is_lo);
        let i = [
            base + local[0],
            base + local[1],
            base + local[2],
            base + local[3],
        ];
        // SAFETY: adjacent group mapping partitions the state into disjoint
        // 4-amplitude groups, and group < state.len() / 4.
        unsafe {
            prepared.apply_group_ptr(ptr.as_f64_ptr(), i);
        }
    }

    #[cfg(feature = "parallel")]
    #[inline(always)]
    fn apply_fused_2q_par(&mut self, q0: usize, q1: usize, mat: &[[Complex64; 4]; 4]) {
        let mask0 = 1usize << q0;
        let mask1 = 1usize << q1;
        let (lo, hi) = if q0 < q1 { (q0, q1) } else { (q1, q0) };
        let n_iter = 1usize << (self.num_qubits - 2);
        let ptr = SendPtr(self.state.as_mut_ptr());
        let prepared = simd::PreparedGate2q::new(mat);

        // SAFETY: insert_zero_bit bijection gives disjoint index groups per iteration.
        (0..n_iter)
            .into_par_iter()
            .with_min_len(MIN_PAR_ITERS)
            .for_each(move |k| {
                let base = insert_zero_bit(insert_zero_bit(k, lo), hi);
                let i = [base, base | mask1, base | mask0, base | mask0 | mask1];
                // SAFETY: disjoint groups, ptr valid for state lifetime.
                unsafe {
                    prepared.apply_group_ptr(ptr.as_f64_ptr(), i);
                }
            });
    }

    #[cfg(feature = "parallel")]
    #[inline(always)]
    fn apply_fused_2q_adjacent_par(&mut self, q0: usize, q1: usize, mat: &[[Complex64; 4]; 4]) {
        let lo = q0.min(q1);
        let stride = 1usize << lo;
        let block_size = stride << 2;
        let q0_is_lo = q0 == lo;
        let n_groups = self.state.len() >> 2;
        let ptr = SendPtr(self.state.as_mut_ptr());
        let prepared = simd::PreparedGate2q::new(mat);

        (0..n_groups)
            .into_par_iter()
            .with_min_len(MIN_PAR_ITERS)
            .for_each(move |group| {
                Self::apply_fused_2q_adjacent_group(
                    ptr, group, lo, stride, block_size, q0_is_lo, &prepared,
                );
            });
    }

    #[inline(always)]
    pub(super) fn apply_batch_phase(&mut self, control: usize, phases: &[(usize, Complex64)]) {
        #[cfg(feature = "parallel")]
        if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
            self.apply_batch_phase_par(control, phases);
            return;
        }

        let groups = self.batch_tables.phase();
        let num_groups = build_batch_phase_tables(phases, groups);

        #[cfg(target_arch = "x86_64")]
        if simd::has_bmi2() && simd::has_fma() {
            // SAFETY: BMI2+FMA availability verified above.
            unsafe {
                apply_batch_phase_bmi2(&mut self.state, control, groups, num_groups);
            }
            return;
        }

        let ctrl_mask = 1usize << control;
        let half = 1usize << (self.num_qubits - 1);
        for k in 0..half {
            let i = insert_zero_bit(k, control) | ctrl_mask;
            let mut combined = ONE;
            for group in groups.iter().take(num_groups) {
                let bits = extract_bits(i, &group.shifts, group.len);
                combined *= group.table[bits];
            }
            self.state[i] *= combined;
        }
    }

    #[cfg(feature = "parallel")]
    #[inline(always)]
    fn apply_batch_phase_par(&mut self, control: usize, phases: &[(usize, Complex64)]) {
        let groups = self.batch_tables.phase();
        let num_groups = build_batch_phase_tables(phases, groups);
        let groups = &*groups;

        #[cfg(target_arch = "x86_64")]
        let use_bmi2 = simd::has_bmi2() && simd::has_fma();

        self.state
            .par_chunks_mut(MIN_PAR_ELEMS)
            .enumerate()
            .for_each(|(tile_idx, tile)| {
                let base = tile_idx * MIN_PAR_ELEMS;

                #[cfg(target_arch = "x86_64")]
                if use_bmi2 {
                    // SAFETY: BMI2+FMA availability verified above.
                    unsafe {
                        batch_phase_tile_bmi2(tile, base, control, groups, num_groups);
                    }
                    return;
                }

                let ctrl_mask = 1usize << control;
                let tile_bits = tile.len().trailing_zeros() as usize;
                if control >= tile_bits {
                    if (base & ctrl_mask) == 0 {
                        return;
                    }
                    for (j, amp) in tile.iter_mut().enumerate() {
                        let i = base + j;
                        let mut combined = ONE;
                        for group in groups.iter().take(num_groups) {
                            let bits = extract_bits(i, &group.shifts, group.len);
                            combined *= group.table[bits];
                        }
                        *amp *= combined;
                    }
                } else {
                    let half = tile.len() >> 1;
                    for k in 0..half {
                        let j = insert_zero_bit(k, control) | ctrl_mask;
                        let i = base + j;
                        let mut combined = ONE;
                        for group in groups.iter().take(num_groups) {
                            let bits = extract_bits(i, &group.shifts, group.len);
                            combined *= group.table[bits];
                        }
                        tile[j] *= combined;
                    }
                }
            });
    }

    /// Apply a QFT on qubits `0..num`, a tiled DIF FFT over every `2^num` amplitude
    /// chunk, in the [`Gate::QftBlock`] form the flags name.
    ///
    /// With `D` the DIF network and `P` the chunk bit reversal, the forward swap form is
    /// `P D` and the forward no-swap form `D`. Conjugate twiddles turn `D` into
    /// `conj(D)`, giving `P conj(D)` for the inverse swap form and `P conj(D) P` for the
    /// inverse no-swap form. Mirroring the labels conjugates each form by `P`.
    pub(super) fn apply_qft_block(
        &mut self,
        num: usize,
        inverse: bool,
        swaps: bool,
        big_endian: bool,
    ) {
        let gather = num >= QFT_GATHER_MIN_QUBITS;
        self.apply_qft_block_with(num, inverse, swaps, big_endian, gather);
    }

    /// [`Self::apply_qft_block`] with the high-stage route chosen: `gather` runs the
    /// stages above the tile as gathered passes when at least [`QFT_GATHER_MIN_STAGES`]
    /// of them sit there, and otherwise they run as full-state radix-4 pairs.
    fn apply_qft_block_with(
        &mut self,
        num: usize,
        inverse: bool,
        swaps: bool,
        big_endian: bool,
        gather: bool,
    ) {
        assert!(num >= 1 && num <= self.num_qubits);
        let reverse_input = big_endian != (inverse && !swaps);
        let reverse_output = big_endian != (swaps || inverse);

        let n = num;
        let inv_sqrt2 = std::f64::consts::FRAC_1_SQRT_2;

        if reverse_input {
            apply_bit_reverse_chunks(&mut self.state, n);
        }

        // Cache-tiled DIF FFT:
        //   - High-stride stages run as gathered passes or full-state passes.
        //   - Low-stride stages run together inside each cache-sized tile.
        let widest = qft_widest_tile_bits_for(
            n,
            qft_max_tile_bits(),
            qft_min_tiles_log2(n, self.num_qubits),
        );
        let tile_bits = if gather && n - widest >= QFT_GATHER_MIN_STAGES {
            qft_gather_high_stages(&mut self.state, n, widest, inverse);
            widest
        } else {
            let tile_bits = qft_tile_bits(n, self.num_qubits);
            self.qft_high_stage_pairs(n, tile_bits, inverse);
            tile_bits
        };

        // Phase 2: low-stride stages on cache-resident tiles.
        let tile_size = 1usize << tile_bits;

        let compact = qft_compact_twiddles(inverse);
        let apply_low_stages_in_tile = |tile: &mut [Complex64]| {
            let mut stage_top = tile_bits;
            let mut stages_left = tile_bits;
            if stages_left % 2 == 1 {
                stage_top -= 1;
                let stride = 1usize << stage_top;
                run_radix2_stage_seq(tile, stride, &compact[stride..stride << 1], 1, inv_sqrt2);
                stages_left -= 1;
            }
            while stages_left > 0 {
                let inner_stride = 1usize << (stage_top - 2);
                fft_stage_pair_in_slice(tile, inner_stride, compact, inv_sqrt2);
                stage_top -= 2;
                stages_left -= 2;
            }
        };

        #[cfg(feature = "parallel")]
        let low_done_parallel =
            if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS && self.state.len() / tile_size >= 4 {
                self.state
                    .par_chunks_mut(tile_size)
                    .for_each(apply_low_stages_in_tile);
                true
            } else {
                false
            };
        #[cfg(not(feature = "parallel"))]
        let low_done_parallel = false;

        if !low_done_parallel {
            for tile in self.state.chunks_mut(tile_size) {
                apply_low_stages_in_tile(tile);
            }
        }

        if reverse_output {
            apply_bit_reverse_chunks(&mut self.state, n);
        }
    }

    /// Run the DIF stages above `tile_bits` as full-state passes: radix-4 pairs, after
    /// one radix-2 stage when the count is odd.
    fn qft_high_stage_pairs(&mut self, n: usize, tile_bits: usize, inverse: bool) {
        let fft_len = 1usize << n;
        let inv_sqrt2 = std::f64::consts::FRAC_1_SQRT_2;
        let twiddles_scaled = qft_twiddles_scaled(n, inverse);
        let twiddles_scaled = twiddles_scaled.as_ref();

        let high_stride_count = n - tile_bits;
        let mut stage_top = n;
        let mut stages_left = high_stride_count;

        let run_radix2_stage = |state: &mut Vec<Complex64>, stage: usize| {
            let stride = 1usize << stage;
            let block_size = stride << 1;
            let twiddle_step = fft_len / block_size;

            #[cfg(feature = "parallel")]
            if state.len() >= 1usize << PARALLEL_THRESHOLD_QUBITS {
                fft_stage_par(
                    state,
                    stride,
                    block_size,
                    twiddle_step,
                    twiddles_scaled,
                    inv_sqrt2,
                );
                return;
            }

            run_radix2_stage_seq(state, stride, twiddles_scaled, twiddle_step, inv_sqrt2);
        };

        if stages_left % 2 == 1 {
            stage_top -= 1;
            run_radix2_stage(&mut self.state, stage_top);
            stages_left -= 1;
        }

        while stages_left > 0 {
            // Fuse pair (stage_top - 1, stage_top - 2): outer stride
            // = 2^(stage_top - 1), inner stride = 2^(stage_top - 2).
            let inner_stride = 1usize << (stage_top - 2);
            fft_stage_pair_par(&mut self.state, inner_stride, twiddles_scaled, inv_sqrt2);
            stage_top -= 2;
            stages_left -= 2;
        }
    }

    #[inline]
    pub(super) fn apply_qft_block_textbook(
        &mut self,
        start: usize,
        num: usize,
        inverse: bool,
        swaps: bool,
        big_endian: bool,
    ) {
        assert!(start + num <= self.num_qubits);
        assert!(num >= 1);

        let h = Gate::H.matrix_2x2();
        for step in qft_textbook_steps(start, num, inverse, swaps, big_endian) {
            match step {
                QftTextbookStep::Hadamard(q) => self.apply_single_gate(q, h),
                QftTextbookStep::CPhase {
                    control,
                    target,
                    theta,
                } => self.apply_cu_phase(control, target, Complex64::from_polar(1.0, theta)),
                QftTextbookStep::Swap(a, b) => self.apply_swap(a, b),
            }
        }
    }

    #[inline(always)]
    pub(super) fn apply_measure(&mut self, qubit: usize, classical_bit: usize) {
        #[cfg(feature = "parallel")]
        if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
            self.apply_measure_par(qubit, classical_bit);
            return;
        }

        let half = 1usize << qubit;
        let block_size = half << 1;
        let norm_sq = self.pending_norm * self.pending_norm;
        let zero = Complex64::new(0.0, 0.0);

        let mut prob_one = 0.0f64;
        for block in self.state.chunks(block_size) {
            for amp in &block[half..] {
                prob_one += amp.norm_sqr();
            }
        }
        prob_one *= norm_sq;

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

        let inv_norm = measurement_inv_norm(outcome, prob_one);

        for block in self.state.chunks_mut(block_size) {
            if outcome {
                for amp in &mut block[..half] {
                    *amp = zero;
                }
            } else {
                for amp in &mut block[half..] {
                    *amp = zero;
                }
            }
        }

        self.pending_norm *= inv_norm;
    }

    /// Project `qubit` onto `outcome` without drawing, given the probability of
    /// one the caller read with [`Self::qubit_probability_one`]: the collapse
    /// half of a measurement, for a caller that chose the outcome itself.
    pub(crate) fn collapse_qubit(&mut self, qubit: usize, outcome: bool, prob_one: f64) {
        let half = 1usize << qubit;
        let block_size = half << 1;

        #[cfg(feature = "parallel")]
        if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
            if self.state.len() / block_size >= 4 {
                self.state
                    .par_chunks_mut(block_size)
                    .with_min_len(chunk_min_len(block_size))
                    .for_each(|chunk| {
                        let (lo, hi) = chunk.split_at_mut(half);
                        simd::zero_slice(if outcome { lo } else { hi });
                    });
            } else {
                for chunk in self.state.chunks_mut(block_size) {
                    let (lo, hi) = chunk.split_at_mut(half);
                    let dropped = if outcome { lo } else { hi };
                    dropped
                        .par_chunks_mut(MIN_PAR_ELEMS)
                        .for_each(simd::zero_slice);
                }
            }
            self.pending_norm *= measurement_inv_norm(outcome, prob_one);
            return;
        }

        for chunk in self.state.chunks_mut(block_size) {
            let (lo, hi) = chunk.split_at_mut(half);
            simd::zero_slice(if outcome { lo } else { hi });
        }
        self.pending_norm *= measurement_inv_norm(outcome, prob_one);
    }

    /// Reset `qubit` along the branch where it measured `outcome`: collapse
    /// onto that outcome and fold it down to `|0>`, as [`Self::apply_reset`]
    /// does after its own draw.
    pub(crate) fn reset_qubit_from(&mut self, qubit: usize, outcome: bool, prob_one: f64) {
        let inv_norm = measurement_inv_norm(outcome, prob_one);

        #[cfg(feature = "parallel")]
        if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
            self.fold_reset_par(qubit, outcome);
            self.pending_norm *= inv_norm;
            return;
        }

        let half = 1usize << qubit;
        for block in self.state.chunks_mut(half << 1) {
            let (lo, hi) = block.split_at_mut(half);
            fold_reset_pair(lo, hi, outcome);
        }
        self.pending_norm *= inv_norm;
    }

    #[inline(always)]
    pub(super) fn qubit_probability_one(&self, qubit: usize) -> f64 {
        #[cfg(feature = "parallel")]
        if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
            return self.qubit_probability_one_par(qubit);
        }

        let half = 1usize << qubit;
        let block_size = half << 1;
        let norm_sq = self.pending_norm * self.pending_norm;

        let mut prob_one = 0.0f64;
        for block in self.state.chunks(block_size) {
            prob_one += simd::norm_sqr_sum(&block[half..]);
        }
        prob_one * norm_sq
    }

    #[cfg(feature = "parallel")]
    fn qubit_probability_one_par(&self, qubit: usize) -> f64 {
        let half = 1usize << qubit;
        let block_size = half << 1;
        let norm_sq = self.pending_norm * self.pending_norm;
        let num_blocks = self.state.len() / block_size;

        let prob_one = if num_blocks >= 4 {
            self.state
                .par_chunks(block_size)
                .with_min_len(chunk_min_len(block_size))
                .map(|chunk| simd::norm_sqr_sum(&chunk[half..]))
                .sum::<f64>()
        } else {
            self.state
                .chunks(block_size)
                .map(|chunk| {
                    chunk[half..]
                        .par_chunks(MIN_PAR_ELEMS)
                        .map(simd::norm_sqr_sum)
                        .sum::<f64>()
                })
                .sum::<f64>()
        };
        prob_one * norm_sq
    }

    #[inline(always)]
    pub(super) fn reduced_density_matrix_one(&self, qubit: usize) -> [[Complex64; 2]; 2] {
        #[cfg(feature = "parallel")]
        if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
            return self.reduced_density_matrix_one_par(qubit);
        }

        let half = 1usize << qubit;
        let block_size = half << 1;
        let norm_sq = self.pending_norm * self.pending_norm;

        let mut p0 = 0.0f64;
        let mut p1 = 0.0f64;
        let mut r = Complex64::new(0.0, 0.0);
        for block in self.state.chunks(block_size) {
            let (b0, b1, br) = super::rdm_block_sums(block, half);
            p0 += b0;
            p1 += b1;
            r += br;
        }

        let scale = Complex64::new(norm_sq, 0.0);
        let r = r * scale;
        [
            [Complex64::new(p0 * norm_sq, 0.0), r.conj()],
            [r, Complex64::new(p1 * norm_sq, 0.0)],
        ]
    }

    /// `rho[t][t'] = sum_e psi[idx(t, e)] * conj(psi[idx(t', e)])`, where `t`
    /// packs `q0` in bit 1 and `q1` in bit 0 and `e` runs over the remaining
    /// `n - 2` qubits. Requires `q0 != q1`; the bit insertion is a bijection
    /// onto the block bases only for distinct positions.
    pub(super) fn reduced_density_matrix_two(&self, q0: usize, q1: usize) -> [[Complex64; 4]; 4] {
        debug_assert_ne!(q0, q1);
        let offsets = [
            0usize,
            1usize << q1,
            1usize << q0,
            (1usize << q0) | (1usize << q1),
        ];
        let (lo, hi) = (q0.min(q1), q0.max(q1));
        let state = self.state.as_slice();
        let groups = state.len() >> 2;
        let block = |g: usize| {
            let base = insert_zero_bit(insert_zero_bit(g, lo), hi);
            let mut amps = [Complex64::new(0.0, 0.0); 4];
            for (t, off) in offsets.iter().enumerate() {
                amps[t] = state[base | off];
            }
            let mut out = [[Complex64::new(0.0, 0.0); 4]; 4];
            for (t, row) in out.iter_mut().enumerate() {
                for (tp, entry) in row.iter_mut().enumerate() {
                    *entry = amps[t] * amps[tp].conj();
                }
            }
            out
        };
        let add = |mut a: [[Complex64; 4]; 4], b: [[Complex64; 4]; 4]| {
            for (row, brow) in a.iter_mut().zip(b.iter()) {
                for (entry, bentry) in row.iter_mut().zip(brow.iter()) {
                    *entry += bentry;
                }
            }
            a
        };
        let zero = || [[Complex64::new(0.0, 0.0); 4]; 4];

        #[cfg(feature = "parallel")]
        let mut rho = if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
            (0..groups)
                .into_par_iter()
                .with_min_len(MIN_PAR_ITERS)
                .map(block)
                .reduce(zero, add)
        } else {
            (0..groups).map(block).fold(zero(), add)
        };
        #[cfg(not(feature = "parallel"))]
        let mut rho = (0..groups).map(block).fold(zero(), add);

        let scale = self.pending_norm * self.pending_norm;
        for row in rho.iter_mut() {
            for entry in row.iter_mut() {
                *entry *= scale;
            }
        }
        rho
    }

    #[cfg(feature = "parallel")]
    fn reduced_density_matrix_one_par(&self, qubit: usize) -> [[Complex64; 2]; 2] {
        let half = 1usize << qubit;
        let block_size = half << 1;
        let norm_sq = self.pending_norm * self.pending_norm;
        let num_blocks = self.state.len() / block_size;
        let zero_sums = || (0.0, 0.0, Complex64::new(0.0, 0.0));

        let (p0, p1, r) = if num_blocks >= 4 {
            self.state
                .par_chunks(block_size)
                .with_min_len(chunk_min_len(block_size))
                .map(|block| super::rdm_block_sums(block, half))
                .reduce(zero_sums, super::rdm_sum_add)
        } else {
            self.state
                .chunks(block_size)
                .map(|block| {
                    let (lo, hi) = block.split_at(half);
                    lo.par_chunks(MIN_PAR_ELEMS)
                        .zip(hi.par_chunks(MIN_PAR_ELEMS))
                        .map(|(lo_t, hi_t)| super::rdm_pair_sums(lo_t, hi_t))
                        .reduce(zero_sums, super::rdm_sum_add)
                })
                .fold(zero_sums(), super::rdm_sum_add)
        };

        let scale = Complex64::new(norm_sq, 0.0);
        let r = r * scale;
        [
            [Complex64::new(p0 * norm_sq, 0.0), r.conj()],
            [r, Complex64::new(p1 * norm_sq, 0.0)],
        ]
    }

    #[inline(always)]
    /// Reset `qubit` to `|0>` as one trajectory of the reset channel: sample
    /// the measurement outcome, collapse onto it, then fold the `|1>` branch
    /// down to `|0>`. See the [`Backend::reset`](crate::backend::Backend::reset)
    /// contract for why a projection onto `|0>` is not equivalent.
    pub(super) fn apply_reset(&mut self, qubit: usize) {
        let prob_one = self.qubit_probability_one(qubit);
        let outcome = self.rng.random::<f64>() < prob_one;
        self.reset_qubit_from(qubit, outcome, prob_one);
    }

    /// Apply multiple single-qubit gates in a multi-tier tiled pass.
    ///
    /// Three tiers based on gate target qubit:
    /// - L2 tier (target 0..13): 256KB tiles, all applied per tile in L2 cache
    /// - L3 tier (target 14..16): 2MB tiles, applied per tile in L3 cache where
    ///   one tile per worker fits it, otherwise folded into the shared traversal
    /// - Individual (target 17+): one shared traversal, or a single
    ///   full-state pass when only one gate lands there
    ///
    /// Where that costs more passes, the targets above the L2 tier are gathered into
    /// subcube tiles as `Multi2q` gathers its high qubits (see [`subcube_plan`]), up
    /// to the same high budget per pass, together with the low targets when they all
    /// fit one tile.
    #[inline(always)]
    pub(crate) fn apply_multi_1q(&mut self, gates: &[(usize, [[Complex64; 2]; 2])]) {
        if gates.is_empty() {
            return;
        }
        if gates.len() == 1 {
            self.apply_single_gate(gates[0].0, gates[0].1);
            return;
        }

        #[cfg(feature = "parallel")]
        if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
            apply_multi_1q_par(&mut self.state, gates);
            return;
        }

        let max_l2_target = multi_gate_max_l2_target();
        let (route, _) = multi_1q_route(gates, self.num_qubits, max_l2_target, true, 1);
        let mut small_gates = SmallVec::new();
        let mut medium_gates = SmallVec::new();
        let mut large_gates = SmallVec::new();
        split_multi_1q_tiers(
            gates,
            max_l2_target,
            matches!(route, Multi1qRoute::Tiered),
            &mut small_gates,
            &mut medium_gates,
            &mut large_gates,
        );

        if let Multi1qRoute::Whole(plan) = &route {
            let tiled = whole_subcube_gates(plan, small_gates, &large_gates);
            apply_multi_1q_subcubes(&mut self.state, plan, &tiled);
            return;
        }

        if !small_gates.is_empty() {
            let tile_size = multi_gate_l2_tile();
            for tile in self.state.chunks_mut(tile_size) {
                for &(target, ref prepared) in &small_gates {
                    prepared.apply_tiled(tile, target);
                }
            }
        }

        if !medium_gates.is_empty() {
            let outer_block = 1usize << (MULTI_GATE_MAX_L3_TARGET + 1);
            let tile_size = MULTI_GATE_L3_TILE.max(outer_block);
            for tile in self.state.chunks_mut(tile_size) {
                for &(target, ref prepared) in &medium_gates {
                    prepared.apply_tiled(tile, target);
                }
            }
        }

        if let Multi1qRoute::Split(runs) = &route {
            apply_multi_1q_runs(&mut self.state, &large_gates, runs);
            return;
        }

        for (target, mat) in large_gates {
            self.apply_single_gate(target, mat);
        }
    }

    pub(super) fn apply_multi_1q_diagonal(&mut self, gates: &[(usize, [[Complex64; 2]; 2])]) {
        if gates.is_empty() {
            return;
        }
        if gates.len() == 1 {
            self.apply_diagonal_gate(gates[0].0, gates[0].1[0][0], gates[0].1[1][1]);
            return;
        }

        #[cfg(feature = "parallel")]
        if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
            self.apply_multi_1q_diagonal_par(gates);
            return;
        }

        let mut small_gates = SmallVec::new();
        let mut medium_gates = SmallVec::new();
        let mut large_gates = SmallVec::new();
        split_multi_1q_diagonal_tiers(
            gates,
            multi_gate_max_l2_target(),
            &mut small_gates,
            &mut medium_gates,
            &mut large_gates,
        );

        if !small_gates.is_empty() {
            let tile_size = multi_gate_l2_tile();
            for tile in self.state.chunks_mut(tile_size) {
                for &(target, d0, d1, skip_lo) in &small_gates {
                    simd::apply_diagonal_sequential(tile, target, d0, d1, skip_lo);
                }
            }
        }

        if !medium_gates.is_empty() {
            let outer_block = 1usize << (MULTI_GATE_MAX_L3_TARGET + 1);
            let tile_size = MULTI_GATE_L3_TILE.max(outer_block);
            for tile in self.state.chunks_mut(tile_size) {
                for &(target, d0, d1, skip_lo) in &medium_gates {
                    simd::apply_diagonal_sequential(tile, target, d0, d1, skip_lo);
                }
            }
        }

        for (target, d0, d1) in large_gates {
            let skip_lo = is_phase_one(d0);
            simd::apply_diagonal_sequential(&mut self.state, target, d0, d1, skip_lo);
        }
    }

    #[cfg(feature = "parallel")]
    #[inline(always)]
    fn apply_multi_1q_diagonal_par(&mut self, gates: &[(usize, [[Complex64; 2]; 2])]) {
        let mut small_gates = SmallVec::new();
        let mut medium_gates = SmallVec::new();
        let mut large_gates = SmallVec::new();
        let tile_size = multi_gate_l2_tile_par(self.state.len());
        split_multi_1q_diagonal_tiers(
            gates,
            max_target_for_tile(tile_size),
            &mut small_gates,
            &mut medium_gates,
            &mut large_gates,
        );

        if !small_gates.is_empty() {
            self.state
                .par_chunks_mut(tile_size)
                .with_min_len(chunk_min_len(tile_size))
                .for_each(|tile| {
                    for &(target, d0, d1, skip_lo) in &small_gates {
                        simd::apply_diagonal_sequential(tile, target, d0, d1, skip_lo);
                    }
                });
        }

        if !medium_gates.is_empty() {
            let outer_block = 1usize << (MULTI_GATE_MAX_L3_TARGET + 1);
            let tile_size = MULTI_GATE_L3_TILE.max(outer_block);
            self.state
                .par_chunks_mut(tile_size)
                .with_min_len(chunk_min_len(tile_size))
                .for_each(|tile| {
                    for &(target, d0, d1, skip_lo) in &medium_gates {
                        simd::apply_diagonal_sequential(tile, target, d0, d1, skip_lo);
                    }
                });
        }

        for (target, d0, d1) in large_gates {
            let skip_lo = is_phase_one(d0);
            let block_size = 1usize << (target + 1);
            let num_blocks = self.state.len() / block_size;
            if num_blocks < 4 {
                apply_diagonal_high_target_par(&mut self.state, target, d0, d1, skip_lo);
                continue;
            }
            let tile_size = block_size.max(MIN_PAR_ELEMS);
            self.state.par_chunks_mut(tile_size).for_each(|tile| {
                simd::apply_diagonal_sequential(tile, target, d0, d1, skip_lo);
            });
        }
    }

    fn apply_multi_2q_subcube(
        &mut self,
        gates: &[(usize, usize, [[Complex64; 4]; 4])],
        plan: &SubcubePlan,
    ) {
        let prepared = prepare_2q(gates);
        let rest = plan.rest(self.num_qubits);
        let state = self.state.as_mut_ptr();
        for outer in 0..1usize << (self.num_qubits - plan.tile_bits) {
            // SAFETY: `state` holds 2^num_qubits amplitudes and each `outer`
            // addresses its own subcube, applied one at a time here.
            unsafe { apply_subcube(state, plan, &rest, outer, &prepared) };
        }
    }

    #[cfg(feature = "parallel")]
    fn apply_multi_2q_subcube_par(
        &mut self,
        gates: &[(usize, usize, [[Complex64; 4]; 4])],
        plan: &SubcubePlan,
    ) {
        let prepared = prepare_2q(gates);
        let rest = plan.rest(self.num_qubits);
        let ptr = SendPtr(self.state.as_mut_ptr());
        (0..1usize << (self.num_qubits - plan.tile_bits))
            .into_par_iter()
            .for_each(|outer| {
                // SAFETY: `ptr` holds 2^num_qubits amplitudes and distinct
                // `outer` values address disjoint subcubes, so no two tasks
                // touch the same amplitude.
                unsafe { apply_subcube(ptr.as_complex_ptr(), plan, &rest, outer, &prepared) };
            });
    }

    /// Apply multiple two-qubit gates in a cache-tiled pass.
    ///
    /// A batch whose qubits reach past the lowest tile bits takes the subcube
    /// path (see [`subcube_plan`]): each subcube (256 KB by default, set by
    /// `PRISM_MULTI_2Q_TILE_BITS`) is gathered into a per-thread tile, every
    /// gate is applied there in order, and the tile is scattered back, so the
    /// whole batch costs one pass over the state plus in-cache arithmetic.
    /// A batch that sits entirely below the tile bits (14 by default) tiles the
    /// state in place. Any other list falls back to the tiered pass by
    /// `max(q0, q1)`: L2 tiles up to qubit 13, L3 tiles up to 16, one full pass
    /// per gate above that. That pass runs tier by tier, so it keeps application
    /// order only within a tier.
    #[inline(always)]
    pub(crate) fn apply_multi_2q(&mut self, gates: &[(usize, usize, [[Complex64; 4]; 4])]) {
        if gates.is_empty() {
            return;
        }
        if gates.len() == 1 {
            self.apply_fused_2q(gates[0].0, gates[0].1, &gates[0].2);
            return;
        }
        if let Some(plan) = subcube_plan(gates, self.num_qubits) {
            #[cfg(feature = "parallel")]
            if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
                self.apply_multi_2q_subcube_par(gates, &plan);
                return;
            }
            self.apply_multi_2q_subcube(gates, &plan);
            return;
        }

        #[cfg(feature = "parallel")]
        if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
            self.apply_multi_2q_par(gates);
            return;
        }

        let tile_size = multi_2q_small_tile(self.num_qubits);
        let mut small_gates = SmallVec::new();
        let mut medium_gates = SmallVec::new();
        split_multi_2q_tiers(
            gates,
            max_target_for_tile(tile_size),
            &mut small_gates,
            &mut medium_gates,
        );

        if !small_gates.is_empty() {
            let tile_qubits = tile_size.trailing_zeros() as usize;
            for tile in self.state.chunks_mut(tile_size) {
                let n = tile.len().trailing_zeros() as usize;
                for &(q0, q1, ref prepared) in &small_gates {
                    prepared.apply_tiled(tile, n.min(tile_qubits), q0, q1);
                }
            }
        }

        if !medium_gates.is_empty() {
            let tile_size = MULTI_GATE_L3_TILE;
            let tile_qubits = tile_size.trailing_zeros() as usize;
            for tile in self.state.chunks_mut(tile_size) {
                let n = tile.len().trailing_zeros() as usize;
                for &(q0, q1, ref prepared) in &medium_gates {
                    prepared.apply_tiled(tile, n.min(tile_qubits), q0, q1);
                }
            }
        }

        self.apply_multi_2q_untiled(gates);
    }

    /// Apply, one full pass each, the gates whose high qubit is above the L3 tile.
    #[inline(always)]
    fn apply_multi_2q_untiled(&mut self, gates: &[(usize, usize, [[Complex64; 4]; 4])]) {
        for &(q0, q1, ref mat) in gates {
            if q0.max(q1) > MULTI_GATE_MAX_L3_TARGET {
                self.apply_fused_2q(q0, q1, mat);
            }
        }
    }

    /// The tiered pass with the in-place tile cut for the pool. The L3 tier stays one
    /// sweep even where the state is a single L3 tile: at 17 qubits per-gate parallel
    /// passes in its place read 28% slower on `statevector/random_d10/17`.
    #[cfg(feature = "parallel")]
    #[inline(always)]
    fn apply_multi_2q_par(&mut self, gates: &[(usize, usize, [[Complex64; 4]; 4])]) {
        let tile_size = multi_2q_small_tile(self.num_qubits);
        let mut small_gates = SmallVec::new();
        let mut medium_gates = SmallVec::new();
        split_multi_2q_tiers(
            gates,
            max_target_for_tile(tile_size),
            &mut small_gates,
            &mut medium_gates,
        );

        if !small_gates.is_empty() {
            let tile_qubits = tile_size.trailing_zeros() as usize;
            self.state
                .par_chunks_mut(tile_size)
                .with_min_len(chunk_min_len(tile_size))
                .for_each(|tile| {
                    let n = tile.len().trailing_zeros() as usize;
                    for &(q0, q1, ref prepared) in &small_gates {
                        prepared.apply_tiled(tile, n.min(tile_qubits), q0, q1);
                    }
                });
        }

        if !medium_gates.is_empty() {
            let tile_size = MULTI_GATE_L3_TILE;
            let tile_qubits = tile_size.trailing_zeros() as usize;
            self.state
                .par_chunks_mut(tile_size)
                .with_min_len(chunk_min_len(tile_size))
                .for_each(|tile| {
                    let n = tile.len().trailing_zeros() as usize;
                    for &(q0, q1, ref prepared) in &medium_gates {
                        prepared.apply_tiled(tile, n.min(tile_qubits), q0, q1);
                    }
                });
        }

        self.apply_multi_2q_untiled(gates);
    }

    #[cfg(feature = "parallel")]
    #[inline(always)]
    fn apply_measure_par(&mut self, qubit: usize, classical_bit: usize) {
        let prob_one = self.qubit_probability_one_par(qubit);
        let outcome = self.rng.random::<f64>() < prob_one;
        self.classical_bits[classical_bit] = outcome;
        self.collapse_qubit(qubit, outcome, prob_one);
    }

    #[cfg(feature = "parallel")]
    fn fold_reset_par(&mut self, qubit: usize, outcome: bool) {
        let half = 1usize << qubit;
        let block_size = half << 1;
        let num_blocks = self.state.len() / block_size;

        if num_blocks >= 4 {
            self.state
                .par_chunks_mut(block_size)
                .with_min_len(chunk_min_len(block_size))
                .for_each(|chunk| {
                    let (lo, hi) = chunk.split_at_mut(half);
                    fold_reset_pair(lo, hi, outcome);
                });
            return;
        }

        for chunk in self.state.chunks_mut(block_size) {
            let (lo, hi) = chunk.split_at_mut(half);
            lo.par_chunks_mut(MIN_PAR_ELEMS)
                .zip(hi.par_chunks_mut(MIN_PAR_ELEMS))
                .for_each(|(lo_t, hi_t)| fold_reset_pair(lo_t, hi_t, outcome));
        }
    }

    pub(super) fn apply_diagonal_batch(&mut self, entries: &[DiagEntry]) {
        #[cfg(feature = "parallel")]
        if self.num_qubits >= PARALLEL_THRESHOLD_QUBITS {
            self.apply_diagonal_batch_par(entries);
            return;
        }

        if !fill_diagonal_batch_tables(entries, self.batch_tables.diag()) {
            self.apply_diagonal_batch_fallback(entries);
            return;
        }
        let built = &*self.batch_tables.diag();
        if built.num_groups == 0 {
            return;
        }

        #[cfg(target_arch = "x86_64")]
        if simd::has_bmi2() && simd::has_fma() {
            // SAFETY: has_bmi2() + has_fma() confirmed above
            unsafe {
                diagonal_batch_tile_bmi2(
                    &mut self.state,
                    0,
                    &built.tables,
                    &built.group_pext_masks,
                    built.num_groups,
                );
            }
            return;
        }

        apply_diagonal_batch_scalar(
            &mut self.state,
            &built.tables,
            &built.unique_qubits,
            &built.group_sizes,
            built.num_groups,
        );
    }

    #[cfg(feature = "parallel")]
    fn apply_diagonal_batch_par(&mut self, entries: &[DiagEntry]) {
        if !fill_diagonal_batch_tables(entries, self.batch_tables.diag()) {
            self.state
                .par_chunks_mut(MIN_PAR_ELEMS)
                .enumerate()
                .for_each(|(chunk_idx, chunk)| {
                    let base = chunk_idx * MIN_PAR_ELEMS;
                    for (j, amp) in chunk.iter_mut().enumerate() {
                        *amp *= diag_entries_phase(base + j, entries);
                    }
                });
            return;
        }
        let built = &*self.batch_tables.diag();
        if built.num_groups == 0 {
            return;
        }

        #[cfg(target_arch = "x86_64")]
        if simd::has_bmi2() && simd::has_fma() {
            self.state
                .par_chunks_mut(MIN_PAR_ELEMS)
                .enumerate()
                .for_each(|(chunk_idx, chunk)| {
                    let base = chunk_idx * MIN_PAR_ELEMS;
                    // SAFETY: has_bmi2() + has_fma() confirmed above
                    unsafe {
                        diagonal_batch_tile_bmi2(
                            chunk,
                            base,
                            &built.tables,
                            &built.group_pext_masks,
                            built.num_groups,
                        );
                    }
                });
            return;
        }

        let group_shifts = diag_batch_group_shifts(&built.unique_qubits, &built.group_sizes);
        self.state
            .par_chunks_mut(MIN_PAR_ELEMS)
            .enumerate()
            .for_each(|(chunk_idx, chunk)| {
                diagonal_batch_tile_scalar(
                    chunk,
                    chunk_idx * MIN_PAR_ELEMS,
                    &built.tables,
                    &group_shifts,
                    &built.group_sizes,
                    built.num_groups,
                );
            });
    }

    fn apply_diagonal_batch_fallback(&mut self, entries: &[DiagEntry]) {
        for (i, amp) in self.state.iter_mut().enumerate() {
            *amp *= diag_entries_phase(i, entries);
        }
    }
}

#[cfg(target_arch = "x86_64")]
#[target_feature(enable = "bmi2,fma")]
unsafe fn diagonal_batch_tile_bmi2(
    tile: &mut [Complex64],
    base_idx: usize,
    tables: &[[Complex64; DIAG_BATCH_TABLE_SIZE]; MAX_DIAG_BATCH_GROUPS],
    pext_masks: &[u64; MAX_DIAG_BATCH_GROUPS],
    num_groups: usize,
) {
    use std::arch::x86_64::_pext_u64;
    let one = Complex64::new(1.0, 0.0);
    for (j, amp) in tile.iter_mut().enumerate() {
        let i = base_idx + j;
        let mut combined = one;
        for g in 0..num_groups {
            let bits = _pext_u64(i as u64, pext_masks[g]) as usize;
            combined *= tables[g][bits];
        }
        *amp *= combined;
    }
}

/// Per-group qubit positions, the non-x86_64 stand-in for the PEXT masks.
///
/// `unique_qubits` holds the groups concatenated in order, and a group can be shorter than
/// `DIAG_BATCH_MAX_QUBITS_PER_GROUP`, so `group_sizes` drives the split.
fn diag_batch_group_shifts(
    unique_qubits: &[usize],
    group_sizes: &[usize; MAX_DIAG_BATCH_GROUPS],
) -> [[usize; DIAG_BATCH_MAX_QUBITS_PER_GROUP]; MAX_DIAG_BATCH_GROUPS] {
    let mut group_shifts = [[0usize; DIAG_BATCH_MAX_QUBITS_PER_GROUP]; MAX_DIAG_BATCH_GROUPS];
    let mut flat = 0;
    for (shifts, &size) in group_shifts.iter_mut().zip(group_sizes) {
        shifts[..size].copy_from_slice(&unique_qubits[flat..flat + size]);
        flat += size;
    }
    group_shifts
}

fn diagonal_batch_tile_scalar(
    tile: &mut [Complex64],
    base_idx: usize,
    tables: &[[Complex64; DIAG_BATCH_TABLE_SIZE]; MAX_DIAG_BATCH_GROUPS],
    group_shifts: &[[usize; DIAG_BATCH_MAX_QUBITS_PER_GROUP]; MAX_DIAG_BATCH_GROUPS],
    group_sizes: &[usize; MAX_DIAG_BATCH_GROUPS],
    num_groups: usize,
) {
    let one = Complex64::new(1.0, 0.0);
    for (j, amp) in tile.iter_mut().enumerate() {
        let i = base_idx + j;
        let mut combined = one;
        for (g, shifts) in group_shifts.iter().enumerate().take(num_groups) {
            let mut bits = 0usize;
            for (p, &shift) in shifts[..group_sizes[g]].iter().enumerate() {
                bits |= ((i >> shift) & 1) << p;
            }
            combined *= tables[g][bits];
        }
        *amp *= combined;
    }
}

fn apply_diagonal_batch_scalar(
    state: &mut [Complex64],
    tables: &[[Complex64; DIAG_BATCH_TABLE_SIZE]; MAX_DIAG_BATCH_GROUPS],
    unique_qubits: &[usize],
    group_sizes: &[usize; MAX_DIAG_BATCH_GROUPS],
    num_groups: usize,
) {
    let group_shifts = diag_batch_group_shifts(unique_qubits, group_sizes);
    diagonal_batch_tile_scalar(state, 0, tables, &group_shifts, group_sizes, num_groups);
}

#[cfg(all(test, target_arch = "x86_64"))]
mod pext_agreement_tests {
    use super::*;
    use rand::SeedableRng;
    use rand_chacha::ChaCha8Rng;

    const NUM_QUBITS: usize = 12;

    // The scalar and BMI2 kernels index their tables differently (shift loop vs
    // PEXT over a mask), so a table shape that indexes one bit per unique qubit
    // silently diverges from a producer that emits a repeated qubit. Both batch
    // families are checked against the same amplitude for every basis index.
    fn combined_phase_scalar_batch_phase(
        idx: usize,
        groups: &[BatchPhaseGroup],
        n: usize,
    ) -> Complex64 {
        let mut combined = Complex64::new(1.0, 0.0);
        for group in groups.iter().take(n) {
            combined *= group.table[extract_bits(idx, &group.shifts, group.len)];
        }
        combined
    }

    #[test]
    fn batch_phase_pext_matches_scalar_lookup() {
        if !std::is_x86_feature_detected!("bmi2") {
            panic!("bmi2 is required to check the PEXT path against the scalar path");
        }
        let mut rng = ChaCha8Rng::seed_from_u64(42);
        let one = Complex64::new(1.0, 0.0);

        for case in 0..64 {
            let len = 1 + case % BatchPhaseData::MAX_PHASES.min(NUM_QUBITS - 1);
            let mut targets: Vec<usize> = (1..NUM_QUBITS).collect();
            for i in (1..targets.len()).rev() {
                targets.swap(i, rng.random_range(0..=i));
            }
            let phases: Vec<(usize, Complex64)> = targets[..len]
                .iter()
                .map(|&q| (q, Complex64::from_polar(1.0, rng.random_range(-3.0..3.0))))
                .collect();

            let mut groups = [BatchPhaseGroup {
                table: [one; BATCH_PHASE_TABLE_SIZE],
                shifts: [0; BATCH_PHASE_GROUP_SIZE],
                len: 0,
                pext_mask: 0,
            }; MAX_BATCH_PHASE_GROUPS];
            let num_groups = build_batch_phase_tables(&phases, &mut groups);

            for idx in 0..(1usize << NUM_QUBITS) {
                let scalar = combined_phase_scalar_batch_phase(idx, &groups, num_groups);
                let mut pext = one;
                for group in groups.iter().take(num_groups) {
                    // SAFETY: BMI2 checked above.
                    let bits = unsafe { std::arch::x86_64::_pext_u64(idx as u64, group.pext_mask) };
                    pext *= group.table[bits as usize];
                }
                assert!(
                    (scalar - pext).norm() < 1e-12,
                    "case {case} idx {idx}: scalar {scalar} vs pext {pext}"
                );
            }
        }
    }

    // Pair-aware grouping can leave every group short of
    // `DIAG_BATCH_MAX_QUBITS_PER_GROUP`, so the scalar path can no longer split
    // `unique_qubits` at a fixed stride. Two six-qubit chains pin that shape: they
    // cannot share a group, and neither fills one.
    #[test]
    fn diagonal_batch_pext_matches_scalar_lookup_with_short_groups() {
        if !std::is_x86_feature_detected!("bmi2") || !std::is_x86_feature_detected!("fma") {
            panic!("bmi2 and fma are required to check the PEXT path against the scalar path");
        }
        let mut rng = ChaCha8Rng::seed_from_u64(42);
        let mut entries = Vec::new();
        for chain in [0usize, 6] {
            for q in chain..chain + 5 {
                entries.push(DiagEntry::Phase2q {
                    q0: q,
                    q1: q + 1,
                    phase: Complex64::from_polar(1.0, rng.random_range(-3.0..3.0)),
                });
            }
            entries.push(DiagEntry::Parity2q {
                q0: chain,
                q1: chain + 3,
                same: Complex64::from_polar(1.0, rng.random_range(-3.0..3.0)),
                diff: Complex64::from_polar(1.0, rng.random_range(-3.0..3.0)),
            });
            entries.push(DiagEntry::Phase1q {
                qubit: chain + 2,
                d0: Complex64::from_polar(1.0, rng.random_range(-3.0..3.0)),
                d1: Complex64::from_polar(1.0, rng.random_range(-3.0..3.0)),
            });
        }

        let built = build_diagonal_batch_tables(&entries).expect("groupable");
        assert_eq!(built.num_groups, 2);
        assert_eq!(built.group_sizes[..2], [6, 6]);

        let dim = 1usize << NUM_QUBITS;
        let mut scalar: Vec<Complex64> = (0..dim)
            .map(|i| Complex64::from_polar(1.0, i as f64 * 0.017))
            .collect();
        let mut pext = scalar.clone();
        let shifts = diag_batch_group_shifts(&built.unique_qubits, &built.group_sizes);
        diagonal_batch_tile_scalar(
            &mut scalar,
            0,
            &built.tables,
            &shifts,
            &built.group_sizes,
            built.num_groups,
        );
        // SAFETY: BMI2 and FMA checked above.
        unsafe {
            diagonal_batch_tile_bmi2(
                &mut pext,
                0,
                &built.tables,
                &built.group_pext_masks,
                built.num_groups,
            );
        }
        for (idx, (a, b)) in scalar.iter().zip(&pext).enumerate() {
            assert!((a - b).norm() < 1e-12, "idx {idx}: scalar {a} vs pext {b}");
        }
    }

    // A 24-qubit ring plus four pairs of a matching: runs of eight consecutive
    // edges would put the pairs in a sixteen-qubit group, past the PEXT table.
    #[test]
    fn batch_rzz_pext_packs_a_ring_with_a_matching() {
        if !std::is_x86_feature_detected!("bmi2") {
            panic!("bmi2 is required to check the PEXT path against the scalar path");
        }
        let one = Complex64::new(1.0, 0.0);
        let mut edges: Vec<(usize, usize, f64)> = (0..24)
            .map(|q| (q, (q + 1) % 24, 0.3 + 0.01 * q as f64))
            .collect();
        edges.extend([(3, 17, 0.7), (5, 11, -0.4), (20, 8, 1.1), (1, 14, 0.2)]);
        edges.extend([(6, 22, 0.9), (9, 2, -1.3), (13, 19, 0.5), (4, 15, 0.8)]);

        let mut bmi2_groups = [BatchRzzBmi2Group {
            table: [one; BATCH_RZZ_BMI2_TABLE_SIZE],
            pext_mask: 0,
        }; BATCH_RZZ_BMI2_MAX_GROUPS];
        let num_groups = build_batch_rzz_bmi2_tables(&edges, &mut bmi2_groups)
            .expect("the ring and matching fit the PEXT groups");

        for idx in (0..1usize << 24).step_by(9973) {
            let mut angle = 0.0;
            for &(q0, q1, theta) in &edges {
                let parity = ((idx >> q0) ^ (idx >> q1)) & 1;
                angle += if parity == 0 {
                    -theta / 2.0
                } else {
                    theta / 2.0
                };
            }
            let expected = Complex64::from_polar(1.0, angle);
            let mut pext = one;
            for group in bmi2_groups.iter().take(num_groups) {
                // SAFETY: BMI2 checked above.
                let bits = unsafe { std::arch::x86_64::_pext_u64(idx as u64, group.pext_mask) };
                pext *= group.table[bits as usize];
            }
            assert!(
                (expected - pext).norm() < 1e-12,
                "idx {idx}: {expected} vs {pext}"
            );
        }
    }

    #[test]
    fn batch_rzz_pext_matches_scalar_lookup() {
        if !std::is_x86_feature_detected!("bmi2") {
            panic!("bmi2 is required to check the PEXT path against the scalar path");
        }
        let mut rng = ChaCha8Rng::seed_from_u64(42);
        let one = Complex64::new(1.0, 0.0);

        for case in 0..64 {
            let num_edges = 1 + case % BatchRzzData::MAX_EDGES;
            let edges: Vec<(usize, usize, f64)> = (0..num_edges)
                .map(|_| {
                    let q0 = rng.random_range(0..NUM_QUBITS);
                    let mut q1 = rng.random_range(0..NUM_QUBITS);
                    while q1 == q0 {
                        q1 = rng.random_range(0..NUM_QUBITS);
                    }
                    (q0, q1, rng.random_range(-3.0..3.0))
                })
                .collect();

            let mut scalar_groups = [BatchRzzGroup {
                table: [one; BATCH_RZZ_TABLE_SIZE],
                q0s: [0; BATCH_RZZ_GROUP_SIZE],
                q1s: [0; BATCH_RZZ_GROUP_SIZE],
                len: 0,
            }; MAX_BATCH_RZZ_GROUPS];
            let num_groups = build_batch_rzz_tables(&edges, &mut scalar_groups);

            let mut bmi2_groups = [BatchRzzBmi2Group {
                table: [one; BATCH_RZZ_BMI2_TABLE_SIZE],
                pext_mask: 0,
            }; BATCH_RZZ_BMI2_MAX_GROUPS];
            let Some(bmi2_num_groups) = build_batch_rzz_bmi2_tables(&edges, &mut bmi2_groups)
            else {
                continue;
            };

            for idx in 0..(1usize << NUM_QUBITS) {
                let mut scalar = one;
                for group in scalar_groups.iter().take(num_groups) {
                    scalar *= group.table[extract_rzz_bits(idx, group)];
                }
                let mut pext = one;
                for group in bmi2_groups.iter().take(bmi2_num_groups) {
                    // SAFETY: BMI2 checked above.
                    let bits = unsafe { std::arch::x86_64::_pext_u64(idx as u64, group.pext_mask) };
                    pext *= group.table[bits as usize];
                }
                assert!(
                    (scalar - pext).norm() < 1e-12,
                    "case {case} idx {idx}: scalar {scalar} vs pext {pext}"
                );
            }
        }
    }
}

#[cfg(test)]
mod qft_twiddle_cache_tests {
    use super::*;

    fn cached_sizes(cache: &RwLock<QftTwiddleCache>) -> Vec<usize> {
        let guard = cache.read().unwrap();
        (0..guard.entries.len())
            .filter(|&n| guard.entries[n].is_some())
            .collect()
    }

    // A table for size n holds 2^(n-1) entries, so 2^(n+3) bytes: n=1 is 16,
    // n=2 is 32, n=3 is 64, n=5 is 256.
    #[test]
    fn evicts_least_recently_used_over_the_cap() {
        let cache = RwLock::new(QftTwiddleCache::new(96));
        cached_qft_twiddles(&cache, 2);
        cached_qft_twiddles(&cache, 3);
        assert_eq!(cached_sizes(&cache), vec![2, 3]);
        cached_qft_twiddles(&cache, 2);
        cached_qft_twiddles(&cache, 1);
        assert_eq!(cached_sizes(&cache), vec![1, 2]);
        assert_eq!(cache.read().unwrap().bytes, 48);
    }

    #[test]
    fn oversized_table_sits_alone_until_the_next_insert() {
        let cache = RwLock::new(QftTwiddleCache::new(96));
        cached_qft_twiddles(&cache, 2);
        cached_qft_twiddles(&cache, 5);
        assert_eq!(cached_sizes(&cache), vec![5]);
        assert_eq!(cache.read().unwrap().bytes, 256);
        cached_qft_twiddles(&cache, 1);
        assert_eq!(cached_sizes(&cache), vec![1]);
        assert_eq!(cache.read().unwrap().bytes, 16);
    }

    #[test]
    fn concurrent_misses_share_one_table() {
        let cache = RwLock::new(QftTwiddleCache::new(
            QFT_TWIDDLE_CACHE_DEFAULT_LIMIT_MB * 1024 * 1024,
        ));
        let tables: Vec<QftTwiddleTable> = std::thread::scope(|scope| {
            let handles: Vec<_> = (0..8)
                .map(|_| scope.spawn(|| cached_qft_twiddles(&cache, 6)))
                .collect();
            handles.into_iter().map(|h| h.join().unwrap()).collect()
        });
        assert_eq!(cached_sizes(&cache), vec![6]);
        let guard = cache.read().unwrap();
        let cached = &guard.entries[6].as_ref().unwrap().table;
        assert!(tables.iter().all(|t| Arc::ptr_eq(t, cached)));
        assert_eq!(guard.bytes, qft_twiddle_table_bytes(cached.len()));
    }
}
#[cfg(test)]
mod qft_layout_tests {
    use super::*;

    #[test]
    fn compact_twiddles_match_whole_state_table_bitwise() {
        // The whole-state table has to reach the compact table's widest stride.
        let n = qft_max_tile_bits().max(16);
        let total = 1usize << n;
        for inverse in [false, true] {
            let whole = qft_twiddles_scaled(n, inverse);
            let compact = qft_compact_twiddles(inverse);
            let mut s = 1;
            while s < compact.len() {
                let step = total / (2 * s);
                for k in 0..s {
                    assert_eq!(
                        compact[s + k],
                        whole[k * step],
                        "inverse {inverse}, stride {s}, k {k}"
                    );
                }
                s <<= 1;
            }
        }
    }

    #[test]
    fn qft_tiles_leave_paired_stages_and_eight_tiles() {
        let widths: Vec<usize> = (10..=26).map(|n| qft_tile_bits_for(n, 15, 3)).collect();
        assert_eq!(
            widths,
            [
                10, 11, 12, 13, 13, 13, 13, 13, 14, 15, 14, 15, 14, 15, 14, 15, 14
            ]
        );
    }

    #[test]
    fn tiled_bit_reversal_matches_index_reversal() {
        for bits in [3, 9, 10, 11, 12, 15, 16, 17] {
            let mut state: Vec<Complex64> = (0..1usize << bits)
                .map(|i| Complex64::new(i as f64, -(i as f64)))
                .collect();
            apply_bit_reverse_permutation(&mut state, bits);
            for (i, amp) in state.iter().enumerate() {
                let src = reverse_low_bits(i, bits) as f64;
                assert_eq!(*amp, Complex64::new(src, -src), "bits {bits}, index {i}");
            }
        }
    }
}

#[cfg(test)]
mod qft_gather_tests {
    use super::*;
    use rand::SeedableRng;
    use rand_chacha::ChaCha8Rng;

    const TOL: f64 = 1e-12;

    fn random_state(bits: usize, seed: u64) -> Vec<Complex64> {
        let mut rng = ChaCha8Rng::seed_from_u64(seed);
        let mut state: Vec<Complex64> = (0..1usize << bits)
            .map(|_| Complex64::new(rng.random::<f64>() - 0.5, rng.random::<f64>() - 0.5))
            .collect();
        let norm = state.iter().map(|a| a.norm_sqr()).sum::<f64>().sqrt();
        state.iter_mut().for_each(|a| *a /= norm);
        state
    }

    fn max_diff(a: &[Complex64], b: &[Complex64]) -> f64 {
        a.iter()
            .zip(b)
            .map(|(x, y)| (x - y).norm())
            .fold(0.0, f64::max)
    }

    #[derive(Clone, Copy, Debug)]
    enum Route {
        Gathered,
        Paired,
        Textbook,
    }

    fn run(
        state: &[Complex64],
        num: usize,
        form: (bool, bool, bool),
        route: Route,
    ) -> Vec<Complex64> {
        let (inverse, swaps, big_endian) = form;
        let mut backend = StatevectorBackend::new(42);
        backend.num_qubits = state.len().ilog2() as usize;
        backend.state = state.to_vec();
        match route {
            Route::Gathered => backend.apply_qft_block_with(num, inverse, swaps, big_endian, true),
            Route::Paired => backend.apply_qft_block_with(num, inverse, swaps, big_endian, false),
            Route::Textbook => backend.apply_qft_block_textbook(0, num, inverse, swaps, big_endian),
        }
        backend.state
    }

    fn forms() -> impl Iterator<Item = (bool, bool, bool)> {
        (0..8).map(|f| (f & 1 != 0, f & 2 != 0, f & 4 != 0))
    }

    fn assert_routes_agree(state_bits: usize, num: usize, textbook: bool) {
        let state = random_state(state_bits, 42 + (state_bits * 64 + num) as u64);
        for form in forms() {
            let gathered = run(&state, num, form, Route::Gathered);
            let mut against = vec![Route::Paired];
            if textbook {
                against.push(Route::Textbook);
            }
            for route in against {
                let diff = max_diff(&gathered, &run(&state, num, form, route));
                assert!(
                    diff < TOL,
                    "{num} of {state_bits} qubits, (inverse, swaps, big_endian) {form:?}: \
                     gathered vs {route:?} differ by {diff:e}"
                );
            }
        }
    }

    // Reference DIF stages straight from the definition: stage `t` twiddles position
    // `k` of its block by `exp(2 pi i k / 2^(t + 1))`, read from the whole-state table.
    fn dif_stages(
        state: &mut [Complex64],
        n: usize,
        stages: std::ops::Range<usize>,
        inverse: bool,
    ) {
        let table = qft_twiddles_scaled(n, inverse);
        for chunk in state.chunks_mut(1 << n) {
            for t in stages.clone().rev() {
                let stride = 1usize << t;
                let step = (1usize << n) / (2 * stride);
                run_radix2_stage_seq(chunk, stride, &table, step, std::f64::consts::FRAC_1_SQRT_2);
            }
        }
    }

    #[test]
    fn gathered_passes_match_the_dif_definition_for_any_split() {
        // Splits cover odd and even pass widths, rows clamped to the offsets below a
        // narrow pass, and several 2^n chunks in one state.
        let cases: &[(usize, usize, &[usize])] = &[
            (10, 10, &[3, 3, 3]),
            (12, 12, &[1, 2, 4]),
            (12, 10, &[5, 4]),
            (14, 13, &[2, 7]),
            (16, 16, &[9, 6]),
            (17, 15, &[6]),
        ];
        for &(state_bits, n, split) in cases {
            for inverse in [false, true] {
                let mut expected = random_state(state_bits, 42 + state_bits as u64);
                let mut actual = expected.clone();
                let bottom = n - split.iter().sum::<usize>();
                dif_stages(&mut expected, n, bottom..n, inverse);
                let mut top = n;
                for &bits in split {
                    qft_gather_pass(&mut actual, QftGatherPass::new(top - bits, bits, inverse));
                    top -= bits;
                }
                let diff = max_diff(&actual, &expected);
                assert!(
                    diff < TOL,
                    "{n} of {state_bits} qubits, split {split:?}, inverse {inverse}: {diff:e}"
                );
            }
        }
    }

    // The textbook circuit costs one full pass per gate, so the widest registers check
    // against the paired route alone, which the narrower ones tie to the textbook.
    #[test]
    fn gathered_route_matches_paired_and_textbook_routes() {
        for n in QFT_GATHER_MIN_QUBITS..=22 {
            assert_routes_agree(n, n, n <= 20);
        }
    }

    #[test]
    fn gathered_route_matches_on_sub_register_chunks() {
        for (state_bits, num) in [(19, 16), (20, 18), (21, 19), (22, 20)] {
            assert_routes_agree(state_bits, num, state_bits <= 20);
        }
    }

    #[cfg(target_arch = "x86_64")]
    #[test]
    fn row_kernels_match_scalar_butterflies() {
        if !simd::has_avx2_fma() {
            return;
        }
        let inv_sqrt2 = std::f64::consts::FRAC_1_SQRT_2;
        let rows: Vec<Vec<Complex64>> = (0..4).map(|r| random_state(4, 7 + r)).collect();
        let w = [0.3, 1.7, -2.4].map(|angle| Complex64::from_polar(inv_sqrt2, angle));

        let (mut lo, mut hi) = (rows[0].clone(), rows[1].clone());
        // SAFETY: AVX2 detected above; the rows share one even length.
        unsafe { radix2_rows_avx2fma(&mut lo, &mut hi, w[0]) };
        for j in 0..lo.len() {
            let (elo, ehi) = radix2_butterfly_values(rows[0][j], rows[1][j], w[0], inv_sqrt2);
            assert!((lo[j] - elo).norm() < TOL && (hi[j] - ehi).norm() < TOL);
        }

        let mut q: [Vec<Complex64>; 4] = std::array::from_fn(|r| rows[r].clone());
        let [qa, qb, qc, qd] = q.each_mut().map(Vec::as_mut_slice);
        // SAFETY: AVX2 detected above; the rows share one even length.
        unsafe { radix4_rows_avx2fma(qa, qb, qc, qd, w[0], w[1], w[2]) };
        for j in 0..rows[0].len() {
            let expected = radix4_butterfly_values(
                rows[0][j], rows[1][j], rows[2][j], rows[3][j], w[0], w[1], w[2], inv_sqrt2,
            );
            let expected = [expected.0, expected.1, expected.2, expected.3];
            for (row, e) in q.iter().zip(expected) {
                assert!((row[j] - e).norm() < TOL);
            }
        }

        let phases = qft_gather_phases(16, true);
        let (anchor, step) = (phases.phase(12345), 37);
        let mut out = vec![Complex64::new(0.0, 0.0); rows[0].len()];
        // SAFETY: AVX2 detected above; `out` and the row share one even length.
        unsafe { qft_scatter_row_avx2fma(&mut out, &rows[0], anchor, &phases.low, step) };
        for (j, o) in out.iter().enumerate() {
            let expected = rows[0][j] * (anchor * phases.low[j * step]);
            assert!((o - expected).norm() < TOL);
        }
    }

    #[test]
    #[cfg_attr(debug_assertions, ignore = "256 MiB states, run with --release")]
    fn gathered_route_matches_paired_route_at_24_qubits() {
        assert_routes_agree(24, 24, false);
    }
}

#[cfg(test)]
mod multi_1q_tier_tests {
    use super::*;
    use crate::backend::Backend;
    use rand::SeedableRng;
    use rand_chacha::ChaCha8Rng;

    fn random_unitary(rng: &mut ChaCha8Rng) -> [[Complex64; 2]; 2] {
        let mut angle = || rng.random::<f64>() * std::f64::consts::TAU;
        let (theta, phi, lambda) = (angle() / 2.0, angle(), angle());
        let (c, s) = (theta.cos(), theta.sin());
        [
            [Complex64::new(c, 0.0), -Complex64::from_polar(s, lambda)],
            [
                Complex64::from_polar(s, phi),
                Complex64::from_polar(c, phi + lambda),
            ],
        ]
    }

    fn random_state(num_qubits: usize, rng: &mut ChaCha8Rng) -> Vec<Complex64> {
        (0..1usize << num_qubits)
            .map(|_| Complex64::new(rng.random::<f64>() - 0.5, rng.random::<f64>() - 0.5))
            .collect()
    }

    fn apply_scalar(state: &mut [Complex64], target: usize, mat: &[[Complex64; 2]; 2]) {
        let bit = 1usize << target;
        for i in (0..state.len()).filter(|i| i & bit == 0) {
            let (a, b) = (state[i], state[i | bit]);
            state[i] = mat[0][0] * a + mat[0][1] * b;
            state[i | bit] = mat[1][0] * a + mat[1][1] * b;
        }
    }

    fn assert_matches(got: &[Complex64], want: &[Complex64], label: &str) {
        for (i, (got, want)) in got.iter().zip(want).enumerate() {
            assert!(
                (got - want).norm() < 1e-12,
                "{label}, amplitude {i}: {got} vs {want}"
            );
        }
    }

    // Targets span the L2, L3 and untiled tiers in no particular order. The route
    // depends on the thread count, so every route must agree bitwise.
    #[cfg(feature = "parallel")]
    #[test]
    fn medium_tier_routes_match_per_gate_application() {
        const NUM_QUBITS: usize = 18;
        let mut rng = ChaCha8Rng::seed_from_u64(42);
        let start = random_state(NUM_QUBITS, &mut rng);
        let gates: Vec<TargetedGate1q> = [17, 0, 14, 5, 16, 13, 15]
            .into_iter()
            .map(|target| (target, random_unitary(&mut rng)))
            .collect();

        let mut expected = start.clone();
        for (target, mat) in &gates {
            apply_scalar(&mut expected, *target, mat);
        }
        let (gathered, _) = multi_1q_route(
            &gates,
            NUM_QUBITS,
            max_target_for_tile(multi_gate_l2_tile_par(start.len())),
            false,
            MAX_SHARED_1Q_TARGETS,
        );
        assert!(!matches!(gathered, Multi1qRoute::Tiered));
        let routes = [
            (true, Multi1qRoute::Tiered),
            (false, Multi1qRoute::Tiered),
            (false, gathered),
        ]
        .map(|(medium_tiled, route)| {
            let mut state = start.clone();
            apply_multi_1q_par_routed(&mut state, &gates, medium_tiled, &route);
            state
        });
        assert_matches(&routes[0], &expected, "L3 tier");
        for route in &routes[1..] {
            assert!(
                route.iter().zip(&routes[0]).all(|(a, b)| a == b),
                "the routes must agree bitwise, or results depend on the thread count"
            );
        }
    }

    // Every shape the router separates: targets in the tile only, a full layer, and
    // one to ten targets at the top of the register, alone and with low targets.
    // Under `parallel` the routed result must also equal the tiers bitwise.
    #[test]
    fn gathered_routes_match_per_gate_application() {
        let mut rng = ChaCha8Rng::seed_from_u64(42);
        for num_qubits in [14usize, 16, 18, 20] {
            let start = random_state(num_qubits, &mut rng);
            let mut cases: Vec<Vec<usize>> = vec![
                (0..num_qubits.min(12)).rev().collect(),
                (0..num_qubits).collect(),
            ];
            for k in 1..=10 {
                let high: Vec<usize> = (num_qubits - k..num_qubits).rev().collect();
                let mut mixed: Vec<usize> = [0, 3, 7, 11]
                    .into_iter()
                    .filter(|&q| q < num_qubits - k)
                    .collect();
                mixed.extend(&high);
                cases.push(high);
                cases.push(mixed);
            }
            for targets in cases {
                let gates: Vec<TargetedGate1q> = targets
                    .iter()
                    .map(|&target| (target, random_unitary(&mut rng)))
                    .collect();
                let mut expected = start.clone();
                for (target, mat) in &gates {
                    apply_scalar(&mut expected, *target, mat);
                }
                let mut backend = StatevectorBackend::new(42);
                backend.init(num_qubits, 0).unwrap();
                backend.state.copy_from_slice(&start);
                backend.apply_multi_1q(&gates);
                let label = format!("{num_qubits}q targets {targets:?}");
                assert_matches(&backend.state, &expected, &label);

                #[cfg(feature = "parallel")]
                if num_qubits >= PARALLEL_THRESHOLD_QUBITS {
                    let mut tiered = start.clone();
                    let medium_tiled = keep_l3_tier(tiered.len());
                    apply_multi_1q_par_routed(
                        &mut tiered,
                        &gates,
                        medium_tiled,
                        &Multi1qRoute::Tiered,
                    );
                    assert!(tiered == backend.state, "{label}: routes differ bitwise");
                }
            }
        }
    }
}