prism-q 0.35.0

use super::kernels::{rowmul_words, xor_words};
use super::*;
use crate::backend::Backend;
use crate::circuit::{Circuit, Instruction};
use crate::error::PrismError;
use crate::gates::Gate;
use crate::sim;
use num_complex::Complex64;

#[test]
fn test_init_tableau() {
    let mut b = StabilizerBackend::new(42);
    b.init(3, 0).unwrap();
    assert_eq!(b.n, 3);
    let stride = b.stride();
    let nw = b.num_words;
    // Destabilizer 0 = X_0
    assert_eq!(b.xz[0] & 1, 1);
    assert_eq!((b.xz[0] >> 1) & 1, 0);
    assert_eq!(b.xz[nw] & 1, 0);
    // Stabilizer 0 (row 3) = Z_0
    assert_eq!(b.xz[3 * stride] & 1, 0);
    assert_eq!(b.xz[3 * stride + nw] & 1, 1);
}

#[test]
fn test_x_flips() {
    let mut c = Circuit::new(1, 1);
    c.add_gate(Gate::X, &[0]);
    c.add_measure(0, 0);
    let mut b = StabilizerBackend::new(42);
    sim::run_on(&mut b, &c).unwrap();
    assert_eq!(b.classical_results(), &[true]);
}

#[test]
fn test_z_on_zero() {
    let mut c = Circuit::new(1, 1);
    c.add_gate(Gate::Z, &[0]);
    c.add_measure(0, 0);
    let mut b = StabilizerBackend::new(42);
    sim::run_on(&mut b, &c).unwrap();
    assert_eq!(b.classical_results(), &[false]);
}

#[test]
fn test_y_flips() {
    let mut c = Circuit::new(1, 1);
    c.add_gate(Gate::Y, &[0]);
    c.add_measure(0, 0);
    let mut b = StabilizerBackend::new(42);
    sim::run_on(&mut b, &c).unwrap();
    assert_eq!(b.classical_results(), &[true]);
}

#[test]
fn test_hzh_equals_x() {
    let mut c = Circuit::new(1, 1);
    c.add_gate(Gate::H, &[0]);
    c.add_gate(Gate::Z, &[0]);
    c.add_gate(Gate::H, &[0]);
    c.add_measure(0, 0);
    let mut b = StabilizerBackend::new(42);
    sim::run_on(&mut b, &c).unwrap();
    assert_eq!(b.classical_results(), &[true]);
}

#[test]
fn test_s_squared_is_z() {
    let mut c = Circuit::new(1, 1);
    c.add_gate(Gate::H, &[0]);
    c.add_gate(Gate::S, &[0]);
    c.add_gate(Gate::S, &[0]);
    c.add_gate(Gate::H, &[0]);
    c.add_measure(0, 0);
    let mut b = StabilizerBackend::new(42);
    sim::run_on(&mut b, &c).unwrap();
    assert_eq!(b.classical_results(), &[true]);
}

#[test]
fn test_s_sdg_cancel() {
    let mut c = Circuit::new(1, 1);
    c.add_gate(Gate::H, &[0]);
    c.add_gate(Gate::S, &[0]);
    c.add_gate(Gate::Sdg, &[0]);
    c.add_gate(Gate::H, &[0]);
    c.add_measure(0, 0);
    let mut b = StabilizerBackend::new(42);
    sim::run_on(&mut b, &c).unwrap();
    assert_eq!(b.classical_results(), &[false]);
}

#[test]
fn test_h_superposition_deterministic() {
    let mut c = Circuit::new(1, 1);
    c.add_gate(Gate::H, &[0]);
    c.add_measure(0, 0);

    let mut b1 = StabilizerBackend::new(42);
    sim::run_on(&mut b1, &c).unwrap();
    let r1 = b1.classical_results()[0];

    let mut b2 = StabilizerBackend::new(42);
    sim::run_on(&mut b2, &c).unwrap();
    let r2 = b2.classical_results()[0];

    assert_eq!(r1, r2);
}

#[test]
fn test_bell_correlated() {
    let mut c = Circuit::new(2, 2);
    c.add_gate(Gate::H, &[0]);
    c.add_gate(Gate::Cx, &[0, 1]);
    c.add_measure(0, 0);
    c.add_measure(1, 1);

    for seed in [42, 100, 999, 12345] {
        let mut b = StabilizerBackend::new(seed);
        sim::run_on(&mut b, &c).unwrap();
        assert_eq!(
            b.classical_results()[0],
            b.classical_results()[1],
            "Bell state measurements must be equal (seed {seed})"
        );
    }
}

#[test]
fn test_ghz_4_correlated() {
    let mut c = Circuit::new(4, 4);
    c.add_gate(Gate::H, &[0]);
    for i in 0..3 {
        c.add_gate(Gate::Cx, &[i, i + 1]);
    }
    for i in 0..4 {
        c.add_measure(i, i);
    }

    for seed in [42, 100, 999] {
        let mut b = StabilizerBackend::new(seed);
        sim::run_on(&mut b, &c).unwrap();
        let results = b.classical_results();
        assert!(
            results.iter().all(|&x| x == results[0]),
            "GHZ-4 measurements must be equal (seed {seed})"
        );
    }
}

#[test]
fn test_swap() {
    let mut c = Circuit::new(2, 2);
    c.add_gate(Gate::X, &[1]);
    c.add_gate(Gate::Swap, &[0, 1]);
    c.add_measure(0, 0);
    c.add_measure(1, 1);
    let mut b = StabilizerBackend::new(42);
    sim::run_on(&mut b, &c).unwrap();
    assert_eq!(b.classical_results(), &[true, false]);
}

#[test]
fn test_cz_on_11() {
    let mut c = Circuit::new(2, 2);
    c.add_gate(Gate::X, &[0]);
    c.add_gate(Gate::X, &[1]);
    c.add_gate(Gate::Cz, &[0, 1]);
    c.add_measure(0, 0);
    c.add_measure(1, 1);
    let mut b = StabilizerBackend::new(42);
    sim::run_on(&mut b, &c).unwrap();
    assert_eq!(b.classical_results(), &[true, true]);
}

#[test]
fn test_rejects_t_gate() {
    let mut c = Circuit::new(1, 0);
    c.add_gate(Gate::T, &[0]);
    let mut b = StabilizerBackend::new(42);
    b.init(1, 0).unwrap();
    let err = b.apply(&c.instructions[0]).unwrap_err();
    assert!(matches!(err, PrismError::BackendUnsupported { .. }));
}

#[test]
fn test_rejects_rx_gate() {
    let mut c = Circuit::new(1, 0);
    c.add_gate(Gate::Rx(0.5), &[0]);
    let mut b = StabilizerBackend::new(42);
    b.init(1, 0).unwrap();
    let err = b.apply(&c.instructions[0]).unwrap_err();
    assert!(matches!(err, PrismError::BackendUnsupported { .. }));
}

#[test]
fn test_probs_zero_state() {
    let mut b = StabilizerBackend::new(42);
    b.init(1, 0).unwrap();
    let probs = b.probabilities().unwrap();
    assert_eq!(probs, vec![1.0, 0.0]);
}

#[test]
fn test_probs_one_state() {
    let mut c = Circuit::new(1, 0);
    c.add_gate(Gate::X, &[0]);
    let mut b = StabilizerBackend::new(42);
    sim::run_on(&mut b, &c).unwrap();
    let probs = b.probabilities().unwrap();
    assert_eq!(probs, vec![0.0, 1.0]);
}

#[test]
fn test_probs_plus_state() {
    let mut c = Circuit::new(1, 0);
    c.add_gate(Gate::H, &[0]);
    let mut b = StabilizerBackend::new(42);
    sim::run_on(&mut b, &c).unwrap();
    let probs = b.probabilities().unwrap();
    assert!((probs[0] - 0.5).abs() < 1e-12);
    assert!((probs[1] - 0.5).abs() < 1e-12);
}

#[test]
fn test_probs_bell_state() {
    let mut c = Circuit::new(2, 0);
    c.add_gate(Gate::H, &[0]);
    c.add_gate(Gate::Cx, &[0, 1]);
    let mut b = StabilizerBackend::new(42);
    sim::run_on(&mut b, &c).unwrap();
    let probs = b.probabilities().unwrap();
    assert!((probs[0] - 0.5).abs() < 1e-12);
    assert!(probs[1].abs() < 1e-12);
    assert!(probs[2].abs() < 1e-12);
    assert!((probs[3] - 0.5).abs() < 1e-12);
}

