prism-q 0.34.0

//! Clifford propagation passes feeding the compiled sampler: backward
//! Heisenberg propagation of measured Z operators (bit-packed, batched) and
//! column-major forward stabilizer simulation.

use super::accumulator::transpose_64x64_dispatch;
use super::{PauliVec, flip_bit, get_bit, set_bit};
use crate::circuit::{Circuit, Instruction};
use crate::error::{PrismError, Result};
use crate::gates::Gate;

/// Heisenberg-picture backward propagation of a Pauli through a Clifford gate,
/// `P <- U†·P·U`, up to sign: the frame carries no phase, so `HYH = -Y` lands
/// as `Y`.
pub fn propagate_backward(pauli: &mut PauliVec, gate: &Gate, targets: &[usize]) {
    match gate {
        Gate::H => {
            let q = targets[0];
            let xb = get_bit(&pauli.x, q);
            let zb = get_bit(&pauli.z, q);
            set_bit(&mut pauli.x, q, zb);
            set_bit(&mut pauli.z, q, xb);
        }
        Gate::S => {
            let q = targets[0];
            if get_bit(&pauli.x, q) {
                flip_bit(&mut pauli.z, q);
            }
        }
        Gate::Sdg => {
            let q = targets[0];
            if get_bit(&pauli.x, q) {
                flip_bit(&mut pauli.z, q);
            }
        }
        Gate::X => {}
        Gate::Y => {}
        Gate::Z => {}
        Gate::SX => {
            let q = targets[0];
            if get_bit(&pauli.z, q) {
                flip_bit(&mut pauli.x, q);
            }
        }
        Gate::SXdg => {
            let q = targets[0];
            if get_bit(&pauli.z, q) {
                flip_bit(&mut pauli.x, q);
            }
        }
        Gate::Cx => {
            let ctrl = targets[0];
            let tgt = targets[1];
            if get_bit(&pauli.x, ctrl) {
                flip_bit(&mut pauli.x, tgt);
            }
            if get_bit(&pauli.z, tgt) {
                flip_bit(&mut pauli.z, ctrl);
            }
        }
        Gate::Cz => {
            let q0 = targets[0];
            let q1 = targets[1];
            let x0 = get_bit(&pauli.x, q0);
            let x1 = get_bit(&pauli.x, q1);
            if x1 {
                flip_bit(&mut pauli.z, q0);
            }
            if x0 {
                flip_bit(&mut pauli.z, q1);
            }
        }
        Gate::Swap => {
            let q0 = targets[0];
            let q1 = targets[1];
            let x0 = get_bit(&pauli.x, q0);
            let x1 = get_bit(&pauli.x, q1);
            set_bit(&mut pauli.x, q0, x1);
            set_bit(&mut pauli.x, q1, x0);
            let z0 = get_bit(&pauli.z, q0);
            let z1 = get_bit(&pauli.z, q1);
            set_bit(&mut pauli.z, q0, z1);
            set_bit(&mut pauli.z, q1, z0);
        }
        Gate::Id => {}
        _ => {}
    }
}

/// Batch backward-propagate Z_q through all gates for every measurement simultaneously.
///
/// Transposed representation: each qubit's m measurement x/z bits are packed into
/// m/64 u64 words, so gate propagation is bulk bitwise ops (64× fewer operations).
///
/// Returns (PauliVec, classical_bit, sign) per measurement.
pub(super) fn build_measurement_rows(circuit: &Circuit) -> Vec<(PauliVec, usize, bool)> {
    let n = circuit.num_qubits;
    let num_qubit_words = n.div_ceil(64);

    let measurements: Vec<(usize, usize)> = circuit
        .instructions
        .iter()
        .filter_map(|inst| match inst {
            Instruction::Measure {
                qubit,
                classical_bit,
            } => Some((*qubit, *classical_bit)),
            _ => None,
        })
        .collect();

    let m = measurements.len();
    if m == 0 {
        return Vec::new();
    }

    let m_words = m.div_ceil(64);

    let mut x_packed: Vec<Vec<u64>> = vec![vec![0u64; m_words]; n];
    let mut z_packed: Vec<Vec<u64>> = vec![vec![0u64; m_words]; n];
    let mut sign_packed: Vec<u64> = vec![0u64; m_words];

    for (meas_idx, &(qubit, _)) in measurements.iter().enumerate() {
        z_packed[qubit][meas_idx / 64] |= 1u64 << (meas_idx % 64);
    }

    let gate_instructions: Vec<(&Gate, &[usize])> = circuit
        .instructions
        .iter()
        .filter_map(|inst| match inst {
            Instruction::Gate { gate, targets } => Some((gate, targets.as_slice())),
            Instruction::Conditional { gate, targets, .. } => Some((gate, targets.as_slice())),
            _ => None,
        })
        .collect();

    for &(gate, targets) in gate_instructions.iter().rev() {
        batch_propagate_backward(
            &mut x_packed,
            &mut z_packed,
            &mut sign_packed,
            gate,
            targets,
            m_words,
        );
    }

    let mut rows: Vec<(PauliVec, usize, bool)> = Vec::with_capacity(m);
    for (meas_idx, &(_, classical_bit)) in measurements.iter().enumerate() {
        let mut pauli = PauliVec::new(num_qubit_words);
        for q in 0..n {
            if x_packed[q][meas_idx / 64] >> (meas_idx % 64) & 1 != 0 {
                pauli.x[q / 64] |= 1u64 << (q % 64);
            }
            if z_packed[q][meas_idx / 64] >> (meas_idx % 64) & 1 != 0 {
                pauli.z[q / 64] |= 1u64 << (q % 64);
            }
        }
        let sign = (sign_packed[meas_idx / 64] >> (meas_idx % 64)) & 1 != 0;
        rows.push((pauli, classical_bit, sign));
    }

    rows
}

#[cfg(target_arch = "x86_64")]
#[target_feature(enable = "avx2")]
unsafe fn batch_propagate_h_avx2(xq: &mut [u64], zq: &mut [u64], sign: &mut [u64], m_words: usize) {
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        use std::arch::x86_64::*;
        let chunks = m_words / 4;
        let xp = xq.as_mut_ptr() as *mut __m256i;
        let zp = zq.as_mut_ptr() as *mut __m256i;
        let sp = sign.as_mut_ptr() as *mut __m256i;
        for i in 0..chunks {
            let xv = _mm256_loadu_si256(xp.add(i));
            let zv = _mm256_loadu_si256(zp.add(i));
            let sv = _mm256_loadu_si256(sp.add(i));
            _mm256_storeu_si256(sp.add(i), _mm256_xor_si256(sv, _mm256_and_si256(xv, zv)));
            _mm256_storeu_si256(xp.add(i), zv);
            _mm256_storeu_si256(zp.add(i), xv);
        }
        let tail = chunks * 4;
        for w in tail..m_words {
            *sign.get_unchecked_mut(w) ^= *xq.get_unchecked(w) & *zq.get_unchecked(w);
        }
        for w in tail..m_words {
            let tmp = *xq.get_unchecked(w);
            *xq.get_unchecked_mut(w) = *zq.get_unchecked(w);
            *zq.get_unchecked_mut(w) = tmp;
        }
    }
}

#[cfg(target_arch = "x86_64")]
#[target_feature(enable = "avx2")]
unsafe fn batch_propagate_s_avx2(
    xq: &mut [u64],
    zq: &mut [u64],
    sign: &mut [u64],
    m_words: usize,
    negate_z: bool,
) {
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        use std::arch::x86_64::*;
        let chunks = m_words / 4;
        let xp = xq.as_ptr() as *const __m256i;
        let zp = zq.as_mut_ptr() as *mut __m256i;
        let sp = sign.as_mut_ptr() as *mut __m256i;
        if negate_z {
            for i in 0..chunks {
                let xv = _mm256_loadu_si256(xp.add(i));
                let zv = _mm256_loadu_si256(zp.add(i));
                let sv = _mm256_loadu_si256(sp.add(i));
                _mm256_storeu_si256(sp.add(i), _mm256_xor_si256(sv, _mm256_andnot_si256(zv, xv)));
                _mm256_storeu_si256(zp.add(i), _mm256_xor_si256(zv, xv));
            }
            let tail = chunks * 4;
            for w in tail..m_words {
                *sign.get_unchecked_mut(w) ^= *xq.get_unchecked(w) & !*zq.get_unchecked(w);
                *zq.get_unchecked_mut(w) ^= *xq.get_unchecked(w);
            }
        } else {
            for i in 0..chunks {
                let xv = _mm256_loadu_si256(xp.add(i));
                let zv = _mm256_loadu_si256(zp.add(i));
                let sv = _mm256_loadu_si256(sp.add(i));
                _mm256_storeu_si256(sp.add(i), _mm256_xor_si256(sv, _mm256_and_si256(xv, zv)));
                _mm256_storeu_si256(zp.add(i), _mm256_xor_si256(zv, xv));
            }
            let tail = chunks * 4;
            for w in tail..m_words {
                *sign.get_unchecked_mut(w) ^= *xq.get_unchecked(w) & *zq.get_unchecked(w);
                *zq.get_unchecked_mut(w) ^= *xq.get_unchecked(w);
            }
        }
    }
}

