[chore]: unwrap project

22 files changed, 5620 insertions, 0 deletions

diff --git a/src/core/circuit.rs b/src/core/circuit.rs
new file mode 100644
index 0000000..6647f58
--- /dev/null
+++ b/src/core/circuit.rs
@@ -0,0 +1,477 @@
+use super::{CustomGate, QuantumState, Runtime, RuntimeConfig};
+use crate::{format_amplitude, format_probability, Vector};
+use core::fmt;
+use std::sync::Arc;
+
+#[derive(Clone)]
+pub enum GateOp {
+    H(usize),
+    X(usize),
+    Y(usize),
+    Z(usize),
+    S(usize),
+    T(usize),
+    Sdg(usize),
+    Tdg(usize),
+    Sx(usize),
+    Sxdg(usize),
+    Rx(usize, f64),
+    Ry(usize, f64),
+    Rz(usize, f64),
+    P(usize, f64),
+    U1(usize, f64),
+    U2(usize, f64, f64),
+    U3(usize, f64, f64, f64),
+    CNOT(usize, usize),
+    CZ(usize, usize),
+    SWAP(usize, usize),
+    CRx(usize, usize, f64),
+    CRy(usize, usize, f64),
+    CRz(usize, usize, f64),
+    CP(usize, usize, f64),
+    CCNOT(usize, usize, usize),
+    CSWAP(usize, usize, usize),
+    Measure(usize, usize),
+    Custom(Arc<CustomGate>, Vec<usize>),
+}
+
+impl GateOp {
+    pub fn name(&self) -> &str {
+        match self {
+            GateOp::H(_) => "H",
+            GateOp::X(_) => "X",
+            GateOp::Y(_) => "Y",
+            GateOp::Z(_) => "Z",
+            GateOp::S(_) => "S",
+            GateOp::T(_) => "T",
+            GateOp::Sdg(_) => "S†",
+            GateOp::Tdg(_) => "T†",
+            GateOp::Sx(_) => "√X",
+            GateOp::Sxdg(_) => "√X†",
+            GateOp::Rx(_, _) => "Rx",
+            GateOp::Ry(_, _) => "Ry",
+            GateOp::Rz(_, _) => "Rz",
+            GateOp::P(_, _) => "P",
+            GateOp::U1(_, _) => "U1",
+            GateOp::U2(_, _, _) => "U2",
+            GateOp::U3(_, _, _, _) => "U3",
+            GateOp::CRx(_, _, _) => "CRx",
+            GateOp::CRy(_, _, _) => "CRy",
+            GateOp::CRz(_, _, _) => "CRz",
+            GateOp::CP(_, _, _) => "CP",
+            GateOp::CNOT(_, _) => "CNOT",
+            GateOp::CZ(_, _) => "CZ",
+            GateOp::SWAP(_, _) => "SWAP",
+            GateOp::CCNOT(_, _, _) => "CCNOT",
+            GateOp::CSWAP(_, _, _) => "CSWAP",
+            GateOp::Measure(_, _) => "M",
+            GateOp::Custom(gate, _) => &gate.name,
+        }
+    }
+
+    pub fn quantum_targets(&self) -> Vec<usize> {
+        match self {
+            GateOp::H(t)
+            | GateOp::X(t)
+            | GateOp::Y(t)
+            | GateOp::Z(t)
+            | GateOp::S(t)
+            | GateOp::T(t)
+            | GateOp::Sdg(t)
+            | GateOp::Tdg(t)
+            | GateOp::Sx(t)
+            | GateOp::Sxdg(t)
+            | GateOp::Rx(t, _)
+            | GateOp::Ry(t, _)
+            | GateOp::Rz(t, _)
+            | GateOp::P(t, _)
+            | GateOp::U1(t, _)
+            | GateOp::U2(t, _, _)
+            | GateOp::U3(t, _, _, _) => vec![*t],
+            GateOp::CNOT(c, t)
+            | GateOp::CZ(c, t)
+            | GateOp::SWAP(c, t)
+            | GateOp::CRx(c, t, _)
+            | GateOp::CRy(c, t, _)
+            | GateOp::CRz(c, t, _)
+            | GateOp::CP(c, t, _) => vec![*c, *t],
+            GateOp::CCNOT(c1, c2, t) | GateOp::CSWAP(c1, c2, t) => vec![*c1, *c2, *t],
+            GateOp::Measure(q, _) => vec![*q],
+            GateOp::Custom(_, targets) => targets.clone(),
+        }
+    }
+
+    pub fn classical_targets(&self) -> Vec<usize> {
+        match self {
+            GateOp::Measure(_, c) => vec![*c],
+            _ => vec![],
+        }
+    }
+
+    pub fn is_measurement(&self) -> bool {
+        matches!(self, GateOp::Measure(_, _))
+    }
+
+    pub fn is_custom(&self) -> bool {
+        matches!(self, GateOp::Custom(_, _))
+    }
+
+    pub fn is_non_clifford(&self) -> bool {
+        matches!(
+            self,
+            GateOp::T(_)
+                | GateOp::Tdg(_)
+                | GateOp::Sx(_)
+                | GateOp::Sxdg(_)
+                | GateOp::Rx(_, _)
+                | GateOp::Ry(_, _)
+                | GateOp::Rz(_, _)
+                | GateOp::P(_, _)
+                | GateOp::U1(_, _)
+                | GateOp::U2(_, _, _)
+                | GateOp::U3(_, _, _, _)
+                | GateOp::CRx(_, _, _)
+                | GateOp::CRy(_, _, _)
+                | GateOp::CRz(_, _, _)
+                | GateOp::CP(_, _, _)
+        )
+    }
+}
+
+pub struct QuantumCircuit {
+    num_qubits: usize,
+    num_classical: usize,
+    operations: Vec<GateOp>,
+    computed_state: Option<QuantumState>,
+}
+
+impl QuantumCircuit {
+    pub fn new(num_qubits: usize) -> QuantumCircuit {
+        QuantumCircuit {
+            num_qubits,
+            num_classical: 0,
+            operations: Vec::new(),
+            computed_state: None,
+        }
+    }
+
+    pub fn with_classical(num_qubits: usize, num_classical: usize) -> QuantumCircuit {
+        QuantumCircuit {
+            num_qubits,
+            num_classical,
+            operations: Vec::new(),
+            computed_state: None,
+        }
+    }
+
+    pub fn num_qubits(&self) -> usize {
+        self.num_qubits
+    }
+
+    pub fn num_classical(&self) -> usize {
+        self.num_classical
+    }
+
+    pub fn operations(&self) -> &[GateOp] {
+        &self.operations
+    }
+
+    pub fn is_computed(&self) -> bool {
+        self.computed_state.is_some()
+    }
+
+    pub fn compute(&mut self) -> &QuantumState {
+        self.compute_with(Runtime::default())
+    }
+
+    pub fn compute_with(&mut self, runtime: Runtime) -> &QuantumState {
+        if self.computed_state.is_none() {
+            self.computed_state = Some(runtime.compute(self.num_qubits, &self.operations));
+        }
+        self.computed_state.as_ref().unwrap()
+    }
+
+    pub fn compute_with_config(&mut self, config: RuntimeConfig) -> &QuantumState {
+        if self.computed_state.is_none() {
+            self.computed_state = Some(config.compute(self.num_qubits, &self.operations));
+        }
+        self.computed_state.as_ref().unwrap()
+    }
+
+    pub fn state(&mut self) -> &QuantumState {
+        self.compute()
+    }
+
+    pub fn state_with(&mut self, runtime: Runtime) -> &QuantumState {
+        self.compute_with(runtime)
+    }
+
+    pub fn state_with_config(&mut self, config: RuntimeConfig) -> &QuantumState {
+        self.compute_with_config(config)
+    }
+
+    pub fn h(&mut self, target: usize) -> &mut Self {
+        self.operations.push(GateOp::H(target));
+        self.computed_state = None;
+        self
+    }
+
+    pub fn x(&mut self, target: usize) -> &mut Self {
+        self.operations.push(GateOp::X(target));
+        self.computed_state = None;
+        self
+    }
+
+    pub fn y(&mut self, target: usize) -> &mut Self {
+        self.operations.push(GateOp::Y(target));
+        self.computed_state = None;
+        self
+    }
+
+    pub fn z(&mut self, target: usize) -> &mut Self {
+        self.operations.push(GateOp::Z(target));
+        self.computed_state = None;
+        self
+    }
+
+    pub fn s(&mut self, target: usize) -> &mut Self {
+        self.operations.push(GateOp::S(target));
+        self.computed_state = None;
+        self
+    }
+
+    pub fn t(&mut self, target: usize) -> &mut Self {
+        self.operations.push(GateOp::T(target));
+        self.computed_state = None;
+        self
+    }
+
+    pub fn sdg(&mut self, target: usize) -> &mut Self {
+        self.operations.push(GateOp::Sdg(target));
+        self.computed_state = None;
+        self
+    }
+
+    pub fn tdg(&mut self, target: usize) -> &mut Self {
+        self.operations.push(GateOp::Tdg(target));
+        self.computed_state = None;
+        self
+    }
+
+    pub fn sx(&mut self, target: usize) -> &mut Self {
+        self.operations.push(GateOp::Sx(target));
+        self.computed_state = None;
+        self
+    }
+
+    pub fn sxdg(&mut self, target: usize) -> &mut Self {
+        self.operations.push(GateOp::Sxdg(target));
+        self.computed_state = None;
+        self
+    }
+
+    pub fn rx(&mut self, target: usize, theta: f64) -> &mut Self {
+        self.operations.push(GateOp::Rx(target, theta));
+        self.computed_state = None;
+        self
+    }
+
+    pub fn ry(&mut self, target: usize, theta: f64) -> &mut Self {
+        self.operations.push(GateOp::Ry(target, theta));
+        self.computed_state = None;
+        self
+    }
+
+    pub fn rz(&mut self, target: usize, theta: f64) -> &mut Self {
+        self.operations.push(GateOp::Rz(target, theta));
+        self.computed_state = None;
+        self
+    }
+
+    pub fn p(&mut self, target: usize, theta: f64) -> &mut Self {
+        self.operations.push(GateOp::P(target, theta));
+        self.computed_state = None;
+        self
+    }
+
+    pub fn u1(&mut self, target: usize, lambda: f64) -> &mut Self {
+        self.operations.push(GateOp::U1(target, lambda));
+        self.computed_state = None;
+        self
+    }
+
+    pub fn u2(&mut self, target: usize, phi: f64, lambda: f64) -> &mut Self {
+        self.operations.push(GateOp::U2(target, phi, lambda));
+        self.computed_state = None;
+        self
+    }
+
+    pub fn u3(&mut self, target: usize, theta: f64, phi: f64, lambda: f64) -> &mut Self {
+        self.operations.push(GateOp::U3(target, theta, phi, lambda));
+        self.computed_state = None;
+        self
+    }
+
+    pub fn crx(&mut self, control: usize, target: usize, theta: f64) -> &mut Self {
+        self.operations.push(GateOp::CRx(control, target, theta));
+        self.computed_state = None;
+        self
+    }
+
+    pub fn cry(&mut self, control: usize, target: usize, theta: f64) -> &mut Self {
+        self.operations.push(GateOp::CRy(control, target, theta));
+        self.computed_state = None;
+        self
+    }
+
+    pub fn crz(&mut self, control: usize, target: usize, theta: f64) -> &mut Self {
+        self.operations.push(GateOp::CRz(control, target, theta));
+        self.computed_state = None;
+        self
+    }
+
+    pub fn cp(&mut self, control: usize, target: usize, theta: f64) -> &mut Self {
+        self.operations.push(GateOp::CP(control, target, theta));
+        self.computed_state = None;
+        self
+    }
+
+    pub fn cnot(&mut self, control: usize, target: usize) -> &mut Self {
+        self.operations.push(GateOp::CNOT(control, target));
+        self.computed_state = None;
+        self
+    }
+
+    pub fn cx(&mut self, control: usize, target: usize) -> &mut Self {
+        self.cnot(control, target)
+    }
+
+    pub fn cz(&mut self, control: usize, target: usize) -> &mut Self {
+        self.operations.push(GateOp::CZ(control, target));
+        self.computed_state = None;
+        self
+    }
+
+    pub fn swap(&mut self, qubit1: usize, qubit2: usize) -> &mut Self {
+        self.operations.push(GateOp::SWAP(qubit1, qubit2));
+        self.computed_state = None;
+        self
+    }
+
+    pub fn ccnot(&mut self, control1: usize, control2: usize, target: usize) -> &mut Self {
+        self.operations
+            .push(GateOp::CCNOT(control1, control2, target));
+        self.computed_state = None;
+        self
+    }
+
+    pub fn toffoli(&mut self, control1: usize, control2: usize, target: usize) -> &mut Self {
+        self.ccnot(control1, control2, target)
+    }
+
+    pub fn cswap(&mut self, control: usize, target1: usize, target2: usize) -> &mut Self {
+        self.operations
+            .push(GateOp::CSWAP(control, target1, target2));
+        self.computed_state = None;
+        self
+    }
+
+    pub fn fredkin(&mut self, control: usize, target1: usize, target2: usize) -> &mut Self {
+        self.cswap(control, target1, target2)
+    }
+
+    pub fn measure(&mut self, qubit: usize, classical: usize) -> &mut Self {
+        if classical >= self.num_classical {
+            self.num_classical = classical + 1;
+        }
+        self.operations.push(GateOp::Measure(qubit, classical));
+        self
+    }
+
+    pub fn measure_all(&mut self) -> &mut Self {
+        for i in 0..self.num_qubits {
+            self.measure(i, i);
+        }
+        self
+    }
+
+    pub fn custom(&mut self, gate: &Arc<CustomGate>, targets: &[usize]) -> &mut Self {
+        self.operations
+            .push(GateOp::Custom(Arc::clone(gate), targets.to_vec()));
+        self.computed_state = None;
+        self
+    }
+
+    pub fn apply_custom(&mut self, gate: CustomGate, targets: &[usize]) -> &mut Self {
+        self.operations
+            .push(GateOp::Custom(Arc::new(gate), targets.to_vec()));
+        self.computed_state = None;
+        self
+    }
+
+    pub fn reset(&mut self) -> &mut Self {
+        self.operations.clear();
+        self.computed_state = None;
+        self
+    }
+
+    pub fn probability(&mut self, state_index: usize) -> f64 {
+        self.compute();
+        let state = self.computed_state.as_ref().unwrap();
+        let amp = state.get(state_index);
+        amp.norm2()
+    }
+
+    pub fn probabilities(&mut self) -> Vec<f64> {
+        self.compute();
+        let n = 1 << self.num_qubits;
+        let state = self.computed_state.as_ref().unwrap();
+        (0..n).map(|i| state.get(i).norm2()).collect()
+    }
+
+    pub fn print_probabilities(&mut self) {
+        let probs = self.probabilities();
+        let n = self.num_qubits;
+        println!("Probabilities:");
+        for (i, p) in probs.iter().enumerate() {
+            if *p > 1e-10 {
+                let basis: String = format!("{:0width$b}", i, width = n);
+                println!("  |{}⟩: {}", basis, format_probability(*p));
+            }
+        }
+    }
+}
+
+impl fmt::Display for QuantumCircuit {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        writeln!(
+            f,
+            "QuantumCircuit ({} qubits, {} classical)",
+            self.num_qubits, self.num_classical
+        )?;
+        writeln!(f, "Operations:")?;
+        for (i, op) in self.operations.iter().enumerate() {
+            match op {
+                GateOp::Measure(q, c) => writeln!(f, "  {}: {} q{} → c{}", i, op.name(), q, c)?,
+                GateOp::Custom(gate, targets) => {
+                    writeln!(f, "  {}: [{}] on {:?}", i, gate.name, targets)?
+                }
+                _ => writeln!(f, "  {}: {} on {:?}", i, op.name(), op.quantum_targets())?,
+            }
+        }
+        if let Some(state) = &self.computed_state {
+            writeln!(f, "State:")?;
+            let n = 1 << self.num_qubits;
+            for i in 0..n {
+                let amp = state.get(i);
+                if amp.real.abs() > 1e-10 || amp.imaginary.abs() > 1e-10 {
+                    let basis: String = format!("{:0width$b}", i, width = self.num_qubits);
+                    writeln!(f, "  |{}⟩: {}", basis, format_amplitude(&amp))?;
+                }
+            }
+        } else {
+            writeln!(f, "State: (not computed)")?;
+        }
+        Ok(())
+    }
+}
diff --git a/src/core/classical_components.rs b/src/core/classical_components.rs
new file mode 100644
index 0000000..b2018d4
--- /dev/null
+++ b/src/core/classical_components.rs
@@ -0,0 +1,63 @@
+use core::ops;
+
+#[derive(Clone, Copy)]
+pub struct ClassicalBit<'a> {
+    state: bool,
+    name: &'a str,
+}
+
+#[derive(Clone)]
+pub struct ClassicalRegister<'a> {
+    bits: Vec<ClassicalBit<'a>>,
+    name: &'a str,
+}
+
+impl<'a> ClassicalBit<'a> {
+    pub fn new(name: &'a str, state: bool) -> ClassicalBit<'a> {
+        ClassicalBit { name, state }
+    }
+
+    pub fn get_name(&self) -> &'a str {
+        self.name
+    }
+
+    pub fn get_state(&self) -> bool {
+        self.state
+    }
+}
+
+impl<'a> ClassicalRegister<'a> {
+    pub fn new(name: &'a str, names: &'a [&'a str]) -> ClassicalRegister<'a> {
+        let mut bits: Vec<ClassicalBit<'a>> = Vec::new();
+        for &name in names {
+            bits.push(ClassicalBit::new(name, false));
+        }
+        ClassicalRegister { name, bits }
+    }
+
+    pub fn set_bits(&mut self, bits: Vec<ClassicalBit<'a>>) {
+        self.bits = bits;
+    }
+
+    pub fn get_bits(&self) -> Vec<ClassicalBit<'a>> {
+        self.bits.clone()
+    }
+
+    pub fn get_name(&self) -> &'a str {
+        self.name
+    }
+}
+
+impl<'a> ops::Index<usize> for ClassicalRegister<'a> {
+    type Output = ClassicalBit<'a>;
+
+    fn index(&self, index: usize) -> &Self::Output {
+        &self.bits[index]
+    }
+}
+
+impl<'a> ops::IndexMut<usize> for ClassicalRegister<'a> {
+    fn index_mut(&mut self, index: usize) -> &mut Self::Output {
+        &mut self.bits[index]
+    }
+}
diff --git a/src/core/custom_gate.rs b/src/core/custom_gate.rs
new file mode 100644
index 0000000..51d3a56
--- /dev/null
+++ b/src/core/custom_gate.rs
@@ -0,0 +1,245 @@
+use crate::{Complex, Matrix, QuantumGate};
+
+#[derive(Clone)]
+pub enum CustomGateDefinition {
+    Matrix(Matrix<Complex<f64>>),
+    Composite(Vec<(CompositeOp, Vec<usize>)>),
+}
+
+#[derive(Clone, Copy)]
+pub enum CompositeOp {
+    H,
+    X,
+    Y,
+    Z,
+    S,
+    T,
+    CNOT,
+    CZ,
+    SWAP,
+    CCNOT,
+    CSWAP,
+}
+
+#[derive(Clone)]
+pub struct CustomGate {
+    pub name: String,
+    pub num_qubits: usize,
+    pub definition: CustomGateDefinition,
+}
+
+impl CustomGate {
+    pub fn from_matrix(name: &str, matrix: Matrix<Complex<f64>>) -> Self {
+        let dim = matrix.rows;
+        let num_qubits = (dim as f64).log2() as usize;
+        assert_eq!(
+            1 << num_qubits,
+            dim,
+            "Matrix dimension must be a power of 2"
+        );
+        assert_eq!(matrix.rows, matrix.cols, "Matrix must be square");
+
+        CustomGate {
+            name: String::from(name),
+            num_qubits,
+            definition: CustomGateDefinition::Matrix(matrix),
+        }
+    }
+
+    pub fn from_composite(
+        name: &str,
+        num_qubits: usize,
+        ops: Vec<(CompositeOp, Vec<usize>)>,
+    ) -> Self {
+        CustomGate {
+            name: String::from(name),
+            num_qubits,
+            definition: CustomGateDefinition::Composite(ops),
+        }
+    }
+
+    pub fn to_quantum_gate(&self) -> QuantumGate<'static> {
+        match &self.definition {
+            CustomGateDefinition::Matrix(matrix) => {
+                let name: &'static str = Box::leak(self.name.clone().into_boxed_str());
+                QuantumGate {
+                    name,
+                    matrix: matrix.clone(),
+                    num_qubits: self.num_qubits,
+                }
+            }
+            CustomGateDefinition::Composite(ops) => {
+                let matrix = self.compute_composite_matrix(ops);
+                let name: &'static str = Box::leak(self.name.clone().into_boxed_str());
+                QuantumGate {
+                    name,
+                    matrix,
+                    num_qubits: self.num_qubits,
+                }
+            }
+        }
+    }
+
+    fn compute_composite_matrix(&self, ops: &[(CompositeOp, Vec<usize>)]) -> Matrix<Complex<f64>> {
+        use crate::gates::*;
+        use crate::Complex;
+
+        let dim = 1 << self.num_qubits;
+        let mut result = Matrix::new(dim, dim, vec![Complex::new(0.0, 0.0); dim * dim]);
+        for i in 0..dim {
+            result.data[i * dim + i] = Complex::new(1.0, 0.0);
+        }
+
+        for (op, targets) in ops {
+            let gate: &QuantumGate = match op {
+                CompositeOp::H => &HADAMARD,
+                CompositeOp::X => &PAULI_X,
+                CompositeOp::Y => &PAULI_Y,
+                CompositeOp::Z => &PAULI_Z,
+                CompositeOp::S => &S_GATE,
+                CompositeOp::T => &T_GATE,
+                CompositeOp::CNOT => &CNOT,
+                CompositeOp::CZ => &CZ,
+                CompositeOp::SWAP => &SWAP,
+                CompositeOp::CCNOT => &TOFFOLI,
+                CompositeOp::CSWAP => &FREDKIN,
+            };
+
+            let full_gate = build_full_operator(&gate.matrix, targets, self.num_qubits);
+            result = matrix_multiply(&full_gate, &result);
+        }
+
+        result
+    }
+}
+
+fn build_full_operator(
+    gate_matrix: &Matrix<Complex<f64>>,
+    targets: &[usize],
+    total_qubits: usize,
+) -> Matrix<Complex<f64>> {
+    let dim = 1 << total_qubits;
+    let gate_dim = gate_matrix.rows;
+    let num_gate_qubits = targets.len();
+
+    let mut result = Matrix::new(dim, dim, vec![Complex::new(0.0, 0.0); dim * dim]);
+
+    for i in 0..dim {
+        for j in 0..dim {
+            let mut gate_i = 0usize;
+            let mut gate_j = 0usize;
+            let mut match_non_targets = true;
+
+            for q in 0..total_qubits {
+                let bit_i = (i >> (total_qubits - 1 - q)) & 1;
+                let bit_j = (j >> (total_qubits - 1 - q)) & 1;
+
+                if let Some(pos) = targets.iter().position(|&t| t == q) {
+                    gate_i |= bit_i << (num_gate_qubits - 1 - pos);
+                    gate_j |= bit_j << (num_gate_qubits - 1 - pos);
+                } else if bit_i != bit_j {
+                    match_non_targets = false;
+                    break;
+                }
+            }
+
+            if match_non_targets {
+                result.data[i * dim + j] = gate_matrix.data[gate_i * gate_dim + gate_j];
+            }
+        }
+    }
+
+    result
+}
+
+fn matrix_multiply(a: &Matrix<Complex<f64>>, b: &Matrix<Complex<f64>>) -> Matrix<Complex<f64>> {
+    let n = a.rows;
+    let mut result = Matrix::new(n, n, vec![Complex::new(0.0, 0.0); n * n]);
+
+    for i in 0..n {
+        for j in 0..n {
+            let mut sum = Complex::new(0.0, 0.0);
+            for k in 0..n {
+                sum += a.data[i * n + k] * b.data[k * n + j];
+            }
+            result.data[i * n + j] = sum;
+        }
+    }
+
+    result
+}
+
+pub struct CustomGateBuilder {
+    name: String,
+    num_qubits: usize,
+    ops: Vec<(CompositeOp, Vec<usize>)>,
+}
+
+impl CustomGateBuilder {
+    pub fn new(name: &str, num_qubits: usize) -> Self {
+        CustomGateBuilder {
+            name: String::from(name),
+            num_qubits,
+            ops: Vec::new(),
+        }
+    }
+
+    pub fn h(mut self, target: usize) -> Self {
+        self.ops.push((CompositeOp::H, vec![target]));
+        self
+    }
+
+    pub fn x(mut self, target: usize) -> Self {
+        self.ops.push((CompositeOp::X, vec![target]));
+        self
+    }
+
+    pub fn y(mut self, target: usize) -> Self {
+        self.ops.push((CompositeOp::Y, vec![target]));
+        self
+    }
+
+    pub fn z(mut self, target: usize) -> Self {
+        self.ops.push((CompositeOp::Z, vec![target]));
+        self
+    }
+
+    pub fn s(mut self, target: usize) -> Self {
+        self.ops.push((CompositeOp::S, vec![target]));
+        self
+    }
+
+    pub fn t(mut self, target: usize) -> Self {
+        self.ops.push((CompositeOp::T, vec![target]));
+        self
+    }
+
+    pub fn cnot(mut self, control: usize, target: usize) -> Self {
+        self.ops.push((CompositeOp::CNOT, vec![control, target]));
+        self
+    }
+
+    pub fn cz(mut self, control: usize, target: usize) -> Self {
+        self.ops.push((CompositeOp::CZ, vec![control, target]));
+        self
+    }
+
+    pub fn swap(mut self, a: usize, b: usize) -> Self {
+        self.ops.push((CompositeOp::SWAP, vec![a, b]));
+        self
+    }
+
+    pub fn ccnot(mut self, c1: usize, c2: usize, target: usize) -> Self {
+        self.ops.push((CompositeOp::CCNOT, vec![c1, c2, target]));
+        self
+    }
+
+    pub fn cswap(mut self, control: usize, t1: usize, t2: usize) -> Self {
+        self.ops.push((CompositeOp::CSWAP, vec![control, t1, t2]));
+        self
+    }
+
+    pub fn build(self) -> CustomGate {
+        CustomGate::from_composite(&self.name, self.num_qubits, self.ops)
+    }
+}
diff --git a/src/core/gates.rs b/src/core/gates.rs
new file mode 100644
index 0000000..cee4ebc
--- /dev/null
+++ b/src/core/gates.rs
@@ -0,0 +1,253 @@
+use crate::{complex, matrix, Complex, Matrix, QuantumGate};
+use std::f64::consts::FRAC_1_SQRT_2;
+
+pub fn rx_matrix(theta: f64) -> Matrix<Complex<f64>> {
+    let cos = (theta / 2.0).cos();
+    let sin = (theta / 2.0).sin();
+    matrix!(
+        [complex!(cos, 0.0), complex!(0.0, -sin)];
+        [complex!(0.0, -sin), complex!(cos, 0.0)]
+    )
+}
+
+pub fn ry_matrix(theta: f64) -> Matrix<Complex<f64>> {
+    let cos = (theta / 2.0).cos();
+    let sin = (theta / 2.0).sin();
+    matrix!(
+        [complex!(cos, 0.0), complex!(-sin, 0.0)];
+        [complex!(sin, 0.0), complex!(cos, 0.0)]
+    )
+}
+
+pub fn rz_matrix(theta: f64) -> Matrix<Complex<f64>> {
+    let half = theta / 2.0;
+    matrix!(
+        [complex!(half.cos(), -half.sin()), complex!(0.0, 0.0)];
+        [complex!(0.0, 0.0), complex!(half.cos(), half.sin())]
+    )
+}
+
+pub fn p_matrix(theta: f64) -> Matrix<Complex<f64>> {
+    matrix!(
+        [complex!(1.0, 0.0), complex!(0.0, 0.0)];
+        [complex!(0.0, 0.0), complex!(theta.cos(), theta.sin())]
+    )
+}
+
+pub fn u1_matrix(lambda: f64) -> Matrix<Complex<f64>> {
+    p_matrix(lambda)
+}
+
+pub fn u2_matrix(phi: f64, lambda: f64) -> Matrix<Complex<f64>> {
+    let inv_sqrt2 = FRAC_1_SQRT_2;
+    matrix!(
+        [complex!(inv_sqrt2, 0.0), complex!(-inv_sqrt2 * lambda.cos(), -inv_sqrt2 * lambda.sin())];
+        [complex!(inv_sqrt2 * phi.cos(), inv_sqrt2 * phi.sin()), complex!((phi + lambda).cos() * inv_sqrt2, (phi + lambda).sin() * inv_sqrt2)]
+    )
+}
+
+pub fn u3_matrix(theta: f64, phi: f64, lambda: f64) -> Matrix<Complex<f64>> {
+    let cos = (theta / 2.0).cos();
+    let sin = (theta / 2.0).sin();
+    matrix!(
+        [complex!(cos, 0.0), complex!(-sin * lambda.cos(), -sin * lambda.sin())];
+        [complex!(sin * phi.cos(), sin * phi.sin()), complex!(cos * (phi + lambda).cos(), cos * (phi + lambda).sin())]
+    )
+}
+
+pub fn crx_matrix(theta: f64) -> Matrix<Complex<f64>> {
+    let cos = (theta / 2.0).cos();
+    let sin = (theta / 2.0).sin();
+    matrix!(
+        [complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
+        [complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
+        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(cos, 0.0), complex!(0.0, -sin)];
+        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, -sin), complex!(cos, 0.0)]
+    )
+}
+
+pub fn cry_matrix(theta: f64) -> Matrix<Complex<f64>> {
+    let cos = (theta / 2.0).cos();
+    let sin = (theta / 2.0).sin();
+    matrix!(
+        [complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
+        [complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
+        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(cos, 0.0), complex!(-sin, 0.0)];
+        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(sin, 0.0), complex!(cos, 0.0)]
+    )
+}
+
+pub fn crz_matrix(theta: f64) -> Matrix<Complex<f64>> {
+    let half = theta / 2.0;
+    matrix!(
+        [complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
+        [complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
+        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(half.cos(), -half.sin()), complex!(0.0, 0.0)];
+        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(half.cos(), half.sin())]
+    )
+}
+
+pub fn cp_matrix(theta: f64) -> Matrix<Complex<f64>> {
+    matrix!(
+        [complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
+        [complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
+        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0)];
+        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(theta.cos(), theta.sin())]
+    )
+}
+
+#[rustfmt::skip]
+lazy_static::lazy_static! {
+    pub static ref HADAMARD: QuantumGate<'static> = QuantumGate {
+        name: "H",
+        matrix: matrix!([complex!(1.0, 0.0), complex!( 1.0, 0.0)];
+                        [complex!(1.0, 0.0), complex!(-1.0, 0.0)]) *
+                complex!(1.0/2.0_f64.sqrt(), 0.0),
+        num_qubits: 1,
+    };
+
+    pub static ref PAULI_X: QuantumGate<'static> = QuantumGate {
+        name: "X",
+        matrix: matrix!([complex!(0.0, 0.0), complex!(1.0, 0.0)];
+                        [complex!(1.0, 0.0), complex!(0.0, 0.0)]),
+        num_qubits: 1,
+    };
+
+    pub static ref PAULI_Y: QuantumGate<'static> = QuantumGate {
+        name: "Y", 
+        matrix: matrix!([complex!(0.0, 0.0), complex!(0.0, -1.0)];
+                        [complex!(0.0, 1.0), complex!(0.0,  0.0)]),
+        num_qubits: 1,
+    };
+
+    pub static ref PAULI_Z: QuantumGate<'static> = QuantumGate {
+        name: "Z", 
+        matrix: matrix!([complex!(1.0, 0.0), complex!( 0.0, 0.0)];
+                        [complex!(0.0, 0.0), complex!(-1.0, 0.0)]),
+        num_qubits: 1,
+    };
+    
+    pub static ref S_GATE: QuantumGate<'static> = QuantumGate {
+        name: "S",
+        matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0)];
+                        [complex!(0.0, 0.0), complex!(0.0, 1.0)]),
+        num_qubits: 1,
+    };
+    
+    pub static ref T_GATE: QuantumGate<'static> = QuantumGate {
+        name: "T",
+        matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0)];
+                        [complex!(0.0, 0.0), complex!(core::f64::consts::FRAC_1_SQRT_2, core::f64::consts::FRAC_1_SQRT_2)]),
+        num_qubits: 1,
+    };
+    
+    pub static ref SDG_GATE: QuantumGate<'static> = QuantumGate {
+        name: "S†",
+        matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0)];
+                        [complex!(0.0, 0.0), complex!(0.0, -1.0)]),
+        num_qubits: 1,
+    };
+    
+    pub static ref TDG_GATE: QuantumGate<'static> = QuantumGate {
+        name: "T†",
+        matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0)];
+                        [complex!(0.0, 0.0), complex!(core::f64::consts::FRAC_1_SQRT_2, -core::f64::consts::FRAC_1_SQRT_2)]),
+        num_qubits: 1,
+    };
+    
+    pub static ref SX_GATE: QuantumGate<'static> = QuantumGate {
+        name: "√X",
+        matrix: matrix!([complex!(0.5, 0.5), complex!(0.5, -0.5)];
+                        [complex!(0.5, -0.5), complex!(0.5, 0.5)]),
+        num_qubits: 1,
+    };
+    
+    pub static ref SXDG_GATE: QuantumGate<'static> = QuantumGate {
+        name: "√X†",
+        matrix: matrix!([complex!(0.5, -0.5), complex!(0.5, 0.5)];
+                        [complex!(0.5, 0.5), complex!(0.5, -0.5)]),
+        num_qubits: 1,
+    };
+    
+    pub static ref IDENTITY: QuantumGate<'static> = QuantumGate {
+        name: "I",
+        matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0)];
+                        [complex!(0.0, 0.0), complex!(1.0, 0.0)]),
+        num_qubits: 1,
+    };
+
+    pub static ref CNOT: QuantumGate<'static> = QuantumGate {
+        name: "CNOT", 
+        matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
+                        [complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
+                        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0)];
+                        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0)]),
+        num_qubits: 2,
+    };
+    
+    pub static ref CZ: QuantumGate<'static> = QuantumGate {
+        name: "CZ",
+        matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!( 0.0, 0.0)];
+                        [complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!( 0.0, 0.0)];
+                        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!( 0.0, 0.0)];
+                        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(-1.0, 0.0)]),
+        num_qubits: 2,
+    };
+    
+    pub static ref SWAP: QuantumGate<'static> = QuantumGate {
+        name: "SWAP",
+        matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
+                        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0)];
+                        [complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
+                        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0)]),
+        num_qubits: 2,
+    };
+    
+    pub static ref ISWAP: QuantumGate<'static> = QuantumGate {
+        name: "iSWAP",
+        matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
+                        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 1.0), complex!(0.0, 0.0)];
+                        [complex!(0.0, 0.0), complex!(0.0, 1.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
+                        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0)]),
+        num_qubits: 2,
+    };
+    
+    pub static ref SQRT_SWAP: QuantumGate<'static> = QuantumGate {
+        name: "√SWAP",
+        matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
+                        [complex!(0.0, 0.0), complex!(0.5, 0.5), complex!(0.5, -0.5), complex!(0.0, 0.0)];
+                        [complex!(0.0, 0.0), complex!(0.5, -0.5), complex!(0.5, 0.5), complex!(0.0, 0.0)];
+                        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0)]),
+        num_qubits: 2,
+    };
+
+    pub static ref TOFFOLI: QuantumGate<'static> = QuantumGate {
+        name: "CCNOT",
+        matrix: matrix!(
+            [complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
+            [complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
+            [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
+            [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
+            [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
+            [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
+            [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0)];
+            [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0)]
+        ),
+        num_qubits: 3,
+    };
+    
+    pub static ref FREDKIN: QuantumGate<'static> = QuantumGate {
+        name: "CSWAP",
+        matrix: matrix!(
+            [complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
+            [complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
+            [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
+            [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
+            [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
+            [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0)];
+            [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
+            [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0)]
+        ),
+        num_qubits: 3,
+    };
+}
diff --git a/src/core/kernel.rs b/src/core/kernel.rs
new file mode 100644
index 0000000..7f73458
--- /dev/null
+++ b/src/core/kernel.rs
@@ -0,0 +1,621 @@
+use crate::maths::simd::{
+    apply_single_qubit_gate_simd, apply_single_qubit_gate_simd_parallel, SimdCapability,
+};
+use crate::{complex, Complex, Matrix};
+use rayon::prelude::*;
+use std::collections::HashSet;
+
+#[derive(Clone, Copy, Debug, PartialEq, Eq)]
+pub enum GateType {
+    Diagonal,
+    NonDiagonal,
+    Controlled,
+}
+
+#[derive(Clone)]
+pub struct Kernel {
+    pub matrix: Matrix<Complex<f64>>,
+    pub targets: Vec<usize>,
+    pub name: String,
+    pub gate_type: GateType,
+}
+
+impl Kernel {
+    pub fn new(name: &str, matrix: Matrix<Complex<f64>>, targets: Vec<usize>) -> Self {
+        let gate_type = Self::detect_gate_type(name, &matrix);
+        Self {
+            matrix,
+            targets,
+            name: name.to_string(),
+            gate_type,
+        }
+    }
+
+    fn detect_gate_type(name: &str, matrix: &Matrix<Complex<f64>>) -> GateType {
+        let diagonal_gates = [
+            "Z", "S", "T", "Sdg", "Tdg", "Rz", "P", "U1", "CZ", "CP", "CRz",
+        ];
+        if diagonal_gates.iter().any(|&g| name.starts_with(g)) {
+            return GateType::Diagonal;
+        }
+
+        let controlled_gates = [
+            "CNOT", "CZ", "SWAP", "CRx", "CRy", "CRz", "CP", "CCNOT", "CSWAP",
+        ];
+        if controlled_gates.iter().any(|&g| name.starts_with(g)) {
+            return GateType::Controlled;
+        }
+
+        if matrix.rows == 2 && matrix.cols == 2 {
+            let is_diag = matrix.data[1].real.abs() < 1e-10
+                && matrix.data[1].imaginary.abs() < 1e-10
+                && matrix.data[2].real.abs() < 1e-10
+                && matrix.data[2].imaginary.abs() < 1e-10;
+            if is_diag {
+                return GateType::Diagonal;
+            }
+        }
+
+        GateType::NonDiagonal
+    }
+
+    pub fn num_qubits(&self) -> usize {
+        self.targets.len()
+    }
+
+    pub fn target_set(&self) -> HashSet<usize> {
+        self.targets.iter().cloned().collect()
+    }
+
+    pub fn shares_qubits(&self, other: &Kernel) -> bool {
+        self.targets.iter().any(|t| other.targets.contains(t))
+    }
+
+    pub fn commutes_with(&self, other: &Kernel) -> bool {
+        if !self.shares_qubits(other) {
+            return true;
+        }
+
+        if self.gate_type == GateType::Diagonal && other.gate_type == GateType::Diagonal
+            && self.targets == other.targets {
+                return true;
+            }
+
+        false
+    }
+
+    pub fn can_fuse_with(&self, other: &Kernel) -> bool {
+        if self.targets.len() != 1 || other.targets.len() != 1 {
+            return false;
+        }
+        self.targets[0] == other.targets[0]
+    }
+
+    pub fn fuse(&self, other: &Kernel) -> Option<Kernel> {
+        if !self.can_fuse_with(other) {
+            return None;
+        }
+        let fused_matrix = other.matrix.dot(&self.matrix)?;
+        let new_type =
+            if self.gate_type == GateType::Diagonal && other.gate_type == GateType::Diagonal {
+                GateType::Diagonal
+            } else {
+                GateType::NonDiagonal
+            };
+        Some(Kernel {
+            matrix: fused_matrix,
+            targets: self.targets.clone(),
+            name: format!("{}+{}", self.name, other.name),
+            gate_type: new_type,
+        })
+    }
+}
+
+pub struct KernelBatch {
+    kernels: Vec<Kernel>,
+    num_qubits: usize,
+}
+
+impl KernelBatch {
+    pub fn new(num_qubits: usize) -> Self {
+        Self {
+            kernels: Vec::new(),
+            num_qubits,
+        }
+    }
+
+    pub fn add(&mut self, kernel: Kernel) {
+        self.kernels.push(kernel);
+    }
+
+    pub fn len(&self) -> usize {
+        self.kernels.len()
+    }
+
+    pub fn is_empty(&self) -> bool {
+        self.kernels.is_empty()
+    }
+
+    pub fn kernels(&self) -> &[Kernel] {
+        &self.kernels
+    }
+
+    pub fn optimize(&mut self) {
+        if self.kernels.len() < 2 {
+            return;
+        }
+
+        let mut optimized: Vec<Kernel> = Vec::with_capacity(self.kernels.len());
+        let mut i = 0;
+
+        while i < self.kernels.len() {
+            let current = &self.kernels[i];
+
+            if i + 1 < self.kernels.len() {
+                let next = &self.kernels[i + 1];
+                if let Some(fused) = current.fuse(next) {
+                    optimized.push(fused);
+                    i += 2;
+                    continue;
+                }
+            }
+
+            optimized.push(current.clone());
+            i += 1;
+        }
+
+        self.kernels = optimized;
+    }
+
+    pub fn execute(&self, state: &mut Vec<Complex<f64>>) {
+        for kernel in &self.kernels {
+            *state = apply_kernel(state, kernel, self.num_qubits);
+        }
+    }
+
+    pub fn execute_parallel(&self, state: &mut Vec<Complex<f64>>) {
+        for kernel in &self.kernels {
+            *state = apply_kernel_parallel(state, kernel, self.num_qubits);
+        }
+    }
+
+    pub fn execute_simd(&self, state: &mut Vec<Complex<f64>>) {
+        for kernel in &self.kernels {
+            if kernel.targets.len() == 1 {
+                let gate = matrix_to_2x2(&kernel.matrix);
+                apply_single_qubit_gate_simd(state, &gate, kernel.targets[0], self.num_qubits);
+            } else {
+                *state = apply_kernel(state, kernel, self.num_qubits);
+            }
+        }
+    }
+
+    pub fn execute_simd_parallel(&self, state: &mut Vec<Complex<f64>>) {
+        for kernel in &self.kernels {
+            if kernel.targets.len() == 1 && self.num_qubits >= 10 {
+                let gate = matrix_to_2x2(&kernel.matrix);
+                apply_single_qubit_gate_simd_parallel(
+                    state,
+                    &gate,
+                    kernel.targets[0],
+                    self.num_qubits,
+                );
+            } else if kernel.targets.len() == 1 {
+                let gate = matrix_to_2x2(&kernel.matrix);
+                apply_single_qubit_gate_simd(state, &gate, kernel.targets[0], self.num_qubits);
+            } else {
+                *state = apply_kernel_parallel(state, kernel, self.num_qubits);
+            }
+        }
+    }
+
+    pub fn simd_capability(&self) -> SimdCapability {
+        SimdCapability::detect()
+    }
+}
+
+fn matrix_to_2x2(matrix: &Matrix<Complex<f64>>) -> [[Complex<f64>; 2]; 2] {
+    [
+        [matrix.data[0], matrix.data[1]],
+        [matrix.data[2], matrix.data[3]],
+    ]
+}
+
+fn apply_kernel(state: &[Complex<f64>], kernel: &Kernel, num_qubits: usize) -> Vec<Complex<f64>> {
+    let dim = 1 << num_qubits;
+    let g = kernel.targets.len();
+    let gate_dim = 1 << g;
+
+    let target_bits: Vec<usize> = kernel.targets.iter().map(|&t| num_qubits - 1 - t).collect();
+
+    let mut non_target_mask: usize = (1 << num_qubits) - 1;
+    for &pos in &target_bits {
+        non_target_mask &= !(1 << pos);
+    }
+
+    let mut new_state = vec![complex!(0.0, 0.0); dim];
+
+    for (i, new_val) in new_state.iter_mut().enumerate() {
+        let mut target_idx = 0usize;
+        for (k, &pos) in target_bits.iter().enumerate() {
+            if (i >> pos) & 1 == 1 {
+                target_idx |= 1 << (g - 1 - k);
+            }
+        }
+
+        let mut sum = complex!(0.0, 0.0);
+
+        for j in 0..gate_dim {
+            let gate_elem = kernel.matrix.data[target_idx * gate_dim + j];
+
+            if gate_elem.real.abs() < 1e-15 && gate_elem.imaginary.abs() < 1e-15 {
+                continue;
+            }
+
+            let mut source_idx = i & non_target_mask;
+            for (k, &pos) in target_bits.iter().enumerate() {
+                if (j >> (g - 1 - k)) & 1 == 1 {
+                    source_idx |= 1 << pos;
+                }
+            }
+
+            sum += gate_elem * state[source_idx];
+        }
+
+        *new_val = sum;
+    }
+
+    new_state
+}
+
+fn apply_kernel_parallel(
+    state: &[Complex<f64>],
+    kernel: &Kernel,
+    num_qubits: usize,
+) -> Vec<Complex<f64>> {
+    let dim = 1 << num_qubits;
+    let g = kernel.targets.len();
+    let gate_dim = 1 << g;
+
+    let target_bits: Vec<usize> = kernel.targets.iter().map(|&t| num_qubits - 1 - t).collect();
+
+    let mut non_target_mask: usize = (1 << num_qubits) - 1;
+    for &pos in &target_bits {
+        non_target_mask &= !(1 << pos);
+    }
+
+    (0..dim)
+        .into_par_iter()
+        .map(|i| {
+            let mut target_idx = 0usize;
+            for (k, &pos) in target_bits.iter().enumerate() {
+                if (i >> pos) & 1 == 1 {
+                    target_idx |= 1 << (g - 1 - k);
+                }
+            }
+
+            let mut sum = complex!(0.0, 0.0);
+
+            for j in 0..gate_dim {
+                let gate_elem = kernel.matrix.data[target_idx * gate_dim + j];
+
+                if gate_elem.real.abs() < 1e-15 && gate_elem.imaginary.abs() < 1e-15 {
+                    continue;
+                }
+
+                let mut source_idx = i & non_target_mask;
+                for (k, &pos) in target_bits.iter().enumerate() {
+                    if (j >> (g - 1 - k)) & 1 == 1 {
+                        source_idx |= 1 << pos;
+                    }
+                }
+
+                sum += gate_elem * state[source_idx];
+            }
+
+            sum
+        })
+        .collect()
+}
+
+pub struct KernelBuilder {
+    num_qubits: usize,
+}
+
+impl KernelBuilder {
+    pub fn new(num_qubits: usize) -> Self {
+        Self { num_qubits }
+    }
+
+    pub fn num_qubits(&self) -> usize {
+        self.num_qubits
+    }
+}
+
+#[derive(Clone)]
+pub struct ExecutionLayer {
+    pub kernels: Vec<Kernel>,
+}
+
+impl ExecutionLayer {
+    pub fn new() -> Self {
+        Self {
+            kernels: Vec::new(),
+        }
+    }
+
+    pub fn can_add(&self, kernel: &Kernel) -> bool {
+        !self.kernels.iter().any(|k| k.shares_qubits(kernel))
+    }
+
+    pub fn add(&mut self, kernel: Kernel) {
+        self.kernels.push(kernel);
+    }
+
+    pub fn affected_qubits(&self) -> HashSet<usize> {
+        self.kernels
+            .iter()
+            .flat_map(|k| k.targets.iter().cloned())
+            .collect()
+    }
+}
+
+impl Default for ExecutionLayer {
+    fn default() -> Self {
+        Self::new()
+    }
+}
+
+pub struct StructureAwareKernelBatch {
+    kernels: Vec<Kernel>,
+    layers: Vec<ExecutionLayer>,
+    num_qubits: usize,
+    optimised: bool,
+}
+
+impl StructureAwareKernelBatch {
+    pub fn new(num_qubits: usize) -> Self {
+        Self {
+            kernels: Vec::new(),
+            layers: Vec::new(),
+            num_qubits,
+            optimised: false,
+        }
+    }
+
+    pub fn add(&mut self, kernel: Kernel) {
+        self.kernels.push(kernel);
+        self.optimised = false;
+    }
+
+    pub fn len(&self) -> usize {
+        self.kernels.len()
+    }
+
+    pub fn is_empty(&self) -> bool {
+        self.kernels.is_empty()
+    }
+
+    pub fn kernels(&self) -> &[Kernel] {
+        &self.kernels
+    }
+
+    pub fn layers(&self) -> &[ExecutionLayer] {
+        &self.layers
+    }
+
+    pub fn num_layers(&self) -> usize {
+        self.layers.len()
+    }
+
+    pub fn optimise(&mut self) {
+        if self.optimised || self.kernels.len() < 2 {
+            return;
+        }
+
+        self.reorder_commuting_gates();
+        self.multi_pass_fusion();
+        self.build_execution_layers();
+        self.optimised = true;
+    }
+
+    fn reorder_commuting_gates(&mut self) {
+        let mut changed = true;
+        let mut iterations = 0;
+        const MAX_ITERATIONS: usize = 100;
+
+        while changed && iterations < MAX_ITERATIONS {
+            changed = false;
+            iterations += 1;
+
+            for i in 0..self.kernels.len().saturating_sub(1) {
+                let current = &self.kernels[i];
+                let next = &self.kernels[i + 1];
+
+                if current.targets.len() == 1
+                    && next.targets.len() == 1
+                    && current.targets[0] != next.targets[0]
+                    && current.commutes_with(next)
+                {
+                    for j in (i + 2)..self.kernels.len() {
+                        let candidate = &self.kernels[j];
+
+                        if candidate.targets.len() == 1
+                            && candidate.targets[0] == current.targets[0]
+                        {
+                            let can_move = (i + 1..j).all(|k| {
+                                let between = &self.kernels[k];
+                                !between.shares_qubits(current) || current.commutes_with(between)
+                            });
+
+                            if can_move && current.can_fuse_with(candidate) {
+                                let kernel_to_move = self.kernels.remove(j);
+                                self.kernels.insert(i + 1, kernel_to_move);
+                                changed = true;
+                                break;
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    fn multi_pass_fusion(&mut self) {
+        let mut changed = true;
+        let mut iterations = 0;
+        const MAX_ITERATIONS: usize = 50;
+
+        while changed && iterations < MAX_ITERATIONS {
+            changed = false;
+            iterations += 1;
+
+            let mut new_kernels: Vec<Kernel> = Vec::with_capacity(self.kernels.len());
+            let mut i = 0;
+
+            while i < self.kernels.len() {
+                if i + 1 < self.kernels.len() {
+                    let current = &self.kernels[i];
+                    let next = &self.kernels[i + 1];
+
+                    if let Some(fused) = current.fuse(next) {
+                        new_kernels.push(fused);
+                        i += 2;
+                        changed = true;
+                        continue;
+                    }
+                }
+
+                new_kernels.push(self.kernels[i].clone());
+                i += 1;
+            }
+
+            self.kernels = new_kernels;
+        }
+    }
+
+    fn build_execution_layers(&mut self) {
+        self.layers.clear();
+
+        for kernel in &self.kernels {
+            let mut placed = false;
+
+            for layer in &mut self.layers {
+                if layer.can_add(kernel) {
+                    layer.add(kernel.clone());
+                    placed = true;
+                    break;
+                }
+            }
+
+            if !placed {
+                let mut new_layer = ExecutionLayer::new();
+                new_layer.add(kernel.clone());
+                self.layers.push(new_layer);
+            }
+        }
+    }
+
+    pub fn execute(&self, state: &mut Vec<Complex<f64>>) {
+        for kernel in &self.kernels {
+            *state = apply_kernel(state, kernel, self.num_qubits);
+        }
+    }
+
+    pub fn execute_parallel(&self, state: &mut Vec<Complex<f64>>) {
+        for kernel in &self.kernels {
+            *state = apply_kernel_parallel(state, kernel, self.num_qubits);
+        }
+    }
+
+    pub fn execute_layered(&self, state: &mut Vec<Complex<f64>>) {
+        for layer in &self.layers {
+            for kernel in &layer.kernels {
+                *state = apply_kernel(state, kernel, self.num_qubits);
+            }
+        }
+    }
+
+    pub fn execute_layered_parallel(&self, state: &mut Vec<Complex<f64>>) {
+        for layer in &self.layers {
+            for kernel in &layer.kernels {
+                *state = apply_kernel_parallel(state, kernel, self.num_qubits);
+            }
+        }
+    }
+
+    pub fn execute_simd(&self, state: &mut Vec<Complex<f64>>) {
+        for kernel in &self.kernels {
+            if kernel.targets.len() == 1 {
+                let gate = matrix_to_2x2(&kernel.matrix);
+                apply_single_qubit_gate_simd(state, &gate, kernel.targets[0], self.num_qubits);
+            } else {
+                *state = apply_kernel(state, kernel, self.num_qubits);
+            }
+        }
+    }
+
+    pub fn execute_simd_parallel(&self, state: &mut Vec<Complex<f64>>) {
+        for kernel in &self.kernels {
+            if kernel.targets.len() == 1 && self.num_qubits >= 10 {
+                let gate = matrix_to_2x2(&kernel.matrix);
+                apply_single_qubit_gate_simd_parallel(
+                    state,
+                    &gate,
+                    kernel.targets[0],
+                    self.num_qubits,
+                );
+            } else if kernel.targets.len() == 1 {
+                let gate = matrix_to_2x2(&kernel.matrix);
+                apply_single_qubit_gate_simd(state, &gate, kernel.targets[0], self.num_qubits);
+            } else {
+                *state = apply_kernel_parallel(state, kernel, self.num_qubits);
+            }
+        }
+    }
+
+    pub fn stats(&self) -> KernelStats {
+        let single_qubit = self.kernels.iter().filter(|k| k.targets.len() == 1).count();
+        let two_qubit = self.kernels.iter().filter(|k| k.targets.len() == 2).count();
+        let multi_qubit = self.kernels.iter().filter(|k| k.targets.len() > 2).count();
+        let diagonal = self
+            .kernels
+            .iter()
+            .filter(|k| k.gate_type == GateType::Diagonal)
+            .count();
+
+        KernelStats {
+            total_kernels: self.kernels.len(),
+            single_qubit,
+            two_qubit,
+            multi_qubit,
+            diagonal,
+            execution_layers: self.layers.len(),
+        }
+    }
+}
+
+#[derive(Debug, Clone)]
+pub struct KernelStats {
+    pub total_kernels: usize,
+    pub single_qubit: usize,
+    pub two_qubit: usize,
+    pub multi_qubit: usize,
+    pub diagonal: usize,
+    pub execution_layers: usize,
+}
+
+impl std::fmt::Display for KernelStats {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        write!(
+            f,
+            "Kernels: {} (1q: {}, 2q: {}, 3q+: {}, diag: {}), Layers: {}",
+            self.total_kernels,
+            self.single_qubit,
+            self.two_qubit,
+            self.multi_qubit,
+            self.diagonal,
+            self.execution_layers
+        )
+    }
+}
diff --git a/src/core/mod.rs b/src/core/mod.rs
new file mode 100644
index 0000000..c38cf38
--- /dev/null
+++ b/src/core/mod.rs
@@ -0,0 +1,17 @@
+pub mod circuit;
+pub mod classical_components;
+pub mod custom_gate;
+pub mod gates;
+pub mod kernel;
+pub mod noise;
+pub mod quantum_components;
+pub mod runtime;
+
+pub use circuit::*;
+pub use classical_components::*;
+pub use custom_gate::*;
+pub use gates::*;
+pub use kernel::*;
+pub use noise::*;
+pub use quantum_components::*;
+pub use runtime::*;
diff --git a/src/core/noise.rs b/src/core/noise.rs
new file mode 100644
index 0000000..8d56953
--- /dev/null
+++ b/src/core/noise.rs
@@ -0,0 +1,559 @@
+use crate::{complex, Complex, Matrix};
+
+#[derive(Clone, Debug)]
+pub struct KrausOperator {
+    pub matrix: Matrix<Complex<f64>>,
+    pub name: String,
+}
+
+impl KrausOperator {
+    pub fn new(name: &str, matrix: Matrix<Complex<f64>>) -> Self {
+        Self {
+            matrix,
+            name: name.to_string(),
+        }
+    }
+}
+
+#[derive(Clone, Debug)]
+pub struct NoiseChannel {
+    pub name: String,
+    pub operators: Vec<KrausOperator>,
+    pub num_qubits: usize,
+}
+
+impl NoiseChannel {
+    pub fn new(name: &str, operators: Vec<KrausOperator>, num_qubits: usize) -> Self {
+        Self {
+            name: name.to_string(),
+            operators,
+            num_qubits,
+        }
+    }
+
+    pub fn depolarising(p: f64) -> Self {
+        let sqrt_1_p = (1.0 - p).sqrt();
+        let sqrt_p3 = (p / 3.0).sqrt();
+
+        let k0 = Matrix::new(
+            2,
+            2,
+            vec![
+                complex!(sqrt_1_p, 0.0),
+                complex!(0.0, 0.0),
+                complex!(0.0, 0.0),
+                complex!(sqrt_1_p, 0.0),
+            ],
+        );
+
+        let k1 = Matrix::new(
+            2,
+            2,
+            vec![
+                complex!(0.0, 0.0),
+                complex!(sqrt_p3, 0.0),
+                complex!(sqrt_p3, 0.0),
+                complex!(0.0, 0.0),
+            ],
+        );
+
+        let k2 = Matrix::new(
+            2,
+            2,
+            vec![
+                complex!(0.0, 0.0),
+                complex!(0.0, -sqrt_p3),
+                complex!(0.0, sqrt_p3),
+                complex!(0.0, 0.0),
+            ],
+        );
+
+        let k3 = Matrix::new(
+            2,
+            2,
+            vec![
+                complex!(sqrt_p3, 0.0),
+                complex!(0.0, 0.0),
+                complex!(0.0, 0.0),
+                complex!(-sqrt_p3, 0.0),
+            ],
+        );
+
+        Self::new(
+            "Depolarising",
+            vec![
+                KrausOperator::new("K0", k0),
+                KrausOperator::new("K1(X)", k1),
+                KrausOperator::new("K2(Y)", k2),
+                KrausOperator::new("K3(Z)", k3),
+            ],
+            1,
+        )
+    }
+
+    pub fn amplitude_damping(gamma: f64) -> Self {
+        let sqrt_gamma = gamma.sqrt();
+        let sqrt_1_gamma = (1.0 - gamma).sqrt();
+
+        let k0 = Matrix::new(
+            2,
+            2,
+            vec![
+                complex!(1.0, 0.0),
+                complex!(0.0, 0.0),
+                complex!(0.0, 0.0),
+                complex!(sqrt_1_gamma, 0.0),
+            ],
+        );
+
+        let k1 = Matrix::new(
+            2,
+            2,
+            vec![
+                complex!(0.0, 0.0),
+                complex!(sqrt_gamma, 0.0),
+                complex!(0.0, 0.0),
+                complex!(0.0, 0.0),
+            ],
+        );
+
+        Self::new(
+            "AmplitudeDamping",
+            vec![
+                KrausOperator::new("K0", k0),
+                KrausOperator::new("K1", k1),
+            ],
+            1,
+        )
+    }
+
+    pub fn phase_damping(gamma: f64) -> Self {
+        let sqrt_gamma = gamma.sqrt();
+        let sqrt_1_gamma = (1.0 - gamma).sqrt();
+
+        let k0 = Matrix::new(
+            2,
+            2,
+            vec![
+                complex!(1.0, 0.0),
+                complex!(0.0, 0.0),
+                complex!(0.0, 0.0),
+                complex!(sqrt_1_gamma, 0.0),
+            ],
+        );
+
+        let k1 = Matrix::new(
+            2,
+            2,
+            vec![
+                complex!(0.0, 0.0),
+                complex!(0.0, 0.0),
+                complex!(0.0, 0.0),
+                complex!(sqrt_gamma, 0.0),
+            ],
+        );
+
+        Self::new(
+            "PhaseDamping",
+            vec![
+                KrausOperator::new("K0", k0),
+                KrausOperator::new("K1", k1),
+            ],
+            1,
+        )
+    }
+
+    pub fn bit_flip(p: f64) -> Self {
+        let sqrt_1_p = (1.0 - p).sqrt();
+        let sqrt_p = p.sqrt();
+
+        let k0 = Matrix::new(
+            2,
+            2,
+            vec![
+                complex!(sqrt_1_p, 0.0),
+                complex!(0.0, 0.0),
+                complex!(0.0, 0.0),
+                complex!(sqrt_1_p, 0.0),
+            ],
+        );
+
+        let k1 = Matrix::new(
+            2,
+            2,
+            vec![
+                complex!(0.0, 0.0),
+                complex!(sqrt_p, 0.0),
+                complex!(sqrt_p, 0.0),
+                complex!(0.0, 0.0),
+            ],
+        );
+
+        Self::new(
+            "BitFlip",
+            vec![
+                KrausOperator::new("K0(I)", k0),
+                KrausOperator::new("K1(X)", k1),
+            ],
+            1,
+        )
+    }
+
+    pub fn phase_flip(p: f64) -> Self {
+        let sqrt_1_p = (1.0 - p).sqrt();
+        let sqrt_p = p.sqrt();
+
+        let k0 = Matrix::new(
+            2,
+            2,
+            vec![
+                complex!(sqrt_1_p, 0.0),
+                complex!(0.0, 0.0),
+                complex!(0.0, 0.0),
+                complex!(sqrt_1_p, 0.0),
+            ],
+        );
+
+        let k1 = Matrix::new(
+            2,
+            2,
+            vec![
+                complex!(sqrt_p, 0.0),
+                complex!(0.0, 0.0),
+                complex!(0.0, 0.0),
+                complex!(-sqrt_p, 0.0),
+            ],
+        );
+
+        Self::new(
+            "PhaseFlip",
+            vec![
+                KrausOperator::new("K0(I)", k0),
+                KrausOperator::new("K1(Z)", k1),
+            ],
+            1,
+        )
+    }
+
+    pub fn bit_phase_flip(p: f64) -> Self {
+        let sqrt_1_p = (1.0 - p).sqrt();
+        let sqrt_p = p.sqrt();
+
+        let k0 = Matrix::new(
+            2,
+            2,
+            vec![
+                complex!(sqrt_1_p, 0.0),
+                complex!(0.0, 0.0),
+                complex!(0.0, 0.0),
+                complex!(sqrt_1_p, 0.0),
+            ],
+        );
+
+        let k1 = Matrix::new(
+            2,
+            2,
+            vec![
+                complex!(0.0, 0.0),
+                complex!(0.0, -sqrt_p),
+                complex!(0.0, sqrt_p),
+                complex!(0.0, 0.0),
+            ],
+        );
+
+        Self::new(
+            "BitPhaseFlip",
+            vec![
+                KrausOperator::new("K0(I)", k0),
+                KrausOperator::new("K1(Y)", k1),
+            ],
+            1,
+        )
+    }
+
+    pub fn generalised_amplitude_damping(p: f64, gamma: f64) -> Self {
+        let sqrt_p = p.sqrt();
+        let sqrt_1_p = (1.0 - p).sqrt();
+        let sqrt_gamma = gamma.sqrt();
+        let sqrt_1_gamma = (1.0 - gamma).sqrt();
+
+        let k0 = Matrix::new(
+            2,
+            2,
+            vec![
+                complex!(sqrt_p, 0.0),
+                complex!(0.0, 0.0),
+                complex!(0.0, 0.0),
+                complex!(sqrt_p * sqrt_1_gamma, 0.0),
+            ],
+        );
+
+        let k1 = Matrix::new(
+            2,
+            2,
+            vec![
+                complex!(0.0, 0.0),
+                complex!(sqrt_p * sqrt_gamma, 0.0),
+                complex!(0.0, 0.0),
+                complex!(0.0, 0.0),
+            ],
+        );
+
+        let k2 = Matrix::new(
+            2,
+            2,
+            vec![
+                complex!(sqrt_1_p * sqrt_1_gamma, 0.0),
+                complex!(0.0, 0.0),
+                complex!(0.0, 0.0),
+                complex!(sqrt_1_p, 0.0),
+            ],
+        );
+
+        let k3 = Matrix::new(
+            2,
+            2,
+            vec![
+                complex!(0.0, 0.0),
+                complex!(0.0, 0.0),
+                complex!(sqrt_1_p * sqrt_gamma, 0.0),
+                complex!(0.0, 0.0),
+            ],
+        );
+
+        Self::new(
+            "GeneralisedAmplitudeDamping",
+            vec![
+                KrausOperator::new("K0", k0),
+                KrausOperator::new("K1", k1),
+                KrausOperator::new("K2", k2),
+                KrausOperator::new("K3", k3),
+            ],
+            1,
+        )
+    }
+}
+
+#[derive(Clone)]
+pub struct DensityMatrix {
+    pub data: Vec<Complex<f64>>,
+    pub dim: usize,
+    pub num_qubits: usize,
+}
+
+impl DensityMatrix {
+    pub fn new(num_qubits: usize) -> Self {
+        let dim = 1 << num_qubits;
+        let mut data = vec![complex!(0.0, 0.0); dim * dim];
+        data[0] = complex!(1.0, 0.0);
+        Self {
+            data,
+            dim,
+            num_qubits,
+        }
+    }
+
+    pub fn from_state_vector(state: &[Complex<f64>]) -> Self {
+        let dim = state.len();
+        let num_qubits = (dim as f64).log2() as usize;
+        let mut data = vec![complex!(0.0, 0.0); dim * dim];
+
+        for i in 0..dim {
+            for j in 0..dim {
+                data[i * dim + j] = state[i] * state[j].get_conjugate();
+            }
+        }
+
+        Self {
+            data,
+            dim,
+            num_qubits,
+        }
+    }
+
+    pub fn get(&self, row: usize, col: usize) -> Complex<f64> {
+        self.data[row * self.dim + col]
+    }
+
+    pub fn set(&mut self, row: usize, col: usize, value: Complex<f64>) {
+        self.data[row * self.dim + col] = value;
+    }
+
+    pub fn trace(&self) -> Complex<f64> {
+        let mut sum = complex!(0.0, 0.0);
+        for i in 0..self.dim {
+            sum += self.get(i, i);
+        }
+        sum
+    }
+
+    pub fn purity(&self) -> f64 {
+        let mut sum = complex!(0.0, 0.0);
+        for i in 0..self.dim {
+            for j in 0..self.dim {
+                let rho_ij = self.get(i, j);
+                let rho_ji = self.get(j, i);
+                sum += rho_ij * rho_ji;
+            }
+        }
+        sum.real
+    }
+
+    pub fn is_pure(&self, tolerance: f64) -> bool {
+        (self.purity() - 1.0).abs() < tolerance
+    }
+
+    pub fn probabilities(&self) -> Vec<f64> {
+        (0..self.dim).map(|i| self.get(i, i).real).collect()
+    }
+
+    pub fn apply_unitary(&mut self, gate: &Matrix<Complex<f64>>, targets: &[usize]) {
+        let g = targets.len();
+        let gate_dim = 1 << g;
+
+        let target_bits: Vec<usize> = targets
+            .iter()
+            .map(|&t| self.num_qubits - 1 - t)
+            .collect();
+
+        let mut non_target_mask: usize = (1 << self.num_qubits) - 1;
+        for &pos in &target_bits {
+            non_target_mask &= !(1 << pos);
+        }
+
+        let mut new_data = vec![complex!(0.0, 0.0); self.dim * self.dim];
+
+        for i in 0..self.dim {
+            for j in 0..self.dim {
+                let mut sum = complex!(0.0, 0.0);
+
+                for k in 0..gate_dim {
+                    for l in 0..gate_dim {
+                        let mut src_i = i & non_target_mask;
+                        let mut src_j = j & non_target_mask;
+
+                        for (idx, &pos) in target_bits.iter().enumerate() {
+                            if (k >> (g - 1 - idx)) & 1 == 1 {
+                                src_i |= 1 << pos;
+                            }
+                            if (l >> (g - 1 - idx)) & 1 == 1 {
+                                src_j |= 1 << pos;
+                            }
+                        }
+
+                        let mut tgt_i = 0usize;
+                        let mut tgt_j = 0usize;
+                        for (idx, &pos) in target_bits.iter().enumerate() {
+                            if (i >> pos) & 1 == 1 {
+                                tgt_i |= 1 << (g - 1 - idx);
+                            }
+                            if (j >> pos) & 1 == 1 {
+                                tgt_j |= 1 << (g - 1 - idx);
+                            }
+                        }
+
+                        let u_ik = gate.data[tgt_i * gate_dim + k];
+                        let u_jl_dag = gate.data[tgt_j * gate_dim + l].get_conjugate();
+                        let rho_kl = self.get(src_i, src_j);
+
+                        sum += u_ik * rho_kl * u_jl_dag;
+                    }
+                }
+
+                new_data[i * self.dim + j] = sum;
+            }
+        }
+
+        self.data = new_data;
+    }
+
+    pub fn apply_noise_channel(&mut self, channel: &NoiseChannel, target: usize) {
+        if channel.num_qubits != 1 {
+            panic!("Only single-qubit noise channels are currently supported");
+        }
+
+        let target_bit = self.num_qubits - 1 - target;
+        let mut new_data = vec![complex!(0.0, 0.0); self.dim * self.dim];
+
+        for kraus in &channel.operators {
+            let k = &kraus.matrix;
+
+            for i in 0..self.dim {
+                for j in 0..self.dim {
+                    let i_target = (i >> target_bit) & 1;
+                    let j_target = (j >> target_bit) & 1;
+
+                    for ki in 0..2 {
+                        for kj in 0..2 {
+                            let src_i = (i & !(1 << target_bit)) | (ki << target_bit);
+                            let src_j = (j & !(1 << target_bit)) | (kj << target_bit);
+
+                            let k_elem = k.data[i_target * 2 + ki];
+                            let k_dag_elem = k.data[j_target * 2 + kj].get_conjugate();
+                            let rho_elem = self.get(src_i, src_j);
+
+                            new_data[i * self.dim + j] += k_elem * rho_elem * k_dag_elem;
+                        }
+                    }
+                }
+            }
+        }
+
+        self.data = new_data;
+    }
+
+    pub fn measure_probability(&self, qubit: usize, outcome: usize) -> f64 {
+        let target_bit = self.num_qubits - 1 - qubit;
+        let mut prob = 0.0;
+
+        for i in 0..self.dim {
+            if (i >> target_bit) & 1 == outcome {
+                prob += self.get(i, i).real;
+            }
+        }
+
+        prob
+    }
+
+    pub fn fidelity_with_pure_state(&self, state: &[Complex<f64>]) -> f64 {
+        let mut sum = complex!(0.0, 0.0);
+
+        for i in 0..self.dim {
+            for j in 0..self.dim {
+                sum += state[i].get_conjugate() * self.get(i, j) * state[j];
+            }
+        }
+
+        sum.real
+    }
+}
+
+impl std::fmt::Display for DensityMatrix {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        writeln!(f, "DensityMatrix ({} qubits, {}×{}):", self.num_qubits, self.dim, self.dim)?;
+        writeln!(f, "  Trace: {:.6}", self.trace().real)?;
+        writeln!(f, "  Purity: {:.6}", self.purity())?;
+        writeln!(f, "  Pure: {}", self.is_pure(1e-10))?;
+        writeln!(f, "  Probabilities: {:?}", self.probabilities())?;
+        Ok(())
+    }
+}
+
+impl std::fmt::Debug for DensityMatrix {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        writeln!(f, "DensityMatrix {}×{}:", self.dim, self.dim)?;
+        for i in 0..self.dim {
+            write!(f, "  [")?;
+            for j in 0..self.dim {
+                let val = self.get(i, j);
+                if j > 0 {
+                    write!(f, ", ")?;
+                }
+                write!(f, "{:.4}+{:.4}i", val.real, val.imaginary)?;
+            }
+            writeln!(f, "]")?;
+        }
+        Ok(())
+    }
+}
+
diff --git a/src/core/quantum_components.rs b/src/core/quantum_components.rs
new file mode 100644
index 0000000..f07b8ef
--- /dev/null
+++ b/src/core/quantum_components.rs
@@ -0,0 +1,318 @@
+use crate::{column_vector, complex, ColumnVector, Complex, Float, Matrix, Vector, VectorMatrix};
+use core::{fmt, ops};
+
+#[macro_export]
+macro_rules! count {
+    () => { 0 };
+    ($head:expr $(,$tail:expr)*) => { 1 + count!($( $tail ),*) };
+}
+
+#[macro_export]
+macro_rules! qubit {
+    ($(($re:expr, $im:expr)),*) => {
+        {
+            let mut vector = Vec::new();
+            $(
+                vector.push(complex!($re, $im));
+            )*
+            QuantumBit::new(vector)
+        }
+    };
+}
+
+#[macro_export]
+macro_rules! quantum_register {
+    ($($bit:expr),*) => {
+        {
+            const N: usize = count!($($bit),*);
+            let mut bits: [QuantumBit; N] = [$($bit),*];
+            QuantumRegister::from(&mut bits)
+        }
+    };
+}
+
+pub type QuantumState = ColumnVector<Complex<f64>>;
+impl QuantumState {
+    pub fn state_0() -> QuantumState {
+        column_vector![complex!(1.0, 0.0), complex!(0.0, 0.0)]
+    }
+
+    pub fn state_1() -> QuantumState {
+        column_vector![complex!(0.0, 0.0), complex!(1.0, 0.0)]
+    }
+}
+
+fn identity_matrix<T: Float>(size: usize) -> Matrix<T> {
+    let mut data = vec![T::zero(); size * size];
+    for i in 0..size {
+        data[i * size + i] = T::one();
+    }
+    Matrix::new(size, size, data)
+}
+
+#[derive(Clone)]
+pub struct QuantumBit<'a> {
+    state: QuantumState,
+    name: &'a str,
+}
+
+#[derive(Clone)]
+pub struct QuantumRegister<'a> {
+    state_vector: QuantumState,
+    name: &'a str,
+    qubits: Vec<QuantumBit<'a>>,
+}
+
+#[derive(Clone)]
+pub struct QuantumGate<'a> {
+    pub name: &'a str,
+    pub matrix: Matrix<Complex<f64>>,
+    pub num_qubits: usize,
+}
+
+impl<'a> QuantumGate<'a> {
+    pub fn new(name: &'a str, matrix: Matrix<Complex<f64>>, num_qubits: usize) -> Self {
+        let expected_dim = 1 << num_qubits;
+        assert_eq!(
+            matrix.rows, expected_dim,
+            "Gate matrix rows must be 2^num_qubits"
+        );
+        assert_eq!(
+            matrix.cols, expected_dim,
+            "Gate matrix cols must be 2^num_qubits"
+        );
+        QuantumGate {
+            name,
+            matrix,
+            num_qubits,
+        }
+    }
+
+    pub fn from_matrix(name: &'a str, matrix: Matrix<Complex<f64>>) -> Self {
+        assert_eq!(matrix.rows, matrix.cols, "Gate matrix must be square");
+        let dim = matrix.rows;
+        assert!(
+            dim > 0 && (dim & (dim - 1)) == 0,
+            "Matrix dimension must be a power of 2"
+        );
+        let num_qubits = (dim as f64).log2() as usize;
+        QuantumGate {
+            name,
+            matrix,
+            num_qubits,
+        }
+    }
+}
+
+impl<'a> QuantumBit<'a> {
+    pub fn new(name: &'a str, state: QuantumState) -> QuantumBit<'a> {
+        QuantumBit { name, state }
+    }
+
+    pub fn get_state(&self) -> QuantumState {
+        self.state.clone()
+    }
+
+    pub fn get_name(&self) -> &'a str {
+        self.name
+    }
+}
+
+impl<'a> QuantumRegister<'a> {
+    pub fn new(name: &'a str, names: &[&'a str]) -> QuantumRegister<'a> {
+        let mut bits: Vec<QuantumBit<'a>> = Vec::new();
+        for &name in names {
+            bits.push(QuantumBit::new(name, QuantumState::state_0()))
+        }
+
+        QuantumRegister::from(name, &mut bits)
+    }
+
+    pub fn from(name: &'a str, bits: &mut [QuantumBit<'a>]) -> QuantumRegister<'a> {
+        let mut register = QuantumRegister {
+            name,
+            qubits: bits.to_vec(),
+            state_vector: ColumnVector::new(vec![]),
+        };
+
+        register.update();
+        register
+    }
+
+    fn update(&mut self) {
+        let matrices: Vec<Matrix<Complex<f64>>> = self
+            .qubits
+            .iter()
+            .map(|qubit| qubit.state.to_matrix())
+            .collect();
+        let mut new_result = matrices[0].clone();
+        for matrix in &matrices[1..] {
+            new_result = new_result.kronecker(matrix);
+        }
+
+        self.state_vector = ColumnVector::from_matrix(&new_result);
+    }
+
+    pub fn get_bits(&self) -> Vec<QuantumBit<'_>> {
+        self.qubits.clone()
+    }
+
+    pub fn get_state(&self) -> QuantumState {
+        self.state_vector.clone()
+    }
+
+    pub fn get_name(&self) -> &'a str {
+        self.name
+    }
+
+    pub fn num_qubits(&self) -> usize {
+        self.qubits.len()
+    }
+
+    pub fn apply_gate(&mut self, gate: &QuantumGate, targets: &[usize]) {
+        let n = self.num_qubits();
+
+        assert_eq!(
+            gate.num_qubits,
+            targets.len(),
+            "Number of target qubits must match gate's qubit count"
+        );
+        for &t in targets {
+            assert!(
+                t < n,
+                "Target qubit index {} out of range for {}-qubit register",
+                t,
+                n
+            );
+        }
+
+        let mut sorted_targets = targets.to_vec();
+        sorted_targets.sort();
+        for i in 1..sorted_targets.len() {
+            assert_ne!(
+                sorted_targets[i],
+                sorted_targets[i - 1],
+                "Duplicate target qubit indices are not allowed"
+            );
+        }
+
+        let full_operator = self.build_full_operator(gate, targets);
+
+        self.state_vector = self
+            .state_vector
+            .mul_matrix(&full_operator)
+            .expect("Matrix multiplication failed during gate application");
+    }
+
+    fn build_full_operator(&self, gate: &QuantumGate, targets: &[usize]) -> Matrix<Complex<f64>> {
+        let n = self.num_qubits();
+        let g = gate.num_qubits;
+        let dim = 1 << n;
+
+        let mut contiguous = true;
+        for i in 1..targets.len() {
+            if targets[i] != targets[i - 1] + 1 {
+                contiguous = false;
+                break;
+            }
+        }
+
+        if contiguous && g == n {
+            return gate.matrix.clone();
+        }
+
+        if contiguous {
+            return self.build_contiguous_operator(gate, targets[0]);
+        }
+
+        let mut result = Matrix::new(dim, dim, vec![complex!(0.0, 0.0); dim * dim]);
+
+        for col in 0..dim {
+            for row in 0..dim {
+                let mut target_row_bits = 0usize;
+                let mut target_col_bits = 0usize;
+
+                for (i, &t) in targets.iter().enumerate() {
+                    let qubit_pos = n - 1 - t;
+                    if (row >> qubit_pos) & 1 == 1 {
+                        target_row_bits |= 1 << (g - 1 - i);
+                    }
+                    if (col >> qubit_pos) & 1 == 1 {
+                        target_col_bits |= 1 << (g - 1 - i);
+                    }
+                }
+
+                let mut non_target_match = true;
+                for q in 0..n {
+                    if !targets.contains(&q) {
+                        let qubit_pos = n - 1 - q;
+                        if ((row >> qubit_pos) & 1) != ((col >> qubit_pos) & 1) {
+                            non_target_match = false;
+                            break;
+                        }
+                    }
+                }
+
+                if non_target_match {
+                    result.set(row, col, gate.matrix.get(target_row_bits, target_col_bits));
+                }
+            }
+        }
+
+        result
+    }
+
+    fn build_contiguous_operator(
+        &self,
+        gate: &QuantumGate,
+        start_idx: usize,
+    ) -> Matrix<Complex<f64>> {
+        let n = self.num_qubits();
+        let g = gate.num_qubits;
+
+        let mut result: Option<Matrix<Complex<f64>>> = None;
+
+        for i in 0..n {
+            let part: Matrix<Complex<f64>> = if i == start_idx {
+                gate.matrix.clone()
+            } else if i > start_idx && i < start_idx + g {
+                continue;
+            } else {
+                identity_matrix(2)
+            };
+
+            result = Some(match result {
+                None => part,
+                Some(r) => r.kronecker(&part),
+            });
+        }
+
+        result.unwrap_or_else(|| identity_matrix(1 << n))
+    }
+
+    pub fn apply_gates(&mut self, operations: &[(&QuantumGate, &[usize])]) {
+        for (gate, targets) in operations {
+            self.apply_gate(gate, targets);
+        }
+    }
+}
+
+impl<'a> ops::Index<usize> for QuantumRegister<'a> {
+    type Output = QuantumBit<'a>;
+
+    fn index(&self, index: usize) -> &Self::Output {
+        &self.qubits[index]
+    }
+}
+
+impl<'a> ops::IndexMut<usize> for QuantumRegister<'a> {
+    fn index_mut(&mut self, index: usize) -> &mut Self::Output {
+        &mut self.qubits[index]
+    }
+}
+
+impl<'a> fmt::Display for QuantumGate<'a> {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        write!(f, "{}", self.name)
+    }
+}
diff --git a/src/core/runtime.rs b/src/core/runtime.rs
new file mode 100644
index 0000000..6ea7b76
--- /dev/null
+++ b/src/core/runtime.rs
@@ -0,0 +1,585 @@
+use super::{
+    GateOp, Kernel, KernelBatch, QuantumGate, QuantumRegister, QuantumState,
+    StructureAwareKernelBatch,
+};
+use crate::gates::{
+    cp_matrix, crx_matrix, cry_matrix, crz_matrix, p_matrix, rx_matrix, ry_matrix, rz_matrix,
+    u1_matrix, u2_matrix, u3_matrix, CNOT, CZ, FREDKIN, HADAMARD, PAULI_X, PAULI_Y, PAULI_Z,
+    SDG_GATE, SWAP, SXDG_GATE, SX_GATE, S_GATE, TDG_GATE, TOFFOLI, T_GATE,
+};
+use crate::maths::simd::{apply_single_qubit_gate_simd, apply_single_qubit_gate_simd_parallel};
+use crate::maths::vector::Vector;
+use crate::{complex, Complex, Matrix};
+use rayon::prelude::*;
+
+const PARALLEL_THRESHOLD: usize = 8;
+
+#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
+pub struct RuntimeConfig {
+    pub parallel: bool,
+    pub simd: bool,
+    pub batched: bool,
+    pub structure_aware: bool,
+    pub parallel_threshold: usize,
+}
+
+impl RuntimeConfig {
+    pub fn new() -> Self {
+        Self {
+            parallel: false,
+            simd: false,
+            batched: false,
+            structure_aware: false,
+            parallel_threshold: PARALLEL_THRESHOLD,
+        }
+    }
+
+    pub fn parallel(mut self) -> Self {
+        self.parallel = true;
+        self
+    }
+
+    pub fn simd(mut self) -> Self {
+        self.simd = true;
+        self
+    }
+
+    pub fn batched(mut self) -> Self {
+        self.batched = true;
+        self
+    }
+
+    pub fn structure_aware(mut self) -> Self {
+        self.structure_aware = true;
+        self
+    }
+
+    pub fn with_threshold(mut self, threshold: usize) -> Self {
+        self.parallel_threshold = threshold;
+        self
+    }
+
+    pub fn optimal() -> Self {
+        Self::new().structure_aware().simd().parallel()
+    }
+
+    pub fn compute(&self, num_qubits: usize, operations: &[GateOp]) -> QuantumState {
+        let dim = 1 << num_qubits;
+        let mut state: Vec<Complex<f64>> = vec![complex!(0.0, 0.0); dim];
+        state[0] = complex!(1.0, 0.0);
+
+        let use_parallel = self.parallel && num_qubits >= self.parallel_threshold;
+
+        if self.structure_aware {
+            let mut batch = Runtime::build_structure_aware_batch(num_qubits, operations);
+            batch.optimise();
+            self.execute_kernels(&mut state, batch.kernels(), num_qubits, use_parallel);
+        } else if self.batched {
+            let mut batch = Runtime::build_kernel_batch(num_qubits, operations);
+            batch.optimize();
+            self.execute_kernels(&mut state, batch.kernels(), num_qubits, use_parallel);
+        } else {
+            let batch = Runtime::build_kernel_batch(num_qubits, operations);
+            self.execute_kernels(&mut state, batch.kernels(), num_qubits, use_parallel);
+        }
+
+        QuantumState::new(state)
+    }
+
+    fn execute_kernels(
+        &self,
+        state: &mut Vec<Complex<f64>>,
+        kernels: &[Kernel],
+        num_qubits: usize,
+        use_parallel: bool,
+    ) {
+        for kernel in kernels {
+            if self.simd && kernel.targets.len() == 1 {
+                let gate = matrix_to_2x2(&kernel.matrix);
+                if use_parallel {
+                    apply_single_qubit_gate_simd_parallel(
+                        state,
+                        &gate,
+                        kernel.targets[0],
+                        num_qubits,
+                    );
+                } else {
+                    apply_single_qubit_gate_simd(state, &gate, kernel.targets[0], num_qubits);
+                }
+            } else if use_parallel {
+                *state = apply_gate_parallel(state, &kernel.matrix, &kernel.targets, num_qubits);
+            } else {
+                *state = apply_kernel_direct(state, kernel, num_qubits);
+            }
+        }
+    }
+}
+
+impl std::fmt::Display for RuntimeConfig {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        let mut features = Vec::new();
+        if self.structure_aware {
+            features.push("structure-aware");
+        }
+        if self.batched && !self.structure_aware {
+            features.push("batched");
+        }
+        if self.simd {
+            features.push("SIMD");
+        }
+        if self.parallel {
+            features.push("parallel");
+        }
+        if features.is_empty() {
+            features.push("basic");
+        }
+        write!(f, "Runtime[{}]", features.join("+"))
+    }
+}
+
+#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
+pub enum Runtime {
+    #[default]
+    BasicRT,
+    BasicRTMT,
+    BatchedRT,
+    BatchedRTMT,
+    SimdRT,
+    SimdRTMT,
+    StructureAwareRT,
+    StructureAwareMT,
+    WFEvolution,
+    WFEvolutionMT,
+    GPUAccelerated,
+    Custom(RuntimeConfig),
+}
+
+impl Runtime {
+    pub fn custom() -> RuntimeConfig {
+        RuntimeConfig::new()
+    }
+
+    pub fn optimal() -> RuntimeConfig {
+        RuntimeConfig::optimal()
+    }
+
+    pub fn to_config(&self) -> RuntimeConfig {
+        match self {
+            Runtime::BasicRT => RuntimeConfig::new(),
+            Runtime::BasicRTMT => RuntimeConfig::new().parallel(),
+            Runtime::BatchedRT => RuntimeConfig::new().batched(),
+            Runtime::BatchedRTMT => RuntimeConfig::new().batched().parallel(),
+            Runtime::SimdRT => RuntimeConfig::new().batched().simd(),
+            Runtime::SimdRTMT => RuntimeConfig::new().batched().simd().parallel(),
+            Runtime::StructureAwareRT => RuntimeConfig::new().structure_aware().simd(),
+            Runtime::StructureAwareMT => RuntimeConfig::new().structure_aware().simd().parallel(),
+            Runtime::Custom(config) => *config,
+            _ => RuntimeConfig::new(),
+        }
+    }
+
+    pub fn compute(&self, num_qubits: usize, operations: &[GateOp]) -> QuantumState {
+        match self {
+            Runtime::BasicRT => Self::compute_basic(num_qubits, operations),
+            Runtime::BasicRTMT => Self::compute_basic_mt(num_qubits, operations),
+            Runtime::Custom(config) => config.compute(num_qubits, operations),
+            Runtime::WFEvolution => {
+                unimplemented!("WFEvolution (Schrödinger equation) runtime not yet implemented")
+            }
+            Runtime::WFEvolutionMT => {
+                unimplemented!(
+                    "WFEvolutionMT (multi-threaded Schrödinger) runtime not yet implemented"
+                )
+            }
+            Runtime::GPUAccelerated => {
+                unimplemented!("GPUAccelerated runtime not yet implemented")
+            }
+            _ => self.to_config().compute(num_qubits, operations),
+        }
+    }
+
+    pub fn build_kernel_batch(num_qubits: usize, operations: &[GateOp]) -> KernelBatch {
+        let mut batch = KernelBatch::new(num_qubits);
+
+        for op in operations {
+            if let Some(kernel) = Self::op_to_kernel(op) {
+                batch.add(kernel);
+            }
+        }
+
+        batch
+    }
+
+    fn op_to_kernel(op: &GateOp) -> Option<Kernel> {
+        let (matrix, targets, name): (Matrix<Complex<f64>>, Vec<usize>, &str) = match op {
+            GateOp::H(t) => (HADAMARD.matrix.clone(), vec![*t], "H"),
+            GateOp::X(t) => (PAULI_X.matrix.clone(), vec![*t], "X"),
+            GateOp::Y(t) => (PAULI_Y.matrix.clone(), vec![*t], "Y"),
+            GateOp::Z(t) => (PAULI_Z.matrix.clone(), vec![*t], "Z"),
+            GateOp::S(t) => (S_GATE.matrix.clone(), vec![*t], "S"),
+            GateOp::T(t) => (T_GATE.matrix.clone(), vec![*t], "T"),
+            GateOp::Sdg(t) => (SDG_GATE.matrix.clone(), vec![*t], "Sdg"),
+            GateOp::Tdg(t) => (TDG_GATE.matrix.clone(), vec![*t], "Tdg"),
+            GateOp::Sx(t) => (SX_GATE.matrix.clone(), vec![*t], "Sx"),
+            GateOp::Sxdg(t) => (SXDG_GATE.matrix.clone(), vec![*t], "Sxdg"),
+            GateOp::Rx(t, theta) => (rx_matrix(*theta), vec![*t], "Rx"),
+            GateOp::Ry(t, theta) => (ry_matrix(*theta), vec![*t], "Ry"),
+            GateOp::Rz(t, theta) => (rz_matrix(*theta), vec![*t], "Rz"),
+            GateOp::P(t, theta) => (p_matrix(*theta), vec![*t], "P"),
+            GateOp::U1(t, lambda) => (u1_matrix(*lambda), vec![*t], "U1"),
+            GateOp::U2(t, phi, lambda) => (u2_matrix(*phi, *lambda), vec![*t], "U2"),
+            GateOp::U3(t, theta, phi, lambda) => (u3_matrix(*theta, *phi, *lambda), vec![*t], "U3"),
+            GateOp::CNOT(c, t) => (CNOT.matrix.clone(), vec![*c, *t], "CNOT"),
+            GateOp::CZ(c, t) => (CZ.matrix.clone(), vec![*c, *t], "CZ"),
+            GateOp::SWAP(a, b) => (SWAP.matrix.clone(), vec![*a, *b], "SWAP"),
+            GateOp::CRx(c, t, theta) => (crx_matrix(*theta), vec![*c, *t], "CRx"),
+            GateOp::CRy(c, t, theta) => (cry_matrix(*theta), vec![*c, *t], "CRy"),
+            GateOp::CRz(c, t, theta) => (crz_matrix(*theta), vec![*c, *t], "CRz"),
+            GateOp::CP(c, t, theta) => (cp_matrix(*theta), vec![*c, *t], "CP"),
+            GateOp::CCNOT(c1, c2, t) => (TOFFOLI.matrix.clone(), vec![*c1, *c2, *t], "CCNOT"),
+            GateOp::CSWAP(c, t1, t2) => (FREDKIN.matrix.clone(), vec![*c, *t1, *t2], "CSWAP"),
+            GateOp::Measure(_, _) => return None,
+            GateOp::Custom(gate, tgts) => {
+                let qg = gate.to_quantum_gate();
+                (qg.matrix, tgts.clone(), "Custom")
+            }
+        };
+
+        Some(Kernel::new(name, matrix, targets))
+    }
+
+    pub fn build_structure_aware_batch(
+        num_qubits: usize,
+        operations: &[GateOp],
+    ) -> StructureAwareKernelBatch {
+        let mut batch = StructureAwareKernelBatch::new(num_qubits);
+
+        for op in operations {
+            if let Some(kernel) = Self::op_to_kernel(op) {
+                batch.add(kernel);
+            }
+        }
+
+        batch
+    }
+
+    fn compute_basic(num_qubits: usize, operations: &[GateOp]) -> QuantumState {
+        let names: Vec<String> = (0..num_qubits).map(|i| format!("q{}", i)).collect();
+        let leaked_names: &'static [String] = Box::leak(names.into_boxed_slice());
+        let name_refs: Vec<&'static str> = leaked_names.iter().map(|s| s.as_str()).collect();
+
+        let mut register = QuantumRegister::new(
+            Box::leak(Box::new("circuit".to_string())).as_str(),
+            &name_refs,
+        );
+
+        for op in operations {
+            match op {
+                // Clifford gates
+                GateOp::H(t) => register.apply_gate(&HADAMARD, &[*t]),
+                GateOp::X(t) => register.apply_gate(&PAULI_X, &[*t]),
+                GateOp::Y(t) => register.apply_gate(&PAULI_Y, &[*t]),
+                GateOp::Z(t) => register.apply_gate(&PAULI_Z, &[*t]),
+                GateOp::S(t) => register.apply_gate(&S_GATE, &[*t]),
+                GateOp::CNOT(c, t) => register.apply_gate(&CNOT, &[*c, *t]),
+                GateOp::CZ(c, t) => register.apply_gate(&CZ, &[*c, *t]),
+                GateOp::SWAP(a, b) => register.apply_gate(&SWAP, &[*a, *b]),
+                GateOp::CCNOT(c1, c2, t) => register.apply_gate(&TOFFOLI, &[*c1, *c2, *t]),
+                GateOp::CSWAP(c, t1, t2) => register.apply_gate(&FREDKIN, &[*c, *t1, *t2]),
+
+                // Non-Clifford fixed gates
+                GateOp::T(t) => register.apply_gate(&T_GATE, &[*t]),
+                GateOp::Sdg(t) => register.apply_gate(&SDG_GATE, &[*t]),
+                GateOp::Tdg(t) => register.apply_gate(&TDG_GATE, &[*t]),
+                GateOp::Sx(t) => register.apply_gate(&SX_GATE, &[*t]),
+                GateOp::Sxdg(t) => register.apply_gate(&SXDG_GATE, &[*t]),
+
+                // Parametric single-qubit gates (non-Clifford for most angles)
+                GateOp::Rx(t, theta) => {
+                    let gate = QuantumGate {
+                        name: "Rx",
+                        matrix: rx_matrix(*theta),
+                        num_qubits: 1,
+                    };
+                    register.apply_gate(&gate, &[*t]);
+                }
+                GateOp::Ry(t, theta) => {
+                    let gate = QuantumGate {
+                        name: "Ry",
+                        matrix: ry_matrix(*theta),
+                        num_qubits: 1,
+                    };
+                    register.apply_gate(&gate, &[*t]);
+                }
+                GateOp::Rz(t, theta) => {
+                    let gate = QuantumGate {
+                        name: "Rz",
+                        matrix: rz_matrix(*theta),
+                        num_qubits: 1,
+                    };
+                    register.apply_gate(&gate, &[*t]);
+                }
+                GateOp::P(t, theta) => {
+                    let gate = QuantumGate {
+                        name: "P",
+                        matrix: p_matrix(*theta),
+                        num_qubits: 1,
+                    };
+                    register.apply_gate(&gate, &[*t]);
+                }
+                GateOp::U1(t, lambda) => {
+                    let gate = QuantumGate {
+                        name: "U1",
+                        matrix: u1_matrix(*lambda),
+                        num_qubits: 1,
+                    };
+                    register.apply_gate(&gate, &[*t]);
+                }
+                GateOp::U2(t, phi, lambda) => {
+                    let gate = QuantumGate {
+                        name: "U2",
+                        matrix: u2_matrix(*phi, *lambda),
+                        num_qubits: 1,
+                    };
+                    register.apply_gate(&gate, &[*t]);
+                }
+                GateOp::U3(t, theta, phi, lambda) => {
+                    let gate = QuantumGate {
+                        name: "U3",
+                        matrix: u3_matrix(*theta, *phi, *lambda),
+                        num_qubits: 1,
+                    };
+                    register.apply_gate(&gate, &[*t]);
+                }
+
+                // Controlled parametric gates
+                GateOp::CRx(c, t, theta) => {
+                    let gate = QuantumGate {
+                        name: "CRx",
+                        matrix: crx_matrix(*theta),
+                        num_qubits: 2,
+                    };
+                    register.apply_gate(&gate, &[*c, *t]);
+                }
+                GateOp::CRy(c, t, theta) => {
+                    let gate = QuantumGate {
+                        name: "CRy",
+                        matrix: cry_matrix(*theta),
+                        num_qubits: 2,
+                    };
+                    register.apply_gate(&gate, &[*c, *t]);
+                }
+                GateOp::CRz(c, t, theta) => {
+                    let gate = QuantumGate {
+                        name: "CRz",
+                        matrix: crz_matrix(*theta),
+                        num_qubits: 2,
+                    };
+                    register.apply_gate(&gate, &[*c, *t]);
+                }
+                GateOp::CP(c, t, theta) => {
+                    let gate = QuantumGate {
+                        name: "CP",
+                        matrix: cp_matrix(*theta),
+                        num_qubits: 2,
+                    };
+                    register.apply_gate(&gate, &[*c, *t]);
+                }
+
+                // Measurement and custom gates
+                GateOp::Measure(_, _) => {}
+                GateOp::Custom(gate, targets) => {
+                    let quantum_gate = gate.to_quantum_gate();
+                    register.apply_gate(&quantum_gate, targets);
+                }
+            }
+        }
+
+        register.get_state()
+    }
+
+    fn compute_basic_mt(num_qubits: usize, operations: &[GateOp]) -> QuantumState {
+        // For small circuits, fall back to single-threaded (overhead not worth it)
+        if num_qubits < PARALLEL_THRESHOLD {
+            return Self::compute_basic(num_qubits, operations);
+        }
+
+        let dim = 1 << num_qubits;
+
+        // Initialize state to |0...0⟩
+        let mut state: Vec<Complex<f64>> = vec![complex!(0.0, 0.0); dim];
+        state[0] = complex!(1.0, 0.0);
+
+        for op in operations {
+            let (gate_matrix, targets): (Matrix<Complex<f64>>, Vec<usize>) = match op {
+                // Clifford gates
+                GateOp::H(t) => (HADAMARD.matrix.clone(), vec![*t]),
+                GateOp::X(t) => (PAULI_X.matrix.clone(), vec![*t]),
+                GateOp::Y(t) => (PAULI_Y.matrix.clone(), vec![*t]),
+                GateOp::Z(t) => (PAULI_Z.matrix.clone(), vec![*t]),
+                GateOp::S(t) => (S_GATE.matrix.clone(), vec![*t]),
+                GateOp::CNOT(c, t) => (CNOT.matrix.clone(), vec![*c, *t]),
+                GateOp::CZ(c, t) => (CZ.matrix.clone(), vec![*c, *t]),
+                GateOp::SWAP(a, b) => (SWAP.matrix.clone(), vec![*a, *b]),
+                GateOp::CCNOT(c1, c2, t) => (TOFFOLI.matrix.clone(), vec![*c1, *c2, *t]),
+                GateOp::CSWAP(c, t1, t2) => (FREDKIN.matrix.clone(), vec![*c, *t1, *t2]),
+
+                // Non-Clifford fixed gates
+                GateOp::T(t) => (T_GATE.matrix.clone(), vec![*t]),
+                GateOp::Sdg(t) => (SDG_GATE.matrix.clone(), vec![*t]),
+                GateOp::Tdg(t) => (TDG_GATE.matrix.clone(), vec![*t]),
+                GateOp::Sx(t) => (SX_GATE.matrix.clone(), vec![*t]),
+                GateOp::Sxdg(t) => (SXDG_GATE.matrix.clone(), vec![*t]),
+
+                // Parametric single-qubit gates
+                GateOp::Rx(t, theta) => (rx_matrix(*theta), vec![*t]),
+                GateOp::Ry(t, theta) => (ry_matrix(*theta), vec![*t]),
+                GateOp::Rz(t, theta) => (rz_matrix(*theta), vec![*t]),
+                GateOp::P(t, theta) => (p_matrix(*theta), vec![*t]),
+                GateOp::U1(t, lambda) => (u1_matrix(*lambda), vec![*t]),
+                GateOp::U2(t, phi, lambda) => (u2_matrix(*phi, *lambda), vec![*t]),
+                GateOp::U3(t, theta, phi, lambda) => (u3_matrix(*theta, *phi, *lambda), vec![*t]),
+
+                // Controlled parametric gates
+                GateOp::CRx(c, t, theta) => (crx_matrix(*theta), vec![*c, *t]),
+                GateOp::CRy(c, t, theta) => (cry_matrix(*theta), vec![*c, *t]),
+                GateOp::CRz(c, t, theta) => (crz_matrix(*theta), vec![*c, *t]),
+                GateOp::CP(c, t, theta) => (cp_matrix(*theta), vec![*c, *t]),
+
+                // Measurement (skip) and custom gates
+                GateOp::Measure(_, _) => continue,
+                GateOp::Custom(custom_gate, tgts) => {
+                    let quantum_gate = custom_gate.to_quantum_gate();
+                    state = apply_gate_parallel(&state, &quantum_gate.matrix, tgts, num_qubits);
+                    continue;
+                }
+            };
+
+            state = apply_gate_parallel(&state, &gate_matrix, &targets, num_qubits);
+        }
+
+        QuantumState::new(state)
+    }
+}
+
+/// Apply a gate to the state vector in parallel using sparse application
+/// This is O(2^n * 2^g) instead of O(2^2n) for full matrix multiplication
+fn apply_gate_parallel(
+    state: &[Complex<f64>],
+    gate_matrix: &Matrix<Complex<f64>>,
+    targets: &[usize],
+    num_qubits: usize,
+) -> Vec<Complex<f64>> {
+    let dim = 1 << num_qubits;
+    let g = targets.len();
+    let gate_dim = 1 << g;
+
+    // Convert target qubit indices to bit positions (from MSB)
+    let target_bits: Vec<usize> = targets.iter().map(|&t| num_qubits - 1 - t).collect();
+
+    // Create a mask for non-target qubits
+    let mut non_target_mask: usize = (1 << num_qubits) - 1;
+    for &pos in &target_bits {
+        non_target_mask &= !(1 << pos);
+    }
+
+    // Parallel computation of new state
+    let new_state: Vec<Complex<f64>> = (0..dim)
+        .into_par_iter()
+        .map(|i| {
+            // Extract the target qubit bits from index i
+            let mut target_idx = 0usize;
+            for (k, &pos) in target_bits.iter().enumerate() {
+                if (i >> pos) & 1 == 1 {
+                    target_idx |= 1 << (g - 1 - k);
+                }
+            }
+
+            // Compute the contribution to state[i]
+            let mut sum = complex!(0.0, 0.0);
+
+            // For each possible input state that could contribute
+            for j in 0..gate_dim {
+                // Get the gate matrix element
+                let gate_elem = gate_matrix.data[target_idx * gate_dim + j];
+
+                // Skip if zero (sparse optimization)
+                if gate_elem.real.abs() < 1e-15 && gate_elem.imaginary.abs() < 1e-15 {
+                    continue;
+                }
+
+                // Compute the source index by replacing target bits in i with bits from j
+                let mut source_idx = i & non_target_mask;
+                for (k, &pos) in target_bits.iter().enumerate() {
+                    if (j >> (g - 1 - k)) & 1 == 1 {
+                        source_idx |= 1 << pos;
+                    }
+                }
+
+                sum += gate_elem * state[source_idx];
+            }
+
+            sum
+        })
+        .collect();
+
+    new_state
+}
+
+fn matrix_to_2x2(matrix: &Matrix<Complex<f64>>) -> [[Complex<f64>; 2]; 2] {
+    [
+        [matrix.data[0], matrix.data[1]],
+        [matrix.data[2], matrix.data[3]],
+    ]
+}
+
+fn apply_kernel_direct(
+    state: &[Complex<f64>],
+    kernel: &Kernel,
+    num_qubits: usize,
+) -> Vec<Complex<f64>> {
+    let dim = 1 << num_qubits;
+    let g = kernel.targets.len();
+    let gate_dim = 1 << g;
+
+    let target_bits: Vec<usize> = kernel.targets.iter().map(|&t| num_qubits - 1 - t).collect();
+
+    let mut non_target_mask: usize = (1 << num_qubits) - 1;
+    for &pos in &target_bits {
+        non_target_mask &= !(1 << pos);
+    }
+
+    let mut new_state = vec![complex!(0.0, 0.0); dim];
+
+    for (i, new_val) in new_state.iter_mut().enumerate() {
+        let mut target_idx = 0usize;
+        for (k, &pos) in target_bits.iter().enumerate() {
+            if (i >> pos) & 1 == 1 {
+                target_idx |= 1 << (g - 1 - k);
+            }
+        }
+
+        let mut sum = complex!(0.0, 0.0);
+
+        for j in 0..gate_dim {
+            let gate_elem = kernel.matrix.data[target_idx * gate_dim + j];
+
+            if gate_elem.real.abs() < 1e-15 && gate_elem.imaginary.abs() < 1e-15 {
+                continue;
+            }
+
+            let mut source_idx = i & non_target_mask;
+            for (k, &pos) in target_bits.iter().enumerate() {
+                if (j >> (g - 1 - k)) & 1 == 1 {
+                    source_idx |= 1 << pos;
+                }
+            }
+
+            sum += gate_elem * state[source_idx];
+        }
+
+        *new_val = sum;
+    }
+
+    new_state
+}
diff --git a/src/lib.rs b/src/lib.rs
new file mode 100644
index 0000000..6f54357
--- /dev/null
+++ b/src/lib.rs
@@ -0,0 +1,21 @@
+pub mod core;
+pub mod maths;
+pub mod visualizer;
+
+pub use maths::complex::*;
+pub use maths::format::*;
+pub use maths::matrix::*;
+pub use maths::numeric::*;
+pub use maths::simd::*;
+pub use maths::vector::*;
+
+pub use core::circuit::*;
+pub use core::classical_components::*;
+pub use core::custom_gate::*;
+pub use core::gates;
+pub use core::kernel::*;
+pub use core::noise::*;
+pub use core::quantum_components::*;
+pub use core::runtime::*;
+
+pub use visualizer::*;
diff --git a/src/maths/complex.rs b/src/maths/complex.rs
new file mode 100644
index 0000000..31eae69
--- /dev/null
+++ b/src/maths/complex.rs
@@ -0,0 +1,180 @@
+use crate::Float;
+use core::{fmt, ops};
+
+#[macro_export]
+macro_rules! complex {
+    ($real:expr, $imaginary:expr) => {
+        $crate::Complex::new($real, $imaginary)
+    };
+}
+
+macro_rules! impl_ops {
+    ($trait:ident, $method:ident, $op:tt) => {
+        impl<T: Float> ops::$trait for Complex<T> {
+            type Output = Complex<T>;
+
+            fn $method(self, other: Complex<T>) -> Complex<T> {
+                Complex {
+                    real: self.real $op other.real,
+                    imaginary: self.imaginary $op other.imaginary,
+                }
+            }
+        }
+    };
+
+    ($trait:ident, $method:ident, $op:tt, real) => {
+        impl<T: Float> ops::$trait<T> for Complex<T> {
+            type Output = Complex<T>;
+
+            fn $method(self, other: T) -> Complex<T> {
+                Complex {
+                    real: self.real $op other,
+                    imaginary: self.imaginary,
+                }
+            }
+        }
+    };
+
+    ($trait_assign:ident, $method_assign:ident, $op:tt, assign) => {
+        impl<T: Float> ops::$trait_assign for Complex<T> {
+            fn $method_assign(&mut self, other: Complex<T>) {
+                self.real = self.real $op other.real;
+                self.imaginary = self.imaginary $op other.imaginary;
+            }
+        }
+    };
+
+    ($trait_assign:ident, $method_assign:ident, $op:tt, assign_real) => {
+        impl<T: Float> ops::$trait_assign<T> for Complex<T> {
+            fn $method_assign(&mut self, other: T) {
+                self.real = self.real $op other;
+            }
+        }
+    };
+}
+
+#[derive(Copy, Clone, PartialOrd, PartialEq)]
+pub struct Complex<T: Float> {
+    pub real: T,
+    pub imaginary: T,
+}
+
+impl<T: Float + fmt::Debug> fmt::Debug for Complex<T> {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        write!(
+            f,
+            "Complex {{ real: {:?}, imaginary: {:?} }}",
+            self.real, self.imaginary
+        )
+    }
+}
+
+impl<T: Float + fmt::Display> fmt::Display for Complex<T> {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        write!(f, "{} + {}i", self.real, self.imaginary)
+    }
+}
+
+impl<T: Float> ops::Neg for Complex<T> {
+    type Output = Complex<T>;
+
+    fn neg(self) -> Complex<T> {
+        Complex {
+            real: -self.real,
+            imaginary: -self.imaginary,
+        }
+    }
+}
+
+impl<T: Float> From<T> for Complex<T> {
+    fn from(real: T) -> Complex<T> {
+        Complex {
+            real,
+            imaginary: T::zero(),
+        }
+    }
+}
+
+impl<T: Float> Complex<T> {
+    pub fn new(real: T, imaginary: T) -> Complex<T> {
+        Complex { real, imaginary }
+    }
+
+    pub fn get_conjugate(&self) -> Complex<T> {
+        Complex {
+            real: self.real,
+            imaginary: -self.imaginary,
+        }
+    }
+
+    pub fn conjugate(&mut self) {
+        self.imaginary = -self.imaginary;
+    }
+
+    pub fn phase(&self) -> T {
+        T::atan2(self.imaginary, self.real)
+    }
+
+    pub fn norm2(&self) -> T {
+        self.real * self.real + self.imaginary * self.imaginary
+    }
+
+    pub fn abs(&self) -> T {
+        T::sqrt(self.norm2())
+    }
+}
+
+impl_ops!(Add, add, +);
+impl_ops!(Sub, sub, -);
+
+impl<T: Float> ops::Mul for Complex<T> {
+    type Output = Complex<T>;
+
+    fn mul(self, other: Complex<T>) -> Complex<T> {
+        // (a + bi) * (c + di) = (ac - bd) + (ad + bc)i
+        Complex {
+            real: self.real * other.real - self.imaginary * other.imaginary,
+            imaginary: self.real * other.imaginary + self.imaginary * other.real,
+        }
+    }
+}
+
+impl<T: Float> ops::Div for Complex<T> {
+    type Output = Complex<T>;
+
+    fn div(self, other: Complex<T>) -> Complex<T> {
+        // (a + bi) / (c + di) = ((ac + bd) + (bc - ad)i) / (c² + d²)
+        let denom = other.real * other.real + other.imaginary * other.imaginary;
+        Complex {
+            real: (self.real * other.real + self.imaginary * other.imaginary) / denom,
+            imaginary: (self.imaginary * other.real - self.real * other.imaginary) / denom,
+        }
+    }
+}
+
+impl_ops!(AddAssign, add_assign, +, assign);
+impl_ops!(SubAssign, sub_assign, -, assign);
+
+impl<T: Float> ops::MulAssign for Complex<T> {
+    fn mul_assign(&mut self, other: Complex<T>) {
+        let new_real = self.real * other.real - self.imaginary * other.imaginary;
+        let new_imag = self.real * other.imaginary + self.imaginary * other.real;
+        self.real = new_real;
+        self.imaginary = new_imag;
+    }
+}
+
+impl<T: Float> ops::DivAssign for Complex<T> {
+    fn div_assign(&mut self, other: Complex<T>) {
+        let denom = other.real * other.real + other.imaginary * other.imaginary;
+        let new_real = (self.real * other.real + self.imaginary * other.imaginary) / denom;
+        let new_imag = (self.imaginary * other.real - self.real * other.imaginary) / denom;
+        self.real = new_real;
+        self.imaginary = new_imag;
+    }
+}
+
+impl_ops!(Add, add, +, real);
+impl_ops!(Sub, sub, -, real);
+impl_ops!(Mul, mul, *, real);
+impl_ops!(Div, div, /, real);
diff --git a/src/maths/format.rs b/src/maths/format.rs
new file mode 100644
index 0000000..b6b0191
--- /dev/null
+++ b/src/maths/format.rs
@@ -0,0 +1,141 @@
+use crate::Complex;
+
+const EPSILON: f64 = 1e-10;
+const SQRT_2: f64 = std::f64::consts::SQRT_2;
+const INV_SQRT_2: f64 = 0.7071067811865475;
+const INV_SQRT_8: f64 = 0.3535533905932738;
+const INV_SQRT_32: f64 = 0.1767766952966369;
+
+fn approx_eq(a: f64, b: f64) -> bool {
+    (a - b).abs() < EPSILON
+}
+
+fn format_real_symbolic(v: f64) -> Option<String> {
+    let abs_v = v.abs();
+    let sign = if v < 0.0 { "-" } else { "" };
+
+    if approx_eq(abs_v, 0.0) {
+        return Some("0".to_string());
+    }
+    if approx_eq(abs_v, 1.0) {
+        return Some(format!("{}1", sign));
+    }
+    if approx_eq(abs_v, 0.5) {
+        return Some(format!("{}½", sign));
+    }
+    if approx_eq(abs_v, 0.25) {
+        return Some(format!("{}¼", sign));
+    }
+    if approx_eq(abs_v, 0.75) {
+        return Some(format!("{}¾", sign));
+    }
+    if approx_eq(abs_v, 0.125) {
+        return Some(format!("{}⅛", sign));
+    }
+    if approx_eq(abs_v, SQRT_2) {
+        return Some(format!("{}√2", sign));
+    }
+    if approx_eq(abs_v, INV_SQRT_2) {
+        return Some(format!("{}¹⁄√2", sign));
+    }
+    if approx_eq(abs_v, INV_SQRT_8) {
+        return Some(format!("{}¹⁄√8", sign));
+    }
+    if approx_eq(abs_v, INV_SQRT_32) {
+        return Some(format!("{}¹⁄√32", sign));
+    }
+    if approx_eq(abs_v, 2.0) {
+        return Some(format!("{}2", sign));
+    }
+    if approx_eq(abs_v, 1.0 / 3.0) {
+        return Some(format!("{}⅓", sign));
+    }
+    if approx_eq(abs_v, 2.0 / 3.0) {
+        return Some(format!("{}⅔", sign));
+    }
+
+    None
+}
+
+pub fn format_amplitude(c: &Complex<f64>) -> String {
+    let re = c.real;
+    let im = c.imaginary;
+
+    let re_zero = approx_eq(re.abs(), 0.0);
+    let im_zero = approx_eq(im.abs(), 0.0);
+
+    if re_zero && im_zero {
+        return "0".to_string();
+    }
+
+    if im_zero {
+        if let Some(s) = format_real_symbolic(re) {
+            return s;
+        }
+        return format!("{:.4}", re);
+    }
+
+    if re_zero {
+        if approx_eq(im.abs(), 1.0) {
+            return if im > 0.0 {
+                "i".to_string()
+            } else {
+                "-i".to_string()
+            };
+        }
+        if let Some(s) = format_real_symbolic(im) {
+            return format!("{}i", s);
+        }
+        return format!("{:.4}i", im);
+    }
+
+    let re_str = format_real_symbolic(re).unwrap_or_else(|| format!("{:.4}", re));
+    let im_str = if approx_eq(im.abs(), 1.0) {
+        if im > 0.0 {
+            "+i".to_string()
+        } else {
+            "-i".to_string()
+        }
+    } else {
+        let im_sym = format_real_symbolic(im.abs());
+        let sign = if im > 0.0 { "+" } else { "-" };
+        match im_sym {
+            Some(s) => format!("{}{}i", sign, s.trim_start_matches('-')),
+            None => format!("{}{:.4}i", sign, im.abs()),
+        }
+    };
+
+    format!("{}{}", re_str, im_str)
+}
+
+pub fn format_probability(p: f64) -> String {
+    if approx_eq(p, 0.0) {
+        return "0".to_string();
+    }
+    if approx_eq(p, 1.0) {
+        return "1".to_string();
+    }
+    if approx_eq(p, 0.5) {
+        return "½".to_string();
+    }
+    if approx_eq(p, 0.25) {
+        return "¼".to_string();
+    }
+    if approx_eq(p, 0.75) {
+        return "¾".to_string();
+    }
+    if approx_eq(p, 0.125) {
+        return "⅛".to_string();
+    }
+    if approx_eq(p, 0.0625) {
+        return "¹⁄₁₆".to_string();
+    }
+    if approx_eq(p, 1.0 / 3.0) {
+        return "⅓".to_string();
+    }
+    if approx_eq(p, 2.0 / 3.0) {
+        return "⅔".to_string();
+    }
+
+    format!("{:.4}", p)
+}
diff --git a/src/maths/matrix.rs b/src/maths/matrix.rs
new file mode 100644
index 0000000..c9d96a4
--- /dev/null
+++ b/src/maths/matrix.rs
@@ -0,0 +1,310 @@
+use super::Float;
+use core::{fmt, ops};
+
+#[macro_export]
+macro_rules! matrix {
+    ( $( $( $x:expr ),* );* ) => {{
+        let mut data = Vec::new();
+        let mut rows = 0;
+        let mut cols = 0;
+
+        $(
+            let row_data = $( $x )*;
+            if cols == 0 {
+                cols = row_data.len();
+            }
+            assert_eq!(cols, row_data.len(), "All rows must have the same number of columns.");
+            data.extend(row_data);
+            rows += 1;
+        )*
+
+        $crate::Matrix::new(rows, cols, data)
+    }};
+}
+
+macro_rules! impl_matrix_ops {
+    ($($trait:ident, $method:ident, $other:ty, $output:ty, $scale_fn:ident),* $(,)?) => {
+        $(
+            impl<T: Float> core::ops::$trait<$other> for Matrix<T> {
+                type Output = $output;
+
+                fn $method(self, other: $other) -> Self::Output {
+                    self.$scale_fn(other)
+                }
+            }
+        )*
+    };
+    ($($trait:ident, $method:ident, $other:ty, $scale_fn:ident),* $(,)?) => {
+        $(
+            impl<T: Float> core::ops::$trait<$other> for Matrix<T> {
+                fn $method(&mut self, other: $other) {
+                    *self = self.$scale_fn(other);
+                }
+            }
+        )*
+    };
+}
+
+#[derive(Clone)]
+pub struct Matrix<T: Float> {
+    pub data: Vec<T>,
+    pub rows: usize,
+    pub cols: usize,
+}
+
+impl<T: Float> Matrix<T> {
+    pub fn new(rows: usize, cols: usize, data: Vec<T>) -> Self {
+        Matrix { data, rows, cols }
+    }
+
+    pub fn get(&self, row: usize, col: usize) -> T {
+        self.data[row * self.cols + col]
+    }
+
+    pub fn set(&mut self, row: usize, col: usize, value: T) {
+        self.data[row * self.cols + col] = value;
+    }
+
+    pub fn dot(&self, other: &Self) -> Option<Matrix<T>> {
+        if self.cols != other.rows {
+            return None;
+        }
+
+        let mut result = Matrix::new(
+            self.rows,
+            other.cols,
+            vec![T::zero(); self.rows * other.cols],
+        );
+        for i in 0..self.rows {
+            for j in 0..other.cols {
+                let mut sum = T::zero();
+                for k in 0..self.cols {
+                    sum += self.get(i, k) * other.get(k, j) ;
+                }
+                result.set(i, j, sum);
+            }
+        }
+        Some(result)
+    }
+
+    pub fn kronecker(&self, other: &Self) -> Matrix<T> {
+        let new_rows = self.rows * other.rows;
+        let new_cols = self.cols * other.cols;
+
+        let mut result = Matrix::new(new_rows, new_cols, vec![T::zero(); new_rows * new_cols]);
+
+        for i in 0..self.rows {
+            for j in 0..self.cols {
+                let self_val = self.get(i, j);
+                for k in 0..other.rows {
+                    for l in 0..other.cols {
+                        let result_row = i * other.rows + k;
+                        let result_col = j * other.cols + l;
+                        result.set(result_row, result_col, self_val * other.get(k, l));
+                    }
+                }
+            }
+        }
+
+        result
+    }
+
+    pub fn transpose(&self) -> Matrix<T> {
+        let mut result = Matrix::new(self.cols, self.rows, vec![T::zero(); self.cols * self.rows]);
+
+        for i in 0..self.rows {
+            for j in 0..self.cols {
+                let value = self.get(i, j);
+                result.set(j, i, value);
+            }
+        }
+
+        result
+    }
+
+    pub fn add_to(&self, other: &Self) -> Option<Matrix<T>> {
+        if self.rows != other.rows || self.cols != other.cols {
+            return None;
+        }
+
+        let mut result = Matrix::new(self.rows, self.cols, vec![T::zero(); self.rows * self.cols]);
+
+        for i in 0..self.rows {
+            for j in 0..self.cols {
+                let sum = self.get(i, j) + other.get(i, j);
+                result.set(i, j, sum);
+            }
+        }
+        Some(result)
+    }
+
+    pub fn subtract(&self, other: &Self) -> Option<Matrix<T>> {
+        if self.rows != other.rows || self.cols != other.cols {
+            return None;
+        }
+
+        let mut result = Matrix::new(self.rows, self.cols, vec![T::zero(); self.rows * self.cols]);
+
+        for i in 0..self.rows {
+            for j in 0..self.cols {
+                let diff = self.get(i, j) - other.get(i, j);
+                result.set(i, j, diff);
+            }
+        }
+        Some(result)
+    }
+
+    pub fn scale(&self, scalar: T) -> Matrix<T> {
+        let mut result = Matrix::new(self.rows, self.cols, vec![T::zero(); self.rows * self.cols]);
+
+        for i in 0..self.rows {
+            for j in 0..self.cols {
+                let scaled_value = self.get(i, j) * scalar;
+                result.set(i, j, scaled_value);
+            }
+        }
+        result
+    }
+}
+
+impl<T: Float> ops::Index<(usize, usize)> for Matrix<T> {
+    type Output = T;
+
+    fn index(&self, index: (usize, usize)) -> &Self::Output {
+        &self.data[index.0 * self.cols + index.1]
+    }
+}
+
+impl<T: Float> ops::IndexMut<(usize, usize)> for Matrix<T> {
+    fn index_mut(&mut self, index: (usize, usize)) -> &mut Self::Output {
+        &mut self.data[index.0 * self.cols + index.1]
+    }
+}
+
+impl<T: Float> ops::AddAssign<&Matrix<T>> for Matrix<T> {
+    fn add_assign(&mut self, other: &Matrix<T>) {
+        if let Some(result) = self.add_to(other) {
+            *self = result;
+        }
+    }
+}
+
+impl<T: Float> ops::SubAssign<&Matrix<T>> for Matrix<T> {
+    fn sub_assign(&mut self, other: &Matrix<T>) {
+        if let Some(result) = self.subtract(other) {
+            *self = result;
+        }
+    }
+}
+
+impl_matrix_ops! {
+    Add, add, &Matrix<T>, Option<Matrix<T>>, add_to,
+    Sub, sub, &Matrix<T>, Option<Matrix<T>>, subtract,
+    Mul, mul, T, Matrix<T>, scale,
+    Div, div, T, Matrix<T>, scale,
+}
+
+impl_matrix_ops! {
+    MulAssign, mul_assign, T, scale,
+    DivAssign, div_assign, T, scale,
+}
+
+impl<T: Float + fmt::Debug> fmt::Debug for Matrix<T> {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        for i in 0..self.rows {
+            for j in 0..self.cols {
+                write!(f, "{:?} ", self.get(i, j))?;
+            }
+            writeln!(f)?;
+        }
+        Ok(())
+    }
+}
+
+impl<T: Float + fmt::Display> fmt::Display for Matrix<T> {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        let elements: Vec<String> = self.data.iter().map(ToString::to_string).collect();
+        let is_complex = elements.iter().any(|element| element.contains("i"));
+
+        let normalized: Vec<(f64, f64)> = self
+            .data
+            .iter()
+            .map(|element| {
+                let element_string = element.to_string();
+
+                if is_complex {
+                    let element_string = element_string.trim_end_matches('i').trim();
+                    let element_split: Vec<&str> = element_string.split_whitespace().collect();
+                    let real = element_split[0].parse::<f64>().unwrap();
+                    let imaginary = element_split
+                        .get(2)
+                        .map_or(0.0, |&s| s.parse::<f64>().unwrap());
+                    (real, imaginary)
+                } else {
+                    (element_string.parse::<f64>().unwrap(), 0.0)
+                }
+            })
+            .collect();
+
+        let max_widths = normalized
+            .iter()
+            .fold((0, 0), |(max_0, max_1), &(real, imag)| {
+                let new_max_0 = max_0.max(format!("{:.2}", real).len());
+                let new_max_1 = if is_complex {
+                    max_1.max(format!("{:.2}", imag.abs()).len())
+                } else {
+                    max_1
+                };
+                (new_max_0, new_max_1)
+            });
+
+        let aligned: Vec<String> = normalized
+            .iter()
+            .map(|&(real, imag)| {
+                if is_complex {
+                    format!(
+                        "{:>rewidth$.2} {} {:>imwidth$.2}i",
+                        real,
+                        if imag > 0.0 { "+" } else { "-" },
+                        imag.abs(),
+                        rewidth = max_widths.0,
+                        imwidth = max_widths.1,
+                    )
+                } else {
+                    format!("{:>width$.2}", real, width = max_widths.0)
+                }
+            })
+            .collect();
+
+        for i in 0..self.rows {
+            if i == 0 {
+                write!(f, "┌")?;
+            } else if i == self.rows - 1 {
+                write!(f, "└")?;
+            } else {
+                write!(f, "│")?;
+            }
+
+            for j in 0..self.cols {
+                write!(f, "{}", aligned[i + j * self.rows])?;
+                if j != self.cols - 1 {
+                    write!(f, ", ")?;
+                }
+            }
+
+            if i == 0 {
+                write!(f, "┐")?;
+            } else if i == self.rows - 1 {
+                write!(f, "┘")?;
+            } else {
+                write!(f, "│")?;
+            }
+
+            if i != self.rows - 1 {
+                writeln!(f)?;
+            }
+        }
+
+        Ok(())
+    }
+}
diff --git a/src/maths/mod.rs b/src/maths/mod.rs
new file mode 100644
index 0000000..85f4872
--- /dev/null
+++ b/src/maths/mod.rs
@@ -0,0 +1,14 @@
+pub mod complex;
+pub mod format;
+pub mod matrix;
+pub mod numeric;
+pub mod simd;
+pub mod vector;
+pub mod vector_ops;
+
+pub use complex::*;
+pub use format::*;
+pub use matrix::*;
+pub use numeric::*;
+pub use simd::*;
+pub use vector::*;
diff --git a/src/maths/numeric.rs b/src/maths/numeric.rs
new file mode 100644
index 0000000..f5a649e
--- /dev/null
+++ b/src/maths/numeric.rs
@@ -0,0 +1,101 @@
+use crate::Complex;
+use core::ops;
+
+macro_rules! impl_numeric {
+    ($($t:ty),*) => {
+        $(
+            impl Numeric for $t {
+                fn zero() -> Self {
+                    0 as $t
+                }
+
+                fn one() -> Self {
+                    1 as $t
+                }
+            }
+        )*
+    };
+}
+
+macro_rules! impl_cnumeric {
+    ($($t:ty),*) => {
+        $(impl Numeric for Complex<$t> {
+            fn zero() -> Self { Complex::new(0.0, 0.0) }
+            fn one() -> Self { Complex::new(1.0, 0.0) }
+        })*
+    };
+}
+
+macro_rules! impl_float {
+    ($($t:ty, $sqrt_fn:path, $atan2_fn:path),*) => {
+        $(
+            impl Float for $t {
+                fn sqrt(self) -> Self {
+                    $sqrt_fn(self)
+                }
+
+                fn atan2(y: Self, x: Self) -> Self {
+                    $atan2_fn(y, x)
+                }
+            }
+        )*
+    };
+}
+
+macro_rules! impl_cfloat {
+    ($($t:ty, $sqrt_fn:path, $atan2_fn:path, $cos_fn:path, $sin_fn:path),*) => {
+        $(
+            impl Float for Complex<$t> {
+                fn sqrt(self) -> Self {
+                    let r = self.abs();
+                    let theta = self.phase();
+
+                    let sqrt_r = $sqrt_fn(r);
+                    let sqrt_theta = theta / 2.0;
+
+                    Complex::new(
+                        sqrt_r * $cos_fn(sqrt_theta),
+                        sqrt_r * $sin_fn(sqrt_theta),
+                    )
+                }
+
+                fn atan2(y: Self, x: Self) -> Self {
+                    Complex::new(
+                        $atan2_fn(y.real, x.real),
+                        $atan2_fn(y.imaginary, x.imaginary),
+                    )
+                }
+            }
+        )*
+    };
+}
+
+pub trait Numeric:
+    Copy
+    + PartialOrd
+    + ops::Add<Output = Self>
+    + ops::Mul<Output = Self>
+    + ops::Sub<Output = Self>
+    + ops::Div<Output = Self>
+    + ops::Neg<Output = Self>
+    + ops::AddAssign
+    + ops::SubAssign
+    + ops::MulAssign
+    + ops::DivAssign
+{
+    fn zero() -> Self;
+    fn one() -> Self;
+}
+
+impl_numeric!(i32, i64, f32, f64);
+impl_cnumeric!(f32, f64);
+impl_float!(f32, libm::sqrtf, libm::atan2f);
+impl_float!(f64, libm::sqrt, libm::atan2);
+impl_cfloat!(f32, libm::sqrtf, libm::atan2f, libm::cosf, libm::sinf);
+impl_cfloat!(f64, libm::sqrt, libm::atan2, libm::cos, libm::sin);
+
+pub trait Integer: Numeric {}
+pub trait Float: Numeric {
+    fn sqrt(self) -> Self;
+    fn atan2(y: Self, x: Self) -> Self;
+}
diff --git a/src/maths/simd.rs b/src/maths/simd.rs
new file mode 100644
index 0000000..de0370c
--- /dev/null
+++ b/src/maths/simd.rs
@@ -0,0 +1,510 @@
+use crate::{complex, Complex};
+
+#[cfg(target_arch = "x86_64")]
+use std::arch::x86_64::*;
+
+#[cfg(target_arch = "aarch64")]
+use std::arch::aarch64::*;
+
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub enum SimdCapability {
+    None,
+    #[cfg(any(target_arch = "x86_64", target_arch = "x86"))]
+    Avx2,
+    #[cfg(any(target_arch = "x86_64", target_arch = "x86"))]
+    Avx512,
+    #[cfg(target_arch = "aarch64")]
+    Neon,
+}
+
+impl SimdCapability {
+    pub fn detect() -> Self {
+        #[cfg(any(target_arch = "x86_64", target_arch = "x86"))]
+        {
+            if is_x86_feature_detected!("avx512f") && is_x86_feature_detected!("avx512dq") {
+                return SimdCapability::Avx512;
+            }
+            if is_x86_feature_detected!("avx2") && is_x86_feature_detected!("fma") {
+                return SimdCapability::Avx2;
+            }
+        }
+
+        #[cfg(target_arch = "aarch64")]
+        {
+            return SimdCapability::Neon;
+        }
+
+        #[allow(unreachable_code)]
+        SimdCapability::None
+    }
+
+    pub fn name(&self) -> &'static str {
+        match self {
+            SimdCapability::None => "Scalar",
+            #[cfg(any(target_arch = "x86_64", target_arch = "x86"))]
+            SimdCapability::Avx2 => "AVX2+FMA",
+            #[cfg(any(target_arch = "x86_64", target_arch = "x86"))]
+            SimdCapability::Avx512 => "AVX-512",
+            #[cfg(target_arch = "aarch64")]
+            SimdCapability::Neon => "NEON",
+        }
+    }
+}
+
+pub fn apply_single_qubit_gate_simd(
+    state: &mut [Complex<f64>],
+    gate: &[[Complex<f64>; 2]; 2],
+    target: usize,
+    num_qubits: usize,
+) {
+    let capability = SimdCapability::detect();
+
+    match capability {
+        #[cfg(target_arch = "x86_64")]
+        SimdCapability::Avx2 => unsafe {
+            apply_single_qubit_avx2(state, gate, target, num_qubits);
+        },
+        #[cfg(target_arch = "x86_64")]
+        SimdCapability::Avx512 => unsafe {
+            apply_single_qubit_avx512(state, gate, target, num_qubits);
+        },
+        #[cfg(target_arch = "aarch64")]
+        SimdCapability::Neon => unsafe {
+            apply_single_qubit_neon(state, gate, target, num_qubits);
+        },
+        _ => {
+            apply_single_qubit_scalar(state, gate, target, num_qubits);
+        }
+    }
+}
+
+#[cfg(target_arch = "x86_64")]
+#[target_feature(enable = "avx2", enable = "fma")]
+unsafe fn apply_single_qubit_avx2(
+    state: &mut [Complex<f64>],
+    gate: &[[Complex<f64>; 2]; 2],
+    target: usize,
+    num_qubits: usize,
+) {
+    let target_bit = num_qubits - 1 - target;
+    let step = 1 << target_bit;
+    let dim = 1 << num_qubits;
+
+    let g00 = gate[0][0];
+    let g01 = gate[0][1];
+    let g10 = gate[1][0];
+    let g11 = gate[1][1];
+
+    let pairs: Vec<(usize, usize)> = (0..dim)
+        .filter(|&i| (i >> target_bit) & 1 == 0)
+        .map(|i| (i, i | step))
+        .collect();
+
+    let chunks = pairs.len() / 2;
+
+    for chunk_idx in 0..chunks {
+        let (i0, j0) = pairs[chunk_idx * 2];
+        let (i1, j1) = pairs[chunk_idx * 2 + 1];
+
+        let s0_re = _mm256_set_pd(
+            state[j1].real,
+            state[i1].real,
+            state[j0].real,
+            state[i0].real,
+        );
+        let s0_im = _mm256_set_pd(
+            state[j1].imaginary,
+            state[i1].imaginary,
+            state[j0].imaginary,
+            state[i0].imaginary,
+        );
+
+        let g_re_0 = _mm256_set_pd(g01.real, g00.real, g01.real, g00.real);
+        let g_im_0 = _mm256_set_pd(g01.imaginary, g00.imaginary, g01.imaginary, g00.imaginary);
+        let g_re_1 = _mm256_set_pd(g11.real, g10.real, g11.real, g10.real);
+        let g_im_1 = _mm256_set_pd(g11.imaginary, g10.imaginary, g11.imaginary, g10.imaginary);
+
+        let prod0_re = _mm256_fmsub_pd(s0_re, g_re_0, _mm256_mul_pd(s0_im, g_im_0));
+        let prod0_im = _mm256_fmadd_pd(s0_re, g_im_0, _mm256_mul_pd(s0_im, g_re_0));
+
+        let prod1_re = _mm256_fmsub_pd(s0_re, g_re_1, _mm256_mul_pd(s0_im, g_im_1));
+        let prod1_im = _mm256_fmadd_pd(s0_re, g_im_1, _mm256_mul_pd(s0_im, g_re_1));
+
+        let mut res0_re = [0.0f64; 4];
+        let mut res0_im = [0.0f64; 4];
+        let mut res1_re = [0.0f64; 4];
+        let mut res1_im = [0.0f64; 4];
+
+        _mm256_storeu_pd(res0_re.as_mut_ptr(), prod0_re);
+        _mm256_storeu_pd(res0_im.as_mut_ptr(), prod0_im);
+        _mm256_storeu_pd(res1_re.as_mut_ptr(), prod1_re);
+        _mm256_storeu_pd(res1_im.as_mut_ptr(), prod1_im);
+
+        state[i0] = complex!(res0_re[0] + res0_re[1], res0_im[0] + res0_im[1]);
+        state[j0] = complex!(res1_re[0] + res1_re[1], res1_im[0] + res1_im[1]);
+        state[i1] = complex!(res0_re[2] + res0_re[3], res0_im[2] + res0_im[3]);
+        state[j1] = complex!(res1_re[2] + res1_re[3], res1_im[2] + res1_im[3]);
+    }
+
+    for &(i, j) in pairs.iter().skip(chunks * 2) {
+        let s0 = state[i];
+        let s1 = state[j];
+
+        let new0 = complex!(
+            s0.real * g00.real - s0.imaginary * g00.imaginary + s1.real * g01.real
+                - s1.imaginary * g01.imaginary,
+            s0.real * g00.imaginary
+                + s0.imaginary * g00.real
+                + s1.real * g01.imaginary
+                + s1.imaginary * g01.real
+        );
+
+        let new1 = complex!(
+            s0.real * g10.real - s0.imaginary * g10.imaginary + s1.real * g11.real
+                - s1.imaginary * g11.imaginary,
+            s0.real * g10.imaginary
+                + s0.imaginary * g10.real
+                + s1.real * g11.imaginary
+                + s1.imaginary * g11.real
+        );
+
+        state[i] = new0;
+        state[j] = new1;
+    }
+}
+
+#[cfg(target_arch = "x86_64")]
+#[target_feature(enable = "avx512f", enable = "avx512dq")]
+unsafe fn apply_single_qubit_avx512(
+    state: &mut [Complex<f64>],
+    gate: &[[Complex<f64>; 2]; 2],
+    target: usize,
+    num_qubits: usize,
+) {
+    let target_bit = num_qubits - 1 - target;
+    let step = 1 << target_bit;
+    let dim = 1 << num_qubits;
+
+    let g00 = gate[0][0];
+    let g01 = gate[0][1];
+    let g10 = gate[1][0];
+    let g11 = gate[1][1];
+
+    let pairs: Vec<(usize, usize)> = (0..dim)
+        .filter(|&i| (i >> target_bit) & 1 == 0)
+        .map(|i| (i, i | step))
+        .collect();
+
+    let chunks = pairs.len() / 4;
+
+    for chunk_idx in 0..chunks {
+        let base = chunk_idx * 4;
+        let (i0, j0) = pairs[base];
+        let (i1, j1) = pairs[base + 1];
+        let (i2, j2) = pairs[base + 2];
+        let (i3, j3) = pairs[base + 3];
+
+        let s0_re = _mm512_set_pd(
+            state[j3].real,
+            state[i3].real,
+            state[j2].real,
+            state[i2].real,
+            state[j1].real,
+            state[i1].real,
+            state[j0].real,
+            state[i0].real,
+        );
+        let s0_im = _mm512_set_pd(
+            state[j3].imaginary,
+            state[i3].imaginary,
+            state[j2].imaginary,
+            state[i2].imaginary,
+            state[j1].imaginary,
+            state[i1].imaginary,
+            state[j0].imaginary,
+            state[i0].imaginary,
+        );
+
+        let g_re_0 = _mm512_set_pd(
+            g01.real, g00.real, g01.real, g00.real, g01.real, g00.real, g01.real, g00.real,
+        );
+        let g_im_0 = _mm512_set_pd(
+            g01.imaginary,
+            g00.imaginary,
+            g01.imaginary,
+            g00.imaginary,
+            g01.imaginary,
+            g00.imaginary,
+            g01.imaginary,
+            g00.imaginary,
+        );
+        let g_re_1 = _mm512_set_pd(
+            g11.real, g10.real, g11.real, g10.real, g11.real, g10.real, g11.real, g10.real,
+        );
+        let g_im_1 = _mm512_set_pd(
+            g11.imaginary,
+            g10.imaginary,
+            g11.imaginary,
+            g10.imaginary,
+            g11.imaginary,
+            g10.imaginary,
+            g11.imaginary,
+            g10.imaginary,
+        );
+
+        let prod0_re = _mm512_fmsub_pd(s0_re, g_re_0, _mm512_mul_pd(s0_im, g_im_0));
+        let prod0_im = _mm512_fmadd_pd(s0_re, g_im_0, _mm512_mul_pd(s0_im, g_re_0));
+        let prod1_re = _mm512_fmsub_pd(s0_re, g_re_1, _mm512_mul_pd(s0_im, g_im_1));
+        let prod1_im = _mm512_fmadd_pd(s0_re, g_im_1, _mm512_mul_pd(s0_im, g_re_1));
+
+        let mut res0_re = [0.0f64; 8];
+        let mut res0_im = [0.0f64; 8];
+        let mut res1_re = [0.0f64; 8];
+        let mut res1_im = [0.0f64; 8];
+
+        _mm512_storeu_pd(res0_re.as_mut_ptr(), prod0_re);
+        _mm512_storeu_pd(res0_im.as_mut_ptr(), prod0_im);
+        _mm512_storeu_pd(res1_re.as_mut_ptr(), prod1_re);
+        _mm512_storeu_pd(res1_im.as_mut_ptr(), prod1_im);
+
+        state[i0] = complex!(res0_re[0] + res0_re[1], res0_im[0] + res0_im[1]);
+        state[j0] = complex!(res1_re[0] + res1_re[1], res1_im[0] + res1_im[1]);
+        state[i1] = complex!(res0_re[2] + res0_re[3], res0_im[2] + res0_im[3]);
+        state[j1] = complex!(res1_re[2] + res1_re[3], res1_im[2] + res1_im[3]);
+        state[i2] = complex!(res0_re[4] + res0_re[5], res0_im[4] + res0_im[5]);
+        state[j2] = complex!(res1_re[4] + res1_re[5], res1_im[4] + res1_im[5]);
+        state[i3] = complex!(res0_re[6] + res0_re[7], res0_im[6] + res0_im[7]);
+        state[j3] = complex!(res1_re[6] + res1_re[7], res1_im[6] + res1_im[7]);
+    }
+
+    for &(i, j) in pairs.iter().skip(chunks * 4) {
+        let s0 = state[i];
+        let s1 = state[j];
+
+        let new0 = complex!(
+            s0.real * g00.real - s0.imaginary * g00.imaginary + s1.real * g01.real
+                - s1.imaginary * g01.imaginary,
+            s0.real * g00.imaginary
+                + s0.imaginary * g00.real
+                + s1.real * g01.imaginary
+                + s1.imaginary * g01.real
+        );
+
+        let new1 = complex!(
+            s0.real * g10.real - s0.imaginary * g10.imaginary + s1.real * g11.real
+                - s1.imaginary * g11.imaginary,
+            s0.real * g10.imaginary
+                + s0.imaginary * g10.real
+                + s1.real * g11.imaginary
+                + s1.imaginary * g11.real
+        );
+
+        state[i] = new0;
+        state[j] = new1;
+    }
+}
+
+#[cfg(target_arch = "aarch64")]
+unsafe fn apply_single_qubit_neon(
+    state: &mut [Complex<f64>],
+    gate: &[[Complex<f64>; 2]; 2],
+    target: usize,
+    num_qubits: usize,
+) {
+    let target_bit = num_qubits - 1 - target;
+    let step = 1 << target_bit;
+    let dim = 1 << num_qubits;
+
+    let g00 = gate[0][0];
+    let g01 = gate[0][1];
+    let g10 = gate[1][0];
+    let g11 = gate[1][1];
+
+    let pairs: Vec<(usize, usize)> = (0..dim)
+        .filter(|&i| (i >> target_bit) & 1 == 0)
+        .map(|i| (i, i | step))
+        .collect();
+
+    let chunks = pairs.len() / 2;
+
+    for chunk_idx in 0..chunks {
+        let (i0, j0) = pairs[chunk_idx * 2];
+        let (i1, j1) = pairs[chunk_idx * 2 + 1];
+
+        let s0_0 = state[i0];
+        let s1_0 = state[j0];
+        let s0_1 = state[i1];
+        let s1_1 = state[j1];
+
+        let s0_re = vld1q_f64([s0_0.real, s0_1.real].as_ptr());
+        let s0_im = vld1q_f64([s0_0.imaginary, s0_1.imaginary].as_ptr());
+        let s1_re = vld1q_f64([s1_0.real, s1_1.real].as_ptr());
+        let s1_im = vld1q_f64([s1_0.imaginary, s1_1.imaginary].as_ptr());
+
+        let g00_re = vdupq_n_f64(g00.real);
+        let g00_im = vdupq_n_f64(g00.imaginary);
+        let g01_re = vdupq_n_f64(g01.real);
+        let g01_im = vdupq_n_f64(g01.imaginary);
+        let g10_re = vdupq_n_f64(g10.real);
+        let g10_im = vdupq_n_f64(g10.imaginary);
+        let g11_re = vdupq_n_f64(g11.real);
+        let g11_im = vdupq_n_f64(g11.imaginary);
+
+        let new0_re = vaddq_f64(
+            vfmsq_f64(vmulq_f64(s0_re, g00_re), s0_im, g00_im),
+            vfmsq_f64(vmulq_f64(s1_re, g01_re), s1_im, g01_im),
+        );
+        let new0_im = vaddq_f64(
+            vfmaq_f64(vmulq_f64(s0_re, g00_im), s0_im, g00_re),
+            vfmaq_f64(vmulq_f64(s1_re, g01_im), s1_im, g01_re),
+        );
+
+        let new1_re = vaddq_f64(
+            vfmsq_f64(vmulq_f64(s0_re, g10_re), s0_im, g10_im),
+            vfmsq_f64(vmulq_f64(s1_re, g11_re), s1_im, g11_im),
+        );
+        let new1_im = vaddq_f64(
+            vfmaq_f64(vmulq_f64(s0_re, g10_im), s0_im, g10_re),
+            vfmaq_f64(vmulq_f64(s1_re, g11_im), s1_im, g11_re),
+        );
+
+        state[i0] = complex!(vgetq_lane_f64(new0_re, 0), vgetq_lane_f64(new0_im, 0));
+        state[j0] = complex!(vgetq_lane_f64(new1_re, 0), vgetq_lane_f64(new1_im, 0));
+        state[i1] = complex!(vgetq_lane_f64(new0_re, 1), vgetq_lane_f64(new0_im, 1));
+        state[j1] = complex!(vgetq_lane_f64(new1_re, 1), vgetq_lane_f64(new1_im, 1));
+    }
+
+    for &(i, j) in pairs.iter().skip(chunks * 2) {
+        let s0 = state[i];
+        let s1 = state[j];
+
+        let new0 = complex!(
+            s0.real * g00.real - s0.imaginary * g00.imaginary + s1.real * g01.real
+                - s1.imaginary * g01.imaginary,
+            s0.real * g00.imaginary
+                + s0.imaginary * g00.real
+                + s1.real * g01.imaginary
+                + s1.imaginary * g01.real
+        );
+
+        let new1 = complex!(
+            s0.real * g10.real - s0.imaginary * g10.imaginary + s1.real * g11.real
+                - s1.imaginary * g11.imaginary,
+            s0.real * g10.imaginary
+                + s0.imaginary * g10.real
+                + s1.real * g11.imaginary
+                + s1.imaginary * g11.real
+        );
+
+        state[i] = new0;
+        state[j] = new1;
+    }
+}
+
+fn apply_single_qubit_scalar(
+    state: &mut [Complex<f64>],
+    gate: &[[Complex<f64>; 2]; 2],
+    target: usize,
+    num_qubits: usize,
+) {
+    let target_bit = num_qubits - 1 - target;
+    let step = 1 << target_bit;
+    let dim = 1 << num_qubits;
+
+    let g00 = gate[0][0];
+    let g01 = gate[0][1];
+    let g10 = gate[1][0];
+    let g11 = gate[1][1];
+
+    for i in 0..dim {
+        if (i >> target_bit) & 1 == 1 {
+            continue;
+        }
+
+        let j = i | step;
+        let s0 = state[i];
+        let s1 = state[j];
+
+        let new0 = complex!(
+            s0.real * g00.real - s0.imaginary * g00.imaginary + s1.real * g01.real
+                - s1.imaginary * g01.imaginary,
+            s0.real * g00.imaginary
+                + s0.imaginary * g00.real
+                + s1.real * g01.imaginary
+                + s1.imaginary * g01.real
+        );
+
+        let new1 = complex!(
+            s0.real * g10.real - s0.imaginary * g10.imaginary + s1.real * g11.real
+                - s1.imaginary * g11.imaginary,
+            s0.real * g10.imaginary
+                + s0.imaginary * g10.real
+                + s1.real * g11.imaginary
+                + s1.imaginary * g11.real
+        );
+
+        state[i] = new0;
+        state[j] = new1;
+    }
+}
+
+pub fn apply_single_qubit_gate_simd_parallel(
+    state: &mut [Complex<f64>],
+    gate: &[[Complex<f64>; 2]; 2],
+    target: usize,
+    num_qubits: usize,
+) {
+    use rayon::prelude::*;
+
+    let target_bit = num_qubits - 1 - target;
+    let step = 1 << target_bit;
+    let dim = 1 << num_qubits;
+
+    let g00 = gate[0][0];
+    let g01 = gate[0][1];
+    let g10 = gate[1][0];
+    let g11 = gate[1][1];
+
+    let pairs: Vec<(usize, usize)> = (0..dim)
+        .filter(|&i| (i >> target_bit) & 1 == 0)
+        .map(|i| (i, i | step))
+        .collect();
+
+    let results: Vec<(usize, usize, Complex<f64>, Complex<f64>)> = pairs
+        .par_iter()
+        .map(|&(i, j)| {
+            let s0 = state[i];
+            let s1 = state[j];
+
+            let new0 = complex!(
+                s0.real * g00.real - s0.imaginary * g00.imaginary + s1.real * g01.real
+                    - s1.imaginary * g01.imaginary,
+                s0.real * g00.imaginary
+                    + s0.imaginary * g00.real
+                    + s1.real * g01.imaginary
+                    + s1.imaginary * g01.real
+            );
+
+            let new1 = complex!(
+                s0.real * g10.real - s0.imaginary * g10.imaginary + s1.real * g11.real
+                    - s1.imaginary * g11.imaginary,
+                s0.real * g10.imaginary
+                    + s0.imaginary * g10.real
+                    + s1.real * g11.imaginary
+                    + s1.imaginary * g11.real
+            );
+
+            (i, j, new0, new1)
+        })
+        .collect();
+
+    for (i, j, new0, new1) in results {
+        state[i] = new0;
+        state[j] = new1;
+    }
+}
+
+pub fn get_simd_info() -> String {
+    let cap = SimdCapability::detect();
+    format!("SIMD: {}", cap.name())
+}
diff --git a/src/maths/vector.rs b/src/maths/vector.rs
new file mode 100644
index 0000000..11f3d29
--- /dev/null
+++ b/src/maths/vector.rs
@@ -0,0 +1,258 @@
+use super::{Float, Matrix};
+use core::{fmt, ops};
+
+#[macro_export]
+macro_rules! row_vector {
+    ($($x:expr),*) => {
+        RowVector::new(vec![$($x),*])
+    };
+    ($($x:expr,)*) => {
+        RowVector::new(vec![$($x),*])
+    };
+}
+
+#[macro_export]
+macro_rules! column_vector {
+    ($($x:expr),*) => {
+        ColumnVector::new(vec![$($x),*])
+    };
+    ($($x:expr,)*) => {
+        ColumnVector::new(vec![$($x),*])
+    };
+}
+
+pub trait Vector<T: Float> {
+    fn new(data: Vec<T>) -> Self;
+    fn get(&self, index: usize) -> T;
+    fn set(&mut self, index: usize, value: T);
+    fn size(&self) -> usize;
+
+    fn dot(&self, other: &Self) -> T;
+    fn norm(&self) -> T;
+
+    fn max(&self) -> T;
+    fn min(&self) -> T;
+    fn sum(&self) -> T;
+
+    fn from_matrix(matrix: &Matrix<T>) -> Self;
+}
+
+pub trait VectorMatrix<T: Float> {
+    fn to_matrix(&self) -> Matrix<T>;
+}
+
+#[derive(Clone)]
+pub struct VectorImpl<T: Float, const ROWS: usize, const COLS: usize>(Vec<T>);
+pub type RowVector<T> = VectorImpl<T, 1, 0>;
+pub type ColumnVector<T> = VectorImpl<T, 0, 1>;
+
+impl<T: Float> ColumnVector<T> {
+    pub fn mul_matrix(&self, matrix: &Matrix<T>) -> Option<ColumnVector<T>> {
+        if matrix.cols != self.size() {
+            return None;
+        }
+
+        let mut result = ColumnVector::new(vec![T::zero(); matrix.rows]);
+
+        for i in 0..matrix.rows {
+            let mut sum = T::zero();
+            for j in 0..matrix.cols {
+                sum += matrix.get(i, j) * self.get(j) ;
+            }
+            result.set(i, sum);
+        }
+
+        Some(result)
+    }
+
+    pub fn transpose(&self) -> RowVector<T> {
+        RowVector::new(self.0.clone())
+    }
+}
+
+impl<T: Float> RowVector<T> {
+    pub fn mul_matrix(&self, matrix: &Matrix<T>) -> Option<RowVector<T>> {
+        if self.size() != matrix.rows {
+            return None;
+        }
+
+        let mut result = RowVector::new(vec![T::zero(); matrix.cols]);
+
+        for j in 0..matrix.cols {
+            let mut sum = T::zero();
+            for i in 0..matrix.rows {
+                sum += self.get(i) * matrix.get(i, j) ;
+            }
+            result.set(j, sum);
+        }
+
+        Some(result)
+    }
+
+    pub fn transpose(&self) -> ColumnVector<T> {
+        ColumnVector::new(self.0.clone())
+    }
+}
+
+impl<T: Float> VectorMatrix<T> for RowVector<T> {
+    fn to_matrix(&self) -> Matrix<T> {
+        Matrix::new(1, self.size(), self.0.clone())
+    }
+}
+
+impl<T: Float> VectorMatrix<T> for ColumnVector<T> {
+    fn to_matrix(&self) -> Matrix<T> {
+        Matrix::new(self.size(), 1, self.0.clone())
+    }
+}
+
+impl<T: Float, const ROWS: usize, const COLS: usize> Vector<T> for VectorImpl<T, ROWS, COLS> {
+    fn from_matrix(matrix: &Matrix<T>) -> Self {
+        Self::new(matrix.data.clone())
+    }
+
+    fn new(data: Vec<T>) -> Self {
+        Self(data)
+    }
+
+    fn get(&self, index: usize) -> T {
+        self.0[index]
+    }
+
+    fn set(&mut self, index: usize, value: T) {
+        self.0[index] = value;
+    }
+
+    fn size(&self) -> usize {
+        self.0.len()
+    }
+
+    fn dot(&self, other: &Self) -> T {
+        self.0
+            .iter()
+            .zip(other.0.iter())
+            .map(|(a, b)| *a * *b)
+            .fold(T::zero(), |acc, x| acc + x)
+    }
+
+    fn norm(&self) -> T {
+        self.0
+            .iter()
+            .map(|x| *x * *x)
+            .fold(T::zero(), |acc, x| acc + x)
+            .sqrt()
+    }
+
+    fn max(&self) -> T {
+        *self
+            .0
+            .iter()
+            .max_by(|a, b| a.partial_cmp(b).unwrap())
+            .unwrap_or(&T::zero())
+    }
+
+    fn min(&self) -> T {
+        *self
+            .0
+            .iter()
+            .min_by(|a, b| a.partial_cmp(b).unwrap())
+            .unwrap_or(&T::zero())
+    }
+
+    fn sum(&self) -> T {
+        self.0.iter().fold(T::zero(), |acc, x| acc + *x)
+    }
+}
+
+impl<T: Float, const ROWS: usize, const COLS: usize> VectorImpl<T, ROWS, COLS> {
+    pub fn add_to(&self, other: &Self) -> Option<VectorImpl<T, ROWS, COLS>> {
+        if self.size() != other.size() {
+            return None;
+        }
+
+        let mut result = VectorImpl::new(vec![T::zero(); ROWS * COLS]);
+
+        for i in 0..self.size() {
+            let sum = self.get(i) + other.get(i);
+            result.set(i, sum);
+        }
+
+        Some(result)
+    }
+
+    pub fn subtract(&self, other: &Self) -> Option<VectorImpl<T, ROWS, COLS>> {
+        if self.size() != other.size() {
+            return None;
+        }
+
+        let mut result = VectorImpl::new(vec![T::zero(); ROWS * COLS]);
+
+        for i in 0..self.size() {
+            let sum = self.get(i) - other.get(i);
+            result.set(i, sum);
+        }
+
+        Some(result)
+    }
+
+    pub fn scale(&self, scalar: T) -> VectorImpl<T, ROWS, COLS> {
+        let mut result = VectorImpl::new(vec![T::zero(); ROWS * COLS]);
+
+        for i in 0..self.size() {
+            let product = self.get(i) * scalar;
+            result.set(i, product);
+        }
+
+        result
+    }
+}
+
+impl<T: Float, const ROWS: usize, const COLS: usize> ops::Index<usize>
+    for VectorImpl<T, ROWS, COLS>
+{
+    type Output = T;
+
+    fn index(&self, index: usize) -> &Self::Output {
+        &self.0[index]
+    }
+}
+
+impl<T: Float, const ROWS: usize, const COLS: usize> ops::IndexMut<usize>
+    for VectorImpl<T, ROWS, COLS>
+{
+    fn index_mut(&mut self, index: usize) -> &mut Self::Output {
+        &mut self.0[index]
+    }
+}
+
+impl<T: Float + fmt::Debug> fmt::Debug for RowVector<T> {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        write!(f, "RowVector({:?})", self.0)
+    }
+}
+
+impl<T: Float + fmt::Debug> fmt::Debug for ColumnVector<T> {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        write!(f, "ColumnVector({:?})", self.0)
+    }
+}
+
+impl<T: Float + fmt::Display> fmt::Display for RowVector<T> {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        write!(
+            f,
+            "[{}]",
+            self.0
+                .iter()
+                .map(|x| x.to_string())
+                .collect::<Vec<String>>()
+                .join(", ")
+        )
+    }
+}
+
+impl<T: Float + fmt::Display> fmt::Display for ColumnVector<T> {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        write!(f, "{}", self.to_matrix())
+    }
+}
diff --git a/src/maths/vector_ops.rs b/src/maths/vector_ops.rs
new file mode 100644
index 0000000..1e7148c
--- /dev/null
+++ b/src/maths/vector_ops.rs
@@ -0,0 +1,107 @@
+use super::{Float, Matrix};
+use crate::{ColumnVector, RowVector, VectorImpl};
+use core::ops;
+
+impl<T: Float, const ROWS: usize, const COLS: usize> ops::Add<&VectorImpl<T, ROWS, COLS>>
+    for VectorImpl<T, ROWS, COLS>
+{
+    type Output = Option<VectorImpl<T, ROWS, COLS>>;
+
+    fn add(self, other: &VectorImpl<T, ROWS, COLS>) -> Self::Output {
+        self.add_to(other)
+    }
+}
+
+impl<T: Float, const ROWS: usize, const COLS: usize> ops::Sub<&VectorImpl<T, ROWS, COLS>>
+    for VectorImpl<T, ROWS, COLS>
+{
+    type Output = Option<VectorImpl<T, ROWS, COLS>>;
+
+    fn sub(self, other: &VectorImpl<T, ROWS, COLS>) -> Self::Output {
+        self.subtract(other)
+    }
+}
+
+impl<T: Float, const ROWS: usize, const COLS: usize> ops::Mul<T> for VectorImpl<T, ROWS, COLS> {
+    type Output = VectorImpl<T, ROWS, COLS>;
+
+    fn mul(self, scalar: T) -> Self::Output {
+        self.scale(scalar)
+    }
+}
+
+impl<T: Float, const ROWS: usize, const COLS: usize> ops::Div<T> for VectorImpl<T, ROWS, COLS> {
+    type Output = VectorImpl<T, ROWS, COLS>;
+
+    fn div(self, scalar: T) -> Self::Output {
+        self.scale(T::one() / scalar)
+    }
+}
+
+impl<T: Float, const ROWS: usize, const COLS: usize> ops::AddAssign<VectorImpl<T, ROWS, COLS>>
+    for VectorImpl<T, ROWS, COLS>
+{
+    fn add_assign(&mut self, other: VectorImpl<T, ROWS, COLS>) {
+        if let Some(result) = self.add_to(&other) {
+            *self = result;
+        }
+    }
+}
+
+impl<T: Float, const ROWS: usize, const COLS: usize> ops::SubAssign<VectorImpl<T, ROWS, COLS>>
+    for VectorImpl<T, ROWS, COLS>
+{
+    fn sub_assign(&mut self, other: VectorImpl<T, ROWS, COLS>) {
+        if let Some(result) = self.subtract(&other) {
+            *self = result;
+        }
+    }
+}
+
+impl<T: Float, const ROWS: usize, const COLS: usize> ops::MulAssign<T>
+    for VectorImpl<T, ROWS, COLS>
+{
+    fn mul_assign(&mut self, scalar: T) {
+        *self = self.scale(scalar);
+    }
+}
+
+impl<T: Float, const ROWS: usize, const COLS: usize> ops::DivAssign<T>
+    for VectorImpl<T, ROWS, COLS>
+{
+    fn div_assign(&mut self, scalar: T) {
+        *self = self.scale(T::one() / scalar);
+    }
+}
+
+impl<T: Float> ops::Mul<&Matrix<T>> for RowVector<T> {
+    type Output = Option<RowVector<T>>;
+
+    fn mul(self, matrix: &Matrix<T>) -> Self::Output {
+        self.mul_matrix(matrix)
+    }
+}
+
+impl<T: Float> ops::Mul<&Matrix<T>> for ColumnVector<T> {
+    type Output = Option<ColumnVector<T>>;
+
+    fn mul(self, matrix: &Matrix<T>) -> Self::Output {
+        self.mul_matrix(matrix)
+    }
+}
+
+impl<T: Float> ops::MulAssign<&Matrix<T>> for RowVector<T> {
+    fn mul_assign(&mut self, matrix: &Matrix<T>) {
+        if let Some(result) = self.mul_matrix(matrix) {
+            *self = result;
+        }
+    }
+}
+
+impl<T: Float> ops::MulAssign<&Matrix<T>> for ColumnVector<T> {
+    fn mul_assign(&mut self, matrix: &Matrix<T>) {
+        if let Some(result) = self.mul_matrix(matrix) {
+            *self = result;
+        }
+    }
+}
diff --git a/src/visualizer/horizontal_cli.rs b/src/visualizer/horizontal_cli.rs
new file mode 100644
index 0000000..28cdb5a
--- /dev/null
+++ b/src/visualizer/horizontal_cli.rs
@@ -0,0 +1,463 @@
+use super::renderer::Visualizer;
+use core::fmt;
+use crate::{GateOp, QuantumCircuit};
+
+pub struct HorizontalRenderer<'a> {
+    circuit: &'a QuantumCircuit,
+}
+
+impl<'a> HorizontalRenderer<'a> {
+    pub fn new(circuit: &'a QuantumCircuit) -> Self {
+        HorizontalRenderer { circuit }
+    }
+}
+
+impl<'a> Visualizer for HorizontalRenderer<'a> {
+    fn export(&self) -> String {
+        format!("{}", self)
+    }
+}
+
+impl<'a> fmt::Display for HorizontalRenderer<'a> {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        let nq = self.circuit.num_qubits();
+        let nc = self.circuit.num_classical();
+        let ops = self.circuit.operations();
+
+        let mut q_lines: Vec<String> = (0..nq).map(|i| format!("q{}: ", i)).collect();
+        let mut c_lines: Vec<String> = (0..nc).map(|i| format!("c{}: ", i)).collect();
+
+        let max_label = q_lines
+            .iter()
+            .chain(c_lines.iter())
+            .map(|s| s.len())
+            .max()
+            .unwrap_or(3);
+
+        for line in &mut q_lines {
+            while line.len() < max_label {
+                line.insert(0, ' ');
+            }
+        }
+        for line in &mut c_lines {
+            while line.len() < max_label {
+                line.insert(0, ' ');
+            }
+        }
+        let mut gap_line = " ".repeat(max_label);
+
+        if ops.is_empty() {
+            for line in &q_lines {
+                writeln!(f, "{}───░", line)?;
+            }
+            if nc > 0 {
+                writeln!(f, "{}   ░", gap_line)?;
+                for line in &c_lines {
+                    writeln!(f, "{}═══░", line)?;
+                }
+            }
+            return Ok(());
+        }
+
+        for op in ops {
+            let q_targets = op.quantum_targets();
+
+            let min_q = q_targets.iter().min().copied().unwrap_or(0);
+            let max_q = q_targets.iter().max().copied().unwrap_or(0);
+
+            match op {
+                GateOp::H(t) => {
+                    for (i, line) in q_lines.iter_mut().enumerate() {
+                        if i == *t {
+                            line.push_str("─[H]─");
+                        } else {
+                            line.push_str("─────");
+                        }
+                    }
+                    for line in c_lines.iter_mut() {
+                        line.push_str("═════");
+                    }
+                    gap_line.push_str("     ");
+                }
+                GateOp::X(t) => {
+                    for (i, line) in q_lines.iter_mut().enumerate() {
+                        if i == *t {
+                            line.push_str("─[X]─");
+                        } else {
+                            line.push_str("─────");
+                        }
+                    }
+                    for line in c_lines.iter_mut() {
+                        line.push_str("═════");
+                    }
+                    gap_line.push_str("     ");
+                }
+                GateOp::Y(t) => {
+                    for (i, line) in q_lines.iter_mut().enumerate() {
+                        if i == *t {
+                            line.push_str("─[Y]─");
+                        } else {
+                            line.push_str("─────");
+                        }
+                    }
+                    for line in c_lines.iter_mut() {
+                        line.push_str("═════");
+                    }
+                    gap_line.push_str("     ");
+                }
+                GateOp::Z(t) => {
+                    for (i, line) in q_lines.iter_mut().enumerate() {
+                        if i == *t {
+                            line.push_str("─[Z]─");
+                        } else {
+                            line.push_str("─────");
+                        }
+                    }
+                    for line in c_lines.iter_mut() {
+                        line.push_str("═════");
+                    }
+                    gap_line.push_str("     ");
+                }
+                GateOp::S(t) => {
+                    for (i, line) in q_lines.iter_mut().enumerate() {
+                        if i == *t {
+                            line.push_str("─[S]─");
+                        } else {
+                            line.push_str("─────");
+                        }
+                    }
+                    for line in c_lines.iter_mut() {
+                        line.push_str("═════");
+                    }
+                    gap_line.push_str("     ");
+                }
+                GateOp::T(t) => {
+                    for (i, line) in q_lines.iter_mut().enumerate() {
+                        if i == *t {
+                            line.push_str("─[T]─");
+                        } else {
+                            line.push_str("─────");
+                        }
+                    }
+                    for line in c_lines.iter_mut() {
+                        line.push_str("═════");
+                    }
+                    gap_line.push_str("     ");
+                }
+                GateOp::Sdg(t) => {
+                    for (i, line) in q_lines.iter_mut().enumerate() {
+                        if i == *t {
+                            line.push_str("─[S†]─");
+                        } else {
+                            line.push_str("──────");
+                        }
+                    }
+                    for line in c_lines.iter_mut() {
+                        line.push_str("══════");
+                    }
+                    gap_line.push_str("      ");
+                }
+                GateOp::Tdg(t) => {
+                    for (i, line) in q_lines.iter_mut().enumerate() {
+                        if i == *t {
+                            line.push_str("─[T†]─");
+                        } else {
+                            line.push_str("──────");
+                        }
+                    }
+                    for line in c_lines.iter_mut() {
+                        line.push_str("══════");
+                    }
+                    gap_line.push_str("      ");
+                }
+                GateOp::Sx(t) => {
+                    for (i, line) in q_lines.iter_mut().enumerate() {
+                        if i == *t {
+                            line.push_str("─[√X]─");
+                        } else {
+                            line.push_str("──────");
+                        }
+                    }
+                    for line in c_lines.iter_mut() {
+                        line.push_str("══════");
+                    }
+                    gap_line.push_str("      ");
+                }
+                GateOp::Sxdg(t) => {
+                    for (i, line) in q_lines.iter_mut().enumerate() {
+                        if i == *t {
+                            line.push_str("─[√X†]─");
+                        } else {
+                            line.push_str("───────");
+                        }
+                    }
+                    for line in c_lines.iter_mut() {
+                        line.push_str("═══════");
+                    }
+                    gap_line.push_str("       ");
+                }
+                GateOp::Rx(t, theta) => {
+                    let label = format!("[Rx({:.2})]", theta);
+                    for (i, line) in q_lines.iter_mut().enumerate() {
+                        if i == *t {
+                            line.push_str(&format!("─{}─", label));
+                        } else {
+                            line.push_str(&format!("─{}─", "─".repeat(label.len())));
+                        }
+                    }
+                    for line in c_lines.iter_mut() {
+                        line.push_str(&format!("═{}═", "═".repeat(label.len())));
+                    }
+                    gap_line.push_str(&format!(" {} ", " ".repeat(label.len())));
+                }
+                GateOp::Ry(t, theta) => {
+                    let label = format!("[Ry({:.2})]", theta);
+                    for (i, line) in q_lines.iter_mut().enumerate() {
+                        if i == *t {
+                            line.push_str(&format!("─{}─", label));
+                        } else {
+                            line.push_str(&format!("─{}─", "─".repeat(label.len())));
+                        }
+                    }
+                    for line in c_lines.iter_mut() {
+                        line.push_str(&format!("═{}═", "═".repeat(label.len())));
+                    }
+                    gap_line.push_str(&format!(" {} ", " ".repeat(label.len())));
+                }
+                GateOp::Rz(t, theta) => {
+                    let label = format!("[Rz({:.2})]", theta);
+                    for (i, line) in q_lines.iter_mut().enumerate() {
+                        if i == *t {
+                            line.push_str(&format!("─{}─", label));
+                        } else {
+                            line.push_str(&format!("─{}─", "─".repeat(label.len())));
+                        }
+                    }
+                    for line in c_lines.iter_mut() {
+                        line.push_str(&format!("═{}═", "═".repeat(label.len())));
+                    }
+                    gap_line.push_str(&format!(" {} ", " ".repeat(label.len())));
+                }
+                GateOp::P(t, theta) => {
+                    let label = format!("[P({:.2})]", theta);
+                    for (i, line) in q_lines.iter_mut().enumerate() {
+                        if i == *t {
+                            line.push_str(&format!("─{}─", label));
+                        } else {
+                            line.push_str(&format!("─{}─", "─".repeat(label.len())));
+                        }
+                    }
+                    for line in c_lines.iter_mut() {
+                        line.push_str(&format!("═{}═", "═".repeat(label.len())));
+                    }
+                    gap_line.push_str(&format!(" {} ", " ".repeat(label.len())));
+                }
+                GateOp::U1(t, lambda) => {
+                    let label = format!("[U1({:.2})]", lambda);
+                    for (i, line) in q_lines.iter_mut().enumerate() {
+                        if i == *t {
+                            line.push_str(&format!("─{}─", label));
+                        } else {
+                            line.push_str(&format!("─{}─", "─".repeat(label.len())));
+                        }
+                    }
+                    for line in c_lines.iter_mut() {
+                        line.push_str(&format!("═{}═", "═".repeat(label.len())));
+                    }
+                    gap_line.push_str(&format!(" {} ", " ".repeat(label.len())));
+                }
+                GateOp::U2(t, _, _) => {
+                    let label = "[U2]";
+                    for (i, line) in q_lines.iter_mut().enumerate() {
+                        if i == *t {
+                            line.push_str(&format!("─{}─", label));
+                        } else {
+                            line.push_str(&format!("─{}─", "─".repeat(label.len())));
+                        }
+                    }
+                    for line in c_lines.iter_mut() {
+                        line.push_str(&format!("═{}═", "═".repeat(label.len())));
+                    }
+                    gap_line.push_str(&format!(" {} ", " ".repeat(label.len())));
+                }
+                GateOp::U3(t, _, _, _) => {
+                    let label = "[U3]";
+                    for (i, line) in q_lines.iter_mut().enumerate() {
+                        if i == *t {
+                            line.push_str(&format!("─{}─", label));
+                        } else {
+                            line.push_str(&format!("─{}─", "─".repeat(label.len())));
+                        }
+                    }
+                    for line in c_lines.iter_mut() {
+                        line.push_str(&format!("═{}═", "═".repeat(label.len())));
+                    }
+                    gap_line.push_str(&format!(" {} ", " ".repeat(label.len())));
+                }
+                GateOp::CRx(c, t, theta) | GateOp::CRy(c, t, theta) | GateOp::CRz(c, t, theta) | GateOp::CP(c, t, theta) => {
+                    let label = match op {
+                        GateOp::CRx(_, _, _) => format!("[CRx({:.2})]", theta),
+                        GateOp::CRy(_, _, _) => format!("[CRy({:.2})]", theta),
+                        GateOp::CRz(_, _, _) => format!("[CRz({:.2})]", theta),
+                        GateOp::CP(_, _, _) => format!("[CP({:.2})]", theta),
+                        _ => unreachable!(),
+                    };
+                    for (i, line) in q_lines.iter_mut().enumerate() {
+                        if i == *c {
+                            line.push_str(&format!("─{}─", "●".to_string() + &"─".repeat(label.len() - 1)));
+                        } else if i == *t {
+                            line.push_str(&format!("─{}─", label));
+                        } else if i > min_q && i < max_q {
+                            line.push_str(&format!("─{}─", "│".to_string() + &"─".repeat(label.len() - 1)));
+                        } else {
+                            line.push_str(&format!("─{}─", "─".repeat(label.len())));
+                        }
+                    }
+                    for line in c_lines.iter_mut() {
+                        line.push_str(&format!("═{}═", "═".repeat(label.len())));
+                    }
+                    gap_line.push_str(&format!(" {} ", " ".repeat(label.len())));
+                }
+                GateOp::CNOT(c, t) => {
+                    for (i, line) in q_lines.iter_mut().enumerate() {
+                        if i == *c {
+                            line.push_str("──●──");
+                        } else if i == *t {
+                            line.push_str("──⊕──");
+                        } else if i > min_q && i < max_q {
+                            line.push_str("──│──");
+                        } else {
+                            line.push_str("─────");
+                        }
+                    }
+                    for line in c_lines.iter_mut() {
+                        line.push_str("═════");
+                    }
+                    gap_line.push_str("     ");
+                }
+                GateOp::CZ(c, t) => {
+                    for (i, line) in q_lines.iter_mut().enumerate() {
+                        if i == *c || i == *t {
+                            line.push_str("──●──");
+                        } else if i > min_q && i < max_q {
+                            line.push_str("──│──");
+                        } else {
+                            line.push_str("─────");
+                        }
+                    }
+                    for line in c_lines.iter_mut() {
+                        line.push_str("═════");
+                    }
+                    gap_line.push_str("     ");
+                }
+                GateOp::SWAP(a, b) => {
+                    for (i, line) in q_lines.iter_mut().enumerate() {
+                        if i == *a || i == *b {
+                            line.push_str("──╳──");
+                        } else if i > min_q && i < max_q {
+                            line.push_str("──│──");
+                        } else {
+                            line.push_str("─────");
+                        }
+                    }
+                    for line in c_lines.iter_mut() {
+                        line.push_str("═════");
+                    }
+                    gap_line.push_str("     ");
+                }
+                GateOp::CCNOT(c1, c2, t) => {
+                    for (i, line) in q_lines.iter_mut().enumerate() {
+                        if i == *c1 || i == *c2 {
+                            line.push_str("──●──");
+                        } else if i == *t {
+                            line.push_str("──⊕──");
+                        } else if i > min_q && i < max_q {
+                            line.push_str("──│──");
+                        } else {
+                            line.push_str("─────");
+                        }
+                    }
+                    for line in c_lines.iter_mut() {
+                        line.push_str("═════");
+                    }
+                    gap_line.push_str("     ");
+                }
+                GateOp::CSWAP(c, t1, t2) => {
+                    for (i, line) in q_lines.iter_mut().enumerate() {
+                        if i == *c {
+                            line.push_str("──●──");
+                        } else if i == *t1 || i == *t2 {
+                            line.push_str("──╳──");
+                        } else if i > min_q && i < max_q {
+                            line.push_str("──│──");
+                        } else {
+                            line.push_str("─────");
+                        }
+                    }
+                    for line in c_lines.iter_mut() {
+                        line.push_str("═════");
+                    }
+                    gap_line.push_str("     ");
+                }
+                GateOp::Measure(q, c) => {
+                    for (i, line) in q_lines.iter_mut().enumerate() {
+                        if i == *q {
+                            line.push_str("─[M]─");
+                        } else if i > *q {
+                            line.push_str("──║──");
+                        } else {
+                            line.push_str("─────");
+                        }
+                    }
+                    for (i, line) in c_lines.iter_mut().enumerate() {
+                        if i == *c {
+                            line.push_str("══╩══");
+                        } else if i < *c {
+                            line.push_str("══║══");
+                        } else {
+                            line.push_str("═════");
+                        }
+                    }
+                    gap_line.push_str("  ║  ");
+                }
+                GateOp::Custom(gate, targets) => {
+                    let name = &gate.name;
+                    let label = format!("[{}]", name);
+
+                    for (i, line) in q_lines.iter_mut().enumerate() {
+                        if targets.contains(&i) {
+                            if i == targets[0] {
+                                line.push_str(&format!("─{}─", label));
+                            } else {
+                                line.push_str(&format!("─{}─", "─".repeat(label.len())));
+                            }
+                        } else if i > min_q && i < max_q {
+                            line.push_str(&format!(
+                                "─{}─",
+                                "│".to_string() + &"─".repeat(label.len() - 1)
+                            ));
+                        } else {
+                            line.push_str(&format!("─{}─", "─".repeat(label.len())));
+                        }
+                    }
+                    for line in c_lines.iter_mut() {
+                        line.push_str(&format!("═{}═", "═".repeat(label.len())));
+                    }
+                    gap_line.push_str(&format!(" {} ", " ".repeat(label.len())));
+                }
+            }
+        }
+
+        for line in &q_lines {
+            writeln!(f, "{}░", line)?;
+        }
+        if nc > 0 {
+            writeln!(f, "{}░", gap_line)?;
+            for line in &c_lines {
+                writeln!(f, "{}░", line)?;
+            }
+        }
+
+        Ok(())
+    }
+}
diff --git a/src/visualizer/mod.rs b/src/visualizer/mod.rs
new file mode 100644
index 0000000..ff53c44
--- /dev/null
+++ b/src/visualizer/mod.rs
@@ -0,0 +1,7 @@
+pub mod horizontal_cli;
+pub mod vertical_cli;
+pub mod renderer;
+
+pub use horizontal_cli::*;
+pub use vertical_cli::*;
+pub use renderer::*;
diff --git a/src/visualizer/renderer.rs b/src/visualizer/renderer.rs
new file mode 100644
index 0000000..33f730f
--- /dev/null
+++ b/src/visualizer/renderer.rs
@@ -0,0 +1,3 @@
+pub trait Visualizer {
+    fn export(&self) -> String;
+}
diff --git a/src/visualizer/vertical_cli.rs b/src/visualizer/vertical_cli.rs
new file mode 100644
index 0000000..e9322bf
--- /dev/null
+++ b/src/visualizer/vertical_cli.rs
@@ -0,0 +1,367 @@
+use super::renderer::Visualizer;
+use core::fmt;
+use crate::{GateOp, QuantumCircuit};
+
+pub struct VerticalRenderer<'a> {
+    circuit: &'a QuantumCircuit,
+}
+
+impl<'a> VerticalRenderer<'a> {
+    pub fn new(circuit: &'a QuantumCircuit) -> Self {
+        VerticalRenderer { circuit }
+    }
+
+    fn gate_label(op: &GateOp) -> String {
+        match op {
+            GateOp::H(_) => "[H]".to_string(),
+            GateOp::X(_) => "[X]".to_string(),
+            GateOp::Y(_) => "[Y]".to_string(),
+            GateOp::Z(_) => "[Z]".to_string(),
+            GateOp::S(_) => "[S]".to_string(),
+            GateOp::T(_) => "[T]".to_string(),
+            GateOp::Sdg(_) => "[S†]".to_string(),
+            GateOp::Tdg(_) => "[T†]".to_string(),
+            GateOp::Sx(_) => "[√X]".to_string(),
+            GateOp::Sxdg(_) => "[√X†]".to_string(),
+            GateOp::Rx(_, theta) => format!("[Rx({:.2})]", theta),
+            GateOp::Ry(_, theta) => format!("[Ry({:.2})]", theta),
+            GateOp::Rz(_, theta) => format!("[Rz({:.2})]", theta),
+            GateOp::P(_, theta) => format!("[P({:.2})]", theta),
+            GateOp::U1(_, lambda) => format!("[U1({:.2})]", lambda),
+            GateOp::U2(_, _, _) => "[U2]".to_string(),
+            GateOp::U3(_, _, _, _) => "[U3]".to_string(),
+            GateOp::CRx(_, _, _) => "[CRx]".to_string(),
+            GateOp::CRy(_, _, _) => "[CRy]".to_string(),
+            GateOp::CRz(_, _, _) => "[CRz]".to_string(),
+            GateOp::CP(_, _, _) => "[CP]".to_string(),
+            GateOp::CNOT(_, _) => "●".to_string(),
+            GateOp::CZ(_, _) => "●".to_string(),
+            GateOp::SWAP(_, _) => "╳".to_string(),
+            GateOp::CCNOT(_, _, _) => "●".to_string(),
+            GateOp::CSWAP(_, _, _) => "●".to_string(),
+            GateOp::Measure(_, _) => "[M]".to_string(),
+            GateOp::Custom(gate, _) => format!("[{}]", gate.name),
+        }
+    }
+
+    fn calculate_col_width(&self) -> usize {
+        let min_width = 3;
+        let mut max_label_len = min_width;
+
+        for op in self.circuit.operations() {
+            let label = Self::gate_label(op);
+            let char_count: usize = label.chars().count();
+            if char_count > max_label_len {
+                max_label_len = char_count;
+            }
+        }
+
+        let width = max_label_len + 2;
+        if width % 2 == 0 {
+            width + 1
+        } else {
+            width
+        }
+    }
+}
+
+impl<'a> Visualizer for VerticalRenderer<'a> {
+    fn export(&self) -> String {
+        format!("{}", self)
+    }
+}
+
+impl<'a> fmt::Display for VerticalRenderer<'a> {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        let nq = self.circuit.num_qubits();
+        let nc = self.circuit.num_classical();
+        let ops = self.circuit.operations();
+
+        let col_width = self.calculate_col_width();
+        let gap_width = 3;
+
+        let q_header: String = (0..nq)
+            .map(|i| format!("{:^width$}", format!("q{}", i), width = col_width))
+            .collect::<Vec<_>>()
+            .join(" ");
+
+        let c_header: String = (0..nc)
+            .map(|i| format!("{:^width$}", format!("c{}", i), width = col_width))
+            .collect::<Vec<_>>()
+            .join(" ");
+
+        if nc > 0 {
+            writeln!(f, "{}{}{}", q_header, " ".repeat(gap_width), c_header)?;
+        } else {
+            writeln!(f, "{}", q_header)?;
+        }
+
+        let q_wires: String = (0..nq)
+            .map(|_| format!("{:^width$}", "│", width = col_width))
+            .collect::<Vec<_>>()
+            .join(" ");
+
+        let c_wires: String = (0..nc)
+            .map(|_| format!("{:^width$}", "║", width = col_width))
+            .collect::<Vec<_>>()
+            .join(" ");
+
+        let full_wires = if nc > 0 {
+            format!("{}{}{}", q_wires, " ".repeat(gap_width), c_wires)
+        } else {
+            q_wires.clone()
+        };
+
+        if ops.is_empty() {
+            writeln!(f, "{}", full_wires)?;
+            return Ok(());
+        }
+
+        let q_total = nq * col_width + (nq - 1);
+        let c_total = if nc > 0 { nc * col_width + (nc - 1) } else { 0 };
+        let total_width = q_total + gap_width + c_total;
+
+        for op in ops {
+            writeln!(f, "{}", full_wires)?;
+
+            let q_targets = op.quantum_targets();
+            let min_q = q_targets.iter().min().copied().unwrap_or(0);
+            let max_q = q_targets.iter().max().copied().unwrap_or(0);
+
+            let label = Self::gate_label(op);
+
+            match op {
+                GateOp::H(t)
+                | GateOp::X(t)
+                | GateOp::Y(t)
+                | GateOp::Z(t)
+                | GateOp::S(t)
+                | GateOp::T(t)
+                | GateOp::Sdg(t)
+                | GateOp::Tdg(t)
+                | GateOp::Sx(t)
+                | GateOp::Sxdg(t)
+                | GateOp::Rx(t, _)
+                | GateOp::Ry(t, _)
+                | GateOp::Rz(t, _)
+                | GateOp::P(t, _)
+                | GateOp::U1(t, _)
+                | GateOp::U2(t, _, _)
+                | GateOp::U3(t, _, _, _) => {
+                    let mut line: Vec<char> = vec![' '; total_width];
+
+                    for i in 0..nq {
+                        let col_start = i * (col_width + 1);
+                        let center = col_start + col_width / 2;
+                        if i == *t {
+                            let label_start = col_start + (col_width - label.chars().count()) / 2;
+                            for (j, ch) in label.chars().enumerate() {
+                                line[label_start + j] = ch;
+                            }
+                        } else {
+                            line[center] = '│';
+                        }
+                    }
+
+                    for i in 0..nc {
+                        let center = q_total + gap_width + i * (col_width + 1) + col_width / 2;
+                        line[center] = '║';
+                    }
+
+                    let gate_line: String = line.into_iter().collect();
+                    writeln!(f, "{}", gate_line)?;
+                }
+                GateOp::CNOT(c, t) | GateOp::CZ(c, t) | GateOp::SWAP(c, t)
+                | GateOp::CRx(c, t, _) | GateOp::CRy(c, t, _) | GateOp::CRz(c, t, _) | GateOp::CP(c, t, _) => {
+                    let (sym1, sym2) = match op {
+                        GateOp::CNOT(_, _) => ('●', '⊕'),
+                        GateOp::CZ(_, _) => ('●', '●'),
+                        GateOp::SWAP(_, _) => ('╳', '╳'),
+                        GateOp::CRx(_, _, _) | GateOp::CRy(_, _, _) | GateOp::CRz(_, _, _) | GateOp::CP(_, _, _) => ('●', '□'),
+                        _ => unreachable!(),
+                    };
+
+                    let mut line: Vec<char> = vec![' '; total_width];
+
+                    for i in 0..nq {
+                        let col_start = i * (col_width + 1);
+                        let center = col_start + col_width / 2;
+                        if i < min_q || i > max_q {
+                            line[center] = '│';
+                        } else if i == *c {
+                            line[center] = sym1;
+                        } else if i == *t {
+                            // For controlled parametric gates, show the gate label on target
+                            if matches!(op, GateOp::CRx(_, _, _) | GateOp::CRy(_, _, _) | GateOp::CRz(_, _, _) | GateOp::CP(_, _, _)) {
+                                let label_start = col_start + (col_width - label.chars().count()) / 2;
+                                for (j, ch) in label.chars().enumerate() {
+                                    if label_start + j < line.len() {
+                                        line[label_start + j] = ch;
+                                    }
+                                }
+                            } else {
+                                line[center] = sym2;
+                            }
+                        }
+                    }
+
+                    let min_center = min_q * (col_width + 1) + col_width / 2;
+                    let max_center = max_q * (col_width + 1) + col_width / 2;
+                    for cell in &mut line[(min_center + 1)..max_center] {
+                        if *cell == ' ' {
+                            *cell = '─';
+                        }
+                    }
+
+                    for i in 0..nc {
+                        let center = q_total + gap_width + i * (col_width + 1) + col_width / 2;
+                        line[center] = '║';
+                    }
+
+                    let gate_line: String = line.into_iter().collect();
+                    writeln!(f, "{}", gate_line)?;
+                }
+                GateOp::CCNOT(c1, c2, t) | GateOp::CSWAP(c1, c2, t) => {
+                    let (sym_c, sym_t) = match op {
+                        GateOp::CCNOT(_, _, _) => ('●', '⊕'),
+                        GateOp::CSWAP(_, _, _) => ('●', '╳'),
+                        _ => unreachable!(),
+                    };
+                    let is_cswap = matches!(op, GateOp::CSWAP(_, _, _));
+
+                    let mut line: Vec<char> = vec![' '; total_width];
+
+                    for i in 0..nq {
+                        let center = i * (col_width + 1) + col_width / 2;
+                        if i < min_q || i > max_q {
+                            line[center] = '│';
+                        } else if i == *c1 {
+                            line[center] = sym_c;
+                        } else if i == *c2 {
+                            line[center] = if is_cswap { sym_t } else { sym_c };
+                        } else if i == *t {
+                            line[center] = sym_t;
+                        }
+                    }
+
+                    let min_center = min_q * (col_width + 1) + col_width / 2;
+                    let max_center = max_q * (col_width + 1) + col_width / 2;
+                    for cell in &mut line[(min_center + 1)..max_center] {
+                        if *cell == ' ' {
+                            *cell = '─';
+                        }
+                    }
+
+                    for i in 0..nc {
+                        let center = q_total + gap_width + i * (col_width + 1) + col_width / 2;
+                        line[center] = '║';
+                    }
+
+                    let gate_line: String = line.into_iter().collect();
+                    writeln!(f, "{}", gate_line)?;
+                }
+                GateOp::Measure(mq, mc) => {
+                    let mut line: Vec<char> = vec![' '; total_width];
+
+                    for i in 0..nq {
+                        let col_start = i * (col_width + 1);
+                        let center = col_start + col_width / 2;
+                        if i < *mq {
+                            line[center] = '│';
+                        } else if i == *mq {
+                            let label_start = col_start + (col_width - label.chars().count()) / 2;
+                            for (j, ch) in label.chars().enumerate() {
+                                line[label_start + j] = ch;
+                            }
+                        }
+                    }
+
+                    let mq_col_start = *mq * (col_width + 1);
+                    let mq_center = mq_col_start + col_width / 2;
+                    let mc_start = q_total + gap_width;
+                    let mc_center = mc_start + *mc * (col_width + 1) + col_width / 2;
+
+                    for cell in &mut line[(mq_center + 2)..=mc_center] {
+                        if *cell == ' ' {
+                            *cell = '═';
+                        }
+                    }
+                    line[mc_center] = '╣';
+
+                    for i in 0..nc {
+                        let center = mc_start + i * (col_width + 1) + col_width / 2;
+                        if i > *mc {
+                            line[center] = '║';
+                        }
+                    }
+
+                    let measure_line: String = line.into_iter().collect();
+                    writeln!(f, "{}", measure_line)?;
+                }
+                GateOp::Custom(_, targets) => {
+                    let mut line: Vec<char> = vec![' '; total_width];
+
+                    if targets.len() == 1 {
+                        for i in 0..nq {
+                            let col_start = i * (col_width + 1);
+                            let center = col_start + col_width / 2;
+                            if i == targets[0] {
+                                let label_start =
+                                    col_start + (col_width - label.chars().count()) / 2;
+                                for (j, ch) in label.chars().enumerate() {
+                                    line[label_start + j] = ch;
+                                }
+                            } else {
+                                line[center] = '│';
+                            }
+                        }
+
+                        for i in 0..nc {
+                            let center = q_total + gap_width + i * (col_width + 1) + col_width / 2;
+                            line[center] = '║';
+                        }
+                    } else {
+                        for i in 0..nq {
+                            let col_start = i * (col_width + 1);
+                            let center = col_start + col_width / 2;
+                            if i < min_q || i > max_q {
+                                line[center] = '│';
+                            } else if i == targets[0] {
+                                let label_start =
+                                    col_start + (col_width - label.chars().count()) / 2;
+                                for (j, ch) in label.chars().enumerate() {
+                                    line[label_start + j] = ch;
+                                }
+                            } else if targets.contains(&i) {
+                                line[center] = '□';
+                            }
+                        }
+
+                        let min_center = min_q * (col_width + 1) + col_width / 2;
+                        let max_center = max_q * (col_width + 1) + col_width / 2;
+                        for cell in &mut line[(min_center + 1)..max_center] {
+                            if *cell == ' ' {
+                                *cell = '─';
+                            }
+                        }
+
+                        for i in 0..nc {
+                            let center = q_total + gap_width + i * (col_width + 1) + col_width / 2;
+                            line[center] = '║';
+                        }
+                    }
+
+                    let gate_line: String = line.into_iter().collect();
+                    writeln!(f, "{}", gate_line)?;
+                }
+            }
+        }
+
+        writeln!(f, "{}", full_wires)?;
+
+        let end_line: String = "░".repeat(total_width);
+        writeln!(f, "{}", end_line)?;
+
+        Ok(())
+    }
+}