#[test]
fn test_probs_ghz_4() {
    let mut c = Circuit::new(4, 0);
    c.add_gate(Gate::H, &[0]);
    for i in 0..3 {
        c.add_gate(Gate::Cx, &[i, i + 1]);
    }
    let mut b = StabilizerBackend::new(42);
    sim::run_on(&mut b, &c).unwrap();
    let probs = b.probabilities().unwrap();
    assert_eq!(probs.len(), 16);
    assert!((probs[0] - 0.5).abs() < 1e-12);
    assert!((probs[15] - 0.5).abs() < 1e-12);
    let rest_sum: f64 = probs[1..15].iter().sum();
    assert!(rest_sum.abs() < 1e-12);
}

#[test]
fn test_probabilities_reject_large_qubit_counts() {
    let mut b = StabilizerBackend::new(42);
    b.init(usize::BITS as usize, 0).unwrap();
    let err = b.probabilities().unwrap_err();
    match err {
        PrismError::BackendUnsupported { operation, .. } => {
            assert!(operation.contains("exceeds addressable memory"));
        }
        other => panic!("unexpected error: {other:?}"),
    }
}

#[test]
fn test_1000_qubit_ghz() {
    let n = 1000;
    let mut c = Circuit::new(n, n);
    c.add_gate(Gate::H, &[0]);
    for i in 0..n - 1 {
        c.add_gate(Gate::Cx, &[i, i + 1]);
    }
    for i in 0..n {
        c.add_measure(i, i);
    }
    let mut b = StabilizerBackend::new(42);
    sim::run_on(&mut b, &c).unwrap();
    let results = b.classical_results();
    assert!(results.iter().all(|&x| x == results[0]));
}

#[test]
fn test_id_no_change() {
    let mut c = Circuit::new(1, 1);
    c.add_gate(Gate::Id, &[0]);
    c.add_measure(0, 0);
    let mut b = StabilizerBackend::new(42);
    sim::run_on(&mut b, &c).unwrap();
    assert_eq!(b.classical_results(), &[false]);
}

#[test]
fn test_double_x_is_identity() {
    let mut c = Circuit::new(1, 1);
    c.add_gate(Gate::X, &[0]);
    c.add_gate(Gate::X, &[0]);
    c.add_measure(0, 0);
    let mut b = StabilizerBackend::new(42);
    sim::run_on(&mut b, &c).unwrap();
    assert_eq!(b.classical_results(), &[false]);
}

#[test]
fn test_supports_fused_gates() {
    let b = StabilizerBackend::new(42);
    assert!(!b.supports_fused_gates());
}

#[test]
fn test_rejects_cu_gate() {
    let h_mat = Gate::H.matrix_2x2();
    let mut c = Circuit::new(2, 0);
    c.add_gate(Gate::Cu(Box::new(h_mat)), &[0, 1]);
    let mut b = StabilizerBackend::new(42);
    b.init(2, 0).unwrap();
    let err = b.apply(&c.instructions[0]).unwrap_err();
    assert!(matches!(err, PrismError::BackendUnsupported { .. }));
}

#[test]
fn test_rejects_mcu_gate() {
    use crate::gates::McuData;
    let x_mat = Gate::X.matrix_2x2();
    let mcu = Gate::Mcu(Box::new(McuData {
        mat: x_mat,
        num_controls: 2,
    }));
    let mut c = Circuit::new(3, 0);
    c.add_gate(mcu, &[0, 1, 2]);
    let mut b = StabilizerBackend::new(42);
    b.init(3, 0).unwrap();
    let err = b.apply(&c.instructions[0]).unwrap_err();
    assert!(matches!(err, PrismError::BackendUnsupported { .. }));
}

#[test]
fn test_export_statevector_zero_state() {
    let mut b = StabilizerBackend::new(42);
    b.init(1, 0).unwrap();
    let sv = b.export_statevector().unwrap();
    assert!((sv[0].re - 1.0).abs() < 1e-10);
    assert!(sv[0].im.abs() < 1e-10);
    assert!(sv[1].norm_sqr() < 1e-10);
}

#[test]
fn test_export_statevector_one_state() {
    let mut c = Circuit::new(1, 0);
    c.add_gate(Gate::X, &[0]);
    let mut b = StabilizerBackend::new(42);
    sim::run_on(&mut b, &c).unwrap();
    let sv = b.export_statevector().unwrap();
    assert!(sv[0].norm_sqr() < 1e-10);
    assert!((sv[1].norm_sqr() - 1.0).abs() < 1e-10);
}

#[test]
fn test_export_statevector_plus_state() {
    let mut c = Circuit::new(1, 0);
    c.add_gate(Gate::H, &[0]);
    let mut b = StabilizerBackend::new(42);
    sim::run_on(&mut b, &c).unwrap();
    let sv = b.export_statevector().unwrap();
    let expected = 1.0 / 2.0_f64.sqrt();
    assert!((sv[0].re - expected).abs() < 1e-10);
    assert!((sv[1].re - expected).abs() < 1e-10);
}

#[test]
fn test_export_statevector_minus_state() {
    // H X |0⟩ = H|1⟩ = (|0⟩ - |1⟩)/√2
    let mut c = Circuit::new(1, 0);
    c.add_gate(Gate::X, &[0]);
    c.add_gate(Gate::H, &[0]);
    let mut b = StabilizerBackend::new(42);
    sim::run_on(&mut b, &c).unwrap();
    let sv = b.export_statevector().unwrap();
    let expected = 1.0 / 2.0_f64.sqrt();
    assert!((sv[0].re - expected).abs() < 1e-10);
    assert!((sv[1].re + expected).abs() < 1e-10);
}

#[test]
fn test_export_statevector_bell_state_matches_sv() {
    use crate::backend::statevector::StatevectorBackend;

    let mut c = Circuit::new(2, 0);
    c.add_gate(Gate::H, &[0]);
    c.add_gate(Gate::Cx, &[0, 1]);

    let mut stab = StabilizerBackend::new(42);
    sim::run_on(&mut stab, &c).unwrap();
    let stab_sv = stab.export_statevector().unwrap();

    let mut sv = StatevectorBackend::new(42);
    sim::run_on(&mut sv, &c).unwrap();
    let sv_ref = sv.state_vector();

    for (s, r) in stab_sv.iter().zip(sv_ref.iter()) {
        assert!(
            (s.norm_sqr() - r.norm_sqr()).abs() < 1e-10,
            "prob mismatch: stab={}, sv={}",
            s.norm_sqr(),
            r.norm_sqr()
        );
    }

    // Global phase: find ratio between first non-zero pair
    let global_phase = find_global_phase(&stab_sv, sv_ref);

    // After removing global phase, amplitudes must match
    for (s, r) in stab_sv.iter().zip(sv_ref.iter()) {
        let adjusted = s * global_phase;
        assert!(
            (adjusted - r).norm() < 1e-10,
            "amplitude mismatch after phase: stab*phase={adjusted:?}, sv={r:?}"
        );
    }
}

#[test]
fn test_export_statevector_ghz3_matches_sv() {
    use crate::backend::statevector::StatevectorBackend;

    let mut c = Circuit::new(3, 0);
    c.add_gate(Gate::H, &[0]);
    c.add_gate(Gate::Cx, &[0, 1]);
    c.add_gate(Gate::Cx, &[1, 2]);

    let mut stab = StabilizerBackend::new(42);
    sim::run_on(&mut stab, &c).unwrap();
    let stab_sv = stab.export_statevector().unwrap();

    let mut sv = StatevectorBackend::new(42);
    sim::run_on(&mut sv, &c).unwrap();
    let sv_ref = sv.state_vector();

    let global_phase = find_global_phase(&stab_sv, sv_ref);
    for (s, r) in stab_sv.iter().zip(sv_ref.iter()) {
        let adjusted = s * global_phase;
        assert!(
            (adjusted - r).norm() < 1e-10,
            "GHZ3 mismatch: stab*phase={adjusted:?}, sv={r:?}"
        );
    }
}