#[cfg(target_arch = "x86_64")]
#[target_feature(enable = "avx2")]
unsafe fn batch_propagate_sign_xor_avx2(dst: &mut [u64], src: &[u64], m_words: usize) {
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        use std::arch::x86_64::*;
        let chunks = m_words / 4;
        let dp = dst.as_mut_ptr() as *mut __m256i;
        let sp = src.as_ptr() as *const __m256i;
        for i in 0..chunks {
            let dv = _mm256_loadu_si256(dp.add(i));
            let sv = _mm256_loadu_si256(sp.add(i));
            _mm256_storeu_si256(dp.add(i), _mm256_xor_si256(dv, sv));
        }
        let tail = chunks * 4;
        for w in tail..m_words {
            *dst.get_unchecked_mut(w) ^= *src.get_unchecked(w);
        }
    }
}

#[cfg(target_arch = "x86_64")]
#[target_feature(enable = "avx2")]
unsafe fn batch_propagate_sign_xor2_avx2(dst: &mut [u64], a: &[u64], b: &[u64], m_words: usize) {
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        use std::arch::x86_64::*;
        let chunks = m_words / 4;
        let dp = dst.as_mut_ptr() as *mut __m256i;
        let ap = a.as_ptr() as *const __m256i;
        let bp = b.as_ptr() as *const __m256i;
        for i in 0..chunks {
            let dv = _mm256_loadu_si256(dp.add(i));
            let av = _mm256_loadu_si256(ap.add(i));
            let bv = _mm256_loadu_si256(bp.add(i));
            _mm256_storeu_si256(dp.add(i), _mm256_xor_si256(dv, _mm256_xor_si256(av, bv)));
        }
        let tail = chunks * 4;
        for w in tail..m_words {
            *dst.get_unchecked_mut(w) ^= *a.get_unchecked(w) ^ *b.get_unchecked(w);
        }
    }
}

#[cfg(target_arch = "x86_64")]
#[target_feature(enable = "avx2")]
unsafe fn batch_propagate_sx_avx2(
    xq: &mut [u64],
    zq: &[u64],
    sign: &mut [u64],
    m_words: usize,
    negate_x: bool,
) {
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        use std::arch::x86_64::*;
        let chunks = m_words / 4;
        let xp = xq.as_mut_ptr() as *mut __m256i;
        let zp = zq.as_ptr() as *const __m256i;
        let sp = sign.as_mut_ptr() as *mut __m256i;
        if negate_x {
            for i in 0..chunks {
                let xv = _mm256_loadu_si256(xp.add(i));
                let zv = _mm256_loadu_si256(zp.add(i));
                let sv = _mm256_loadu_si256(sp.add(i));
                _mm256_storeu_si256(sp.add(i), _mm256_xor_si256(sv, _mm256_andnot_si256(xv, zv)));
                _mm256_storeu_si256(xp.add(i), _mm256_xor_si256(xv, zv));
            }
            let tail = chunks * 4;
            for w in tail..m_words {
                *sign.get_unchecked_mut(w) ^= !*xq.get_unchecked(w) & *zq.get_unchecked(w);
                *xq.get_unchecked_mut(w) ^= *zq.get_unchecked(w);
            }
        } else {
            for i in 0..chunks {
                let xv = _mm256_loadu_si256(xp.add(i));
                let zv = _mm256_loadu_si256(zp.add(i));
                let sv = _mm256_loadu_si256(sp.add(i));
                _mm256_storeu_si256(sp.add(i), _mm256_xor_si256(sv, _mm256_and_si256(xv, zv)));
                _mm256_storeu_si256(xp.add(i), _mm256_xor_si256(xv, zv));
            }
            let tail = chunks * 4;
            for w in tail..m_words {
                *sign.get_unchecked_mut(w) ^= *xq.get_unchecked(w) & *zq.get_unchecked(w);
                *xq.get_unchecked_mut(w) ^= *zq.get_unchecked(w);
            }
        }
    }
}

#[cfg(target_arch = "x86_64")]
#[target_feature(enable = "avx2")]
unsafe fn batch_propagate_cx_avx2(
    x_ctrl: &[u64],
    z_ctrl: &mut [u64],
    x_tgt: &mut [u64],
    z_tgt: &[u64],
    sign: &mut [u64],
    m_words: usize,
) {
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        use std::arch::x86_64::*;
        let chunks = m_words / 4;
        let xcp = x_ctrl.as_ptr() as *const __m256i;
        let zcp = z_ctrl.as_mut_ptr() as *mut __m256i;
        let xtp = x_tgt.as_mut_ptr() as *mut __m256i;
        let ztp = z_tgt.as_ptr() as *const __m256i;
        let sp = sign.as_mut_ptr() as *mut __m256i;
        for i in 0..chunks {
            let xc = _mm256_loadu_si256(xcp.add(i));
            let zc = _mm256_loadu_si256(zcp.add(i));
            let xt = _mm256_loadu_si256(xtp.add(i));
            let zt = _mm256_loadu_si256(ztp.add(i));
            let sv = _mm256_loadu_si256(sp.add(i));
            let xnor = _mm256_andnot_si256(_mm256_xor_si256(zc, xt), _mm256_set1_epi64x(-1));
            let flip = _mm256_and_si256(_mm256_and_si256(xc, zt), xnor);
            _mm256_storeu_si256(sp.add(i), _mm256_xor_si256(sv, flip));
            _mm256_storeu_si256(xtp.add(i), _mm256_xor_si256(xt, xc));
            _mm256_storeu_si256(zcp.add(i), _mm256_xor_si256(zc, zt));
        }
        let tail = chunks * 4;
        for w in tail..m_words {
            let xc = *x_ctrl.get_unchecked(w);
            let zc = *z_ctrl.get_unchecked(w);
            let xt = *x_tgt.get_unchecked(w);
            let zt = *z_tgt.get_unchecked(w);
            *sign.get_unchecked_mut(w) ^= xc & zt & !(zc ^ xt);
            *x_tgt.get_unchecked_mut(w) = xt ^ xc;
            *z_ctrl.get_unchecked_mut(w) = zc ^ zt;
        }
    }
}

#[cfg(target_arch = "x86_64")]
#[target_feature(enable = "avx2")]
unsafe fn batch_propagate_cz_avx2(
    x0: &[u64],
    z0: &mut [u64],
    x1: &[u64],
    z1: &mut [u64],
    sign: &mut [u64],
    m_words: usize,
) {
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        use std::arch::x86_64::*;
        let chunks = m_words / 4;
        let x0p = x0.as_ptr() as *const __m256i;
        let z0p = z0.as_mut_ptr() as *mut __m256i;
        let x1p = x1.as_ptr() as *const __m256i;
        let z1p = z1.as_mut_ptr() as *mut __m256i;
        let sp = sign.as_mut_ptr() as *mut __m256i;
        for i in 0..chunks {
            let xv0 = _mm256_loadu_si256(x0p.add(i));
            let zv0 = _mm256_loadu_si256(z0p.add(i));
            let xv1 = _mm256_loadu_si256(x1p.add(i));
            let zv1 = _mm256_loadu_si256(z1p.add(i));
            let sv = _mm256_loadu_si256(sp.add(i));
            let zxor = _mm256_xor_si256(zv0, zv1);
            let flip = _mm256_and_si256(_mm256_and_si256(xv0, xv1), zxor);
            _mm256_storeu_si256(sp.add(i), _mm256_xor_si256(sv, flip));
            _mm256_storeu_si256(z0p.add(i), _mm256_xor_si256(zv0, xv1));
            _mm256_storeu_si256(z1p.add(i), _mm256_xor_si256(zv1, xv0));
        }
        let tail = chunks * 4;
        for w in tail..m_words {
            let xv0 = *x0.get_unchecked(w);
            let zv0 = *z0.get_unchecked(w);
            let xv1 = *x1.get_unchecked(w);
            let zv1 = *z1.get_unchecked(w);
            *sign.get_unchecked_mut(w) ^= xv0 & xv1 & (zv0 ^ zv1);
            *z0.get_unchecked_mut(w) = zv0 ^ xv1;
            *z1.get_unchecked_mut(w) = zv1 ^ xv0;
        }
    }
}

