From c07ab7911bd50d0060841690a863333a9a3f226d Mon Sep 17 00:00:00 2001 From: hachem Date: Tue, 24 Sep 2024 21:33:51 +0200 Subject: rename project --- libpsi-core/Cargo.toml | 8 + libpsi-core/src/core/component.rs | 141 +++++++++++++++++ libpsi-core/src/core/gates.rs | 22 +++ libpsi-core/src/core/mod.rs | 4 + libpsi-core/src/lib.rs | 10 ++ libpsi-core/src/maths/complex.rs | 82 ++++++++++ libpsi-core/src/maths/complex_ops.rs | 147 ++++++++++++++++++ libpsi-core/src/maths/matrix.rs | 190 ++++++++++++++++++++++ libpsi-core/src/maths/matrix_ops.rs | 62 ++++++++ libpsi-core/src/maths/mod.rs | 16 ++ libpsi-core/src/maths/numeric_float.rs | 65 ++++++++ libpsi-core/src/maths/numeric_int.rs | 44 ++++++ libpsi-core/src/maths/numeric_traits.rs | 24 +++ libpsi-core/src/maths/vector.rs | 268 ++++++++++++++++++++++++++++++++ libpsi-core/src/maths/vector_ops.rs | 107 +++++++++++++ libpsi/Cargo.toml | 8 - libpsi/src/core/component.rs | 141 ----------------- libpsi/src/core/gates.rs | 22 --- libpsi/src/core/mod.rs | 4 - libpsi/src/lib.rs | 10 -- libpsi/src/maths/complex.rs | 82 ---------- libpsi/src/maths/complex_ops.rs | 147 ------------------ libpsi/src/maths/matrix.rs | 190 ---------------------- libpsi/src/maths/matrix_ops.rs | 62 -------- libpsi/src/maths/mod.rs | 16 -- libpsi/src/maths/numeric_float.rs | 65 -------- libpsi/src/maths/numeric_int.rs | 44 ------ libpsi/src/maths/numeric_traits.rs | 24 --- libpsi/src/maths/vector.rs | 268 -------------------------------- libpsi/src/maths/vector_ops.rs | 107 ------------- 30 files changed, 1190 insertions(+), 1190 deletions(-) create mode 100644 libpsi-core/Cargo.toml create mode 100644 libpsi-core/src/core/component.rs create mode 100644 libpsi-core/src/core/gates.rs create mode 100644 libpsi-core/src/core/mod.rs create mode 100644 libpsi-core/src/lib.rs create mode 100644 libpsi-core/src/maths/complex.rs create mode 100644 libpsi-core/src/maths/complex_ops.rs create mode 100644 libpsi-core/src/maths/matrix.rs create mode 100644 libpsi-core/src/maths/matrix_ops.rs create mode 100644 libpsi-core/src/maths/mod.rs create mode 100644 libpsi-core/src/maths/numeric_float.rs create mode 100644 libpsi-core/src/maths/numeric_int.rs create mode 100644 libpsi-core/src/maths/numeric_traits.rs create mode 100644 libpsi-core/src/maths/vector.rs create mode 100644 libpsi-core/src/maths/vector_ops.rs delete mode 100644 libpsi/Cargo.toml delete mode 100644 libpsi/src/core/component.rs delete mode 100644 libpsi/src/core/gates.rs delete mode 100644 libpsi/src/core/mod.rs delete mode 100644 libpsi/src/lib.rs delete mode 100644 libpsi/src/maths/complex.rs delete mode 100644 libpsi/src/maths/complex_ops.rs delete mode 100644 libpsi/src/maths/matrix.rs delete mode 100644 libpsi/src/maths/matrix_ops.rs delete mode 100644 libpsi/src/maths/mod.rs delete mode 100644 libpsi/src/maths/numeric_float.rs delete mode 100644 libpsi/src/maths/numeric_int.rs delete mode 100644 libpsi/src/maths/numeric_traits.rs delete mode 100644 libpsi/src/maths/vector.rs delete mode 100644 libpsi/src/maths/vector_ops.rs 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>; +pub type QuantumGate = Matrix>; + +#[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, +} + +pub struct QuantumRegister { + state: ColumnVector>, + qubits: Vec, +} + +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) { + self.bits = bits; + } +} + +impl QuantumRegister { + fn update(&mut self) { + let matrices: Vec>> = + 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> = self.state.mul_matrix(gate).unwrap(); + self.qubits[index] = result; + self.update(); + } +} + +impl ops::Index for QuantumRegister { + type Output = QuantumBit; + + fn index(&self, index: usize) -> &Self::Output { + &self.qubits[index] + } +} + +impl ops::IndexMut for QuantumRegister { + fn index_mut(&mut self, index: usize) -> &mut Self::Output { + &mut self.qubits[index] + } +} + +impl ops::Index for ClassicalRegister { + type Output = i32; + + fn index(&self, index: usize) -> &Self::Output { + &self.bits[index] + } +} + +impl ops::IndexMut 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 { + pub real: T, + pub imaginary: T, +} + +impl fmt::Debug for Complex { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!( + f, + "Complex {{ real: {:?}, imaginary: {:?} }}", + self.real, self.imaginary + ) + } +} + +impl fmt::Display for Complex { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!(f, "{} + {}i", self.real, self.imaginary) + } +} + +impl Numeric for Complex { + fn zero() -> Self { + Complex::new(0.0, 0.0) + } + + fn one() -> Self { + Complex::new(1.0, 0.0) + } +} + +impl Numeric for Complex { + fn zero() -> Self { + Complex::new(0.0, 0.0) + } + + fn one() -> Self { + Complex::new(1.0, 0.0) + } +} + +impl Complex { + pub fn new(real: T, imaginary: T) -> Complex { + Complex { real, imaginary } + } + + pub fn get_conjugate(&self) -> Complex { + 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 ops::Neg for Complex { + type Output = Complex; + + fn neg(self) -> Complex { + Complex { + real: -self.real, + imaginary: -self.imaginary, + } + } +} + +impl From for Complex { + fn from(real: T) -> Complex { + Complex { + real, + imaginary: T::zero(), + } + } +} + +// Complex-Complex +impl ops::Add for Complex { + type Output = Complex; + + fn add(self, other: Complex) -> Complex { + Complex { + real: self.real - other.real, + imaginary: self.imaginary - other.imaginary, + } + } +} + +impl ops::Sub for Complex { + type Output = Complex; + + fn sub(self, other: Complex) -> Complex { + Complex { + real: self.real - other.real, + imaginary: self.imaginary - other.imaginary, + } + } +} + +impl ops::Mul for Complex { + type Output = Complex; + + fn mul(self, other: Complex) -> Complex { + Complex { + real: self.real * other.real - self.imaginary * other.imaginary, + imaginary: self.real * other.imaginary + self.imaginary * other.real, + } + } +} + +impl ops::Div for Complex { + type Output = Complex; + + fn div(self, other: Complex) -> Complex { + 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 ops::AddAssign for Complex { + fn add_assign(&mut self, other: Complex) { + self.real += other.real; + self.imaginary += other.imaginary; + } +} + +impl ops::SubAssign for Complex { + fn sub_assign(&mut self, other: Complex) { + self.real -= other.real; + self.imaginary -= other.imaginary; + } +} + +impl ops::MulAssign for Complex { + fn mul_assign(&mut self, other: Complex) { + 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 ops::DivAssign for Complex { + fn div_assign(&mut self, other: Complex) { + 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 ops::Add for Complex { + type Output = Complex; + + fn add(self, other: T) -> Complex { + Complex { + real: self.real + other, + imaginary: self.imaginary, + } + } +} + +impl ops::Sub for Complex { + type Output = Complex; + + fn sub(self, other: T) -> Complex { + Complex { + real: self.real - other, + imaginary: self.imaginary, + } + } +} + +impl ops::Mul for Complex { + type Output = Complex; + + fn mul(self, other: T) -> Complex { + Complex { + real: self.real * other, + imaginary: self.imaginary * other, + } + } +} + +impl ops::Div for Complex { + type Output = Complex; + + fn div(self, other: T) -> Complex { + 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 { + pub data: Vec, + pub rows: usize, + pub cols: usize, +} + +impl Matrix { + pub fn new(rows: usize, cols: usize, data: Vec) -> 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> { + 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 { + 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 { + 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> { + 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> { + 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 { + 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 ops::Index<(usize, usize)> for Matrix { + type Output = T; + + fn index(&self, index: (usize, usize)) -> &Self::Output { + &self.data[index.0 * self.cols + index.1] + } +} + +impl ops::IndexMut<(usize, usize)> for Matrix { + fn index_mut(&mut self, index: (usize, usize)) -> &mut Self::Output { + &mut self.data[index.0 * self.cols + index.1] + } +} + +impl fmt::Debug for Matrix { + 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 fmt::Display for Matrix { + 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 ops::Add<&Matrix> for Matrix { + type Output = Option>; + + fn add(self, other: &Matrix) -> Self::Output { + self.add_to(other) + } +} + +impl ops::Sub<&Matrix> for Matrix { + type Output = Option>; + + fn sub(self, other: &Matrix) -> Self::Output { + self.subtract(other) + } +} + +impl ops::Mul for Matrix { + type Output = Matrix; + + fn mul(self, scalar: T) -> Self::Output { + self.scale(scalar) + } +} + +impl ops::Div for Matrix { + type Output = Matrix; + + fn div(self, scalar: T) -> Self::Output { + self.scale(T::one() / scalar) + } +} + +impl ops::AddAssign<&Matrix> for Matrix { + fn add_assign(&mut self, other: &Matrix) { + if let Some(result) = self.add_to(other) { + *self = result; + } + } +} + +impl ops::SubAssign<&Matrix> for Matrix { + fn sub_assign(&mut self, other: &Matrix) { + if let Some(result) = self.subtract(other) { + *self = result; + } + } +} + +impl ops::MulAssign for Matrix { + fn mul_assign(&mut self, scalar: T) { + *self = self.scale(scalar); + } +} + +impl ops::DivAssign for Matrix { + 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 { + 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 { + 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 + + ops::Mul + + ops::Sub + + ops::Div + + ops::Neg + + 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 { + fn new(data: Vec) -> 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 { + fn from_matrix(matrix: &Matrix) -> Self; + fn to_matrix(&self) -> Matrix; +} + +#[derive(Clone)] +pub struct VectorImpl(Vec); +pub type RowVector = VectorImpl; +pub type ColumnVector = VectorImpl; + +impl ColumnVector { + pub fn mul_matrix(&self, matrix: &Matrix) -> Option> { + 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 { + RowVector::new(self.0.clone()) + } +} + +impl RowVector { + pub fn mul_matrix(&self, matrix: &Matrix) -> Option> { + 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 { + ColumnVector::new(self.0.clone()) + } +} + +impl VectorMatrix for RowVector { + fn to_matrix(&self) -> Matrix { + Matrix::new(1, self.size(), self.0.clone()) + } + + fn from_matrix(matrix: &Matrix) -> Self { + Self::new(matrix.data.clone()) + } +} + +impl VectorMatrix for ColumnVector { + fn to_matrix(&self) -> Matrix { + Matrix::new(self.size(), 1, self.0.clone()) + } + + fn from_matrix(matrix: &Matrix) -> Self { + Self::new(matrix.data.clone()) + } +} + +impl Vector for VectorImpl { + fn new(data: Vec) -> 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 VectorImpl { + pub fn add_to(&self, other: &Self) -> Option> { + 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> { + 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 { + 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 ops::Index + for VectorImpl +{ + type Output = T; + + fn index(&self, index: usize) -> &Self::Output { + &self.0[index] + } +} + +impl ops::IndexMut + for VectorImpl +{ + fn index_mut(&mut self, index: usize) -> &mut Self::Output { + &mut self.0[index] + } +} + +impl fmt::Debug for RowVector { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!(f, "RowVector({:?})", self.0) + } +} + +impl fmt::Debug for ColumnVector { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!(f, "ColumnVector({:?})", self.0) + } +} + +impl fmt::Display for RowVector { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!( + f, + "[{}]", + self.0 + .iter() + .map(|x| x.to_string()) + .collect::>() + .join(", ") + ) + } +} + +impl fmt::Display for ColumnVector { + 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 ops::Add<&VectorImpl> + for VectorImpl +{ + type Output = Option>; + + fn add(self, other: &VectorImpl) -> Self::Output { + self.add_to(other) + } +} + +impl ops::Sub<&VectorImpl> + for VectorImpl +{ + type Output = Option>; + + fn sub(self, other: &VectorImpl) -> Self::Output { + self.subtract(other) + } +} + +impl ops::Mul for VectorImpl { + type Output = VectorImpl; + + fn mul(self, scalar: T) -> Self::Output { + self.scale(scalar) + } +} + +impl ops::Div for VectorImpl { + type Output = VectorImpl; + + fn div(self, scalar: T) -> Self::Output { + self.scale(T::one() / scalar) + } +} + +impl ops::AddAssign> + for VectorImpl +{ + fn add_assign(&mut self, other: VectorImpl) { + if let Some(result) = self.add_to(&other) { + *self = result; + } + } +} + +impl ops::SubAssign> + for VectorImpl +{ + fn sub_assign(&mut self, other: VectorImpl) { + if let Some(result) = self.subtract(&other) { + *self = result; + } + } +} + +impl ops::MulAssign + for VectorImpl +{ + fn mul_assign(&mut self, scalar: T) { + *self = self.scale(scalar); + } +} + +impl ops::DivAssign + for VectorImpl +{ + fn div_assign(&mut self, scalar: T) { + *self = self.scale(T::one() / scalar); + } +} + +impl ops::Mul<&Matrix> for RowVector { + type Output = Option>; + + fn mul(self, matrix: &Matrix) -> Self::Output { + self.mul_matrix(matrix) + } +} + +impl ops::Mul<&Matrix> for ColumnVector { + type Output = Option>; + + fn mul(self, matrix: &Matrix) -> Self::Output { + self.mul_matrix(matrix) + } +} + +impl ops::MulAssign<&Matrix> for RowVector { + fn mul_assign(&mut self, matrix: &Matrix) { + if let Some(result) = self.mul_matrix(matrix) { + *self = result; + } + } +} + +impl ops::MulAssign<&Matrix> for ColumnVector { + fn mul_assign(&mut self, matrix: &Matrix) { + if let Some(result) = self.mul_matrix(matrix) { + *self = result; + } + } +} diff --git a/libpsi/Cargo.toml b/libpsi/Cargo.toml deleted file mode 100644 index 96b1881..0000000 --- a/libpsi/Cargo.toml +++ /dev/null @@ -1,8 +0,0 @@ -[package] -name = "libpsi" -version = "0.1.0" -edition = "2021" - -[dependencies] -lazy_static = "1.5.0" -libm = "0.2.8" diff --git a/libpsi/src/core/component.rs b/libpsi/src/core/component.rs deleted file mode 100644 index 41b7cc7..0000000 --- a/libpsi/src/core/component.rs +++ /dev/null @@ -1,141 +0,0 @@ -use crate::{complex, ColumnVector, Complex, Matrix, Vector, VectorMatrix}; -use core::{fmt, ops}; - -pub type QuantumBit = ColumnVector>; -pub type QuantumGate = Matrix>; - -#[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, -} - -pub struct QuantumRegister { - state: ColumnVector>, - qubits: Vec, -} - -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) { - self.bits = bits; - } -} - -impl QuantumRegister { - fn update(&mut self) { - let matrices: Vec>> = - 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> = self.state.mul_matrix(gate).unwrap(); - self.qubits[index] = result; - self.update(); - } -} - -impl ops::Index for QuantumRegister { - type Output = QuantumBit; - - fn index(&self, index: usize) -> &Self::Output { - &self.qubits[index] - } -} - -impl ops::IndexMut for QuantumRegister { - fn index_mut(&mut self, index: usize) -> &mut Self::Output { - &mut self.qubits[index] - } -} - -impl ops::Index for ClassicalRegister { - type Output = i32; - - fn index(&self, index: usize) -> &Self::Output { - &self.bits[index] - } -} - -impl ops::IndexMut for ClassicalRegister { - fn index_mut(&mut self, index: usize) -> &mut Self::Output { - &mut self.bits[index] - } -} diff --git a/libpsi/src/core/gates.rs b/libpsi/src/core/gates.rs deleted file mode 100644 index e9077ad..0000000 --- a/libpsi/src/core/gates.rs +++ /dev/null @@ -1,22 +0,0 @@ -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/src/core/mod.rs b/libpsi/src/core/mod.rs deleted file mode 100644 index 374c555..0000000 --- a/libpsi/src/core/mod.rs +++ /dev/null @@ -1,4 +0,0 @@ -pub mod component; -pub mod gates; - -pub use gates::*; diff --git a/libpsi/src/lib.rs b/libpsi/src/lib.rs deleted file mode 100644 index af50663..0000000 --- a/libpsi/src/lib.rs +++ /dev/null @@ -1,10 +0,0 @@ -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/src/maths/complex.rs b/libpsi/src/maths/complex.rs deleted file mode 100644 index bb19510..0000000 --- a/libpsi/src/maths/complex.rs +++ /dev/null @@ -1,82 +0,0 @@ -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 { - pub real: T, - pub imaginary: T, -} - -impl fmt::Debug for Complex { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!( - f, - "Complex {{ real: {:?}, imaginary: {:?} }}", - self.real, self.imaginary - ) - } -} - -impl fmt::Display for Complex { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!(f, "{} + {}i", self.real, self.imaginary) - } -} - -impl Numeric for Complex { - fn zero() -> Self { - Complex::new(0.0, 0.0) - } - - fn one() -> Self { - Complex::new(1.0, 0.0) - } -} - -impl Numeric for Complex { - fn zero() -> Self { - Complex::new(0.0, 0.0) - } - - fn one() -> Self { - Complex::new(1.0, 0.0) - } -} - -impl Complex { - pub fn new(real: T, imaginary: T) -> Complex { - Complex { real, imaginary } - } - - pub fn get_conjugate(&self) -> Complex { - 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/src/maths/complex_ops.rs b/libpsi/src/maths/complex_ops.rs deleted file mode 100644 index bfcf920..0000000 --- a/libpsi/src/maths/complex_ops.rs +++ /dev/null @@ -1,147 +0,0 @@ -use super::{Complex, Float}; -use core::ops; - -impl ops::Neg for Complex { - type Output = Complex; - - fn neg(self) -> Complex { - Complex { - real: -self.real, - imaginary: -self.imaginary, - } - } -} - -impl From for Complex { - fn from(real: T) -> Complex { - Complex { - real, - imaginary: T::zero(), - } - } -} - -// Complex-Complex -impl ops::Add for Complex { - type Output = Complex; - - fn add(self, other: Complex) -> Complex { - Complex { - real: self.real - other.real, - imaginary: self.imaginary - other.imaginary, - } - } -} - -impl ops::Sub for Complex { - type Output = Complex; - - fn sub(self, other: Complex) -> Complex { - Complex { - real: self.real - other.real, - imaginary: self.imaginary - other.imaginary, - } - } -} - -impl ops::Mul for Complex { - type Output = Complex; - - fn mul(self, other: Complex) -> Complex { - Complex { - real: self.real * other.real - self.imaginary * other.imaginary, - imaginary: self.real * other.imaginary + self.imaginary * other.real, - } - } -} - -impl ops::Div for Complex { - type Output = Complex; - - fn div(self, other: Complex) -> Complex { - 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 ops::AddAssign for Complex { - fn add_assign(&mut self, other: Complex) { - self.real += other.real; - self.imaginary += other.imaginary; - } -} - -impl ops::SubAssign for Complex { - fn sub_assign(&mut self, other: Complex) { - self.real -= other.real; - self.imaginary -= other.imaginary; - } -} - -impl ops::MulAssign for Complex { - fn mul_assign(&mut self, other: Complex) { - 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 ops::DivAssign for Complex { - fn div_assign(&mut self, other: Complex) { - 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 ops::Add for Complex { - type Output = Complex; - - fn add(self, other: T) -> Complex { - Complex { - real: self.real + other, - imaginary: self.imaginary, - } - } -} - -impl ops::Sub for Complex { - type Output = Complex; - - fn sub(self, other: T) -> Complex { - Complex { - real: self.real - other, - imaginary: self.imaginary, - } - } -} - -impl ops::Mul for Complex { - type Output = Complex; - - fn mul(self, other: T) -> Complex { - Complex { - real: self.real * other, - imaginary: self.imaginary * other, - } - } -} - -impl ops::Div for Complex { - type Output = Complex; - - fn div(self, other: T) -> Complex { - Complex { - real: self.real / other, - imaginary: self.imaginary / other, - } - } -} diff --git a/libpsi/src/maths/matrix.rs b/libpsi/src/maths/matrix.rs deleted file mode 100644 index 975c92d..0000000 --- a/libpsi/src/maths/matrix.rs +++ /dev/null @@ -1,190 +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) - }}; -} - -#[derive(Clone)] -pub struct Matrix { - pub data: Vec, - pub rows: usize, - pub cols: usize, -} - -impl Matrix { - pub fn new(rows: usize, cols: usize, data: Vec) -> 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> { - 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 { - 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 { - 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> { - 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> { - 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 { - 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 ops::Index<(usize, usize)> for Matrix { - type Output = T; - - fn index(&self, index: (usize, usize)) -> &Self::Output { - &self.data[index.0 * self.cols + index.1] - } -} - -impl ops::IndexMut<(usize, usize)> for Matrix { - fn index_mut(&mut self, index: (usize, usize)) -> &mut Self::Output { - &mut self.data[index.0 * self.cols + index.1] - } -} - -impl fmt::Debug for Matrix { - 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 fmt::Display for Matrix { - 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/src/maths/matrix_ops.rs b/libpsi/src/maths/matrix_ops.rs deleted file mode 100644 index 15ab28b..0000000 --- a/libpsi/src/maths/matrix_ops.rs +++ /dev/null @@ -1,62 +0,0 @@ -use super::{Float, Matrix}; -use core::ops; - -impl ops::Add<&Matrix> for Matrix { - type Output = Option>; - - fn add(self, other: &Matrix) -> Self::Output { - self.add_to(other) - } -} - -impl ops::Sub<&Matrix> for Matrix { - type Output = Option>; - - fn sub(self, other: &Matrix) -> Self::Output { - self.subtract(other) - } -} - -impl ops::Mul for Matrix { - type Output = Matrix; - - fn mul(self, scalar: T) -> Self::Output { - self.scale(scalar) - } -} - -impl ops::Div for Matrix { - type Output = Matrix; - - fn div(self, scalar: T) -> Self::Output { - self.scale(T::one() / scalar) - } -} - -impl ops::AddAssign<&Matrix> for Matrix { - fn add_assign(&mut self, other: &Matrix) { - if let Some(result) = self.add_to(other) { - *self = result; - } - } -} - -impl ops::SubAssign<&Matrix> for Matrix { - fn sub_assign(&mut self, other: &Matrix) { - if let Some(result) = self.subtract(other) { - *self = result; - } - } -} - -impl ops::MulAssign for Matrix { - fn mul_assign(&mut self, scalar: T) { - *self = self.scale(scalar); - } -} - -impl ops::DivAssign for Matrix { - 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 deleted file mode 100644 index a2b9b84..0000000 --- a/libpsi/src/maths/mod.rs +++ /dev/null @@ -1,16 +0,0 @@ -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/src/maths/numeric_float.rs b/libpsi/src/maths/numeric_float.rs deleted file mode 100644 index 3591f24..0000000 --- a/libpsi/src/maths/numeric_float.rs +++ /dev/null @@ -1,65 +0,0 @@ -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 { - 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 { - 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/src/maths/numeric_int.rs b/libpsi/src/maths/numeric_int.rs deleted file mode 100644 index a9ecec7..0000000 --- a/libpsi/src/maths/numeric_int.rs +++ /dev/null @@ -1,44 +0,0 @@ -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 deleted file mode 100644 index dd7fdec..0000000 --- a/libpsi/src/maths/numeric_traits.rs +++ /dev/null @@ -1,24 +0,0 @@ -use core::ops; - -pub trait Numeric: - Copy - + PartialOrd - + ops::Add - + ops::Mul - + ops::Sub - + ops::Div - + ops::Neg - + 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 deleted file mode 100644 index 07ff7fb..0000000 --- a/libpsi/src/maths/vector.rs +++ /dev/null @@ -1,268 +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 { - fn new(data: Vec) -> 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 { - fn from_matrix(matrix: &Matrix) -> Self; - fn to_matrix(&self) -> Matrix; -} - -#[derive(Clone)] -pub struct VectorImpl(Vec); -pub type RowVector = VectorImpl; -pub type ColumnVector = VectorImpl; - -impl ColumnVector { - pub fn mul_matrix(&self, matrix: &Matrix) -> Option> { - 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 { - RowVector::new(self.0.clone()) - } -} - -impl RowVector { - pub fn mul_matrix(&self, matrix: &Matrix) -> Option> { - 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 { - ColumnVector::new(self.0.clone()) - } -} - -impl VectorMatrix for RowVector { - fn to_matrix(&self) -> Matrix { - Matrix::new(1, self.size(), self.0.clone()) - } - - fn from_matrix(matrix: &Matrix) -> Self { - Self::new(matrix.data.clone()) - } -} - -impl VectorMatrix for ColumnVector { - fn to_matrix(&self) -> Matrix { - Matrix::new(self.size(), 1, self.0.clone()) - } - - fn from_matrix(matrix: &Matrix) -> Self { - Self::new(matrix.data.clone()) - } -} - -impl Vector for VectorImpl { - fn new(data: Vec) -> 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 VectorImpl { - pub fn add_to(&self, other: &Self) -> Option> { - 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> { - 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 { - 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 ops::Index - for VectorImpl -{ - type Output = T; - - fn index(&self, index: usize) -> &Self::Output { - &self.0[index] - } -} - -impl ops::IndexMut - for VectorImpl -{ - fn index_mut(&mut self, index: usize) -> &mut Self::Output { - &mut self.0[index] - } -} - -impl fmt::Debug for RowVector { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!(f, "RowVector({:?})", self.0) - } -} - -impl fmt::Debug for ColumnVector { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!(f, "ColumnVector({:?})", self.0) - } -} - -impl fmt::Display for RowVector { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!( - f, - "[{}]", - self.0 - .iter() - .map(|x| x.to_string()) - .collect::>() - .join(", ") - ) - } -} - -impl fmt::Display for ColumnVector { - 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/src/maths/vector_ops.rs b/libpsi/src/maths/vector_ops.rs deleted file mode 100644 index 1e7148c..0000000 --- a/libpsi/src/maths/vector_ops.rs +++ /dev/null @@ -1,107 +0,0 @@ -use super::{Float, Matrix}; -use crate::{ColumnVector, RowVector, VectorImpl}; -use core::ops; - -impl ops::Add<&VectorImpl> - for VectorImpl -{ - type Output = Option>; - - fn add(self, other: &VectorImpl) -> Self::Output { - self.add_to(other) - } -} - -impl ops::Sub<&VectorImpl> - for VectorImpl -{ - type Output = Option>; - - fn sub(self, other: &VectorImpl) -> Self::Output { - self.subtract(other) - } -} - -impl ops::Mul for VectorImpl { - type Output = VectorImpl; - - fn mul(self, scalar: T) -> Self::Output { - self.scale(scalar) - } -} - -impl ops::Div for VectorImpl { - type Output = VectorImpl; - - fn div(self, scalar: T) -> Self::Output { - self.scale(T::one() / scalar) - } -} - -impl ops::AddAssign> - for VectorImpl -{ - fn add_assign(&mut self, other: VectorImpl) { - if let Some(result) = self.add_to(&other) { - *self = result; - } - } -} - -impl ops::SubAssign> - for VectorImpl -{ - fn sub_assign(&mut self, other: VectorImpl) { - if let Some(result) = self.subtract(&other) { - *self = result; - } - } -} - -impl ops::MulAssign - for VectorImpl -{ - fn mul_assign(&mut self, scalar: T) { - *self = self.scale(scalar); - } -} - -impl ops::DivAssign - for VectorImpl -{ - fn div_assign(&mut self, scalar: T) { - *self = self.scale(T::one() / scalar); - } -} - -impl ops::Mul<&Matrix> for RowVector { - type Output = Option>; - - fn mul(self, matrix: &Matrix) -> Self::Output { - self.mul_matrix(matrix) - } -} - -impl ops::Mul<&Matrix> for ColumnVector { - type Output = Option>; - - fn mul(self, matrix: &Matrix) -> Self::Output { - self.mul_matrix(matrix) - } -} - -impl ops::MulAssign<&Matrix> for RowVector { - fn mul_assign(&mut self, matrix: &Matrix) { - if let Some(result) = self.mul_matrix(matrix) { - *self = result; - } - } -} - -impl ops::MulAssign<&Matrix> for ColumnVector { - fn mul_assign(&mut self, matrix: &Matrix) { - if let Some(result) = self.mul_matrix(matrix) { - *self = result; - } - } -} -- cgit v1.3