#[test]
fn test_export_statevector_complex_clifford_matches_sv() {
    use crate::backend::statevector::StatevectorBackend;

    // Deeper Clifford circuit with S, Sdg, CZ, SX to exercise i-phases
    let mut c = Circuit::new(4, 0);
    c.add_gate(Gate::H, &[0]);
    c.add_gate(Gate::S, &[0]);
    c.add_gate(Gate::Cx, &[0, 1]);
    c.add_gate(Gate::H, &[2]);
    c.add_gate(Gate::Cz, &[1, 2]);
    c.add_gate(Gate::SX, &[3]);
    c.add_gate(Gate::Cx, &[2, 3]);
    c.add_gate(Gate::Sdg, &[1]);
    c.add_gate(Gate::H, &[3]);
    c.add_gate(Gate::S, &[2]);
    c.add_gate(Gate::Swap, &[0, 3]);

    let mut stab = StabilizerBackend::new(42);
    sim::run_on(&mut stab, &c).unwrap();
    let stab_sv = stab.export_statevector().unwrap();

    let mut sv = StatevectorBackend::new(42);
    sim::run_on(&mut sv, &c).unwrap();
    let sv_ref = sv.state_vector();

    let global_phase = find_global_phase(&stab_sv, sv_ref);
    for (i, (s, r)) in stab_sv.iter().zip(sv_ref.iter()).enumerate() {
        let adjusted = s * global_phase;
        assert!(
            (adjusted - r).norm() < 1e-10,
            "4q Clifford mismatch at index {i}: stab*phase={adjusted:?}, sv={r:?}"
        );
    }
}

#[test]
fn test_export_statevector_all_paulis_match_sv() {
    use crate::backend::statevector::StatevectorBackend;

    // Y gate introduces i-phase: Y|0⟩ = i|1⟩
    let mut c = Circuit::new(2, 0);
    c.add_gate(Gate::Y, &[0]);
    c.add_gate(Gate::H, &[1]);
    c.add_gate(Gate::S, &[1]);

    let mut stab = StabilizerBackend::new(42);
    sim::run_on(&mut stab, &c).unwrap();
    let stab_sv = stab.export_statevector().unwrap();

    let mut sv = StatevectorBackend::new(42);
    sim::run_on(&mut sv, &c).unwrap();
    let sv_ref = sv.state_vector();

    let global_phase = find_global_phase(&stab_sv, sv_ref);
    for (i, (s, r)) in stab_sv.iter().zip(sv_ref.iter()).enumerate() {
        let adjusted = s * global_phase;
        assert!(
            (adjusted - r).norm() < 1e-10,
            "Pauli test mismatch at {i}: stab*phase={adjusted:?}, sv={r:?}"
        );
    }
}

/// Find the global phase ratio between two state vectors.
/// Returns the complex scalar c such that stab * c ≈ reference.
fn find_global_phase(stab: &[Complex64], reference: &[Complex64]) -> Complex64 {
    for (s, r) in stab.iter().zip(reference.iter()) {
        if s.norm_sqr() > 1e-10 && r.norm_sqr() > 1e-10 {
            return r / s;
        }
    }
    Complex64::new(1.0, 0.0)
}

#[test]
fn test_rejects_cphase_gate() {
    let mut c = Circuit::new(2, 0);
    c.add_gate(Gate::cphase(std::f64::consts::FRAC_PI_4), &[0, 1]);
    let mut b = StabilizerBackend::new(42);
    b.init(2, 0).unwrap();
    let err = b.apply(&c.instructions[0]).unwrap_err();
    assert!(matches!(err, PrismError::BackendUnsupported { .. }));
}

#[test]
fn test_xor_words_various_lengths() {
    for len in 0..=17 {
        let src: Vec<u64> = (0..len)
            .map(|i| 0xAAAA_BBBB_CCCC_0000u64 | i as u64)
            .collect();
        let original: Vec<u64> = (0..len)
            .map(|i| 0x1111_2222_3333_0000u64 | (i as u64 * 7))
            .collect();

        let mut expected = original.clone();
        for i in 0..len {
            expected[i] ^= src[i];
        }

        let mut actual = original.clone();
        if len > 0 {
            // SAFETY: actual and src are valid for len elements and do not overlap.
            unsafe { xor_words(actual.as_mut_ptr(), src.as_ptr(), len) };
        }

        assert_eq!(actual, expected, "xor_words mismatch at len={len}");
    }
}

#[test]
fn test_xor_words_all_ones_and_zeros() {
    let len = 8;
    let src = vec![u64::MAX; len];
    let mut dst = vec![0u64; len];
    // SAFETY: dst and src are valid for len elements and do not overlap.
    unsafe { xor_words(dst.as_mut_ptr(), src.as_ptr(), len) };
    assert!(dst.iter().all(|&v| v == u64::MAX));

    // SAFETY: dst and src are valid for len elements and do not overlap.
    unsafe { xor_words(dst.as_mut_ptr(), src.as_ptr(), len) };
    assert!(dst.iter().all(|&v| v == 0));
}

#[test]
fn test_rowmul_words_zero_src() {
    let nw = 4;
    let mut dst_x = vec![0xFFu64; nw];
    let mut dst_z = vec![0xAAu64; nw];
    let src_x = vec![0u64; nw];
    let src_z = vec![0u64; nw];
    let sum = rowmul_words(&mut dst_x, &mut dst_z, &src_x, &src_z, 0);
    assert_eq!(dst_x, vec![0xFFu64; nw]);
    assert_eq!(dst_z, vec![0xAAu64; nw]);
    assert_eq!(sum & 3, 0);
}

#[test]
fn test_rowmul_words_matches_manual() {
    let nw = 3;
    let src_x = vec![0b1010u64, 0b1100u64, 0b0011u64];
    let src_z = vec![0b0110u64, 0b1010u64, 0b1001u64];
    let orig_dst_x = vec![0b1100u64, 0b0110u64, 0b1010u64];
    let orig_dst_z = vec![0b0011u64, 0b1001u64, 0b0110u64];

    let mut manual_x = orig_dst_x.clone();
    let mut manual_z = orig_dst_z.clone();
    let mut manual_sum = 4u64;
    for w in 0..nw {
        let x1 = src_x[w];
        let z1 = src_z[w];
        let x2 = manual_x[w];
        let z2 = manual_z[w];
        let new_x = x1 ^ x2;
        let new_z = z1 ^ z2;
        manual_x[w] = new_x;
        manual_z[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);
            manual_sum = manual_sum.wrapping_add(2 * pos.count_ones() as u64);
            manual_sum = manual_sum.wrapping_sub(nonzero.count_ones() as u64);
        }
    }

    let mut fn_x = orig_dst_x.clone();
    let mut fn_z = orig_dst_z.clone();
    let fn_sum = rowmul_words(&mut fn_x, &mut fn_z, &src_x, &src_z, 4);

    assert_eq!(fn_x, manual_x);
    assert_eq!(fn_z, manual_z);
    assert_eq!(fn_sum & 3, manual_sum & 3);
}

#[test]
fn test_rowmul_words_phase_y_times_x() {
    let src_x = vec![1u64];
    let src_z = vec![1u64];
    let mut dst_x = vec![1u64];
    let mut dst_z = vec![0u64];
    let sum = rowmul_words(&mut dst_x, &mut dst_z, &src_x, &src_z, 0);
    assert_eq!(dst_x[0], 0);
    assert_eq!(dst_z[0], 1);
    assert!(
        (sum & 3) >= 2,
        "Y*X should give phase -1 (sum&3={})",
        sum & 3
    );
}