#[cfg(target_arch = "x86_64")]
#[target_feature(enable = "avx2")]
unsafe fn batch_propagate_swap_avx2(
    x0: &mut [u64],
    z0: &mut [u64],
    x1: &mut [u64],
    z1: &mut [u64],
    m_words: usize,
) {
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        use std::arch::x86_64::*;
        let chunks = m_words / 4;
        let x0p = x0.as_mut_ptr() as *mut __m256i;
        let z0p = z0.as_mut_ptr() as *mut __m256i;
        let x1p = x1.as_mut_ptr() as *mut __m256i;
        let z1p = z1.as_mut_ptr() as *mut __m256i;
        for i in 0..chunks {
            let xv0 = _mm256_loadu_si256(x0p.add(i));
            let xv1 = _mm256_loadu_si256(x1p.add(i));
            _mm256_storeu_si256(x0p.add(i), xv1);
            _mm256_storeu_si256(x1p.add(i), xv0);
            let zv0 = _mm256_loadu_si256(z0p.add(i));
            let zv1 = _mm256_loadu_si256(z1p.add(i));
            _mm256_storeu_si256(z0p.add(i), zv1);
            _mm256_storeu_si256(z1p.add(i), zv0);
        }
        let tail = chunks * 4;
        for w in tail..m_words {
            let tmp = *x0.get_unchecked(w);
            *x0.get_unchecked_mut(w) = *x1.get_unchecked(w);
            *x1.get_unchecked_mut(w) = tmp;
            let tmp = *z0.get_unchecked(w);
            *z0.get_unchecked_mut(w) = *z1.get_unchecked(w);
            *z1.get_unchecked_mut(w) = tmp;
        }
    }
}

#[cfg(target_arch = "aarch64")]
#[allow(dead_code)]
unsafe fn batch_propagate_h_neon(xq: &mut [u64], zq: &mut [u64], sign: &mut [u64], m_words: usize) {
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        use std::arch::aarch64::*;
        let chunks = m_words / 2;
        let xp = xq.as_mut_ptr();
        let zp = zq.as_mut_ptr();
        let sp = sign.as_mut_ptr();
        for i in 0..chunks {
            let off = i * 2;
            let xv = vld1q_u64(xp.add(off));
            let zv = vld1q_u64(zp.add(off));
            let sv = vld1q_u64(sp.add(off));
            vst1q_u64(sp.add(off), veorq_u64(sv, vandq_u64(xv, zv)));
            vst1q_u64(xp.add(off), zv);
            vst1q_u64(zp.add(off), xv);
        }
        let tail = chunks * 2;
        for w in tail..m_words {
            *sign.get_unchecked_mut(w) ^= *xq.get_unchecked(w) & *zq.get_unchecked(w);
        }
        for w in tail..m_words {
            let tmp = *xq.get_unchecked(w);
            *xq.get_unchecked_mut(w) = *zq.get_unchecked(w);
            *zq.get_unchecked_mut(w) = tmp;
        }
    }
}

#[cfg(target_arch = "aarch64")]
#[allow(dead_code)]
unsafe fn batch_propagate_s_neon(
    xq: &mut [u64],
    zq: &mut [u64],
    sign: &mut [u64],
    m_words: usize,
    negate_z: bool,
) {
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        use std::arch::aarch64::*;
        let chunks = m_words / 2;
        let xp = xq.as_ptr();
        let zp = zq.as_mut_ptr();
        let sp = sign.as_mut_ptr();
        if negate_z {
            for i in 0..chunks {
                let off = i * 2;
                let xv = vld1q_u64(xp.add(off));
                let zv = vld1q_u64(zp.add(off));
                let sv = vld1q_u64(sp.add(off));
                vst1q_u64(sp.add(off), veorq_u64(sv, vbicq_u64(xv, zv)));
                vst1q_u64(zp.add(off), veorq_u64(zv, xv));
            }
            let tail = chunks * 2;
            for w in tail..m_words {
                *sign.get_unchecked_mut(w) ^= *xq.get_unchecked(w) & !*zq.get_unchecked(w);
                *zq.get_unchecked_mut(w) ^= *xq.get_unchecked(w);
            }
        } else {
            for i in 0..chunks {
                let off = i * 2;
                let xv = vld1q_u64(xp.add(off));
                let zv = vld1q_u64(zp.add(off));
                let sv = vld1q_u64(sp.add(off));
                vst1q_u64(sp.add(off), veorq_u64(sv, vandq_u64(xv, zv)));
                vst1q_u64(zp.add(off), veorq_u64(zv, xv));
            }
            let tail = chunks * 2;
            for w in tail..m_words {
                *sign.get_unchecked_mut(w) ^= *xq.get_unchecked(w) & *zq.get_unchecked(w);
                *zq.get_unchecked_mut(w) ^= *xq.get_unchecked(w);
            }
        }
    }
}

#[cfg(target_arch = "aarch64")]
#[allow(dead_code)]
unsafe fn batch_propagate_sign_xor_neon(dst: &mut [u64], src: &[u64], m_words: usize) {
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        use std::arch::aarch64::*;
        let chunks = m_words / 2;
        let dp = dst.as_mut_ptr();
        let sp = src.as_ptr();
        for i in 0..chunks {
            let off = i * 2;
            let dv = vld1q_u64(dp.add(off));
            let sv = vld1q_u64(sp.add(off));
            vst1q_u64(dp.add(off), veorq_u64(dv, sv));
        }
        let tail = chunks * 2;
        for w in tail..m_words {
            *dst.get_unchecked_mut(w) ^= *src.get_unchecked(w);
        }
    }
}

#[cfg(target_arch = "aarch64")]
#[allow(dead_code)]
unsafe fn batch_propagate_sign_xor2_neon(dst: &mut [u64], a: &[u64], b: &[u64], m_words: usize) {
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        use std::arch::aarch64::*;
        let chunks = m_words / 2;
        let dp = dst.as_mut_ptr();
        let ap = a.as_ptr();
        let bp = b.as_ptr();
        for i in 0..chunks {
            let off = i * 2;
            let dv = vld1q_u64(dp.add(off));
            let av = vld1q_u64(ap.add(off));
            let bv = vld1q_u64(bp.add(off));
            vst1q_u64(dp.add(off), veorq_u64(dv, veorq_u64(av, bv)));
        }
        let tail = chunks * 2;
        for w in tail..m_words {
            *dst.get_unchecked_mut(w) ^= *a.get_unchecked(w) ^ *b.get_unchecked(w);
        }
    }
}

#[cfg(target_arch = "aarch64")]
#[allow(dead_code)]
unsafe fn batch_propagate_sx_neon(
    xq: &mut [u64],
    zq: &[u64],
    sign: &mut [u64],
    m_words: usize,
    negate_x: bool,
) {
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        use std::arch::aarch64::*;
        let chunks = m_words / 2;
        let xp = xq.as_mut_ptr();
        let zp = zq.as_ptr();
        let sp = sign.as_mut_ptr();
        if negate_x {
            for i in 0..chunks {
                let off = i * 2;
                let xv = vld1q_u64(xp.add(off));
                let zv = vld1q_u64(zp.add(off));
                let sv = vld1q_u64(sp.add(off));
                vst1q_u64(sp.add(off), veorq_u64(sv, vbicq_u64(zv, xv)));
                vst1q_u64(xp.add(off), veorq_u64(xv, zv));
            }
            let tail = chunks * 2;
            for w in tail..m_words {
                *sign.get_unchecked_mut(w) ^= !*xq.get_unchecked(w) & *zq.get_unchecked(w);
                *xq.get_unchecked_mut(w) ^= *zq.get_unchecked(w);
            }
        } else {
            for i in 0..chunks {
                let off = i * 2;
                let xv = vld1q_u64(xp.add(off));
                let zv = vld1q_u64(zp.add(off));
                let sv = vld1q_u64(sp.add(off));
                vst1q_u64(sp.add(off), veorq_u64(sv, vandq_u64(xv, zv)));
                vst1q_u64(xp.add(off), veorq_u64(xv, zv));
            }
            let tail = chunks * 2;
            for w in tail..m_words {
                *sign.get_unchecked_mut(w) ^= *xq.get_unchecked(w) & *zq.get_unchecked(w);
                *xq.get_unchecked_mut(w) ^= *zq.get_unchecked(w);
            }
        }
    }
}

#[cfg(target_arch = "aarch64")]
#[allow(dead_code)]
unsafe fn batch_propagate_cx_neon(
    x_ctrl: &[u64],
    z_ctrl: &mut [u64],
    x_tgt: &mut [u64],
    z_tgt: &[u64],
    sign: &mut [u64],
    m_words: usize,
) {
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        use std::arch::aarch64::*;
        let chunks = m_words / 2;
        let xcp = x_ctrl.as_ptr();
        let zcp = z_ctrl.as_mut_ptr();
        let xtp = x_tgt.as_mut_ptr();
        let ztp = z_tgt.as_ptr();
        let sp = sign.as_mut_ptr();
        let ones = vdupq_n_u64(!0u64);
        for i in 0..chunks {
            let off = i * 2;
            let xc = vld1q_u64(xcp.add(off));
            let zc = vld1q_u64(zcp.add(off));
            let xt = vld1q_u64(xtp.add(off));
            let zt = vld1q_u64(ztp.add(off));
            let sv = vld1q_u64(sp.add(off));
            let xnor = veorq_u64(veorq_u64(zc, xt), ones);
            let flip = vandq_u64(vandq_u64(xc, zt), xnor);
            vst1q_u64(sp.add(off), veorq_u64(sv, flip));
            vst1q_u64(xtp.add(off), veorq_u64(xt, xc));
            vst1q_u64(zcp.add(off), veorq_u64(zc, zt));
        }
        let tail = chunks * 2;
        for w in tail..m_words {
            let xc = *x_ctrl.get_unchecked(w);
            let zc = *z_ctrl.get_unchecked(w);
            let xt = *x_tgt.get_unchecked(w);
            let zt = *z_tgt.get_unchecked(w);
            *sign.get_unchecked_mut(w) ^= xc & zt & !(zc ^ xt);
            *x_tgt.get_unchecked_mut(w) = xt ^ xc;
            *z_ctrl.get_unchecked_mut(w) = zc ^ zt;
        }
    }
}

