[chore]: unwrap project

8 files changed, 0 insertions, 1621 deletions

diff --git a/libpsi-core/src/maths/complex.rs b/libpsi-core/src/maths/complex.rs
deleted file mode 100644
index 31eae69..0000000
--- a/libpsi-core/src/maths/complex.rs
+++ /dev/null
@@ -1,180 +0,0 @@
-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/libpsi-core/src/maths/format.rs b/libpsi-core/src/maths/format.rs
deleted file mode 100644
index 957f1b6..0000000
--- a/libpsi-core/src/maths/format.rs
+++ /dev/null
@@ -1,141 +0,0 @@
-use crate::Complex;
-
-const EPSILON: f64 = 1e-10;
-const SQRT_2: f64 = 1.4142135623730951;
-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/libpsi-core/src/maths/matrix.rs b/libpsi-core/src/maths/matrix.rs
deleted file mode 100644
index 58492f9..0000000
--- a/libpsi-core/src/maths/matrix.rs
+++ /dev/null
@@ -1,310 +0,0 @@
-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 = 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.clone() * 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 {
-                write!(f, "\n")?;
-            }
-        }
-
-        Ok(())
-    }
-}
diff --git a/libpsi-core/src/maths/mod.rs b/libpsi-core/src/maths/mod.rs
deleted file mode 100644
index 85f4872..0000000
--- a/libpsi-core/src/maths/mod.rs
+++ /dev/null
@@ -1,14 +0,0 @@
-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/libpsi-core/src/maths/numeric.rs b/libpsi-core/src/maths/numeric.rs
deleted file mode 100644
index f5a649e..0000000
--- a/libpsi-core/src/maths/numeric.rs
+++ /dev/null
@@ -1,101 +0,0 @@
-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/libpsi-core/src/maths/simd.rs b/libpsi-core/src/maths/simd.rs
deleted file mode 100644
index de0370c..0000000
--- a/libpsi-core/src/maths/simd.rs
+++ /dev/null
@@ -1,510 +0,0 @@
-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/libpsi-core/src/maths/vector.rs b/libpsi-core/src/maths/vector.rs
deleted file mode 100644
index 87b5353..0000000
--- a/libpsi-core/src/maths/vector.rs
+++ /dev/null
@@ -1,258 +0,0 @@
-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 = 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 = 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/libpsi-core/src/maths/vector_ops.rs b/libpsi-core/src/maths/vector_ops.rs
deleted file mode 100644
index 1e7148c..0000000
--- a/libpsi-core/src/maths/vector_ops.rs
+++ /dev/null
@@ -1,107 +0,0 @@
-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;
-        }
-    }
-}