#[test]
fn test_rowmul_words_simd_large() {
    for nw in [4, 5, 8, 9, 16, 17] {
        let src_x: Vec<u64> = (0..nw)
            .map(|i| 0xAAAA_BBBB_CCCC_0000u64 | i as u64)
            .collect();
        let src_z: Vec<u64> = (0..nw)
            .map(|i| 0x1111_2222_3333_0000u64 | (i * 7) as u64)
            .collect();
        let orig_x: Vec<u64> = (0..nw)
            .map(|i| 0x5555_6666_7777_0000u64 | (i * 3) as u64)
            .collect();
        let orig_z: Vec<u64> = (0..nw)
            .map(|i| 0x9999_AAAA_BBBB_0000u64 | (i * 5) as u64)
            .collect();

        let mut ref_x = orig_x.clone();
        let mut ref_z = orig_z.clone();
        let mut ref_sum = 2u64;
        for w in 0..nw {
            let x1 = src_x[w];
            let z1 = src_z[w];
            let x2 = ref_x[w];
            let z2 = ref_z[w];
            let new_x = x1 ^ x2;
            let new_z = z1 ^ z2;
            ref_x[w] = new_x;
            ref_z[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);
                ref_sum = ref_sum.wrapping_add(2 * pos.count_ones() as u64);
                ref_sum = ref_sum.wrapping_sub(nonzero.count_ones() as u64);
            }
        }

        let mut fn_x = orig_x.clone();
        let mut fn_z = orig_z.clone();
        let fn_sum = rowmul_words(&mut fn_x, &mut fn_z, &src_x, &src_z, 2);

        assert_eq!(fn_x, ref_x, "XOR mismatch at nw={nw}");
        assert_eq!(fn_z, ref_z, "XOR mismatch at nw={nw}");
        assert_eq!(
            fn_sum & 3,
            ref_sum & 3,
            "Phase mismatch at nw={nw}: fn_sum={fn_sum}, ref_sum={ref_sum}"
        );
    }
}

#[test]
fn test_rowmul_refactor_preserves_ghz_correctness() {
    let n = 10;
    let mut c = Circuit::new(n, n);
    c.add_gate(Gate::H, &[0]);
    for i in 0..n - 1 {
        c.add_gate(Gate::Cx, &[i, i + 1]);
    }
    for i in 0..n {
        c.add_measure(i, i);
    }

    for seed in 0..20u64 {
        let mut b = StabilizerBackend::new(seed);
        sim::run_on(&mut b, &c).unwrap();
        let results = b.classical_results();
        let first = results[0];
        assert!(
            results.iter().all(|&r| r == first),
            "GHZ violation at seed {seed}: {results:?}"
        );
    }
}

#[test]
fn test_rowmul_refactor_preserves_probabilities() {
    let mut c = Circuit::new(3, 0);
    c.add_gate(Gate::H, &[0]);
    c.add_gate(Gate::Cx, &[0, 1]);
    let mut b = StabilizerBackend::new(42);
    sim::run_on(&mut b, &c).unwrap();
    let probs = b.probabilities().unwrap();
    assert!((probs[0] - 0.5).abs() < 1e-10);
    assert!(probs[1].abs() < 1e-10);
    assert!(probs[2].abs() < 1e-10);
    assert!((probs[3] - 0.5).abs() < 1e-10);
    assert!(probs[4].abs() < 1e-10);
    assert!(probs[5].abs() < 1e-10);
    assert!((probs[6]).abs() < 1e-10);
    assert!(probs[7].abs() < 1e-10);
}

#[test]
fn test_rowmul_multi_word_correctness() {
    let n = 65;
    let mut c = Circuit::new(n, n);
    c.add_gate(Gate::H, &[0]);
    c.add_gate(Gate::Cx, &[0, 64]);
    c.add_measure(0, 0);
    c.add_measure(64, 1);

    for seed in 0..10u64 {
        let mut b = StabilizerBackend::new(seed);
        sim::run_on(&mut b, &c).unwrap();
        let results = b.classical_results();
        assert_eq!(
            results[0], results[1],
            "Bell pair violation at seed {seed}: q0={}, q64={}",
            results[0], results[1]
        );
    }
}

#[test]
fn test_sgi_500q_clifford_d10_matches_gate_by_gate() {
    use crate::circuits;
    let n = 500;
    let mut circuit = circuits::clifford_heavy_circuit(n, 10, 42);
    circuit.measure_all();

    let mut b1 = StabilizerBackend::new(42);
    b1.init(circuit.num_qubits, circuit.num_classical_bits)
        .unwrap();
    for instr in &circuit.instructions {
        b1.apply(instr).unwrap();
    }
    let r1 = b1.classical_results().to_vec();

    let mut b2 = StabilizerBackend::new(42);
    sim::run_on(&mut b2, &circuit).unwrap();
    let r2 = b2.classical_results().to_vec();

    assert_eq!(
        r1, r2,
        "SGI 500q Clifford d10: gate-by-gate vs apply_instructions mismatch"
    );
}

// A CNOT wall collapses on every qubit: each measurement finds a pivot, so the
// index is widened n times in a row with the SGI path still engaged. That is
// the shape the widening is there for, so it is pinned against a gate-by-gate
// run of the same circuit.
#[test]
fn test_sgi_cnot_wall_measure_all_matches_gate_by_gate() {
    let n = 400;
    let mut c = Circuit::new(n, n);
    for q in 0..n {
        c.add_gate(Gate::H, &[q]);
    }
    for d in 0..10 {
        let mut q = d % 2;
        while q + 1 < n {
            c.add_gate(Gate::Cx, &[q, q + 1]);
            q += 2;
        }
    }
    c.measure_all();

    let mut b1 = StabilizerBackend::new(42);
    b1.init(c.num_qubits, c.num_classical_bits).unwrap();
    for instr in &c.instructions {
        b1.apply(instr).unwrap();
    }
    let r1 = b1.classical_results().to_vec();

    let mut b2 = StabilizerBackend::new(42);
    sim::run_on(&mut b2, &c).unwrap();
    let r2 = b2.classical_results().to_vec();

    assert_eq!(
        r1, r2,
        "CNOT wall 400q measure-all: gate-by-gate vs apply_instructions mismatch"
    );
}

#[test]
fn test_sgi_300q_ghz_all_agree() {
    let n = 300;
    let mut c = Circuit::new(n, n);
    c.add_gate(Gate::H, &[0]);
    for i in 0..n - 1 {
        c.add_gate(Gate::Cx, &[i, i + 1]);
    }
    for i in 0..n {
        c.add_measure(i, i);
    }
    let mut b = StabilizerBackend::new(42);
    sim::run_on(&mut b, &c).unwrap();
    let results = b.classical_results();
    assert!(
        results.iter().all(|&x| x == results[0]),
        "GHZ 300q: not all qubits agree"
    );
}

#[test]
fn test_sgi_index_consistency() {
    let n = 300;
    let mut circuit = crate::circuits::clifford_heavy_circuit(n, 5, 42);
    circuit.num_classical_bits = 0;

    let mut b = StabilizerBackend::new(42);
    b.init(circuit.num_qubits, circuit.num_classical_bits)
        .unwrap();
    b.apply_instructions(&circuit.instructions).unwrap();

    let stride = b.stride();
    let nw = b.num_words;
    for q in 0..n {
        for &g in &b.qubit_active[q] {
            let g = g as usize;
            let row = &b.xz[g * stride..(g + 1) * stride];
            let word = q / 64;
            let bit_mask = 1u64 << (q % 64);
            let x = row[word] & bit_mask;
            let z = row[nw + word] & bit_mask;
            assert!(
                x != 0 || z != 0,
                "qubit_active[{q}] contains generator {g} which has I on qubit {q}"
            );
        }
    }

    for g in 0..2 * n {
        let row = &b.xz[g * stride..(g + 1) * stride];
        for q in 0..n {
            let word = q / 64;
            let bit_mask = 1u64 << (q % 64);
            let x = row[word] & bit_mask;
            let z = row[nw + word] & bit_mask;
            let active = x != 0 || z != 0;
            let in_index = b.qubit_active[q].contains(&(g as u32));
            assert_eq!(
                active, in_index,
                "generator {g} qubit {q}: active={active} but in_index={in_index}"
            );
        }
    }
}

// The brick layer keeps row weight low so the SGI path stays engaged, and the
// long-range CX fan out of the reset target gives the measurement pivot support
// on qubits that the rows anticommuting with it do not already touch. That is
// what makes the post-reset index widening observable.
fn sgi_reset_circuit(n: usize) -> Circuit {
    let mut c = crate::circuits::clifford_heavy_circuit(n, 2, 42);
    for k in 0..8 {
        c.add_gate(Gate::Cx, &[1, 30 + 17 * k]);
    }
    c.add_gate(Gate::H, &[1]);
    c.add_reset(1);
    for instr in crate::circuits::clifford_heavy_circuit(n, 2, 7).instructions {
        c.instructions.push(instr);
    }
    c.measure_all();
    c
}

