15 files changed, 1190 insertions, 0 deletions
diff --git a/libpsi-core/Cargo.toml b/libpsi-core/Cargo.toml
new file mode 100644
index 0000000..2657d02
--- /dev/null
+++ b/libpsi-core/Cargo.toml
@@ -0,0 +1,8 @@
+[package]
+name = "libpsi-core"
+version = "0.1.0"
+edition = "2021"
+
+[dependencies]
+lazy_static = "1.5.0"
+libm = "0.2.8"
diff --git a/libpsi-core/src/core/component.rs b/libpsi-core/src/core/component.rs
new file mode 100644
index 0000000..41b7cc7
--- /dev/null
+++ b/libpsi-core/src/core/component.rs
@@ -0,0 +1,141 @@
+use crate::{complex, ColumnVector, Complex, Matrix, Vector, VectorMatrix};
+use core::{fmt, ops};
+
+pub type QuantumBit = ColumnVector<Complex<f64>>;
+pub type QuantumGate = Matrix<Complex<f64>>;
+
+#[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 struct ClassicalRegister {
+    bits: Vec<i32>,
+}
+
+pub struct QuantumRegister {
+    state: ColumnVector<Complex<f64>>,
+    qubits: Vec<QuantumBit>,
+}
+
+impl QuantumBit {
+    pub fn measure(&self) -> i32 {
+        let alpha_abs = self[0].abs();
+        let beta_abs = self[1].abs();
+
+        let alpha_norm = alpha_abs * alpha_abs;
+        let beta_norm = beta_abs * beta_abs;
+
+        if alpha_norm > beta_norm {
+            0
+        } else {
+            1
+        }
+    }
+
+    pub fn state_0() -> QuantumBit {
+        QuantumBit::new(vec![complex!(1.0, 0.0), complex!(0.0, 0.0)])
+    }
+
+    pub fn state_1() -> QuantumBit {
+        QuantumBit::new(vec![complex!(0.0, 0.0), complex!(1.0, 0.0)])
+    }
+}
+
+impl ClassicalRegister {
+    pub fn new(count: usize) -> ClassicalRegister {
+        ClassicalRegister {
+            bits: Vec::with_capacity(count),
+        }
+    }
+
+    pub fn set_bits(&mut self, bits: Vec<i32>) {
+        self.bits = bits;
+    }
+}
+
+impl QuantumRegister {
+    fn update(&mut self) {
+        let matrices: Vec<Matrix<Complex<f64>>> =
+            self.qubits.iter().map(|qubit| qubit.to_matrix()).collect();
+        let mut new_result = matrices[0].clone();
+        for matrix in &matrices[1..] {
+            new_result = new_result.kronecker(matrix);
+        }
+
+        self.state = ColumnVector::from_matrix(&new_result);
+    }
+
+    pub fn from(bits: &mut [QuantumBit]) -> QuantumRegister {
+        let mut register = QuantumRegister {
+            qubits: bits.to_vec(),
+            state: ColumnVector::new(vec![]),
+        };
+
+        register.update();
+        register
+    }
+
+    pub fn measure(&self, classical_register: &mut ClassicalRegister) {
+        classical_register.set_bits(self.qubits.iter().map(|qubit| qubit.measure()).collect());
+    }
+
+    pub fn apply(&mut self, gate: &QuantumGate, index: usize) {
+        let result: ColumnVector<Complex<f64>> = self.state.mul_matrix(gate).unwrap();
+        self.qubits[index] = result;
+        self.update();
+    }
+}
+
+impl ops::Index<usize> for QuantumRegister {
+    type Output = QuantumBit;
+
+    fn index(&self, index: usize) -> &Self::Output {
+        &self.qubits[index]
+    }
+}
+
+impl ops::IndexMut<usize> for QuantumRegister {
+    fn index_mut(&mut self, index: usize) -> &mut Self::Output {
+        &mut self.qubits[index]
+    }
+}
+
+impl ops::Index<usize> for ClassicalRegister {
+    type Output = i32;
+
+    fn index(&self, index: usize) -> &Self::Output {
+        &self.bits[index]
+    }
+}
+
+impl ops::IndexMut<usize> for ClassicalRegister {
+    fn index_mut(&mut self, index: usize) -> &mut Self::Output {
+        &mut self.bits[index]
+    }
+}
diff --git a/libpsi-core/src/core/gates.rs b/libpsi-core/src/core/gates.rs
new file mode 100644
index 0000000..e9077ad
--- /dev/null
+++ b/libpsi-core/src/core/gates.rs
@@ -0,0 +1,22 @@
+use crate::{complex, matrix, QuantumGate};
+
+#[rustfmt::skip]
+lazy_static::lazy_static! {
+    pub static ref HADAMARD: QuantumGate = 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);
+
+    pub static ref PAULI_X: QuantumGate = matrix!([complex!(0.0, 0.0), complex!(1.0, 0.0)];
+                                                  [complex!(1.0, 0.0), complex!(0.0, 0.0)]);
+
+    pub static ref PAULI_Y: QuantumGate = matrix!([complex!(0.0, 0.0), complex!(0.0, -1.0)];
+                                                  [complex!(0.0, 1.0), complex!(0.0,  0.0)]);
+
+    pub static ref PAULI_Z: QuantumGate = matrix!([complex!(1.0, 0.0), complex!( 0.0, 0.0)];
+                                                  [complex!(0.0, 0.0), complex!(-1.0, 0.0)]);
+
+    pub static ref CNOT: QuantumGate = 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)]);
+}
diff --git a/libpsi-core/src/core/mod.rs b/libpsi-core/src/core/mod.rs
new file mode 100644
index 0000000..374c555
--- /dev/null
+++ b/libpsi-core/src/core/mod.rs
@@ -0,0 +1,4 @@
+pub mod component;
+pub mod gates;
+
+pub use gates::*;
diff --git a/libpsi-core/src/lib.rs b/libpsi-core/src/lib.rs
new file mode 100644
index 0000000..af50663
--- /dev/null
+++ b/libpsi-core/src/lib.rs
@@ -0,0 +1,10 @@
+pub mod core;
+pub mod maths;
+
+pub use maths::complex::*;
+pub use maths::matrix::*;
+pub use maths::numeric_traits::*;
+pub use maths::vector::*;
+
+pub use core::component::*;
+pub use core::gates;
diff --git a/libpsi-core/src/maths/complex.rs b/libpsi-core/src/maths/complex.rs
new file mode 100644
index 0000000..bb19510
--- /dev/null
+++ b/libpsi-core/src/maths/complex.rs
@@ -0,0 +1,82 @@
+use crate::Float;
+use core::fmt;
+
+use super::Numeric;
+
+#[macro_export]
+macro_rules! complex {
+    ($real:expr, $imaginary:expr) => {
+        $crate::Complex::new($real, $imaginary)
+    };
+}
+
+#[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 Numeric for Complex<f32> {
+    fn zero() -> Self {
+        Complex::new(0.0, 0.0)
+    }
+
+    fn one() -> Self {
+        Complex::new(1.0, 0.0)
+    }
+}
+
+impl Numeric for Complex<f64> {
+    fn zero() -> Self {
+        Complex::new(0.0, 0.0)
+    }
+
+    fn one() -> Self {
+        Complex::new(1.0, 0.0)
+    }
+}
+
+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 norm(&self) -> T {
+        self.real * self.real + self.imaginary * self.imaginary
+    }
+
+    pub fn abs(&self) -> T {
+        T::sqrt(self.norm())
+    }
+}
diff --git a/libpsi-core/src/maths/complex_ops.rs b/libpsi-core/src/maths/complex_ops.rs
new file mode 100644
index 0000000..bfcf920
--- /dev/null
+++ b/libpsi-core/src/maths/complex_ops.rs
@@ -0,0 +1,147 @@
+use super::{Complex, Float};
+use core::ops;
+
+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(),
+        }
+    }
+}
+
+// Complex-Complex
+impl<T: Float> ops::Add for Complex<T> {
+    type Output = Complex<T>;
+
+    fn add(self, other: Complex<T>) -> Complex<T> {
+        Complex {
+            real: self.real - other.real,
+            imaginary: self.imaginary - other.imaginary,
+        }
+    }
+}
+
+impl<T: Float> ops::Sub for Complex<T> {
+    type Output = Complex<T>;
+
+    fn sub(self, other: Complex<T>) -> Complex<T> {
+        Complex {
+            real: self.real - other.real,
+            imaginary: self.imaginary - other.imaginary,
+        }
+    }
+}
+
+impl<T: Float> ops::Mul for Complex<T> {
+    type Output = Complex<T>;
+
+    fn mul(self, other: Complex<T>) -> Complex<T> {
+        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> {
+        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,
+        }
+    }
+}
+
+// Complex-Complex, Assign
+impl<T: Float> ops::AddAssign for Complex<T> {
+    fn add_assign(&mut self, other: Complex<T>) {
+        self.real += other.real;
+        self.imaginary += other.imaginary;
+    }
+}
+
+impl<T: Float> ops::SubAssign for Complex<T> {
+    fn sub_assign(&mut self, other: Complex<T>) {
+        self.real -= other.real;
+        self.imaginary -= other.imaginary;
+    }
+}
+
+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_imaginary = self.real * other.imaginary + self.imaginary * other.real;
+        self.real = new_real;
+        self.imaginary = new_imaginary;
+    }
+}
+
+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_imaginary = (self.imaginary * other.real - self.real * other.imaginary) / denom;
+        self.real = new_real;
+        self.imaginary = new_imaginary;
+    }
+}
+
+// Real-Complex
+impl<T: Float> ops::Add<T> for Complex<T> {
+    type Output = Complex<T>;
+
+    fn add(self, other: T) -> Complex<T> {
+        Complex {
+            real: self.real + other,
+            imaginary: self.imaginary,
+        }
+    }
+}
+
+impl<T: Float> ops::Sub<T> for Complex<T> {
+    type Output = Complex<T>;
+
+    fn sub(self, other: T) -> Complex<T> {
+        Complex {
+            real: self.real - other,
+            imaginary: self.imaginary,
+        }
+    }
+}
+
+impl<T: Float> ops::Mul<T> for Complex<T> {
+    type Output = Complex<T>;
+
+    fn mul(self, other: T) -> Complex<T> {
+        Complex {
+            real: self.real * other,
+            imaginary: self.imaginary * other,
+        }
+    }
+}
+
+impl<T: Float> ops::Div<T> for Complex<T> {
+    type Output = Complex<T>;
+
+    fn div(self, other: T) -> Complex<T> {
+        Complex {
+            real: self.real / other,
+            imaginary: self.imaginary / other,
+        }
+    }
+}
diff --git a/libpsi-core/src/maths/matrix.rs b/libpsi-core/src/maths/matrix.rs
new file mode 100644
index 0000000..975c92d
--- /dev/null
+++ b/libpsi-core/src/maths/matrix.rs
@@ -0,0 +1,190 @@
+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)
+    }};
+}
+
+#[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 + 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 {
+        for i in 0..self.rows {
+            write!(f, "[")?;
+            for j in 0..self.cols {
+                write!(f, "{}", self.get(i, j))?;
+                if j != self.cols - 1 {
+                    write!(f, " ")?;
+                }
+            }
+            write!(f, "]")?;
+            if i != self.rows - 1 {
+                write!(f, "\n")?;
+            }
+        }
+        Ok(())
+    }
+}
diff --git a/libpsi-core/src/maths/matrix_ops.rs b/libpsi-core/src/maths/matrix_ops.rs
new file mode 100644
index 0000000..15ab28b
--- /dev/null
+++ b/libpsi-core/src/maths/matrix_ops.rs
@@ -0,0 +1,62 @@
+use super::{Float, Matrix};
+use core::ops;
+
+impl<T: Float> ops::Add<&Matrix<T>> for Matrix<T> {
+    type Output = Option<Matrix<T>>;
+
+    fn add(self, other: &Matrix<T>) -> Self::Output {
+        self.add_to(other)
+    }
+}
+
+impl<T: Float> ops::Sub<&Matrix<T>> for Matrix<T> {
+    type Output = Option<Matrix<T>>;
+
+    fn sub(self, other: &Matrix<T>) -> Self::Output {
+        self.subtract(other)
+    }
+}
+
+impl<T: Float> ops::Mul<T> for Matrix<T> {
+    type Output = Matrix<T>;
+
+    fn mul(self, scalar: T) -> Self::Output {
+        self.scale(scalar)
+    }
+}
+
+impl<T: Float> ops::Div<T> for Matrix<T> {
+    type Output = Matrix<T>;
+
+    fn div(self, scalar: T) -> Self::Output {
+        self.scale(T::one() / scalar)
+    }
+}
+
+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<T: Float> ops::MulAssign<T> for Matrix<T> {
+    fn mul_assign(&mut self, scalar: T) {
+        *self = self.scale(scalar);
+    }
+}
+
+impl<T: Float> ops::DivAssign<T> for Matrix<T> {
+    fn div_assign(&mut self, scalar: T) {
+        *self = self.scale(T::one() / scalar);
+    }
+}
diff --git a/libpsi-core/src/maths/mod.rs b/libpsi-core/src/maths/mod.rs
new file mode 100644
index 0000000..a2b9b84
--- /dev/null
+++ b/libpsi-core/src/maths/mod.rs
@@ -0,0 +1,16 @@
+pub mod complex;
+pub mod complex_ops;
+
+pub mod matrix;
+pub mod matrix_ops;
+
+pub mod numeric_float;
+pub mod numeric_int;
+pub mod numeric_traits;
+
+pub mod vector;
+pub mod vector_ops;
+
+pub use complex::*;
+pub use matrix::*;
+pub use numeric_traits::*;
diff --git a/libpsi-core/src/maths/numeric_float.rs b/libpsi-core/src/maths/numeric_float.rs
new file mode 100644
index 0000000..3591f24
--- /dev/null
+++ b/libpsi-core/src/maths/numeric_float.rs
@@ -0,0 +1,65 @@
+use super::{Complex, Float};
+
+impl Float for f32 {
+    fn sqrt(self) -> Self {
+        libm::sqrtf(self)
+    }
+
+    fn atan2(y: Self, x: Self) -> Self {
+        libm::atan2f(y, x)
+    }
+}
+
+impl Float for f64 {
+    fn sqrt(self) -> Self {
+        libm::sqrt(self)
+    }
+
+    fn atan2(y: Self, x: Self) -> Self {
+        libm::atan2(y, x)
+    }
+}
+
+impl Float for Complex<f32> {
+    fn sqrt(self) -> Self {
+        let r = self.abs();
+        let theta = self.phase();
+
+        let sqrt_r = libm::sqrtf(r);
+        let sqrt_theta = theta / 2.0;
+
+        Complex::new(
+            sqrt_r * libm::cosf(sqrt_theta),
+            sqrt_r * libm::sinf(sqrt_theta),
+        )
+    }
+
+    fn atan2(y: Self, x: Self) -> Self {
+        Complex::new(
+            libm::atan2f(y.real, x.real),
+            libm::atan2f(y.imaginary, x.imaginary),
+        )
+    }
+}
+
+impl Float for Complex<f64> {
+    fn sqrt(self) -> Self {
+        let r = self.abs();
+        let theta = self.phase();
+
+        let sqrt_r = libm::sqrt(r);
+        let sqrt_theta = theta / 2.0;
+
+        Complex::new(
+            sqrt_r * libm::cos(sqrt_theta),
+            sqrt_r * libm::sin(sqrt_theta),
+        )
+    }
+
+    fn atan2(y: Self, x: Self) -> Self {
+        Complex::new(
+            libm::atan2(y.real, x.real),
+            libm::atan2(y.imaginary, x.imaginary),
+        )
+    }
+}
diff --git a/libpsi-core/src/maths/numeric_int.rs b/libpsi-core/src/maths/numeric_int.rs
new file mode 100644
index 0000000..a9ecec7
--- /dev/null
+++ b/libpsi-core/src/maths/numeric_int.rs
@@ -0,0 +1,44 @@
+use super::{Integer, Numeric};
+
+impl Integer for i64 {}
+impl Integer for i32 {}
+
+impl Numeric for i32 {
+    fn zero() -> Self {
+        0
+    }
+
+    fn one() -> Self {
+        1
+    }
+}
+
+impl Numeric for i64 {
+    fn zero() -> Self {
+        0
+    }
+
+    fn one() -> Self {
+        1
+    }
+}
+
+impl Numeric for f32 {
+    fn zero() -> Self {
+        0.0
+    }
+
+    fn one() -> Self {
+        1.0
+    }
+}
+
+impl Numeric for f64 {
+    fn zero() -> Self {
+        0.0
+    }
+
+    fn one() -> Self {
+        1.0
+    }
+}
diff --git a/libpsi-core/src/maths/numeric_traits.rs b/libpsi-core/src/maths/numeric_traits.rs
new file mode 100644
index 0000000..dd7fdec
--- /dev/null
+++ b/libpsi-core/src/maths/numeric_traits.rs
@@ -0,0 +1,24 @@
+use core::ops;
+
+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;
+}
+
+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/vector.rs b/libpsi-core/src/maths/vector.rs
new file mode 100644
index 0000000..07ff7fb
--- /dev/null
+++ b/libpsi-core/src/maths/vector.rs
@@ -0,0 +1,268 @@
+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;
+}
+
+pub trait VectorMatrix<T: Float> {
+    fn from_matrix(matrix: &Matrix<T>) -> Self;
+    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())
+    }
+
+    fn from_matrix(matrix: &Matrix<T>) -> Self {
+        Self::new(matrix.data.clone())
+    }
+}
+
+impl<T: Float> VectorMatrix<T> for ColumnVector<T> {
+    fn to_matrix(&self) -> Matrix<T> {
+        Matrix::new(self.size(), 1, self.0.clone())
+    }
+
+    fn from_matrix(matrix: &Matrix<T>) -> Self {
+        Self::new(matrix.data.clone())
+    }
+}
+
+impl<T: Float, const ROWS: usize, const COLS: usize> Vector<T> for VectorImpl<T, ROWS, COLS> {
+    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, "[")?;
+        for (i, x) in self.0.iter().enumerate() {
+            if i > 0 {
+                write!(f, ",\n ")?;
+            }
+            write!(f, "{}", x)?;
+        }
+        write!(f, "]")
+    }
+}
diff --git a/libpsi-core/src/maths/vector_ops.rs b/libpsi-core/src/maths/vector_ops.rs
new file mode 100644
index 0000000..1e7148c
--- /dev/null
+++ b/libpsi-core/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;
+        }
+    }
+}