#[cfg(target_arch = "aarch64")]
#[allow(dead_code)]
unsafe fn batch_propagate_cz_neon(
    x0: &[u64],
    z0: &mut [u64],
    x1: &[u64],
    z1: &mut [u64],
    sign: &mut [u64],
    m_words: usize,
) {
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        use std::arch::aarch64::*;
        let chunks = m_words / 2;
        let x0p = x0.as_ptr();
        let z0p = z0.as_mut_ptr();
        let x1p = x1.as_ptr();
        let z1p = z1.as_mut_ptr();
        let sp = sign.as_mut_ptr();
        for i in 0..chunks {
            let off = i * 2;
            let xv0 = vld1q_u64(x0p.add(off));
            let zv0 = vld1q_u64(z0p.add(off));
            let xv1 = vld1q_u64(x1p.add(off));
            let zv1 = vld1q_u64(z1p.add(off));
            let sv = vld1q_u64(sp.add(off));
            let zxor = veorq_u64(zv0, zv1);
            let flip = vandq_u64(vandq_u64(xv0, xv1), zxor);
            vst1q_u64(sp.add(off), veorq_u64(sv, flip));
            vst1q_u64(z0p.add(off), veorq_u64(zv0, xv1));
            vst1q_u64(z1p.add(off), veorq_u64(zv1, xv0));
        }
        let tail = chunks * 2;
        for w in tail..m_words {
            let xv0 = *x0.get_unchecked(w);
            let zv0 = *z0.get_unchecked(w);
            let xv1 = *x1.get_unchecked(w);
            let zv1 = *z1.get_unchecked(w);
            *sign.get_unchecked_mut(w) ^= xv0 & xv1 & (zv0 ^ zv1);
            *z0.get_unchecked_mut(w) = zv0 ^ xv1;
            *z1.get_unchecked_mut(w) = zv1 ^ xv0;
        }
    }
}

#[cfg(target_arch = "aarch64")]
#[allow(dead_code)]
unsafe fn batch_propagate_swap_neon(
    x0: &mut [u64],
    z0: &mut [u64],
    x1: &mut [u64],
    z1: &mut [u64],
    m_words: usize,
) {
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        use std::arch::aarch64::*;
        let chunks = m_words / 2;
        let x0p = x0.as_mut_ptr();
        let z0p = z0.as_mut_ptr();
        let x1p = x1.as_mut_ptr();
        let z1p = z1.as_mut_ptr();
        for i in 0..chunks {
            let off = i * 2;
            let xv0 = vld1q_u64(x0p.add(off));
            let xv1 = vld1q_u64(x1p.add(off));
            vst1q_u64(x0p.add(off), xv1);
            vst1q_u64(x1p.add(off), xv0);
            let zv0 = vld1q_u64(z0p.add(off));
            let zv1 = vld1q_u64(z1p.add(off));
            vst1q_u64(z0p.add(off), zv1);
            vst1q_u64(z1p.add(off), zv0);
        }
        let tail = chunks * 2;
        for w in tail..m_words {
            let tmp = *x0.get_unchecked(w);
            *x0.get_unchecked_mut(w) = *x1.get_unchecked(w);
            *x1.get_unchecked_mut(w) = tmp;
            let tmp = *z0.get_unchecked(w);
            *z0.get_unchecked_mut(w) = *z1.get_unchecked(w);
            *z1.get_unchecked_mut(w) = tmp;
        }
    }
}

#[inline(always)]
pub(crate) fn batch_propagate_backward(
    x: &mut [Vec<u64>],
    z: &mut [Vec<u64>],
    sign: &mut [u64],
    gate: &Gate,
    targets: &[usize],
    m_words: usize,
) {
    match *targets {
        [q] => propagate_rows_1q(gate, &mut x[q], &mut z[q], sign, m_words),
        [a, b] => {
            let (x0, x1) = pair_mut(x, a, b);
            let (z0, z1) = pair_mut(z, a, b);
            propagate_rows_2q(gate, x0, z0, x1, z1, sign, m_words);
        }
        _ => {}
    }
}

/// [`batch_propagate_backward`] over rows packed end to end: slot `s` holds words
/// `s * words..(s + 1) * words` of `x` and `z`.
#[inline(always)]
pub(crate) fn batch_propagate_backward_flat(
    x: &mut [u64],
    z: &mut [u64],
    sign: &mut [u64],
    words: usize,
    gate: &Gate,
    slots: &[usize],
) {
    match *slots {
        [q] => propagate_rows_1q(
            gate,
            &mut x[q * words..(q + 1) * words],
            &mut z[q * words..(q + 1) * words],
            sign,
            words,
        ),
        [a, b] => {
            let (x0, x1) = pair_rows_mut(x, a, b, words);
            let (z0, z1) = pair_rows_mut(z, a, b, words);
            propagate_rows_2q(gate, x0, z0, x1, z1, sign, words);
        }
        _ => {}
    }
}

fn pair_mut<T>(v: &mut [T], a: usize, b: usize) -> (&mut T, &mut T) {
    debug_assert_ne!(a, b);
    if a < b {
        let (lo, hi) = v.split_at_mut(b);
        (&mut lo[a], &mut hi[0])
    } else {
        let (lo, hi) = v.split_at_mut(a);
        (&mut hi[0], &mut lo[b])
    }
}

/// Rows `a` and `b` of `v`, `words` words each, as two disjoint slices.
pub(crate) fn pair_rows_mut(
    v: &mut [u64],
    a: usize,
    b: usize,
    words: usize,
) -> (&mut [u64], &mut [u64]) {
    debug_assert_ne!(a, b);
    if a < b {
        let (lo, hi) = v.split_at_mut(b * words);
        (&mut lo[a * words..(a + 1) * words], &mut hi[..words])
    } else {
        let (lo, hi) = v.split_at_mut(a * words);
        (&mut hi[..words], &mut lo[b * words..(b + 1) * words])
    }
}

