10 files changed, 438 insertions, 434 deletions

diff --git a/libpsi/src/maths/complex.rs b/libpsi/src/maths/complex.rs
index 28a3121..bb19510 100644
--- a/libpsi/src/maths/complex.rs
+++ b/libpsi/src/maths/complex.rs
@@ -1,5 +1,7 @@
 use crate::Float;
-use core::{fmt, ops};
+use core::fmt;
+
+use super::Numeric;
 
 #[macro_export]
 macro_rules! complex {
@@ -30,6 +32,26 @@ impl<T: Float + fmt::Display> fmt::Display for Complex<T> {
     }
 }
 
+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 }
@@ -58,162 +80,3 @@ impl<T: Float> Complex<T> {
         T::sqrt(self.norm())
     }
 }
-
-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> 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,
-        }
-    }
-}
-
-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,
-        }
-    }
-}
-
-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;
-    }
-}
-
-impl<T: Float> ops::AddAssign<T> for Complex<T> {
-    fn add_assign(&mut self, other: T) {
-        self.real += other;
-    }
-}
-
-impl<T: Float> ops::SubAssign<T> for Complex<T> {
-    fn sub_assign(&mut self, other: T) {
-        self.real -= other;
-    }
-}
-
-impl<T: Float> ops::MulAssign<T> for Complex<T> {
-    fn mul_assign(&mut self, other: T) {
-        self.real *= other;
-        self.imaginary *= other;
-    }
-}
-
-impl<T: Float> ops::DivAssign<T> for Complex<T> {
-    fn div_assign(&mut self, other: T) {
-        self.real /= other;
-        self.imaginary /= other;
-    }
-}
diff --git a/libpsi/src/maths/complex_ops.rs b/libpsi/src/maths/complex_ops.rs
new file mode 100644
index 0000000..bfcf920
--- /dev/null
+++ b/libpsi/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/src/maths/matrix.rs b/libpsi/src/maths/matrix.rs
index 388f858..795b37f 100644
--- a/libpsi/src/maths/matrix.rs
+++ b/libpsi/src/maths/matrix.rs
@@ -183,63 +183,3 @@ impl<T: Float + fmt::Display> fmt::Display for Matrix<T> {
         Ok(())
     }
 }
-
-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/src/maths/matrix_ops.rs b/libpsi/src/maths/matrix_ops.rs
new file mode 100644
index 0000000..15ab28b
--- /dev/null
+++ b/libpsi/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/src/maths/mod.rs b/libpsi/src/maths/mod.rs
index 734cfb4..a2b9b84 100644
--- a/libpsi/src/maths/mod.rs
+++ b/libpsi/src/maths/mod.rs
@@ -1,8 +1,16 @@
 pub mod complex;
+pub mod complex_ops;
+
 pub mod matrix;
-pub mod numeric_types;
+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_types::*;
+pub use numeric_traits::*;
diff --git a/libpsi/src/maths/numeric_types.rs b/libpsi/src/maths/numeric_float.rs
index f1355de..3591f24 100644
--- a/libpsi/src/maths/numeric_types.rs
+++ b/libpsi/src/maths/numeric_float.rs
@@ -1,28 +1,4 @@
-use super::Complex;
-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;
-}
+use super::{Complex, Float};
 
 impl Float for f32 {
     fn sqrt(self) -> Self {
@@ -87,66 +63,3 @@ impl Float for Complex<f64> {
         )
     }
 }
-
-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
-    }
-}
-
-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)
-    }
-}
diff --git a/libpsi/src/maths/numeric_int.rs b/libpsi/src/maths/numeric_int.rs
new file mode 100644
index 0000000..a9ecec7
--- /dev/null
+++ b/libpsi/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/src/maths/numeric_traits.rs b/libpsi/src/maths/numeric_traits.rs
new file mode 100644
index 0000000..dd7fdec
--- /dev/null
+++ b/libpsi/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/src/maths/vector.rs b/libpsi/src/maths/vector.rs
index 243061d..a21cf36 100644
--- a/libpsi/src/maths/vector.rs
+++ b/libpsi/src/maths/vector.rs
@@ -96,6 +96,26 @@ impl<T: Float> RowVector<T> {
     }
 }
 
+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)
@@ -211,130 +231,6 @@ impl<T: Float, const ROWS: usize, const COLS: usize> ops::IndexMut<usize>
     }
 }
 
-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> 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> 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;
-        }
-    }
-}
-
 impl<T: Float + fmt::Debug> fmt::Debug for RowVector<T> {
     fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
         write!(f, "RowVector({:?})", self.0)
diff --git a/libpsi/src/maths/vector_ops.rs b/libpsi/src/maths/vector_ops.rs
new file mode 100644
index 0000000..1e7148c
--- /dev/null
+++ b/libpsi/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;
+        }
+    }
+}