// Guards the two tests below against passing vacuously: they only exercise the
// defect if the SGI path is still engaged when the reset is reached.
fn assert_sgi_engaged_at_reset(circuit: &Circuit) {
    let reset_at = circuit
        .instructions
        .iter()
        .position(|i| matches!(i, Instruction::Reset { .. }))
        .expect("circuit must contain a reset");

    let mut probe = StabilizerBackend::new(42);
    probe
        .init(circuit.num_qubits, circuit.num_classical_bits)
        .unwrap();
    probe
        .apply_instructions(&circuit.instructions[..reset_at])
        .unwrap();
    assert!(
        probe.sgi_enabled(),
        "SGI path not engaged at the reset; the test would pass vacuously"
    );
}

#[test]
fn test_sgi_reset_random_outcome_matches_gate_by_gate() {
    let n = 300;
    let circuit = sgi_reset_circuit(n);
    assert_sgi_engaged_at_reset(&circuit);

    let mut b1 = StabilizerBackend::new(42);
    b1.init(circuit.num_qubits, circuit.num_classical_bits)
        .unwrap();
    for instr in &circuit.instructions {
        b1.apply(instr).unwrap();
    }
    let r1 = b1.classical_results().to_vec();

    let mut b2 = StabilizerBackend::new(42);
    sim::run_on(&mut b2, &circuit).unwrap();
    let r2 = b2.classical_results().to_vec();

    assert_eq!(
        r1, r2,
        "SGI mid-circuit reset: gate-by-gate vs apply_instructions mismatch"
    );
}

// A collapse widens the index by the pivot's support rather than rebuilding it,
// so past a measurement or reset the lists are a superset of the true support.
// That is the invariant every reader needs: an entry with no support is a
// wasted visit, an unlisted active row is a dropped update.
#[test]
fn test_sgi_reset_leaves_index_covering() {
    let n = 300;
    let circuit = sgi_reset_circuit(n);
    assert_sgi_engaged_at_reset(&circuit);

    // Stop at the reset. Past it the run drops out of the SGI path, and the
    // word-batch route does not maintain the index at all, so a whole-circuit
    // check would assert an invariant that only holds while SGI is engaged.
    let through_reset = circuit
        .instructions
        .iter()
        .position(|i| matches!(i, Instruction::Reset { .. }))
        .expect("circuit must contain a reset")
        + 1;

    let mut b = StabilizerBackend::new(42);
    b.init(circuit.num_qubits, circuit.num_classical_bits)
        .unwrap();
    b.apply_instructions(&circuit.instructions[..through_reset])
        .unwrap();

    let stride = b.stride();
    let nw = b.num_words;
    for q in 0..n {
        let list = &b.qubit_active[q];
        assert!(
            list.windows(2).all(|w| w[0] < w[1]),
            "qubit_active[{q}] is not ascending and duplicate free"
        );
    }
    for g in 0..2 * n {
        let row = &b.xz[g * stride..(g + 1) * stride];
        for q in 0..n {
            let word = q / 64;
            let bit_mask = 1u64 << (q % 64);
            let active = row[word] & bit_mask != 0 || row[nw + word] & bit_mask != 0;
            assert!(
                !active || b.qubit_active[q].contains(&(g as u32)),
                "after reset: generator {g} is active on qubit {q} but the index omits it"
            );
        }
    }
}

// Per-instruction `apply` never maintains `qubit_active`, and the init-time
// counters stay light, so a later bulk entry re-checks `sgi_enabled`, sees light
// counters, and runs SGI against an index that no longer describes the rows. A
// guarded region is the in-repo path that reaches `apply_instructions` after a
// per-instruction prefix.
//
// The region opens with gates because an SGI gate iterates `qubit_active[q]`
// and updates only the rows it lists, which is where a stale index silently
// drops the rest.
//
// The reference drives the same circuit entirely per instruction, so it never
// reaches a bulk entry and never takes the SGI path at all.
#[test]
fn test_sgi_per_instruction_prefix_then_region_matches_gate_by_gate() {
    use crate::circuit::{ClassicalCondition, GuardedRegion};

    let n = 300;
    let partner = 137;

    let mut prefix = Circuit::new(n, 2);
    prefix.add_gate(Gate::H, &[0]);
    for q in 1..n {
        prefix.add_gate(Gate::Cx, &[0, q]);
    }

    let body = vec![
        Instruction::Gate {
            gate: Gate::H,
            targets: smallvec::smallvec![0],
        },
        Instruction::Gate {
            gate: Gate::S,
            targets: smallvec::smallvec![partner],
        },
        Instruction::Gate {
            gate: Gate::Cx,
            targets: smallvec::smallvec![0, partner],
        },
        Instruction::Measure {
            qubit: 0,
            classical_bit: 0,
        },
        Instruction::Measure {
            qubit: partner,
            classical_bit: 1,
        },
    ];
    let region = Instruction::Region(Box::new(GuardedRegion::new(
        ClassicalCondition::BitIsZero(0),
        body.clone(),
    )));

    let mut broken = Vec::new();
    for seed in 0..40u64 {
        let mut bulk = StabilizerBackend::new(seed);
        bulk.init(n, 2).unwrap();
        for instr in &prefix.instructions {
            bulk.apply(instr).unwrap();
        }
        bulk.apply(&region).unwrap();

        let mut dense = StabilizerBackend::new(seed);
        dense.init(n, 2).unwrap();
        for instr in prefix.instructions.iter().chain(body.iter()) {
            dense.apply(instr).unwrap();
        }

        if bulk.classical_results() != dense.classical_results() {
            broken.push(seed);
        }
    }

    assert!(
        broken.is_empty(),
        "the region disagreed with gate-by-gate driving on {} of 40 seeds: {:?}",
        broken.len(),
        broken
    );
}

#[test]
fn test_pcc_random_pairs_matches_gate_by_gate() {
    use crate::circuits;
    let n = 500;
    let mut circuit = circuits::clifford_random_pairs(n, 10, 42);
    circuit.measure_all();

    let mut b1 = StabilizerBackend::new(42);
    b1.init(circuit.num_qubits, circuit.num_classical_bits)
        .unwrap();
    for instr in &circuit.instructions {
        b1.apply(instr).unwrap();
    }
    let r1 = b1.classical_results().to_vec();

    let mut b2 = StabilizerBackend::new(42);
    sim::run_on(&mut b2, &circuit).unwrap();
    let r2 = b2.classical_results().to_vec();

    assert_eq!(
        r1, r2,
        "PCC 500q random-pairs d10: gate-by-gate vs apply_instructions mismatch"
    );
}

// A CX chain over a tableau whose rows carry X and Z bits everywhere, crossing
// seven word boundaries and starting and ending mid-word, so every phase term
// the chain kernel folds is live. Compared against the per-gate path on the
// eager run, where the batch driver owns the chain, and on the lazy run, where
// the kernel sweeps the stabilizer half alone.
#[test]
fn test_cx_chain_matches_gate_by_gate_across_words() {
    use crate::circuits;
    let n = 500;
    for seed in 0..4u64 {
        let mut circuit = circuits::clifford_random_pairs(n, 4, seed);
        for k in 3..n - 2 {
            circuit.add_gate(Gate::Cx, &[k, k + 1]);
        }
        let mut gates_only = circuit.clone();
        circuit.measure_all();

        let mut b1 = StabilizerBackend::new(seed);
        b1.init(n, circuit.num_classical_bits).unwrap();
        for instr in &circuit.instructions {
            b1.apply(instr).unwrap();
        }
        let mut b2 = StabilizerBackend::new(seed);
        sim::run_on(&mut b2, &circuit).unwrap();
        assert_eq!(
            b1.classical_results(),
            b2.classical_results(),
            "seed {seed}: eager chain outcomes"
        );

        gates_only.num_classical_bits = 0;
        let mut l1 = StabilizerBackend::new_lazy(seed);
        l1.init(n, 0).unwrap();
        for instr in &gates_only.instructions {
            l1.apply(instr).unwrap();
        }
        let mut l2 = StabilizerBackend::new_lazy(seed);
        sim::run_on(&mut l2, &gates_only).unwrap();
        let (xz1, ph1) = l1.raw_tableau();
        let (xz2, ph2) = l2.raw_tableau();
        let stride = 2 * n.div_ceil(64);
        assert_eq!(
            &xz1[n * stride..],
            &xz2[n * stride..],
            "seed {seed}: lazy chain rows"
        );
        assert_eq!(&ph1[n..], &ph2[n..], "seed {seed}: lazy chain phases");
    }
}