#[inline(always)]
fn propagate_rows_1q(
    gate: &Gate,
    xq: &mut [u64],
    zq: &mut [u64],
    sign: &mut [u64],
    m_words: usize,
) {
    #[cfg(target_arch = "x86_64")]
    let use_avx2 = m_words >= 4 && is_x86_feature_detected!("avx2");
    #[cfg(target_arch = "aarch64")]
    let use_neon = m_words >= 2;

    match gate {
        Gate::H => {
            #[cfg(target_arch = "x86_64")]
            if use_avx2 {
                // SAFETY: AVX2 detected, slices are valid with m_words elements
                unsafe { batch_propagate_h_avx2(xq, zq, sign, m_words) };
                return;
            }
            #[cfg(target_arch = "aarch64")]
            if use_neon {
                // SAFETY: NEON is baseline on aarch64, slices are valid with m_words elements
                unsafe { batch_propagate_h_neon(xq, zq, sign, m_words) };
                return;
            }
            for w in 0..m_words {
                sign[w] ^= xq[w] & zq[w];
            }
            xq[..m_words].swap_with_slice(&mut zq[..m_words]);
        }
        Gate::S => {
            #[cfg(target_arch = "x86_64")]
            if use_avx2 {
                // SAFETY: AVX2 detected, slices are valid with m_words elements
                unsafe { batch_propagate_s_avx2(xq, zq, sign, m_words, true) };
                return;
            }
            #[cfg(target_arch = "aarch64")]
            if use_neon {
                // SAFETY: NEON is baseline on aarch64, slices are valid with m_words elements
                unsafe { batch_propagate_s_neon(xq, zq, sign, m_words, true) };
                return;
            }
            for w in 0..m_words {
                sign[w] ^= xq[w] & !zq[w];
                zq[w] ^= xq[w];
            }
        }
        Gate::Sdg => {
            #[cfg(target_arch = "x86_64")]
            if use_avx2 {
                // SAFETY: AVX2 detected, slices are valid with m_words elements
                unsafe { batch_propagate_s_avx2(xq, zq, sign, m_words, false) };
                return;
            }
            #[cfg(target_arch = "aarch64")]
            if use_neon {
                // SAFETY: NEON is baseline on aarch64, slices are valid with m_words elements
                unsafe { batch_propagate_s_neon(xq, zq, sign, m_words, false) };
                return;
            }
            for w in 0..m_words {
                sign[w] ^= xq[w] & zq[w];
                zq[w] ^= xq[w];
            }
        }
        Gate::X => {
            #[cfg(target_arch = "x86_64")]
            if use_avx2 {
                // SAFETY: AVX2 detected, slices are valid with m_words elements
                unsafe { batch_propagate_sign_xor_avx2(sign, zq, m_words) };
                return;
            }
            #[cfg(target_arch = "aarch64")]
            if use_neon {
                // SAFETY: NEON is baseline on aarch64, slices are valid with m_words elements
                unsafe { batch_propagate_sign_xor_neon(sign, zq, m_words) };
                return;
            }
            for w in 0..m_words {
                sign[w] ^= zq[w];
            }
        }
        Gate::Y => {
            #[cfg(target_arch = "x86_64")]
            if use_avx2 {
                // SAFETY: AVX2 detected, slices are valid with m_words elements
                unsafe { batch_propagate_sign_xor2_avx2(sign, xq, zq, m_words) };
                return;
            }
            #[cfg(target_arch = "aarch64")]
            if use_neon {
                // SAFETY: NEON is baseline on aarch64, slices are valid with m_words elements
                unsafe { batch_propagate_sign_xor2_neon(sign, xq, zq, m_words) };
                return;
            }
            for w in 0..m_words {
                sign[w] ^= xq[w] ^ zq[w];
            }
        }
        Gate::Z => {
            #[cfg(target_arch = "x86_64")]
            if use_avx2 {
                // SAFETY: AVX2 detected, slices are valid with m_words elements
                unsafe { batch_propagate_sign_xor_avx2(sign, xq, m_words) };
                return;
            }
            #[cfg(target_arch = "aarch64")]
            if use_neon {
                // SAFETY: NEON is baseline on aarch64, slices are valid with m_words elements
                unsafe { batch_propagate_sign_xor_neon(sign, xq, m_words) };
                return;
            }
            for w in 0..m_words {
                sign[w] ^= xq[w];
            }
        }
        Gate::SX => {
            #[cfg(target_arch = "x86_64")]
            if use_avx2 {
                // SAFETY: AVX2 detected, slices are valid with m_words elements
                unsafe { batch_propagate_sx_avx2(xq, zq, sign, m_words, false) };
                return;
            }
            #[cfg(target_arch = "aarch64")]
            if use_neon {
                // SAFETY: NEON is baseline on aarch64, slices are valid with m_words elements
                unsafe { batch_propagate_sx_neon(xq, zq, sign, m_words, false) };
                return;
            }
            for w in 0..m_words {
                sign[w] ^= xq[w] & zq[w];
                xq[w] ^= zq[w];
            }
        }
        Gate::SXdg => {
            #[cfg(target_arch = "x86_64")]
            if use_avx2 {
                // SAFETY: AVX2 detected, slices are valid with m_words elements
                unsafe { batch_propagate_sx_avx2(xq, zq, sign, m_words, true) };
                return;
            }
            #[cfg(target_arch = "aarch64")]
            if use_neon {
                // SAFETY: NEON is baseline on aarch64, slices are valid with m_words elements
                unsafe { batch_propagate_sx_neon(xq, zq, sign, m_words, true) };
                return;
            }
            for w in 0..m_words {
                sign[w] ^= !xq[w] & zq[w];
                xq[w] ^= zq[w];
            }
        }
        _ => {}
    }
}

#[inline(always)]
fn propagate_rows_2q(
    gate: &Gate,
    x0: &mut [u64],
    z0: &mut [u64],
    x1: &mut [u64],
    z1: &mut [u64],
    sign: &mut [u64],
    m_words: usize,
) {
    #[cfg(target_arch = "x86_64")]
    let use_avx2 = m_words >= 4 && is_x86_feature_detected!("avx2");
    #[cfg(target_arch = "aarch64")]
    let use_neon = m_words >= 2;

    match gate {
        Gate::Cx => {
            #[cfg(target_arch = "x86_64")]
            if use_avx2 {
                // SAFETY: AVX2 detected, slices are valid, ctrl != tgt
                unsafe { batch_propagate_cx_avx2(x0, z0, x1, z1, sign, m_words) };
                return;
            }
            #[cfg(target_arch = "aarch64")]
            if use_neon {
                // SAFETY: NEON is baseline on aarch64, slices are valid, ctrl != tgt
                unsafe { batch_propagate_cx_neon(x0, z0, x1, z1, sign, m_words) };
                return;
            }
            for w in 0..m_words {
                sign[w] ^= x0[w] & z1[w] & !(z0[w] ^ x1[w]);
                x1[w] ^= x0[w];
                z0[w] ^= z1[w];
            }
        }
        Gate::Cz => {
            #[cfg(target_arch = "x86_64")]
            if use_avx2 {
                // SAFETY: AVX2 detected, slices are valid, q0 != q1
                unsafe { batch_propagate_cz_avx2(x0, z0, x1, z1, sign, m_words) };
                return;
            }
            #[cfg(target_arch = "aarch64")]
            if use_neon {
                // SAFETY: NEON is baseline on aarch64, slices are valid, q0 != q1
                unsafe { batch_propagate_cz_neon(x0, z0, x1, z1, sign, m_words) };
                return;
            }
            for w in 0..m_words {
                sign[w] ^= x0[w] & x1[w] & (z0[w] ^ z1[w]);
                z0[w] ^= x1[w];
                z1[w] ^= x0[w];
            }
        }
        Gate::Swap => {
            #[cfg(target_arch = "x86_64")]
            if use_avx2 {
                // SAFETY: AVX2 detected, slices are valid, q0 != q1
                unsafe { batch_propagate_swap_avx2(x0, z0, x1, z1, m_words) };
                return;
            }
            #[cfg(target_arch = "aarch64")]
            if use_neon {
                // SAFETY: NEON is baseline on aarch64, slices are valid, q0 != q1
                unsafe { batch_propagate_swap_neon(x0, z0, x1, z1, m_words) };
                return;
            }
            x0[..m_words].swap_with_slice(&mut x1[..m_words]);
            z0[..m_words].swap_with_slice(&mut z1[..m_words]);
        }
        _ => {}
    }
}

/// AVX2 rowmul: src (read-only) multiplied into dst (modified in-place).
/// Both buffers are laid out as [x_words(nw) | z_words(nw)].
/// Returns the accumulated Pauli phase sum.
#[cfg(target_arch = "x86_64")]
#[target_feature(enable = "avx2")]
unsafe fn rowmul_avx2(src_ptr: *const u64, dst_ptr: *mut u64, nw: usize, initial_sum: u64) -> u64 {
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        use std::arch::x86_64::*;
        let chunks = nw / 4;
        let mut sum = initial_sum;
        let src_x = src_ptr;
        let src_z = src_ptr.add(nw);
        let dst_x = dst_ptr;
        let dst_z = dst_ptr.add(nw);

        for i in 0..chunks {
            let off = i * 4;
            let x1 = _mm256_loadu_si256(src_x.add(off) as *const __m256i);
            let z1 = _mm256_loadu_si256(src_z.add(off) as *const __m256i);
            let x2 = _mm256_loadu_si256(dst_x.add(off) as *const __m256i);
            let z2 = _mm256_loadu_si256(dst_z.add(off) as *const __m256i);

            let new_x = _mm256_xor_si256(x1, x2);
            let new_z = _mm256_xor_si256(z1, z2);
            _mm256_storeu_si256(dst_x.add(off) as *mut __m256i, new_x);
            _mm256_storeu_si256(dst_z.add(off) as *mut __m256i, new_z);

            let any = _mm256_or_si256(_mm256_or_si256(x1, z1), _mm256_or_si256(x2, z2));
            if _mm256_testz_si256(any, any) == 0 {
                let x1z1 = _mm256_and_si256(x1, z1);
                let x2z2 = _mm256_and_si256(x2, z2);

                let nonzero = _mm256_and_si256(
                    _mm256_and_si256(_mm256_or_si256(new_x, new_z), _mm256_or_si256(x1, z1)),
                    _mm256_or_si256(x2, z2),
                );
                // pos = (x1&z1&!x2&z2) | (x1&!z1&x2&z2) | (!x1&z1&x2&!z2)
                let pos = _mm256_or_si256(
                    _mm256_or_si256(
                        _mm256_and_si256(x1z1, _mm256_andnot_si256(x2, z2)),
                        _mm256_and_si256(_mm256_andnot_si256(z1, x1), x2z2),
                    ),
                    _mm256_and_si256(_mm256_andnot_si256(x1, z1), _mm256_andnot_si256(z2, x2)),
                );

                let nz_arr: [u64; 4] = std::mem::transmute(nonzero);
                let pos_arr: [u64; 4] = std::mem::transmute(pos);
                let mut nz_count = 0u64;
                let mut pos_count = 0u64;
                for j in 0..4 {
                    nz_count += nz_arr[j].count_ones() as u64;
                    pos_count += pos_arr[j].count_ones() as u64;
                }
                sum = sum.wrapping_add(2 * pos_count);
                sum = sum.wrapping_sub(nz_count);
            }
        }

        let tail = chunks * 4;
        for w in tail..nw {
            let x1 = *src_x.add(w);
            let z1 = *src_z.add(w);
            let x2 = *dst_x.add(w);
            let z2 = *dst_z.add(w);
            let new_x = x1 ^ x2;
            let new_z = z1 ^ z2;
            *dst_x.add(w) = new_x;
            *dst_z.add(w) = new_z;
            if (x1 | z1 | x2 | z2) != 0 {
                let nonzero = (new_x | new_z) & (x1 | z1) & (x2 | z2);
                let pos = (x1 & z1 & !x2 & z2) | (x1 & !z1 & x2 & z2) | (!x1 & z1 & x2 & !z2);
                sum = sum.wrapping_add(2 * pos.count_ones() as u64);
                sum = sum.wrapping_sub(nonzero.count_ones() as u64);
            }
        }
        sum
    }
}

