add: complex numbers

10 files changed, 1615 insertions, 0 deletions

diff --git a/tester/benchmarks.rs b/tester/benchmarks.rs
new file mode 100644
index 0000000..7aea535
--- /dev/null
+++ b/tester/benchmarks.rs
@@ -0,0 +1,90 @@
+use crate::common::{benchmark_circuit, print_section, BenchmarkResult};
+use psi::QuantumCircuit;
+use psi::HorizontalRenderer;
+
+pub fn run_all(results: &mut Vec<BenchmarkResult>) {
+    println!("═══════════════════════════════════════════════════════════════");
+    println!("                    BENCHMARK CIRCUITS");
+    println!("═══════════════════════════════════════════════════════════════\n");
+
+    test_8_qubit(results);
+    test_10_qubit(results);
+    test_12_qubit(results);
+    test_14_qubit(results);
+}
+
+pub fn test_8_qubit(results: &mut Vec<BenchmarkResult>) {
+    print_section("8-qubit Entangled Circuit");
+
+    let builder = || {
+        let mut circuit = QuantumCircuit::new(8);
+        for i in 0..8 {
+            circuit.h(i);
+        }
+        for i in 0..7 {
+            circuit.cnot(i, i + 1);
+        }
+        circuit
+    };
+
+    println!("{}", HorizontalRenderer::new(&builder()));
+    results.push(benchmark_circuit("8-qubit entangled", builder));
+}
+
+pub fn test_10_qubit(results: &mut Vec<BenchmarkResult>) {
+    print_section("10-qubit Entangled Circuit");
+
+    let builder = || {
+        let mut circuit = QuantumCircuit::new(10);
+        for i in 0..10 {
+            circuit.h(i);
+        }
+        for i in 0..9 {
+            circuit.cnot(i, i + 1);
+        }
+        circuit.cz(0, 9);
+        circuit
+    };
+
+    println!("{}", HorizontalRenderer::new(&builder()));
+    results.push(benchmark_circuit("10-qubit entangled", builder));
+}
+
+pub fn test_12_qubit(results: &mut Vec<BenchmarkResult>) {
+    print_section("12-qubit Entangled Circuit");
+
+    let builder = || {
+        let mut circuit = QuantumCircuit::new(12);
+        for i in 0..12 {
+            circuit.h(i);
+        }
+        for i in 0..11 {
+            circuit.cnot(i, i + 1);
+        }
+        circuit.cz(0, 11);
+        circuit.swap(5, 6);
+        circuit
+    };
+
+    println!("{}", HorizontalRenderer::new(&builder()));
+    results.push(benchmark_circuit("12-qubit entangled", builder));
+}
+
+pub fn test_14_qubit(results: &mut Vec<BenchmarkResult>) {
+    print_section("14-qubit Entangled Circuit");
+
+    let builder = || {
+        let mut circuit = QuantumCircuit::new(14);
+        for i in 0..14 {
+            circuit.h(i);
+        }
+        for i in 0..13 {
+            circuit.cnot(i, i + 1);
+        }
+        circuit
+    };
+
+    println!("{}", HorizontalRenderer::new(&builder()));
+    results.push(benchmark_circuit("14-qubit entangled", builder));
+}
+
diff --git a/tester/clifford.rs b/tester/clifford.rs
new file mode 100644
index 0000000..932a918
--- /dev/null
+++ b/tester/clifford.rs
@@ -0,0 +1,126 @@
+use crate::common::{benchmark_circuit, print_circuit, print_section, BenchmarkResult};
+use psi::QuantumCircuit;
+
+pub fn run_all(results: &mut Vec<BenchmarkResult>) {
+    println!("═══════════════════════════════════════════════════════════════");
+    println!("                     CLIFFORD GATES TESTS");
+    println!("═══════════════════════════════════════════════════════════════\n");
+
+    test_bell_state(results);
+    test_ghz_state(results);
+    test_swap_via_cnots(results);
+    test_toffoli(results);
+    test_hadamard_measure(results);
+    test_complex_circuit(results);
+}
+
+pub fn test_bell_state(results: &mut Vec<BenchmarkResult>) {
+    print_section("Bell State with Measurement");
+
+    let builder = || {
+        let mut circuit = QuantumCircuit::with_classical(2, 2);
+        circuit.h(0).cnot(0, 1).measure(0, 0).measure(1, 1);
+        circuit
+    };
+
+    print_circuit(&builder());
+    results.push(benchmark_circuit("Bell State (2 qubits)", builder));
+
+    let mut display = builder();
+    display.compute();
+    println!("{}\n", display);
+}
+
+pub fn test_ghz_state(results: &mut Vec<BenchmarkResult>) {
+    print_section("GHZ State");
+
+    let builder = || {
+        let mut circuit = QuantumCircuit::new(3);
+        circuit.h(0).cnot(0, 1).cnot(0, 2);
+        circuit
+    };
+
+    print_circuit(&builder());
+    results.push(benchmark_circuit("GHZ State (3 qubits)", builder));
+
+    let mut display = builder();
+    display.compute();
+    println!("{}\n", display);
+}
+
+pub fn test_swap_via_cnots(results: &mut Vec<BenchmarkResult>) {
+    print_section("SWAP via 3 CNOTs");
+
+    let builder = || {
+        let mut circuit = QuantumCircuit::new(2);
+        circuit.x(0).cnot(0, 1).cnot(1, 0).cnot(0, 1);
+        circuit
+    };
+
+    print_circuit(&builder());
+    results.push(benchmark_circuit("SWAP via CNOTs (2 qubits)", builder));
+
+    let mut display = builder();
+    display.compute();
+    println!("{}\n", display);
+}
+
+pub fn test_toffoli(results: &mut Vec<BenchmarkResult>) {
+    print_section("Toffoli Gate");
+
+    let builder = || {
+        let mut circuit = QuantumCircuit::new(3);
+        circuit.x(0).x(1).toffoli(0, 1, 2);
+        circuit
+    };
+
+    print_circuit(&builder());
+    results.push(benchmark_circuit("Toffoli (3 qubits)", builder));
+
+    let mut display = builder();
+    display.compute();
+    println!("{}\n", display);
+}
+
+pub fn test_hadamard_measure(results: &mut Vec<BenchmarkResult>) {
+    print_section("Full Circuit with Measurements");
+
+    let builder = || {
+        let mut circuit = QuantumCircuit::with_classical(3, 3);
+        circuit.h(0).h(1).h(2).measure_all();
+        circuit
+    };
+
+    print_circuit(&builder());
+    results.push(benchmark_circuit("3-qubit Hadamard + Measure", builder));
+
+    let mut display = builder();
+    display.compute();
+    println!("{}\n", display);
+}
+
+pub fn test_complex_circuit(results: &mut Vec<BenchmarkResult>) {
+    print_section("Complex Circuit");
+
+    let builder = || {
+        let mut circuit = QuantumCircuit::with_classical(4, 2);
+        circuit
+            .h(0)
+            .h(1)
+            .cnot(0, 2)
+            .cnot(1, 3)
+            .cz(2, 3)
+            .swap(0, 1)
+            .measure(0, 0)
+            .measure(1, 1);
+        circuit
+    };
+
+    print_circuit(&builder());
+    results.push(benchmark_circuit("Complex (4 qubits)", builder));
+
+    let mut display = builder();
+    display.compute();
+    println!("{}\n", display);
+}
+
diff --git a/tester/common.rs b/tester/common.rs
new file mode 100644
index 0000000..9e576da
--- /dev/null
+++ b/tester/common.rs
@@ -0,0 +1,215 @@
+use psi::{QuantumCircuit, QuantumState, Runtime, Vector};
+use psi::{HorizontalRenderer, VerticalRenderer};
+use std::time::{Duration, Instant};
+
+/// A named list of circuit builders used by the benchmark/test suites.
+pub type CircuitCases = Vec<(&'static str, Box<dyn Fn() -> QuantumCircuit>)>;
+
+pub struct BenchmarkResult {
+    pub name: String,
+    pub basic_time: Duration,
+    pub mt_time: Duration,
+    pub results_match: bool,
+}
+
+pub fn benchmark_circuit<F>(name: &str, circuit_builder: F) -> BenchmarkResult
+where
+    F: Fn() -> QuantumCircuit,
+{
+    let mut circuit_st = circuit_builder();
+    let mut circuit_mt = circuit_builder();
+
+    let start_st = Instant::now();
+    circuit_st.compute_with(Runtime::BasicRT);
+    let basic_time = start_st.elapsed();
+
+    let start_mt = Instant::now();
+    circuit_mt.compute_with(Runtime::BasicRTMT);
+    let mt_time = start_mt.elapsed();
+
+    let state_st = circuit_st.state();
+    let state_mt = circuit_mt.state();
+
+    let results_match = states_equal(state_st, state_mt);
+
+    BenchmarkResult {
+        name: name.to_string(),
+        basic_time,
+        mt_time,
+        results_match,
+    }
+}
+
+pub fn states_equal(a: &QuantumState, b: &QuantumState) -> bool {
+    if a.size() != b.size() {
+        return false;
+    }
+    for i in 0..a.size() {
+        let amp_a = a.get(i);
+        let amp_b = b.get(i);
+        let diff_real = (amp_a.real - amp_b.real).abs();
+        let diff_imag = (amp_a.imaginary - amp_b.imaginary).abs();
+        if diff_real > 1e-10 || diff_imag > 1e-10 {
+            return false;
+        }
+    }
+    true
+}
+
+pub fn format_duration(d: Duration) -> String {
+    if d.as_secs() > 0 {
+        format!("{:.3}s", d.as_secs_f64())
+    } else if d.as_millis() > 0 {
+        format!("{:.3}ms", d.as_secs_f64() * 1000.0)
+    } else {
+        format!("{:.3}μs", d.as_secs_f64() * 1_000_000.0)
+    }
+}
+
+pub fn print_section(title: &str) {
+    let width = 61;
+    let padding = width - title.len() - 2;
+    println!("┌{}┐", "─".repeat(width));
+    println!("│ {}{} │", title, " ".repeat(padding));
+    println!("└{}┘\n", "─".repeat(width));
+}
+
+pub fn print_circuit(circuit: &QuantumCircuit) {
+    println!("Horizontal:\n{}", HorizontalRenderer::new(circuit));
+    println!("Vertical:\n{}", VerticalRenderer::new(circuit));
+}
+
+pub fn print_benchmark_table(results: &[BenchmarkResult]) {
+    if results.is_empty() {
+        return;
+    }
+
+    let headers = ["Circuit", "BasicRT", "BasicRTMT", "Speedup", "Match"];
+
+    let formatted: Vec<(String, String, String, String, String)> = results
+        .iter()
+        .map(|r| {
+            let speedup = r.basic_time.as_secs_f64() / r.mt_time.as_secs_f64();
+            (
+                r.name.clone(),
+                format_duration(r.basic_time),
+                format_duration(r.mt_time),
+                if speedup.is_finite() {
+                    format!("{:.2}x", speedup)
+                } else {
+                    "N/A".to_string()
+                },
+                if r.results_match { "✓" } else { "✗" }.to_string(),
+            )
+        })
+        .collect();
+
+    let c1 = formatted
+        .iter()
+        .map(|r| r.0.len())
+        .max()
+        .unwrap()
+        .max(headers[0].len());
+    let c2 = formatted
+        .iter()
+        .map(|r| r.1.len())
+        .max()
+        .unwrap()
+        .max(headers[1].len());
+    let c3 = formatted
+        .iter()
+        .map(|r| r.2.len())
+        .max()
+        .unwrap()
+        .max(headers[2].len());
+    let c4 = formatted
+        .iter()
+        .map(|r| r.3.len())
+        .max()
+        .unwrap()
+        .max(headers[3].len());
+    let c5 = formatted
+        .iter()
+        .map(|r| r.4.chars().count())
+        .max()
+        .unwrap()
+        .max(headers[4].len());
+
+    let top = format!(
+        "╔{}═{}═{}═{}═{}╗",
+        "═".repeat(c1 + 2),
+        "═".repeat(c2 + 2),
+        "═".repeat(c3 + 2),
+        "═".repeat(c4 + 2),
+        "═".repeat(c5 + 2)
+    );
+    let title_sep = format!(
+        "╠{}╤{}╤{}╤{}╤{}╣",
+        "═".repeat(c1 + 2),
+        "═".repeat(c2 + 2),
+        "═".repeat(c3 + 2),
+        "═".repeat(c4 + 2),
+        "═".repeat(c5 + 2)
+    );
+    let header_sep = format!(
+        "╠{}╪{}╪{}╪{}╪{}╣",
+        "═".repeat(c1 + 2),
+        "═".repeat(c2 + 2),
+        "═".repeat(c3 + 2),
+        "═".repeat(c4 + 2),
+        "═".repeat(c5 + 2)
+    );
+    let bottom = format!(
+        "╚{}╧{}╧{}╧{}╧{}╝",
+        "═".repeat(c1 + 2),
+        "═".repeat(c2 + 2),
+        "═".repeat(c3 + 2),
+        "═".repeat(c4 + 2),
+        "═".repeat(c5 + 2)
+    );
+
+    let total_width = c1 + c2 + c3 + c4 + c5 + 14;
+
+    println!("\n{}", top);
+    println!(
+        "║{:^width$}║",
+        "RUNTIME BENCHMARK RESULTS",
+        width = total_width
+    );
+    println!("{}", title_sep);
+    println!(
+        "║ {:<c1$} │ {:^c2$} │ {:^c3$} │ {:^c4$} │ {:^c5$} ║",
+        headers[0], headers[1], headers[2], headers[3], headers[4],
+    );
+    println!("{}", header_sep);
+
+    for (name, basic, mt, speedup, matched) in &formatted {
+        println!(
+            "║ {:<c1$} │ {:>c2$} │ {:>c3$} │ {:>c4$} │ {:^c5$} ║",
+            name, basic, mt, speedup, matched,
+        );
+    }
+
+    println!("{}", bottom);
+}
+
+pub fn print_summary(results: &[BenchmarkResult]) {
+    let all_match = results.iter().all(|r| r.results_match);
+    println!("\n");
+    if all_match {
+        println!("✓ All circuits produced identical results with both runtimes!");
+    } else {
+        println!("✗ WARNING: Some circuits produced different results!");
+    }
+
+    let total_basic: Duration = results.iter().map(|r| r.basic_time).sum();
+    let total_mt: Duration = results.iter().map(|r| r.mt_time).sum();
+    let overall_speedup = total_basic.as_secs_f64() / total_mt.as_secs_f64();
+
+    println!(
+        "\nTotal time - BasicRT: {} | BasicRTMT: {} | Overall speedup: {:.2}x",
+        format_duration(total_basic),
+        format_duration(total_mt),
+        overall_speedup
+    );
+}
diff --git a/tester/custom_gates.rs b/tester/custom_gates.rs
new file mode 100644
index 0000000..ffe5c01
--- /dev/null
+++ b/tester/custom_gates.rs
@@ -0,0 +1,121 @@
+use crate::common::{benchmark_circuit, print_circuit, print_section, BenchmarkResult};
+use psi::{complex, matrix, CustomGate, CustomGateBuilder, QuantumCircuit};
+
+pub fn run_all(results: &mut Vec<BenchmarkResult>) {
+    println!("═══════════════════════════════════════════════════════════════");
+    println!("                    CUSTOM GATES TESTS");
+    println!("═══════════════════════════════════════════════════════════════\n");
+
+    test_bell_gate(results);
+    test_swap_gate(results);
+    test_sqrt_x_gate(results);
+}
+
+pub fn test_bell_gate(results: &mut Vec<BenchmarkResult>) {
+    print_section("Custom Gate: Bell Pair Creator");
+
+    let bell_gate = CustomGateBuilder::new("BELL", 2).h(0).cnot(0, 1).build();
+    let gate_clone = bell_gate.clone();
+
+    let builder = move || {
+        let mut circuit = QuantumCircuit::new(4);
+        circuit
+            .apply_custom(gate_clone.clone(), &[0, 1])
+            .apply_custom(gate_clone.clone(), &[2, 3]);
+        circuit
+    };
+
+    let display_circuit = {
+        let mut circuit = QuantumCircuit::new(4);
+        circuit
+            .apply_custom(bell_gate.clone(), &[0, 1])
+            .apply_custom(bell_gate.clone(), &[2, 3]);
+        circuit
+    };
+    print_circuit(&display_circuit);
+    results.push(benchmark_circuit("Custom BELL (4 qubits)", builder));
+
+    let mut display = {
+        let mut circuit = QuantumCircuit::new(4);
+        circuit
+            .apply_custom(bell_gate.clone(), &[0, 1])
+            .apply_custom(bell_gate.clone(), &[2, 3]);
+        circuit
+    };
+    display.compute();
+    println!("{}\n", display);
+}
+
+pub fn test_swap_gate(results: &mut Vec<BenchmarkResult>) {
+    print_section("Custom Gate: Swap via CNOTs");
+
+    let swap_gate = CustomGateBuilder::new("MYSWAP", 2)
+        .cnot(0, 1)
+        .cnot(1, 0)
+        .cnot(0, 1)
+        .build();
+    let gate_clone = swap_gate.clone();
+
+    let builder = move || {
+        let mut circuit = QuantumCircuit::new(2);
+        circuit.x(0).apply_custom(gate_clone.clone(), &[0, 1]);
+        circuit
+    };
+
+    let display_circuit = {
+        let mut circuit = QuantumCircuit::new(2);
+        circuit.x(0).apply_custom(swap_gate.clone(), &[0, 1]);
+        circuit
+    };
+    print_circuit(&display_circuit);
+    results.push(benchmark_circuit("Custom SWAP (2 qubits)", builder));
+
+    let mut display = {
+        let mut circuit = QuantumCircuit::new(2);
+        circuit.x(0).apply_custom(swap_gate.clone(), &[0, 1]);
+        circuit
+    };
+    display.compute();
+    println!("{}\n", display);
+}
+
+pub fn test_sqrt_x_gate(results: &mut Vec<BenchmarkResult>) {
+    print_section("Custom Gate: Matrix-defined √X gate");
+
+    let sqrt_x_matrix = matrix!(
+        [complex!(0.5, 0.5), complex!(0.5, -0.5)];
+        [complex!(0.5, -0.5), complex!(0.5, 0.5)]
+    );
+    let sqrt_x = CustomGate::from_matrix("√X", sqrt_x_matrix);
+    let gate_clone = sqrt_x.clone();
+
+    let builder = move || {
+        let mut circuit = QuantumCircuit::new(1);
+        circuit
+            .apply_custom(gate_clone.clone(), &[0])
+            .apply_custom(gate_clone.clone(), &[0]);
+        circuit
+    };
+
+    let display_circuit = {
+        let mut circuit = QuantumCircuit::new(1);
+        circuit
+            .apply_custom(sqrt_x.clone(), &[0])
+            .apply_custom(sqrt_x.clone(), &[0]);
+        circuit
+    };
+    print_circuit(&display_circuit);
+    results.push(benchmark_circuit("√X gate (1 qubit)", builder));
+
+    let mut display = {
+        let mut circuit = QuantumCircuit::new(1);
+        circuit
+            .apply_custom(sqrt_x.clone(), &[0])
+            .apply_custom(sqrt_x.clone(), &[0]);
+        circuit
+    };
+    display.compute();
+    println!("{}", display);
+    println!("(Two √X gates should equal X, so |0⟩ becomes |1⟩)\n");
+}
+
diff --git a/tester/kernels.rs b/tester/kernels.rs
new file mode 100644
index 0000000..b42b181
--- /dev/null
+++ b/tester/kernels.rs
@@ -0,0 +1,420 @@
+use crate::common::{print_section, states_equal, BenchmarkResult, CircuitCases};
+use psi::{QuantumCircuit, Runtime, RuntimeConfig};
+use std::f64::consts::PI;
+use std::time::Instant;
+
+pub fn run_all(results: &mut Vec<BenchmarkResult>) {
+    println!("═══════════════════════════════════════════════════════════════");
+    println!("                    KERNEL BATCHING TESTS");
+    println!("═══════════════════════════════════════════════════════════════\n");
+
+    test_kernel_fusion(results);
+    test_batched_vs_basic(results);
+    test_batched_large_circuits(results);
+    test_structure_aware(results);
+    test_composable_runtime(results);
+}
+
+pub fn test_kernel_fusion(results: &mut Vec<BenchmarkResult>) {
+    print_section("Kernel Fusion Test");
+
+    let builder = || {
+        let mut circuit = QuantumCircuit::new(2);
+        circuit.h(0).t(0).s(0).x(0).h(1).z(1);
+        circuit
+    };
+
+    let circuit = builder();
+    let batch = Runtime::build_kernel_batch(2, circuit.operations());
+    let original_count = batch.len();
+    println!("Original kernels: {}", original_count);
+    for (i, k) in batch.kernels().iter().enumerate() {
+        println!("  {}: {} on {:?}", i, k.name, k.targets);
+    }
+
+    let mut optimized_batch = Runtime::build_kernel_batch(2, circuit.operations());
+    optimized_batch.optimize();
+    let optimized_count = optimized_batch.len();
+    println!("\nOptimized kernels: {}", optimized_count);
+    for (i, k) in optimized_batch.kernels().iter().enumerate() {
+        println!("  {}: {} on {:?}", i, k.name, k.targets);
+    }
+
+    let reduction = ((original_count - optimized_count) as f64 / original_count as f64) * 100.0;
+    println!(
+        "\nKernel reduction: {} → {} ({:.0}% fewer)",
+        original_count, optimized_count, reduction
+    );
+
+    let mut basic = builder();
+    let start = Instant::now();
+    basic.compute_with(Runtime::BasicRT);
+    let basic_time = start.elapsed();
+
+    let mut batched = builder();
+    let start = Instant::now();
+    batched.compute_with(Runtime::BatchedRT);
+    let batched_time = start.elapsed();
+
+    let match_result = states_equal(basic.state(), batched.state());
+    println!("Results match: {}\n", if match_result { "✓" } else { "✗" });
+
+    results.push(BenchmarkResult {
+        name: format!("Fusion ({}→{} kernels)", original_count, optimized_count),
+        basic_time,
+        mt_time: batched_time,
+        results_match: match_result,
+    });
+
+    let fusion_heavy = || {
+        let mut circuit = QuantumCircuit::new(1);
+        circuit.h(0).t(0).s(0).x(0).y(0).z(0).h(0).t(0);
+        circuit
+    };
+
+    let circuit2 = fusion_heavy();
+    let batch2 = Runtime::build_kernel_batch(1, circuit2.operations());
+    let orig2 = batch2.len();
+    let mut opt_batch2 = Runtime::build_kernel_batch(1, circuit2.operations());
+    opt_batch2.optimize();
+    let opt2 = opt_batch2.len();
+
+    let mut basic2 = fusion_heavy();
+    let start = Instant::now();
+    basic2.compute_with(Runtime::BasicRT);
+    let basic_time2 = start.elapsed();
+
+    let mut batched2 = fusion_heavy();
+    let start = Instant::now();
+    batched2.compute_with(Runtime::BatchedRT);
+    let batched_time2 = start.elapsed();
+
+    let match2 = states_equal(basic2.state(), batched2.state());
+
+    results.push(BenchmarkResult {
+        name: format!("Heavy fusion ({}→{} kernels)", orig2, opt2),
+        basic_time: basic_time2,
+        mt_time: batched_time2,
+        results_match: match2,
+    });
+}
+
+pub fn test_batched_vs_basic(results: &mut Vec<BenchmarkResult>) {
+    print_section("Batched vs Basic Runtime Comparison");
+
+    let test_cases: CircuitCases = vec![
+        (
+            "Bell State",
+            Box::new(|| {
+                let mut c = QuantumCircuit::new(2);
+                c.h(0).cnot(0, 1);
+                c
+            }),
+        ),
+        (
+            "GHZ State",
+            Box::new(|| {
+                let mut c = QuantumCircuit::new(3);
+                c.h(0).cnot(0, 1).cnot(0, 2);
+                c
+            }),
+        ),
+        (
+            "Rotation Chain",
+            Box::new(|| {
+                let mut c = QuantumCircuit::new(3);
+                c.rx(0, PI / 4.0)
+                    .ry(0, PI / 4.0)
+                    .rz(0, PI / 4.0)
+                    .rx(1, PI / 3.0)
+                    .ry(1, PI / 3.0);
+                c
+            }),
+        ),
+        (
+            "Mixed Gates",
+            Box::new(|| {
+                let mut c = QuantumCircuit::new(4);
+                c.h(0).h(1).h(2).h(3).cnot(0, 1).cnot(2, 3).cz(1, 2);
+                c
+            }),
+        ),
+    ];
+
+    for (name, builder) in test_cases {
+        let mut basic = builder();
+        let start = Instant::now();
+        basic.compute_with(Runtime::BasicRT);
+        let basic_time = start.elapsed();
+
+        let mut batched = builder();
+        let start = Instant::now();
+        batched.compute_with(Runtime::BatchedRT);
+        let batched_time = start.elapsed();
+
+        let match_result = states_equal(basic.state(), batched.state());
+
+        println!(
+            "{}: Basic={:.2}μs, Batched={:.2}μs, Match={}",
+            name,
+            basic_time.as_secs_f64() * 1_000_000.0,
+            batched_time.as_secs_f64() * 1_000_000.0,
+            if match_result { "✓" } else { "✗" }
+        );
+
+        results.push(BenchmarkResult {
+            name: format!("Batched: {}", name),
+            basic_time,
+            mt_time: batched_time,
+            results_match: match_result,
+        });
+    }
+    println!();
+}
+
+pub fn test_batched_large_circuits(results: &mut Vec<BenchmarkResult>) {
+    print_section("Batched Runtime on Large Circuits");
+
+    let sizes = [8, 10, 12];
+
+    for &n in &sizes {
+        let builder = || {
+            let mut circuit = QuantumCircuit::new(n);
+            for i in 0..n {
+                circuit.h(i);
+            }
+            for i in 0..(n - 1) {
+                circuit.cnot(i, i + 1);
+            }
+            circuit
+        };
+
+        let mut basic_mt = builder();
+        let start = Instant::now();
+        basic_mt.compute_with(Runtime::BasicRTMT);
+        let basic_mt_time = start.elapsed();
+
+        let mut batched_mt = builder();
+        let start = Instant::now();
+        batched_mt.compute_with(Runtime::BatchedRTMT);
+        let batched_mt_time = start.elapsed();
+
+        let match_result = states_equal(basic_mt.state(), batched_mt.state());
+
+        println!(
+            "{}-qubit: BasicRTMT={:.3}ms, BatchedRTMT={:.3}ms, Match={}",
+            n,
+            basic_mt_time.as_secs_f64() * 1000.0,
+            batched_mt_time.as_secs_f64() * 1000.0,
+            if match_result { "✓" } else { "✗" }
+        );
+
+        results.push(BenchmarkResult {
+            name: format!("{}-qubit batched", n),
+            basic_time: basic_mt_time,
+            mt_time: batched_mt_time,
+            results_match: match_result,
+        });
+    }
+    println!();
+}
+
+pub fn test_structure_aware(results: &mut Vec<BenchmarkResult>) {
+    print_section("Structure-Aware Kernel Optimisation");
+
+    let commute_test = || {
+        let mut c = QuantumCircuit::new(3);
+        c.t(0).h(1).t(0).h(2).s(0).t(1).rz(0, PI / 4.0);
+        c
+    };
+
+    let circuit = commute_test();
+    let mut batch = Runtime::build_structure_aware_batch(3, circuit.operations());
+    let original = batch.len();
+    println!("Original operations: {}", original);
+    for (i, k) in batch.kernels().iter().enumerate() {
+        println!("  {}: {} on {:?} ({:?})", i, k.name, k.targets, k.gate_type);
+    }
+
+    batch.optimise();
+    let optimised = batch.len();
+    println!("\nAfter optimisation: {}", optimised);
+    for (i, k) in batch.kernels().iter().enumerate() {
+        println!("  {}: {} on {:?}", i, k.name, k.targets);
+    }
+
+    println!("\nExecution layers: {}", batch.num_layers());
+    for (i, layer) in batch.layers().iter().enumerate() {
+        let names: Vec<_> = layer.kernels.iter().map(|k| k.name.as_str()).collect();
+        println!("  Layer {}: {:?}", i, names);
+    }
+
+    let stats = batch.stats();
+    println!("\nStats: {}", stats);
+
+    let mut basic = commute_test();
+    let start = Instant::now();
+    basic.compute_with(Runtime::BasicRT);
+    let basic_time = start.elapsed();
+
+    let mut sa = commute_test();
+    let start = Instant::now();
+    sa.compute_with(Runtime::StructureAwareRT);
+    let sa_time = start.elapsed();
+
+    let match_result = states_equal(basic.state(), sa.state());
+    println!(
+        "\nBasic={:.2}μs, StructureAware={:.2}μs, Match={}",
+        basic_time.as_secs_f64() * 1_000_000.0,
+        sa_time.as_secs_f64() * 1_000_000.0,
+        if match_result { "✓" } else { "✗" }
+    );
+
+    results.push(BenchmarkResult {
+        name: format!("SA: Commuting ({}→{})", original, optimised),
+        basic_time,
+        mt_time: sa_time,
+        results_match: match_result,
+    });
+
+    println!();
+    print_section("Structure-Aware vs Other Runtimes");
+
+    let test_cases: CircuitCases = vec![
+        (
+            "Diagonal-heavy (5q)",
+            Box::new(|| {
+                let mut c = QuantumCircuit::new(5);
+                for q in 0..5 {
+                    c.t(q).s(q).rz(q, PI / 4.0).t(q);
+                }
+                c
+            }),
+        ),
+        (
+            "Interleaved (4q)",
+            Box::new(|| {
+                let mut c = QuantumCircuit::new(4);
+                c.h(0).h(1).h(2).h(3);
+                c.t(0).t(1).t(2).t(3);
+                c.cnot(0, 1).cnot(2, 3);
+                c.s(0).s(1).s(2).s(3);
+                c
+            }),
+        ),
+        (
+            "Deep rotation (3q)",
+            Box::new(|| {
+                let mut c = QuantumCircuit::new(3);
+                for _ in 0..5 {
+                    for q in 0..3 {
+                        c.rx(q, PI / 8.0).ry(q, PI / 8.0).rz(q, PI / 8.0);
+                    }
+                }
+                c
+            }),
+        ),
+    ];
+
+    for (name, builder) in test_cases {
+        let mut batched = builder();
+        let start = Instant::now();
+        batched.compute_with(Runtime::BatchedRT);
+        let batched_time = start.elapsed();
+
+        let mut sa = builder();
+        let start = Instant::now();
+        sa.compute_with(Runtime::StructureAwareRT);
+        let sa_time = start.elapsed();
+
+        let match_result = states_equal(batched.state(), sa.state());
+
+        let speedup = batched_time.as_secs_f64() / sa_time.as_secs_f64();
+        println!(
+            "{}: Batched={:.2}μs, SA={:.2}μs, Speedup={:.2}x, Match={}",
+            name,
+            batched_time.as_secs_f64() * 1_000_000.0,
+            sa_time.as_secs_f64() * 1_000_000.0,
+            speedup,
+            if match_result { "✓" } else { "✗" }
+        );
+
+        results.push(BenchmarkResult {
+            name: format!("SA: {}", name),
+            basic_time: batched_time,
+            mt_time: sa_time,
+            results_match: match_result,
+        });
+    }
+    println!();
+}
+
+pub fn test_composable_runtime(results: &mut Vec<BenchmarkResult>) {
+    print_section("Composable Runtime Configurations");
+
+    let builder = || {
+        let mut c = QuantumCircuit::new(6);
+        for q in 0..6 {
+            c.h(q).t(q).s(q);
+        }
+        for q in 0..5 {
+            c.cnot(q, q + 1);
+        }
+        for q in 0..6 {
+            c.rx(q, PI / 4.0).rz(q, PI / 4.0);
+        }
+        c
+    };
+
+    let configs: Vec<(&str, RuntimeConfig)> = vec![
+        ("Basic", RuntimeConfig::new()),
+        ("Batched", RuntimeConfig::new().batched()),
+        ("SIMD", RuntimeConfig::new().simd()),
+        ("Batched+SIMD", RuntimeConfig::new().batched().simd()),
+        ("SA+SIMD", RuntimeConfig::new().structure_aware().simd()),
+        (
+            "SA+SIMD+Parallel",
+            RuntimeConfig::new().structure_aware().simd().parallel(),
+        ),
+        ("Optimal", Runtime::optimal()),
+    ];
+
+    let mut reference = builder();
+    reference.compute_with(Runtime::BasicRT);
+    let ref_state = reference.state().clone();
+
+    println!("Testing 6-qubit circuit with different runtime configurations:\n");
+
+    for (name, config) in &configs {
+        let mut circuit = builder();
+        let start = Instant::now();
+        circuit.compute_with_config(*config);
+        let time = start.elapsed();
+
+        let match_result = states_equal(&ref_state, circuit.state());
+
+        println!(
+            "{:20} : {:.2}μs, Match={}",
+            name,
+            time.as_secs_f64() * 1_000_000.0,
+            if match_result { "✓" } else { "✗" }
+        );
+
+        results.push(BenchmarkResult {
+            name: format!("Config: {}", name),
+            basic_time: time,
+            mt_time: time,
+            results_match: match_result,
+        });
+    }
+
+    println!("\nConfiguration Display Examples:");
+    println!("  {}", RuntimeConfig::new());
+    println!("  {}", RuntimeConfig::new().batched().simd());
+    println!(
+        "  {}",
+        RuntimeConfig::new().structure_aware().simd().parallel()
+    );
+    println!("  {}", Runtime::optimal());
+    println!();
+}
diff --git a/tester/main.c b/tester/main.c
new file mode 100644
index 0000000..eb629fc
--- /dev/null
+++ b/tester/main.c
@@ -0,0 +1,23 @@
+#include <stdio.h>
+
+#include "maths/complex.h"
+#include "psi.h"
+
+int main(void)
+{
+	printf("psi %s\n", psi_version());
+
+	struct PsiComplex a = psi_new_complex(1.0, 2.0);
+	struct PsiComplex b = psi_new_complex(3.0, -1.0);
+
+	struct PsiComplex sum = psi_add_complex(a, b);
+	struct PsiComplex product = psi_mul_complex(a, b);
+
+	printf("a       = %g + %gi\n", a.real, a.imaginary);
+	printf("b       = %g + %gi\n", b.real, b.imaginary);
+	printf("a + b   = %g + %gi\n", sum.real, sum.imaginary);
+	printf("a * b   = %g + %gi\n", product.real, product.imaginary);
+	printf("|a|     = %g\n", psi_abs_complex(a));
+
+	return 0;
+}
diff --git a/tester/main.rs b/tester/main.rs
new file mode 100644
index 0000000..ecebd34
--- /dev/null
+++ b/tester/main.rs
@@ -0,0 +1,101 @@
+mod benchmarks;
+mod clifford;
+mod common;
+mod custom_gates;
+mod kernels;
+mod noise;
+mod non_clifford;
+mod simd;
+
+use common::{print_benchmark_table, print_summary, BenchmarkResult};
+use std::env;
+
+fn print_header() {
+    println!("═══════════════════════════════════════════════════════════════");
+    println!("                    PSI Quantum Simulator");
+    println!("═══════════════════════════════════════════════════════════════\n");
+}
+
+fn print_usage() {
+    println!("Usage: tester [OPTIONS]");
+    println!();
+    println!("Options:");
+    println!("  all          Run all tests (default)");
+    println!("  clifford     Run Clifford gate tests only");
+    println!("  non-clifford Run non-Clifford gate tests only");
+    println!("  custom       Run custom gate tests only");
+    println!("  kernels      Run kernel batching tests only");
+    println!("  simd         Run SIMD acceleration tests only");
+    println!("  noise        Run noise channel tests only");
+    println!("  bench        Run benchmark tests only");
+    println!("  help         Show this help message");
+    println!();
+    println!("Examples:");
+    println!("  tester                   # Run all tests");
+    println!("  tester clifford          # Run only Clifford gate tests");
+    println!("  tester non-clifford      # Run only rotation/parametric gate tests");
+    println!("  tester kernels           # Run only kernel batching tests");
+    println!("  tester simd              # Run only SIMD tests");
+    println!("  tester noise             # Run only noise channel tests");
+    println!("  tester custom bench      # Run custom gates and benchmarks");
+}
+
+fn main() {
+    let args: Vec<String> = env::args().skip(1).collect();
+
+    if args
+        .iter()
+        .any(|a| a == "help" || a == "--help" || a == "-h")
+    {
+        print_usage();
+        return;
+    }
+
+    print_header();
+
+    let mut results: Vec<BenchmarkResult> = Vec::new();
+
+    let run_all = args.is_empty() || args.iter().any(|a| a == "all");
+    let run_clifford = run_all || args.iter().any(|a| a == "clifford");
+    let run_non_clifford = run_all || args.iter().any(|a| a == "non-clifford");
+    let run_custom = run_all || args.iter().any(|a| a == "custom");
+    let run_kernels = run_all || args.iter().any(|a| a == "kernels");
+    let run_simd = run_all || args.iter().any(|a| a == "simd");
+    let run_noise = run_all || args.iter().any(|a| a == "noise");
+    let run_bench = run_all || args.iter().any(|a| a == "bench");
+
+    if run_clifford {
+        clifford::run_all(&mut results);
+    }
+    
+
+    if run_non_clifford {
+        non_clifford::run_all(&mut results);
+    }
+
+    if run_custom {
+        custom_gates::run_all(&mut results);
+    }
+
+    if run_kernels {
+        kernels::run_all(&mut results);
+    }
+
+    if run_simd {
+        simd::run_all(&mut results);
+    }
+
+    if run_noise {
+        noise::run_all(&mut results);
+    }
+
+    if run_bench {
+        benchmarks::run_all(&mut results);
+    }
+    
+
+    if !results.is_empty() {
+        print_benchmark_table(&results);
+        print_summary(&results);
+    }
+}
diff --git a/tester/noise.rs b/tester/noise.rs
new file mode 100644
index 0000000..02a4ba3
--- /dev/null
+++ b/tester/noise.rs
@@ -0,0 +1,172 @@
+use crate::common::{print_section, BenchmarkResult};
+use psi::{
+    complex, DensityMatrix, NoiseChannel, QuantumCircuit, Runtime, Vector,
+};
+use std::time::Instant;
+
+pub fn run_all(results: &mut Vec<BenchmarkResult>) {
+    println!("═══════════════════════════════════════════════════════════════");
+    println!("                    NOISE CHANNEL TESTS");
+    println!("═══════════════════════════════════════════════════════════════\n");
+
+    test_density_matrix_basics(results);
+    test_noise_channels(results);
+    test_noisy_circuit(results);
+}
+
+pub fn test_density_matrix_basics(results: &mut Vec<BenchmarkResult>) {
+    print_section("Density Matrix Basics");
+
+    let dm = DensityMatrix::new(2);
+    println!("Initial |00⟩ state:");
+    println!("{}", dm);
+
+    let mut circuit = QuantumCircuit::new(2);
+    circuit.h(0).cnot(0, 1);
+    circuit.compute_with(Runtime::BasicRT);
+    let state = circuit.state();
+
+    let state_vec: Vec<_> = (0..state.size())
+        .map(|i| state.get(i))
+        .collect();
+
+    let dm_bell = DensityMatrix::from_state_vector(&state_vec);
+    println!("Bell state |Φ+⟩:");
+    println!("{}", dm_bell);
+    println!("Full matrix:");
+    println!("{:?}", dm_bell);
+
+    let is_pure = dm_bell.is_pure(1e-10);
+    println!("Purity check: {}\n", if is_pure { "✓ Pure" } else { "✗ Mixed" });
+
+    results.push(BenchmarkResult {
+        name: "DM: Bell state".to_string(),
+        basic_time: std::time::Duration::from_micros(0),
+        mt_time: std::time::Duration::from_micros(0),
+        results_match: is_pure,
+    });
+}
+
+pub fn test_noise_channels(results: &mut Vec<BenchmarkResult>) {
+    print_section("Noise Channel Effects");
+
+    let channels: Vec<(&str, NoiseChannel)> = vec![
+        ("Depolarising (p=0.1)", NoiseChannel::depolarising(0.1)),
+        ("Amplitude Damping (γ=0.2)", NoiseChannel::amplitude_damping(0.2)),
+        ("Phase Damping (γ=0.2)", NoiseChannel::phase_damping(0.2)),
+        ("Bit Flip (p=0.1)", NoiseChannel::bit_flip(0.1)),
+        ("Phase Flip (p=0.1)", NoiseChannel::phase_flip(0.1)),
+        ("Bit-Phase Flip (p=0.1)", NoiseChannel::bit_phase_flip(0.1)),
+    ];
+
+    let plus_state = vec![
+        complex!(1.0 / 2.0_f64.sqrt(), 0.0),
+        complex!(1.0 / 2.0_f64.sqrt(), 0.0),
+    ];
+
+    println!("Starting with |+⟩ state: (|0⟩ + |1⟩)/√2\n");
+
+    for (name, channel) in channels {
+        let mut dm = DensityMatrix::from_state_vector(&plus_state);
+        let initial_purity = dm.purity();
+
+        let start = Instant::now();
+        dm.apply_noise_channel(&channel, 0);
+        let elapsed = start.elapsed();
+
+        let final_purity = dm.purity();
+        let fidelity = dm.fidelity_with_pure_state(&plus_state);
+
+        println!("{:30}", name);
+        println!("  Purity: {:.4} → {:.4}", initial_purity, final_purity);
+        println!("  Fidelity with |+⟩: {:.4}", fidelity);
+        println!("  Probabilities: {:?}", dm.probabilities());
+        println!("  Time: {:.2}μs\n", elapsed.as_secs_f64() * 1_000_000.0);
+
+        let purity_decreased = final_purity <= initial_purity + 1e-10;
+
+        results.push(BenchmarkResult {
+            name: format!("Noise: {}", name),
+            basic_time: elapsed,
+            mt_time: elapsed,
+            results_match: purity_decreased,
+        });
+    }
+}
+
+pub fn test_noisy_circuit(results: &mut Vec<BenchmarkResult>) {
+    print_section("Noisy Circuit Simulation");
+
+    let mut circuit = QuantumCircuit::new(2);
+    circuit.h(0).cnot(0, 1);
+    circuit.compute_with(Runtime::BasicRT);
+    let state = circuit.state();
+    let state_vec: Vec<_> = (0..state.size()).map(|i| state.get(i)).collect();
+
+    let mut dm = DensityMatrix::from_state_vector(&state_vec);
+    println!("Bell state before noise:");
+    println!("{}", dm);
+
+    let depol = NoiseChannel::depolarising(0.05);
+
+    let start = Instant::now();
+    dm.apply_noise_channel(&depol, 0);
+    dm.apply_noise_channel(&depol, 1);
+    let elapsed = start.elapsed();
+
+    println!("Bell state after 5% depolarising on both qubits:");
+    println!("{}", dm);
+
+    let fidelity = dm.fidelity_with_pure_state(&state_vec);
+    println!("Fidelity with ideal Bell state: {:.4}", fidelity);
+    println!("Time: {:.2}μs\n", elapsed.as_secs_f64() * 1_000_000.0);
+
+    let mut dm2 = DensityMatrix::from_state_vector(&state_vec);
+    let amp_damp = NoiseChannel::amplitude_damping(0.1);
+
+    dm2.apply_noise_channel(&amp_damp, 0);
+    dm2.apply_noise_channel(&amp_damp, 1);
+
+    println!("Bell state after 10% amplitude damping on both qubits:");
+    println!("{}", dm2);
+    println!("Probabilities show decay towards |00⟩: {:?}", dm2.probabilities());
+
+    results.push(BenchmarkResult {
+        name: "Noisy Bell circuit".to_string(),
+        basic_time: elapsed,
+        mt_time: elapsed,
+        results_match: fidelity > 0.8 && fidelity < 1.0,
+    });
+
+    println!();
+    print_section("T1/T2 Relaxation Simulation");
+
+    let one_state = vec![complex!(0.0, 0.0), complex!(1.0, 0.0)];
+    let mut dm_t1 = DensityMatrix::from_state_vector(&one_state);
+
+    println!("Simulating T1 decay of |1⟩ state:");
+    println!("  Initial: P(0)={:.4}, P(1)={:.4}", dm_t1.probabilities()[0], dm_t1.probabilities()[1]);
+
+    let t1_channel = NoiseChannel::amplitude_damping(0.3);
+    for step in 1..=5 {
+        dm_t1.apply_noise_channel(&t1_channel, 0);
+        println!(
+            "  Step {}: P(0)={:.4}, P(1)={:.4}, Purity={:.4}",
+            step,
+            dm_t1.probabilities()[0],
+            dm_t1.probabilities()[1],
+            dm_t1.purity()
+        );
+    }
+
+    let decayed = dm_t1.probabilities()[0] > 0.8;
+    println!("  Decay complete: {}\n", if decayed { "✓" } else { "✗" });
+
+    results.push(BenchmarkResult {
+        name: "T1 decay simulation".to_string(),
+        basic_time: std::time::Duration::from_micros(0),
+        mt_time: std::time::Duration::from_micros(0),
+        results_match: decayed,
+    });
+}
+
diff --git a/tester/non_clifford.rs b/tester/non_clifford.rs
new file mode 100644
index 0000000..2146fb2
--- /dev/null
+++ b/tester/non_clifford.rs
@@ -0,0 +1,137 @@
+use crate::common::{benchmark_circuit, print_circuit, print_section, BenchmarkResult};
+use psi::QuantumCircuit;
+use std::f64::consts::PI;
+
+pub fn run_all(results: &mut Vec<BenchmarkResult>) {
+    println!("═══════════════════════════════════════════════════════════════");
+    println!("                   NON-CLIFFORD GATES TESTS");
+    println!("═══════════════════════════════════════════════════════════════\n");
+
+    test_fixed_gates(results);
+    test_rotation_gates(results);
+    test_phase_gates(results);
+    test_general_unitaries(results);
+    test_controlled_rotations(results);
+    test_variational_circuit(results);
+}
+
+pub fn test_fixed_gates(results: &mut Vec<BenchmarkResult>) {
+    print_section("Non-Clifford Gates: T, T†, √X, S†");
+
+    let builder = || {
+        let mut circuit = QuantumCircuit::new(2);
+        circuit.h(0).t(0).tdg(0).sx(1).sxdg(1).h(0).s(0).sdg(0);
+        circuit
+    };
+
+    print_circuit(&builder());
+    results.push(benchmark_circuit("Non-Clifford fixed gates", builder));
+
+    let mut display = builder();
+    display.compute();
+    println!("{}\n", display);
+}
+
+pub fn test_rotation_gates(results: &mut Vec<BenchmarkResult>) {
+    print_section("Rotation Gates: Rx, Ry, Rz");
+
+    let builder = || {
+        let mut circuit = QuantumCircuit::new(3);
+        circuit
+            .rx(0, PI / 4.0)
+            .ry(1, PI / 2.0)
+            .rz(2, PI)
+            .rx(0, -PI / 4.0);
+        circuit
+    };
+
+    print_circuit(&builder());
+    results.push(benchmark_circuit("Rotation gates (3 qubits)", builder));
+
+    let mut display = builder();
+    display.compute();
+    println!("{}\n", display);
+}
+
+pub fn test_phase_gates(results: &mut Vec<BenchmarkResult>) {
+    print_section("Phase Gate: P(θ)");
+
+    let builder = || {
+        let mut circuit = QuantumCircuit::new(2);
+        circuit.h(0).p(0, PI / 4.0).h(1).p(1, PI / 2.0);
+        circuit
+    };
+
+    print_circuit(&builder());
+    results.push(benchmark_circuit("Phase gates (2 qubits)", builder));
+
+    let mut display = builder();
+    display.compute();
+    println!("{}\n", display);
+}
+
+pub fn test_general_unitaries(results: &mut Vec<BenchmarkResult>) {
+    print_section("General Unitaries: U1, U2, U3");
+
+    let builder = || {
+        let mut circuit = QuantumCircuit::new(3);
+        circuit
+            .u1(0, PI / 4.0)
+            .u2(1, 0.0, PI)
+            .u3(2, PI / 2.0, 0.0, PI);
+        circuit
+    };
+
+    print_circuit(&builder());
+    results.push(benchmark_circuit("General unitaries (3 qubits)", builder));
+
+    let mut display = builder();
+    display.compute();
+    println!("{}\n", display);
+}
+
+pub fn test_controlled_rotations(results: &mut Vec<BenchmarkResult>) {
+    print_section("Controlled Rotation Gates: CRx, CRy, CRz, CP");
+
+    let builder = || {
+        let mut circuit = QuantumCircuit::new(4);
+        circuit
+            .x(0)
+            .crx(0, 1, PI / 2.0)
+            .x(2)
+            .cry(2, 3, PI / 4.0)
+            .crz(0, 2, PI)
+            .cp(1, 3, PI / 2.0);
+        circuit
+    };
+
+    print_circuit(&builder());
+    results.push(benchmark_circuit(
+        "Controlled rotations (4 qubits)",
+        builder,
+    ));
+
+    let mut display = builder();
+    display.compute();
+    println!("{}\n", display);
+}
+
+pub fn test_variational_circuit(results: &mut Vec<BenchmarkResult>) {
+    print_section("Variational Circuit (VQE-like)");
+
+    let builder = || {
+        let mut circuit = QuantumCircuit::new(3);
+        circuit.ry(0, 0.5).ry(1, 0.3).ry(2, 0.7);
+        circuit.cnot(0, 1).cnot(1, 2);
+        circuit.rx(0, 0.2).rx(1, 0.4).rx(2, 0.6);
+        circuit.cz(0, 2);
+        circuit
+    };
+
+    print_circuit(&builder());
+    results.push(benchmark_circuit("Variational circuit (3 qubits)", builder));
+
+    let mut display = builder();
+    display.compute();
+    println!("{}\n", display);
+}
diff --git a/tester/simd.rs b/tester/simd.rs
new file mode 100644
index 0000000..e91dc6e
--- /dev/null
+++ b/tester/simd.rs
@@ -0,0 +1,210 @@
+use crate::common::{print_section, states_equal, BenchmarkResult, CircuitCases};
+use psi::{get_simd_info, QuantumCircuit, Runtime};
+use std::f64::consts::PI;
+use std::time::Instant;
+
+pub fn run_all(results: &mut Vec<BenchmarkResult>) {
+    println!("═══════════════════════════════════════════════════════════════");
+    println!("                    SIMD ACCELERATION TESTS");
+    println!("═══════════════════════════════════════════════════════════════\n");
+
+    println!("Detected: {}\n", get_simd_info());
+
+    test_simd_correctness(results);
+    test_simd_vs_batched(results);
+    test_simd_large_circuits(results);
+}
+
+pub fn test_simd_correctness(results: &mut Vec<BenchmarkResult>) {
+    print_section("SIMD Correctness Verification");
+
+    let test_cases: CircuitCases = vec![
+        (
+            "Bell State",
+            Box::new(|| {
+                let mut c = QuantumCircuit::new(2);
+                c.h(0).cnot(0, 1);
+                c
+            }),
+        ),
+        (
+            "GHZ-3",
+            Box::new(|| {
+                let mut c = QuantumCircuit::new(3);
+                c.h(0).cnot(0, 1).cnot(0, 2);
+                c
+            }),
+        ),
+        (
+            "Rotation Chain",
+            Box::new(|| {
+                let mut c = QuantumCircuit::new(3);
+                c.rx(0, PI / 4.0)
+                    .ry(0, PI / 4.0)
+                    .rz(0, PI / 4.0)
+                    .rx(1, PI / 3.0)
+                    .ry(1, PI / 3.0);
+                c
+            }),
+        ),
+        (
+            "Mixed Single-Qubit",
+            Box::new(|| {
+                let mut c = QuantumCircuit::new(4);
+                c.h(0).t(0).s(0).x(0).h(1).y(1).z(1).h(2).t(2).h(3).s(3);
+                c
+            }),
+        ),
+    ];
+
+    for (name, builder) in test_cases {
+        let mut basic = builder();
+        basic.compute_with(Runtime::BasicRT);
+
+        let mut simd = builder();
+        simd.compute_with(Runtime::SimdRT);
+
+        let match_result = states_equal(basic.state(), simd.state());
+
+        println!(
+            "{}: {}",
+            name,
+            if match_result {
+                "✓ Match"
+            } else {
+                "✗ MISMATCH"
+            }
+        );
+
+        results.push(BenchmarkResult {
+            name: format!("SIMD verify: {}", name),
+            basic_time: std::time::Duration::from_micros(0),
+            mt_time: std::time::Duration::from_micros(0),
+            results_match: match_result,
+        });
+    }
+    println!();
+}
+
+pub fn test_simd_vs_batched(results: &mut Vec<BenchmarkResult>) {
+    print_section("SIMD vs Batched Runtime Comparison");
+
+    let test_cases: CircuitCases = vec![
+        (
+            "Single-Qubit Heavy (6q)",
+            Box::new(|| {
+                let mut c = QuantumCircuit::new(6);
+                for q in 0..6 {
+                    c.h(q).t(q).s(q).x(q).y(q).z(q);
+                }
+                c
+            }),
+        ),
+        (
+            "Rotation Circuit (5q)",
+            Box::new(|| {
+                let mut c = QuantumCircuit::new(5);
+                for q in 0..5 {
+                    c.rx(q, PI / 4.0).ry(q, PI / 3.0).rz(q, PI / 6.0);
+                }
+                c
+            }),
+        ),
+        (
+            "Deep Single-Qubit (4q)",
+            Box::new(|| {
+                let mut c = QuantumCircuit::new(4);
+                for _ in 0..10 {
+                    for q in 0..4 {
+                        c.h(q).t(q);
+                    }
+                }
+                c
+            }),
+        ),
+    ];
+
+    for (name, builder) in test_cases {
+        let mut batched = builder();
+        let start = Instant::now();
+        batched.compute_with(Runtime::BatchedRT);
+        let batched_time = start.elapsed();
+
+        let mut simd = builder();
+        let start = Instant::now();
+        simd.compute_with(Runtime::SimdRT);
+        let simd_time = start.elapsed();
+
+        let match_result = states_equal(batched.state(), simd.state());
+
+        let speedup = batched_time.as_secs_f64() / simd_time.as_secs_f64();
+        println!(
+            "{}: Batched={:.2}μs, SIMD={:.2}μs, Speedup={:.2}x, Match={}",
+            name,
+            batched_time.as_secs_f64() * 1_000_000.0,
+            simd_time.as_secs_f64() * 1_000_000.0,
+            speedup,
+            if match_result { "✓" } else { "✗" }
+        );
+
+        results.push(BenchmarkResult {
+            name: format!("SIMD: {}", name),
+            basic_time: batched_time,
+            mt_time: simd_time,
+            results_match: match_result,
+        });
+    }
+    println!();
+}
+
+pub fn test_simd_large_circuits(results: &mut Vec<BenchmarkResult>) {
+    print_section("SIMD on Large Circuits (Multi-threaded)");
+
+    let sizes = [8, 10, 12];
+
+    for &n in &sizes {
+        let builder = || {
+            let mut circuit = QuantumCircuit::new(n);
+            for i in 0..n {
+                circuit.h(i);
+            }
+            for i in 0..(n - 1) {
+                circuit.cnot(i, i + 1);
+            }
+            for i in 0..n {
+                circuit.t(i).s(i);
+            }
+            circuit
+        };
+
+        let mut batched_mt = builder();
+        let start = Instant::now();
+        batched_mt.compute_with(Runtime::BatchedRTMT);
+        let batched_time = start.elapsed();
+
+        let mut simd_mt = builder();
+        let start = Instant::now();
+        simd_mt.compute_with(Runtime::SimdRTMT);
+        let simd_time = start.elapsed();
+
+        let match_result = states_equal(batched_mt.state(), simd_mt.state());
+
+        let speedup = batched_time.as_secs_f64() / simd_time.as_secs_f64();
+        println!(
+            "{}-qubit: BatchedMT={:.3}ms, SIMD_MT={:.3}ms, Speedup={:.2}x, Match={}",
+            n,
+            batched_time.as_secs_f64() * 1000.0,
+            simd_time.as_secs_f64() * 1000.0,
+            speedup,
+            if match_result { "✓" } else { "✗" }
+        );
+
+        results.push(BenchmarkResult {
+            name: format!("{}-qubit SIMD", n),
+            basic_time: batched_time,
+            mt_time: simd_time,
+            results_match: match_result,
+        });
+    }
+    println!();
+}