// Runs shorter than the chain kernel's floor still fuse inside a word group,
// including one that ends on bit 63, and a lone chain gate stays a plain gate.
#[test]
fn test_short_cx_chains_fuse_inside_word_groups() {
    use crate::circuits;
    let n = 320;
    for seed in 0..4u64 {
        let mut circuit = circuits::clifford_random_pairs(n, 3, seed);
        for (start, len) in [
            (1usize, 2usize),
            (10, 7),
            (49, 14),
            (64, 15),
            (130, 3),
            (200, 1),
        ] {
            for k in start..start + len {
                circuit.add_gate(Gate::Cx, &[k, k + 1]);
            }
            circuit.add_gate(Gate::S, &[start]);
        }
        circuit.measure_all();

        let mut b1 = StabilizerBackend::new(seed);
        b1.init(n, circuit.num_classical_bits).unwrap();
        for instr in &circuit.instructions {
            b1.apply(instr).unwrap();
        }
        let mut b2 = StabilizerBackend::new(seed);
        sim::run_on(&mut b2, &circuit).unwrap();
        assert_eq!(
            b1.classical_results(),
            b2.classical_results(),
            "seed {seed}"
        );
    }
}

// The per-instruction path must stay lazy across gates, reconstruct on the
// first measurement, and produce the exact outcome sequence of an eager run:
// gates update the stabilizer half identically in both modes, so pivot
// classification and RNG draws align, and deterministic outcomes are unique.
#[test]
fn lazy_destab_apply_path_matches_eager_exactly() {
    for n in [3usize, 5, 8, 16, 80, 300] {
        for seed in 0..20u64 {
            let mut circuit = crate::circuits::clifford_heavy_circuit(n, 10, seed);
            circuit.measure_all();

            let mut eager = StabilizerBackend::new(seed);
            eager.init(n, n).unwrap();
            for inst in &circuit.instructions {
                eager.apply(inst).unwrap();
            }

            let mut lazy = StabilizerBackend::new(seed);
            lazy.init(n, n).unwrap();
            lazy.enable_lazy_destab();
            let mut seen_lazy_gates = false;
            for inst in &circuit.instructions {
                if matches!(inst, Instruction::Measure { .. }) && !seen_lazy_gates {
                    assert!(
                        lazy.lazy_destab,
                        "n={n} seed={seed}: gates materialized the destabilizers"
                    );
                    seen_lazy_gates = true;
                }
                lazy.apply(inst).unwrap();
            }
            assert!(
                !lazy.lazy_destab,
                "n={n} seed={seed}: measurement left the tableau lazy"
            );

            assert_eq!(
                eager.classical_results(),
                lazy.classical_results(),
                "n={n} seed={seed}: outcome sequence diverged"
            );
        }
    }
}

// A materialization on the bulk route must not re-arm the SGI guard over an
// index the dense paths do not maintain. The shape needs every ingredient: a
// CX fan-out from qubit 0 trips the sgi_max_active high-water mark and hands
// the stream to the word-batch path, the uncompute leaves the actual supports
// light so a counter rebuild at the measurement would re-arm the guard, the
// random Bell measurement migrates the pivot row into a destabilizer slot the
// stale index does not list, and the guarded regions re-enter instruction
// dispatch where the guard is re-checked. A stale index there breaks Bell
// correlation.
#[test]
fn lazy_destab_bulk_region_after_measurement_keeps_correlation() {
    use crate::circuit::{ClassicalCondition, GuardedRegion};
    let n = 300;
    for seed in 0..40u64 {
        let mut c = Circuit::new(n, 2);
        for q in 1..=40 {
            c.add_gate(Gate::Cx, &[0, q]);
        }
        for q in 1..=40 {
            c.add_gate(Gate::Cx, &[0, q]);
        }
        c.add_gate(Gate::H, &[0]);
        c.add_gate(Gate::Cx, &[0, 1]);
        c.instructions.push(Instruction::Measure {
            qubit: 0,
            classical_bit: 0,
        });
        for condition in [
            ClassicalCondition::BitIsZero(0),
            ClassicalCondition::BitIsOne(0),
        ] {
            c.instructions
                .push(Instruction::Region(Box::new(GuardedRegion::new(
                    condition,
                    vec![Instruction::Measure {
                        qubit: 1,
                        classical_bit: 1,
                    }],
                ))));
        }

        let mut lazy = StabilizerBackend::new_lazy(seed);
        lazy.init(n, 2).unwrap();
        lazy.apply_instructions(&c.instructions).unwrap();
        let bits = lazy.classical_results();
        assert_eq!(
            bits[0], bits[1],
            "seed={seed}: Bell correlation broken on the lazy bulk route"
        );
    }
}

// Symplectic pairing after Gaussian-elimination reconstruction: destabilizer i
// anticommutes with stabilizer i and commutes with every other generator.
#[test]
fn lazy_destab_reconstruction_restores_symplectic_pairing() {
    fn anticommutes(b: &StabilizerBackend, row_a: usize, row_b: usize) -> bool {
        let nw = b.num_words;
        let stride = 2 * nw;
        let (a, bb) = (&b.xz[row_a * stride..], &b.xz[row_b * stride..]);
        let mut acc = 0u32;
        for w in 0..nw {
            acc ^= (a[w] & bb[nw + w]).count_ones() ^ (a[nw + w] & bb[w]).count_ones();
        }
        acc % 2 == 1
    }

    for n in [3usize, 5, 8, 16, 80, 300] {
        let circuit = crate::circuits::clifford_heavy_circuit(n, 10, 42);
        let mut b = StabilizerBackend::new(42);
        b.init(n, 0).unwrap();
        b.enable_lazy_destab();
        for inst in &circuit.instructions {
            b.apply(inst).unwrap();
        }
        assert!(
            b.lazy_destab,
            "n={n}: gate application must not materialize"
        );
        b.ensure_destabilizers();

        for i in 0..n {
            for j in 0..n {
                assert_eq!(
                    anticommutes(&b, i, n + j),
                    i == j,
                    "n={n}: destab {i} vs stab {j} pairing broken"
                );
                assert!(
                    !anticommutes(&b, i, j),
                    "n={n}: destabs {i},{j} must commute"
                );
                assert!(
                    !anticommutes(&b, n + i, n + j),
                    "n={n}: stabs {i},{j} must commute"
                );
            }
        }
    }
}

#[test]
fn lazy_destab_matches_eager() {
    for n in [3, 5, 10] {
        let circuit = crate::circuits::clifford_heavy_circuit(n, 10, 42);

        let mut eager = StabilizerBackend::new(42);
        eager.init(n, 0).unwrap();
        eager.apply_gates_only(&circuit.instructions).unwrap();
        let eager_probs = eager.probabilities().unwrap();

        let mut lazy = StabilizerBackend::new(42);
        lazy.init(n, 0).unwrap();
        lazy.enable_lazy_destab();
        lazy.apply_gates_only(&circuit.instructions).unwrap();
        let lazy_probs = lazy.probabilities().unwrap();

        for (i, (&e, &l)) in eager_probs.iter().zip(lazy_probs.iter()).enumerate() {
            assert!(
                (e - l).abs() < 1e-10,
                "n={n} prob mismatch at {i}: eager={e}, lazy={l}"
            );
        }
    }
}

#[test]
fn lazy_destab_measure_matches_eager() {
    fn run_measurement_shots(
        cliff: &Circuit,
        n: usize,
        num_shots: usize,
        seed: u64,
        lazy_destab: bool,
    ) -> Vec<Vec<bool>> {
        let mut shots = Vec::with_capacity(num_shots);
        for i in 0..num_shots {
            let mut backend = StabilizerBackend::new(seed.wrapping_add(i as u64));
            backend.init(n, n).unwrap();
            if lazy_destab {
                backend.enable_lazy_destab();
            }
            for inst in &cliff.instructions {
                backend.apply(inst).unwrap();
            }
            for q in 0..n {
                backend
                    .apply(&Instruction::Measure {
                        qubit: q,
                        classical_bit: q,
                    })
                    .unwrap();
            }
            shots.push(backend.classical_results().to_vec());
        }
        shots
    }

    for n in [3, 5, 8] {
        let cliff = crate::circuits::clifford_heavy_circuit(n, 10, 42);
        let num_shots = 2_000;
        let eager_shots = run_measurement_shots(&cliff, n, num_shots, 42, false);
        let lazy_shots = run_measurement_shots(&cliff, n, num_shots, 42, true);

        for q in 0..n {
            let eager_ones = eager_shots.iter().filter(|shot| shot[q]).count();
            let lazy_ones = lazy_shots.iter().filter(|shot| shot[q]).count();
            let delta =
                (eager_ones as isize - lazy_ones as isize).unsigned_abs() as f64 / num_shots as f64;
            assert!(
                delta < 0.08,
                "n={n} q={q}: lazy/eager marginal mismatch too large ({delta:.3})"
            );
        }
    }
}