#[cfg(target_arch = "aarch64")]
#[allow(dead_code)]
unsafe fn rowmul_neon(src_ptr: *const u64, dst_ptr: *mut u64, nw: usize, initial_sum: u64) -> u64 {
    // SAFETY: same contract as the enclosing unsafe fn.
    unsafe {
        use std::arch::aarch64::*;
        let chunks = nw / 2;
        let mut sum = initial_sum;
        let src_x = src_ptr;
        let src_z = src_ptr.add(nw);
        let dst_x = dst_ptr;
        let dst_z = dst_ptr.add(nw);

        for i in 0..chunks {
            let off = i * 2;
            let x1 = vld1q_u64(src_x.add(off));
            let z1 = vld1q_u64(src_z.add(off));
            let x2 = vld1q_u64(dst_x.add(off));
            let z2 = vld1q_u64(dst_z.add(off));

            let new_x = veorq_u64(x1, x2);
            let new_z = veorq_u64(z1, z2);
            vst1q_u64(dst_x.add(off), new_x);
            vst1q_u64(dst_z.add(off), new_z);

            let any = vorrq_u64(vorrq_u64(x1, z1), vorrq_u64(x2, z2));
            let any_arr: [u64; 2] = std::mem::transmute(any);
            if (any_arr[0] | any_arr[1]) != 0 {
                let x1z1 = vandq_u64(x1, z1);
                let x2z2 = vandq_u64(x2, z2);

                let nonzero = vandq_u64(
                    vandq_u64(vorrq_u64(new_x, new_z), vorrq_u64(x1, z1)),
                    vorrq_u64(x2, z2),
                );
                // pos = (x1&z1&!x2&z2) | (x1&!z1&x2&z2) | (!x1&z1&x2&!z2)
                let pos = vorrq_u64(
                    vorrq_u64(
                        vandq_u64(x1z1, vbicq_u64(z2, x2)),
                        vandq_u64(vbicq_u64(x1, z1), x2z2),
                    ),
                    vandq_u64(vbicq_u64(z1, x1), vbicq_u64(x2, z2)),
                );

                let nz_arr: [u64; 2] = std::mem::transmute(nonzero);
                let pos_arr: [u64; 2] = std::mem::transmute(pos);
                let mut nz_count = 0u64;
                let mut pos_count = 0u64;
                for j in 0..2 {
                    nz_count += nz_arr[j].count_ones() as u64;
                    pos_count += pos_arr[j].count_ones() as u64;
                }
                sum = sum.wrapping_add(2 * pos_count);
                sum = sum.wrapping_sub(nz_count);
            }
        }

        let tail = chunks * 2;
        for w in tail..nw {
            let x1 = *src_x.add(w);
            let z1 = *src_z.add(w);
            let x2 = *dst_x.add(w);
            let z2 = *dst_z.add(w);
            let new_x = x1 ^ x2;
            let new_z = z1 ^ z2;
            *dst_x.add(w) = new_x;
            *dst_z.add(w) = new_z;
            if (x1 | z1 | x2 | z2) != 0 {
                let nonzero = (new_x | new_z) & (x1 | z1) & (x2 | z2);
                let pos = (x1 & z1 & !x2 & z2) | (x1 & !z1 & x2 & z2) | (!x1 & z1 & x2 & !z2);
                sum = sum.wrapping_add(2 * pos.count_ones() as u64);
                sum = sum.wrapping_sub(nonzero.count_ones() as u64);
            }
        }
        sum
    }
}

/// Rowmul: multiply src into xz[r_base..], returning the resulting phase.
#[inline(always)]
pub(super) fn rowmul_phase(
    src: &[u64],
    xz: &mut [u64],
    r_base: usize,
    nw: usize,
    src_phase: bool,
    dst_phase: bool,
) -> bool {
    let initial_sum = if src_phase { 2u64 } else { 0 } + if dst_phase { 2u64 } else { 0 };

    #[cfg(target_arch = "x86_64")]
    if nw >= 4 && is_x86_feature_detected!("avx2") {
        // SAFETY: AVX2 detected, src has 2*nw elements, xz[r_base..] has 2*nw elements
        let sum =
            unsafe { rowmul_avx2(src.as_ptr(), xz.as_mut_ptr().add(r_base), nw, initial_sum) };
        return (sum & 3) >= 2;
    }

    #[cfg(target_arch = "aarch64")]
    if nw >= 2 {
        // SAFETY: NEON is baseline on aarch64, src has 2*nw elements, xz[r_base..] has 2*nw elements
        let sum =
            unsafe { rowmul_neon(src.as_ptr(), xz.as_mut_ptr().add(r_base), nw, initial_sum) };
        return (sum & 3) >= 2;
    }

    let mut sum = initial_sum;
    for w in 0..nw {
        let x1 = src[w];
        let z1 = src[nw + w];
        let x2 = xz[r_base + w];
        let z2 = xz[r_base + nw + w];
        let new_x = x1 ^ x2;
        let new_z = z1 ^ z2;
        xz[r_base + w] = new_x;
        xz[r_base + nw + w] = new_z;
        if (x1 | z1 | x2 | z2) != 0 {
            let nonzero = (new_x | new_z) & (x1 | z1) & (x2 | z2);
            let pos = (x1 & z1 & !x2 & z2) | (x1 & !z1 & x2 & z2) | (!x1 & z1 & x2 & !z2);
            sum = sum.wrapping_add(2 * pos.count_ones() as u64);
            sum = sum.wrapping_sub(nonzero.count_ones() as u64);
        }
    }
    (sum & 3) >= 2
}

/// Rowmul variant for scratch buffer: multiply src (from xz) into scratch.
#[inline(always)]
pub(super) fn rowmul_phase_into(
    xz: &[u64],
    s_base: usize,
    scratch: &mut [u64],
    nw: usize,
    src_phase: bool,
    dst_phase: bool,
) -> bool {
    let initial_sum = if src_phase { 2u64 } else { 0 } + if dst_phase { 2u64 } else { 0 };

    #[cfg(target_arch = "x86_64")]
    if nw >= 4 && is_x86_feature_detected!("avx2") {
        // SAFETY: AVX2 detected, xz[s_base..] is src, scratch is dst
        let sum = unsafe {
            rowmul_avx2(
                xz.as_ptr().add(s_base),
                scratch.as_mut_ptr(),
                nw,
                initial_sum,
            )
        };
        return (sum & 3) >= 2;
    }

    #[cfg(target_arch = "aarch64")]
    if nw >= 2 {
        // SAFETY: NEON is baseline on aarch64, xz[s_base..] is src, scratch is dst
        let sum = unsafe {
            rowmul_neon(
                xz.as_ptr().add(s_base),
                scratch.as_mut_ptr(),
                nw,
                initial_sum,
            )
        };
        return (sum & 3) >= 2;
    }

    let mut sum = initial_sum;
    for w in 0..nw {
        let x1 = xz[s_base + w];
        let z1 = xz[s_base + nw + w];
        let x2 = scratch[w];
        let z2 = scratch[nw + w];
        let new_x = x1 ^ x2;
        let new_z = z1 ^ z2;
        scratch[w] = new_x;
        scratch[nw + w] = new_z;
        if (x1 | z1 | x2 | z2) != 0 {
            let nonzero = (new_x | new_z) & (x1 | z1) & (x2 | z2);
            let pos = (x1 & z1 & !x2 & z2) | (x1 & !z1 & x2 & z2) | (!x1 & z1 & x2 & !z2);
            sum = sum.wrapping_add(2 * pos.count_ones() as u64);
            sum = sum.wrapping_sub(nonzero.count_ones() as u64);
        }
    }
    (sum & 3) >= 2
}