#[cfg(feature = "gpu")]
mod gpu_scaffold {
    use super::*;
    use crate::gpu::GpuContext;

    // With the stub context, `GpuTableau::new` returns `BackendUnsupported` at
    // `alloc_zeros` time. `StabilizerBackend::init` must surface that error
    // cleanly rather than panicking.
    #[test]
    fn with_gpu_init_on_stub_returns_unsupported() {
        let ctx = GpuContext::stub_for_tests();
        let mut backend = StabilizerBackend::new(42).with_gpu(ctx);
        let err = backend.init(4, 0).unwrap_err();
        assert!(matches!(err, PrismError::BackendUnsupported { .. }));
    }

    // `with_gpu_auto` is the soft (fail-closed) mode. On the stub context,
    // `GpuTableau::new` fails at `init`, but instead of surfacing the error the
    // backend must degrade to the host tableau and run correctly. A GHZ chain
    // then measures deterministically-correlated bits, proving the fallback
    // produced a working CPU stabilizer rather than a corrupt half-state.
    #[test]
    fn with_gpu_auto_falls_back_to_host_on_stub() {
        use crate::gates::Gate;

        let ctx = GpuContext::stub_for_tests();
        let mut backend = StabilizerBackend::new(42).with_gpu_auto(ctx);
        backend
            .init(4, 4)
            .expect("soft mode must fall back to host, not error");
        assert_eq!(backend.num_qubits(), 4);

        let mut circuit = Circuit::new(4, 4);
        circuit.add_gate(Gate::H, &[0]);
        for i in 0..3 {
            circuit.add_gate(Gate::Cx, &[i, i + 1]);
        }
        for i in 0..4 {
            circuit.add_measure(i, i);
        }
        for inst in &circuit.instructions {
            backend.apply(inst).unwrap();
        }
        let bits = backend.classical_results();
        assert!(
            bits.iter().all(|&b| b == bits[0]),
            "GHZ measurement bits must all agree: {bits:?}"
        );
    }

    // Constructing with `with_gpu(ctx)` but then not initialising must leave
    // the backend in a usable (context-only) state.
    #[test]
    fn with_gpu_stores_context_without_panicking() {
        let ctx = GpuContext::stub_for_tests();
        let backend = StabilizerBackend::new(42).with_gpu(ctx);
        assert_eq!(backend.name(), "stabilizer");
    }

    // Transactional init: if `GpuTableau::new` fails, the backend must retain
    // its pre-init state so subsequent calls do not touch partially-updated
    // host buffers. Here the backend is fresh (`n == 0`), and the invariant
    // is that it stays that way after the failure.
    #[test]
    fn failed_gpu_init_leaves_backend_uninitialised() {
        let ctx = GpuContext::stub_for_tests();
        let mut backend = StabilizerBackend::new(42).with_gpu(ctx);
        let _ = backend.init(4, 0).unwrap_err();
        assert_eq!(backend.num_qubits(), 0);
    }

    // Cloning a CPU-only backend preserves state; this covers the common
    // `stabilizer_rank` path and confirms the GPU-gated Clone impl still
    // works for non-GPU callers.
    #[test]
    fn clone_cpu_only_preserves_tableau() {
        let mut backend = StabilizerBackend::new(42);
        backend.init(3, 0).unwrap();
        let cloned = backend.clone();
        assert_eq!(cloned.n, 3);
        assert_eq!(cloned.xz, backend.xz);
        assert_eq!(cloned.phase, backend.phase);
    }
}

fn statevector_expectations(
    prefix: &Circuit,
    project: Option<(usize, bool)>,
    suffix: &[Instruction],
    observables: &[Vec<crate::sim::unified_pauli::PauliTerm>],
) -> Vec<f64> {
    use crate::backend::statevector::StatevectorBackend;
    let mut sv = StatevectorBackend::new(42);
    sv.init(prefix.num_qubits, prefix.num_classical_bits)
        .unwrap();
    for inst in &prefix.instructions {
        sv.apply(inst).unwrap();
    }
    if let Some((qubit, outcome)) = project {
        let one = Complex64::new(1.0, 0.0);
        let zero = Complex64::new(0.0, 0.0);
        let projector = if outcome {
            [[zero, zero], [zero, one]]
        } else {
            [[one, zero], [zero, zero]]
        };
        sv.apply_1q_matrix(qubit, &projector).unwrap();
    }
    for inst in suffix {
        sv.apply(inst).unwrap();
    }
    let masks: Vec<_> = observables
        .iter()
        .map(|o| crate::sim::pauli_masks(o, prefix.num_qubits).unwrap())
        .collect();
    crate::sim::pauli_expectations_from_masks(sv.state_vector(), &masks).0
}

fn tableau_expectations(
    backend: &mut StabilizerBackend,
    circuit: &Circuit,
    observables: &[Vec<crate::sim::unified_pauli::PauliTerm>],
) -> Vec<f64> {
    sim::run_on(backend, circuit).unwrap();
    backend.pauli_expectations(observables).unwrap()
}

/// Pauli strings read off the stabilizer half of an eager tableau: each
/// generator, then the product of each adjacent pair. Every one is in the
/// group, so its expectation is +1 or -1.
fn group_strings(backend: &StabilizerBackend) -> Vec<Vec<crate::sim::unified_pauli::PauliTerm>> {
    let n = backend.num_qubits();
    let nw = n.div_ceil(64);
    let (xz, _) = backend.raw_tableau();
    let string_of = |rows: &[usize]| -> Vec<crate::sim::unified_pauli::PauliTerm> {
        (0..n)
            .filter_map(|q| {
                let bit =
                    |r: usize, half: usize| (xz[r * 2 * nw + half * nw + q / 64] >> (q % 64)) & 1;
                let x = rows.iter().fold(0, |acc, &r| acc ^ bit(r, 0)) == 1;
                let z = rows.iter().fold(0, |acc, &r| acc ^ bit(r, 1)) == 1;
                match (x, z) {
                    (true, false) => Some(crate::sim::unified_pauli::PauliTerm::x(q)),
                    (false, true) => Some(crate::sim::unified_pauli::PauliTerm::z(q)),
                    (true, true) => Some(crate::sim::unified_pauli::PauliTerm::y(q)),
                    (false, false) => None,
                }
            })
            .collect()
    };
    let mut strings: Vec<Vec<crate::sim::unified_pauli::PauliTerm>> =
        (n..2 * n).map(|r| string_of(&[r])).collect();
    strings.extend((n..2 * n - 1).map(|r| string_of(&[r, r + 1])));
    strings
}

fn assert_expectations_close(got: &[f64], want: &[f64], label: &str) {
    assert_eq!(got.len(), want.len(), "{label}");
    for (i, (g, w)) in got.iter().zip(want).enumerate() {
        assert!((g - w).abs() < 1e-12, "{label} term {i}: got {g}, want {w}");
    }
}

#[test]
fn pauli_expectations_on_basis_and_ghz_states_match_statevector() {
    use crate::sim::unified_pauli::PauliTerm;
    let zero = Circuit::new(2, 0);
    let mut plus = Circuit::new(2, 0);
    plus.add_gate(Gate::H, &[0]);
    plus.add_gate(Gate::H, &[1]);
    plus.add_gate(Gate::S, &[1]);
    let ghz = crate::circuits::ghz_circuit(4);

    let observables = vec![
        vec![],
        vec![PauliTerm::z(0)],
        vec![PauliTerm::x(0)],
        vec![PauliTerm::y(0)],
        vec![PauliTerm::y(1)],
        vec![PauliTerm::z(0), PauliTerm::z(1)],
        vec![PauliTerm::x(0), PauliTerm::x(1)],
        vec![PauliTerm::y(0), PauliTerm::y(1)],
        vec![PauliTerm::x(0), PauliTerm::y(1)],
    ];
    for (circuit, label) in [(&zero, "zero"), (&plus, "plus")] {
        let want = statevector_expectations(circuit, None, &[], &observables);
        let got = tableau_expectations(&mut StabilizerBackend::new(42), circuit, &observables);
        assert_expectations_close(&got, &want, label);
    }

    let ghz_observables = vec![
        vec![],
        vec![PauliTerm::z(0)],
        vec![PauliTerm::z(0), PauliTerm::z(3)],
        vec![PauliTerm::z(0), PauliTerm::z(1), PauliTerm::z(2)],
        vec![
            PauliTerm::x(0),
            PauliTerm::x(1),
            PauliTerm::x(2),
            PauliTerm::x(3),
        ],
        vec![
            PauliTerm::y(0),
            PauliTerm::y(1),
            PauliTerm::x(2),
            PauliTerm::x(3),
        ],
        vec![
            PauliTerm::y(0),
            PauliTerm::y(1),
            PauliTerm::y(2),
            PauliTerm::y(3),
        ],
        vec![PauliTerm::x(0), PauliTerm::x(1)],
    ];
    let want = statevector_expectations(&ghz, None, &[], &ghz_observables);
    assert_eq!(want[5], -1.0);
    let got = tableau_expectations(&mut StabilizerBackend::new(42), &ghz, &ghz_observables);
    assert_expectations_close(&got, &want, "ghz");
    let lazy = tableau_expectations(&mut StabilizerBackend::new_lazy(42), &ghz, &ghz_observables);
    assert_expectations_close(&lazy, &want, "ghz lazy");
}

#[test]
fn pauli_expectations_on_random_clifford_match_statevector() {
    use crate::sim::unified_pauli::{PauliAxis, PauliTerm};
    use rand::RngExt;
    let n = 10;
    let circuit = crate::circuits::clifford_random_pairs(n, 6, 0xDEAD_BEEF);
    let mut rng = ChaCha8Rng::seed_from_u64(7);
    let axes = [PauliAxis::X, PauliAxis::Y, PauliAxis::Z];
    let observables: Vec<Vec<PauliTerm>> = (0..40)
        .map(|_| {
            (0..n)
                .filter_map(|q| {
                    axes.get(rng.random_range(0..4))
                        .map(|&axis| PauliTerm::new(q, axis))
                })
                .collect()
        })
        .collect();
    let mut eager = StabilizerBackend::new(42);
    sim::run_on(&mut eager, &circuit).unwrap();
    let mut observables = observables;
    observables.extend(group_strings(&eager));
    let want = statevector_expectations(&circuit, None, &[], &observables);
    assert!(want[40..].iter().all(|w| w.abs() == 1.0));
    assert!(want[..40].contains(&0.0));
    assert!(want[40..].contains(&-1.0));
    let got = eager.pauli_expectations(&observables).unwrap();
    assert_expectations_close(&got, &want, "eager");
    let lazy = tableau_expectations(&mut StabilizerBackend::new_lazy(42), &circuit, &observables);
    assert_expectations_close(&lazy, &want, "lazy");
}

#[test]
fn pauli_expectation_of_a_term_anticommuting_with_a_generator_is_zero() {
    use crate::sim::unified_pauli::PauliTerm;
    let mut plus = Circuit::new(1, 0);
    plus.add_gate(Gate::H, &[0]);
    let got = tableau_expectations(
        &mut StabilizerBackend::new(42),
        &plus,
        &[vec![PauliTerm::z(0)], vec![PauliTerm::y(0)]],
    );
    assert_eq!(got, vec![0.0, 0.0]);

    let ghz = crate::circuits::ghz_circuit(3);
    let got = tableau_expectations(
        &mut StabilizerBackend::new(42),
        &ghz,
        &[
            vec![PauliTerm::x(0)],
            vec![PauliTerm::z(1), PauliTerm::y(2)],
        ],
    );
    assert_eq!(got, vec![0.0, 0.0]);
}

#[test]
fn pauli_expectations_after_mid_circuit_measurement_match_statevector() {
    use crate::sim::unified_pauli::PauliTerm;
    let prefix = crate::circuits::ghz_circuit(3);
    let mut circuit = Circuit::new(3, 1);
    for inst in &prefix.instructions {
        circuit.instructions.push(inst.clone());
    }
    circuit.add_measure(0, 0);
    circuit.add_gate(Gate::H, &[1]);
    circuit.add_gate(Gate::Cx, &[1, 2]);
    let suffix = circuit.instructions[prefix.instructions.len() + 1..].to_vec();

    let observables = vec![
        vec![PauliTerm::z(0)],
        vec![PauliTerm::z(1)],
        vec![PauliTerm::x(1)],
        vec![PauliTerm::x(1), PauliTerm::x(2)],
        vec![PauliTerm::z(1), PauliTerm::z(2)],
        vec![PauliTerm::z(0), PauliTerm::x(1), PauliTerm::x(2)],
    ];
    let mut stab = StabilizerBackend::new(42);
    let got = tableau_expectations(&mut stab, &circuit, &observables);
    let outcome = stab.classical_results()[0];
    let want = statevector_expectations(&prefix, Some((0, outcome)), &suffix, &observables);
    assert_expectations_close(&got, &want, "post-measurement");
    assert_eq!(got[0], if outcome { -1.0 } else { 1.0 });
}

#[test]
fn pauli_expectations_span_words_on_a_wide_register() {
    use crate::sim::unified_pauli::PauliTerm;
    let n = 130;
    let ghz = crate::circuits::ghz_circuit(n);
    let mut y_pair: Vec<PauliTerm> = (2..n).map(PauliTerm::x).collect();
    y_pair.push(PauliTerm::y(0));
    y_pair.push(PauliTerm::y(1));
    let observables = vec![
        vec![PauliTerm::z(0), PauliTerm::z(n - 1)],
        (0..n).map(PauliTerm::x).collect(),
        y_pair,
        vec![PauliTerm::z(64)],
        vec![PauliTerm::z(0), PauliTerm::z(65), PauliTerm::z(129)],
    ];
    let got = tableau_expectations(&mut StabilizerBackend::new(42), &ghz, &observables);
    assert_eq!(got, vec![1.0, 1.0, -1.0, 0.0, 0.0]);
}

#[test]
fn pauli_expectations_reject_out_of_range_and_duplicate_factors() {
    use crate::sim::unified_pauli::PauliTerm;
    let mut b = StabilizerBackend::new(42);
    b.init(3, 0).unwrap();
    assert!(matches!(
        b.pauli_expectations(&[vec![PauliTerm::z(3)]]),
        Err(PrismError::InvalidQubit {
            index: 3,
            register_size: 3
        })
    ));
    assert!(matches!(
        b.pauli_expectations(&[vec![PauliTerm::z(1), PauliTerm::x(1)]]),
        Err(PrismError::InvalidParameter { .. })
    ));
}

// `stabilizer/random_pairs` prices the batched cross-word kernel, which lives
// on the word-batch route. The sparse gate index covers a run while the rows
// stay light and hands the rest of the circuit over once they do not, so which
// widths reach that kernel is set by this predicate and by nothing a result
// carries. The bench sweep stops at 1000 for the same reason, and the wider
// widths are checked here so the omission does not rest on an untested claim.
#[test]
fn random_pairs_bench_widths_split_on_the_sparse_gate_index() {
    // Lazy destabilizers, as the plan picks for a circuit that never measures.
    fn index_covers(n: usize) -> bool {
        let circuit = crate::circuits::clifford_random_pairs(n, 10, 0xDEAD_BEEF);
        let mut backend = StabilizerBackend::new_lazy(42);
        backend.init(n, 0).unwrap();
        backend.apply_instructions(&circuit.instructions).unwrap();
        backend.sgi_ready()
    }

    for n in [500usize, 1000] {
        assert!(
            !index_covers(n),
            "random_pairs/{n}: the sparse gate index still covers the run at \
             its last gate, so the row never reaches the word-batch route"
        );
    }
    for n in [2000usize, 5000] {
        assert!(
            index_covers(n),
            "random_pairs/{n}: the index now yields at this width, so the row \
             the bench sweep omits would reach the cross-word kernel after all"
        );
    }
}