/// Compute reference measurement outcomes by simulating measurements of the
/// propagated Paulis on |0⟩^n using the Aaronson-Gottesman measurement protocol.
/// O(m x n^2/64), avoids re-simulating all T circuit gates through the stabilizer.
#[allow(clippy::needless_range_loop)]
pub(super) fn compute_reference_bits(
    measurement_rows: &[(PauliVec, usize, bool)],
    n: usize,
) -> Vec<bool> {
    let m = measurement_rows.len();
    let nw = n.div_ceil(64);
    let stride = 2 * nw;

    let mut xz = vec![0u64; 2 * n * stride];
    let mut phase = vec![false; 2 * n];

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

    let mut ref_bits = vec![false; m];
    let mut p_data = vec![0u64; stride];
    let mut scratch = vec![0u64; stride];

    for (meas_idx, (pauli, _, _)) in measurement_rows.iter().enumerate() {
        let mut anti_idx = None;
        for g in n..2 * n {
            let base = g * stride;
            let mut inner = 0u64;
            for w in 0..nw {
                inner ^= (pauli.x[w] & xz[base + nw + w]) ^ (pauli.z[w] & xz[base + w]);
            }
            if inner.count_ones() % 2 == 1 {
                anti_idx = Some(g);
                break;
            }
        }

        match anti_idx {
            Some(p) => {
                p_data.copy_from_slice(&xz[p * stride..(p + 1) * stride]);
                let p_phase = phase[p];

                for r in 0..2 * n {
                    if r == p {
                        continue;
                    }
                    let r_base = r * stride;
                    let mut inner = 0u64;
                    for w in 0..nw {
                        inner ^= (pauli.x[w] & xz[r_base + nw + w]) ^ (pauli.z[w] & xz[r_base + w]);
                    }
                    if inner.count_ones() % 2 == 1 {
                        phase[r] = rowmul_phase(&p_data, &mut xz, r_base, nw, p_phase, phase[r]);
                    }
                }

                let dest_idx = p - n;
                let dest_base = dest_idx * stride;
                xz.copy_within(p * stride..(p + 1) * stride, dest_base);
                phase[dest_idx] = p_phase;

                let p_base = p * stride;
                for w in 0..nw {
                    xz[p_base + w] = pauli.x[w];
                    xz[p_base + nw + w] = pauli.z[w];
                }
                phase[p] = false;

                ref_bits[meas_idx] = false;
            }
            None => {
                scratch.fill(0);
                let mut scratch_phase = false;

                for g in 0..n {
                    let d_base = g * stride;
                    let mut inner = 0u64;
                    for w in 0..nw {
                        inner ^= (pauli.x[w] & xz[d_base + nw + w]) ^ (pauli.z[w] & xz[d_base + w]);
                    }
                    if inner.count_ones() % 2 == 1 {
                        let s_base = (g + n) * stride;
                        let s_phase = phase[g + n];
                        scratch_phase = rowmul_phase_into(
                            &xz,
                            s_base,
                            &mut scratch,
                            nw,
                            s_phase,
                            scratch_phase,
                        );
                    }
                }

                ref_bits[meas_idx] = scratch_phase;
            }
        }
    }

    ref_bits
}

/// Forward conjugation `G P G†` is the backward step by `G†`.
pub(super) fn forward_as_backward_gate(gate: &Gate) -> Result<&Gate> {
    match gate {
        Gate::S => Ok(&Gate::Sdg),
        Gate::Sdg => Ok(&Gate::S),
        Gate::SX => Ok(&Gate::SXdg),
        Gate::SXdg => Ok(&Gate::SX),
        Gate::H | Gate::X | Gate::Y | Gate::Z | Gate::Id | Gate::Cx | Gate::Cz | Gate::Swap => {
            Ok(gate)
        }
        _ => Err(PrismError::IncompatibleBackend {
            backend: "CompiledSampler".to_string(),
            reason: format!("unsupported gate {:?} in column-major forward sim", gate),
        }),
    }
}

/// Column-major forward stabilizer simulation.
///
/// Stores x_cols[qubit][row_word] and z_cols[qubit][row_word] so gate
/// application is O(row_words) vector XOR instead of O(2n) per-row bit ops.
///
/// Returns (xz_rowmajor, phase, nw) transposed back to row-major.
pub(super) fn colmajor_forward_sim(
    n: usize,
    instructions: &[Instruction],
) -> Result<(Vec<u64>, Vec<bool>, usize)> {
    let total_rows = 2 * n;
    let nw = n.div_ceil(64);
    let row_words = total_rows.div_ceil(64);

    let mut x_cols: Vec<Vec<u64>> = vec![vec![0u64; row_words]; n];
    let mut z_cols: Vec<Vec<u64>> = vec![vec![0u64; row_words]; n];
    let mut phase: Vec<u64> = vec![0u64; row_words];

    for q in 0..n {
        x_cols[q][q / 64] |= 1u64 << (q % 64);
        let stab_row = n + q;
        z_cols[q][stab_row / 64] |= 1u64 << (stab_row % 64);
    }

    for inst in instructions {
        let (gate, targets) = match inst {
            Instruction::Gate { gate, targets } => (gate, targets.as_slice()),
            Instruction::Conditional { gate, targets, .. } => (gate, targets.as_slice()),
            _ => continue,
        };

        batch_propagate_backward(
            &mut x_cols,
            &mut z_cols,
            &mut phase,
            forward_as_backward_gate(gate)?,
            targets,
            row_words,
        );
    }

    let stride = 2 * nw;
    let mut xz = vec![0u64; total_rows * stride];
    let mut phase_vec = vec![false; total_rows];

    let mut block = [0u64; 64];
    for (cols, offset) in [(&x_cols, 0), (&z_cols, nw)] {
        for qw in 0..nw {
            let qubits = &cols[qw * 64..n.min(qw * 64 + 64)];
            for rw in 0..row_words {
                let mut any = 0u64;
                for (b, col) in qubits.iter().enumerate() {
                    block[b] = col[rw];
                    any |= col[rw];
                }
                if any == 0 {
                    continue;
                }
                block[qubits.len()..].fill(0);
                transpose_64x64_dispatch(&mut block);
                let rows = (total_rows - rw * 64).min(64);
                for (i, &word) in block[..rows].iter().enumerate() {
                    xz[(rw * 64 + i) * stride + offset + qw] = word;
                }
            }
        }
    }

    for (rw, &pw) in phase.iter().enumerate().take(row_words) {
        let mut pbits = pw;
        while pbits != 0 {
            let bit = pbits.trailing_zeros() as usize;
            let row = rw * 64 + bit;
            if row < total_rows {
                phase_vec[row] = true;
            }
            pbits &= pbits - 1;
        }
    }

    Ok((xz, phase_vec, nw))
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::CircuitBuilder;
    use crate::sim::compiled::PauliVec;

    fn pauli(n: usize, x_bits: &[usize], z_bits: &[usize]) -> PauliVec {
        let mut p = PauliVec::new(n.div_ceil(64));
        for &q in x_bits {
            p.x[q / 64] |= 1u64 << (q % 64);
        }
        for &q in z_bits {
            p.z[q / 64] |= 1u64 << (q % 64);
        }
        p
    }

    #[test]
    fn propagate_h_swaps_x_z() {
        let mut p = pauli(4, &[1], &[]);
        propagate_backward(&mut p, &Gate::H, &[1]);
        assert_eq!(p.x[0] & 0b10, 0);
        assert_eq!(p.z[0] & 0b10, 0b10);
    }

    #[test]
    fn propagate_s_and_sdg_add_z_when_x_set() {
        let mut p = pauli(2, &[0], &[]);
        propagate_backward(&mut p, &Gate::S, &[0]);
        assert_eq!(p.z[0] & 1, 1);
        let mut q = pauli(2, &[0], &[]);
        propagate_backward(&mut q, &Gate::Sdg, &[0]);
        assert_eq!(q.z[0] & 1, 1);
    }

    #[test]
    fn propagate_sx_and_sxdg_add_x_when_z_set() {
        let mut p = pauli(2, &[], &[0]);
        propagate_backward(&mut p, &Gate::SX, &[0]);
        assert_eq!(p.x[0] & 1, 1);
        let mut q = pauli(2, &[], &[0]);
        propagate_backward(&mut q, &Gate::SXdg, &[0]);
        assert_eq!(q.x[0] & 1, 1);
    }

    #[test]
    fn propagate_pauli_gates_noop() {
        let p0 = pauli(2, &[0], &[1]);
        for gate in [Gate::X, Gate::Y, Gate::Z, Gate::Id] {
            let mut p = p0.clone();
            propagate_backward(&mut p, &gate, &[0]);
            assert_eq!(p.x, p0.x);
            assert_eq!(p.z, p0.z);
        }
    }

    #[test]
    fn propagate_cx_flips_target_x_and_control_z() {
        let mut p = pauli(2, &[0], &[1]);
        propagate_backward(&mut p, &Gate::Cx, &[0, 1]);
        assert_eq!(p.x[0] & 0b11, 0b11);
        assert_eq!(p.z[0] & 0b11, 0b11);
    }

    #[test]
    fn propagate_cz_adds_z_on_partner_when_x_set() {
        let mut p = pauli(2, &[0], &[]);
        propagate_backward(&mut p, &Gate::Cz, &[0, 1]);
        assert_eq!(p.z[0] & 0b10, 0b10);
        assert_eq!(p.x[0] & 0b01, 0b01);
    }

    #[test]
    fn propagate_swap_exchanges_pauli_bits() {
        let mut p = pauli(2, &[0], &[1]);
        propagate_backward(&mut p, &Gate::Swap, &[0, 1]);
        assert_eq!(p.x[0] & 0b11, 0b10);
        assert_eq!(p.z[0] & 0b11, 0b01);
    }

    #[test]
    fn propagate_unsupported_gate_is_noop() {
        let p0 = pauli(1, &[], &[0]);
        let mut p = p0.clone();
        propagate_backward(&mut p, &Gate::T, &[0]);
        assert_eq!(p.x, p0.x);
        assert_eq!(p.z, p0.z);
    }

    #[test]
    fn build_measurement_rows_empty_when_no_measurements() {
        let circuit = CircuitBuilder::new(1).h(0).build();
        let rows = build_measurement_rows(&circuit);
        assert!(rows.is_empty());
    }

    #[test]
    fn build_measurement_rows_z_for_unentangled_meas() {
        let circuit = CircuitBuilder::new_with_classical(2, 2)
            .measure(0, 0)
            .measure(1, 1)
            .build();
        let rows = build_measurement_rows(&circuit);
        assert_eq!(rows.len(), 2);
        assert_eq!(rows[0].1, 0);
        assert_eq!(rows[1].1, 1);
        assert_eq!(rows[0].0.x[0], 0);
        assert_eq!(rows[1].0.x[0], 0);
        assert_eq!(rows[0].0.z[0] & 1, 1);
        assert_eq!(rows[1].0.z[0] & 0b10, 0b10);
    }

    #[test]
    fn build_measurement_rows_h_swaps_to_x() {
        let circuit = CircuitBuilder::new_with_classical(1, 1)
            .h(0)
            .measure(0, 0)
            .build();
        let rows = build_measurement_rows(&circuit);
        assert_eq!(rows.len(), 1);
        assert_eq!(rows[0].0.x[0] & 1, 1);
        assert_eq!(rows[0].0.z[0] & 1, 0);
    }

    #[test]
    fn propagate_cz_flips_q0_when_x1_set() {
        let mut p = pauli(2, &[1], &[]);
        propagate_backward(&mut p, &Gate::Cz, &[0, 1]);
        assert_eq!(p.z[0] & 0b01, 0b01);
        assert_eq!(p.x[0] & 0b10, 0b10);
    }

    #[test]
    fn build_measurement_rows_includes_conditional_gates() {
        use crate::circuit::{ClassicalCondition, Instruction, SmallVec};
        use smallvec::smallvec;
        let mut c = crate::circuit::Circuit::new(1, 1);
        c.instructions.push(Instruction::Conditional {
            condition: ClassicalCondition::BitIsOne(0),
            gate: Gate::H,
            targets: smallvec![0],
        });
        c.instructions.push(Instruction::Measure {
            qubit: 0,
            classical_bit: 0,
        });
        let _: SmallVec<[usize; 4]> = smallvec![];
        let rows = build_measurement_rows(&c);
        assert_eq!(rows.len(), 1);
        assert_eq!(rows[0].0.x[0] & 1, 1);
    }

    fn build_avx_pauli(q: usize, m_words: usize) -> Vec<Vec<u64>> {
        let mut v = vec![vec![0u64; m_words]; q.max(1)];
        for word in v[0].iter_mut() {
            *word = 0xA5A5_A5A5_A5A5_A5A5;
        }
        v
    }

    #[test]
    fn batch_h_avx2_path_193_meas() {
        let m_words = 4;
        let mut x = build_avx_pauli(2, m_words);
        let mut z = build_avx_pauli(2, m_words);
        for word in z[0].iter_mut() {
            *word = 0x5A5A_5A5A_5A5A_5A5A;
        }
        let mut sign = vec![0u64; m_words];
        let x_before = x[0].clone();
        let z_before = z[0].clone();
        batch_propagate_backward(&mut x, &mut z, &mut sign, &Gate::H, &[0], m_words);
        assert_eq!(x[0], z_before);
        assert_eq!(z[0], x_before);
    }

    #[test]
    fn batch_s_and_sdg_avx2_path() {
        let m_words = 4;
        let mut x = build_avx_pauli(2, m_words);
        let mut z = build_avx_pauli(2, m_words);
        let mut sign = vec![0u64; m_words];
        batch_propagate_backward(&mut x, &mut z, &mut sign, &Gate::S, &[0], m_words);
        batch_propagate_backward(&mut x, &mut z, &mut sign, &Gate::Sdg, &[0], m_words);
        for word in z[0].iter().take(m_words) {
            assert_eq!(*word & 0xA5A5_A5A5_A5A5_A5A5, 0xA5A5_A5A5_A5A5_A5A5);
        }
    }

    #[test]
    fn batch_x_y_z_avx2_paths() {
        let m_words = 4;
        let mut x = build_avx_pauli(2, m_words);
        let mut z = build_avx_pauli(2, m_words);
        let mut sign = vec![0u64; m_words];
        let z_initial = z[0].clone();
        batch_propagate_backward(&mut x, &mut z, &mut sign, &Gate::X, &[0], m_words);
        for w in 0..m_words {
            assert_eq!(sign[w], z_initial[w]);
        }
        let mut sign2 = vec![0u64; m_words];
        batch_propagate_backward(&mut x, &mut z, &mut sign2, &Gate::Z, &[0], m_words);
        let mut sign3 = vec![0u64; m_words];
        batch_propagate_backward(&mut x, &mut z, &mut sign3, &Gate::Y, &[0], m_words);
        for w in 0..m_words {
            assert_eq!(sign3[w], x[0][w] ^ z[0][w]);
        }
    }

    #[test]
    fn batch_sx_avx2_paths() {
        let m_words = 4;
        let mut x = build_avx_pauli(2, m_words);
        let mut z = build_avx_pauli(2, m_words);
        let mut sign = vec![0u64; m_words];
        batch_propagate_backward(&mut x, &mut z, &mut sign, &Gate::SX, &[0], m_words);
        batch_propagate_backward(&mut x, &mut z, &mut sign, &Gate::SXdg, &[0], m_words);
    }

    #[test]
    fn batch_cx_cz_swap_avx2_paths() {
        let m_words = 4;
        let mut x = build_avx_pauli(2, m_words);
        let mut z = build_avx_pauli(2, m_words);
        let mut sign = vec![0u64; m_words];
        batch_propagate_backward(&mut x, &mut z, &mut sign, &Gate::Cx, &[0, 1], m_words);
        batch_propagate_backward(&mut x, &mut z, &mut sign, &Gate::Cz, &[0, 1], m_words);
        batch_propagate_backward(&mut x, &mut z, &mut sign, &Gate::Swap, &[0, 1], m_words);
        assert_eq!(x.len(), 2);
    }

    #[test]
    fn batch_path_with_tail_words() {
        let m_words = 5;
        let mut x = vec![vec![0xFFu64; m_words]; 2];
        let mut z = vec![vec![0xAAu64; m_words]; 2];
        let mut sign = vec![0u64; m_words];
        batch_propagate_backward(&mut x, &mut z, &mut sign, &Gate::H, &[0], m_words);
        batch_propagate_backward(&mut x, &mut z, &mut sign, &Gate::S, &[1], m_words);
        batch_propagate_backward(&mut x, &mut z, &mut sign, &Gate::Cx, &[0, 1], m_words);
    }

    #[test]
    fn colmajor_forward_sim_matches_stabilizer_tableau() {
        use crate::backend::Backend;
        use crate::backend::stabilizer::StabilizerBackend;
        use rand::{RngExt, SeedableRng};

        let one_qubit = [
            Gate::H,
            Gate::S,
            Gate::Sdg,
            Gate::X,
            Gate::Y,
            Gate::Z,
            Gate::Id,
            Gate::SX,
            Gate::SXdg,
        ];
        let two_qubit = [Gate::Cx, Gate::Cz, Gate::Swap];

        // Row words are 2n/64: 97 and 130 reach the AVX2 threshold of four, 130 and 300
        // add tail words past the vector chunks.
        for n in [5, 64, 97, 130, 300] {
            let mut rng = rand_chacha::ChaCha8Rng::seed_from_u64(42);
            let mut c = Circuit::new(n, 0);
            for _ in 0..12 * n {
                if rng.random_range(0..3) == 0 {
                    let a = rng.random_range(0..n);
                    let b = (a + rng.random_range(1..n)) % n;
                    c.add_gate(
                        two_qubit[rng.random_range(0..two_qubit.len())].clone(),
                        &[a, b],
                    );
                } else {
                    let q = rng.random_range(0..n);
                    c.add_gate(
                        one_qubit[rng.random_range(0..one_qubit.len())].clone(),
                        &[q],
                    );
                }
            }

            let (xz, phase, nw) = colmajor_forward_sim(n, &c.instructions).unwrap();

            let mut reference = StabilizerBackend::new(42);
            reference.init(n, 0).unwrap();
            reference.apply_gates_only(&c.instructions).unwrap();
            let (ref_xz, ref_phase) = reference.export_tableau().unwrap();

            assert_eq!(
                xz[..],
                ref_xz[..2 * n * 2 * nw],
                "tableau bits differ at n={n}"
            );
            assert_eq!(
                phase[..],
                ref_phase[..2 * n],
                "tableau signs differ at n={n}"
            );
        }
    }
}