From 6edb6a95b3b6a42be16bf7dde40204601c22b6cd Mon Sep 17 00:00:00 2001 From: Walter Oggioni Date: Thu, 28 Dec 2023 14:07:00 +0800 Subject: [PATCH] initial commit --- .gitignore | 1 + Cargo.lock | 186 +++++++ Cargo.toml | 28 + examples/benchmark.rs | 42 ++ examples/matrix_example.rs | 51 ++ src/err.rs | 50 ++ src/hludecomposition.rs | 215 ++++++++ src/hmatrix.rs | 1057 ++++++++++++++++++++++++++++++++++++ src/hpivot.rs | 87 +++ src/lib.rs | 97 ++++ src/matrix.rs | 489 +++++++++++++++++ src/misc.rs | 49 ++ src/rational.rs | 791 +++++++++++++++++++++++++++ src/sludecomposition.rs | 208 +++++++ src/smatrix.rs | 956 ++++++++++++++++++++++++++++++++ src/spivot.rs | 94 ++++ 16 files changed, 4401 insertions(+) create mode 100644 .gitignore create mode 100644 Cargo.lock create mode 100644 Cargo.toml create mode 100644 examples/benchmark.rs create mode 100644 examples/matrix_example.rs create mode 100644 src/err.rs create mode 100644 src/hludecomposition.rs create mode 100644 src/hmatrix.rs create mode 100644 src/hpivot.rs create mode 100644 src/lib.rs create mode 100644 src/matrix.rs create mode 100644 src/misc.rs create mode 100644 src/rational.rs create mode 100644 src/sludecomposition.rs create mode 100644 src/smatrix.rs create mode 100644 src/spivot.rs diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..ea8c4bf --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +/target diff --git a/Cargo.lock b/Cargo.lock new file mode 100644 index 0000000..63dcbb8 --- /dev/null +++ b/Cargo.lock @@ -0,0 +1,186 @@ +# This file is automatically @generated by Cargo. +# It is not intended for manual editing. +version = 3 + +[[package]] +name = "autocfg" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d468802bab17cbc0cc575e9b053f41e72aa36bfa6b7f55e3529ffa43161b97fa" + +[[package]] +name = "cfg-if" +version = "1.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "baf1de4339761588bc0619e3cbc0120ee582ebb74b53b4efbf79117bd2da40fd" + +[[package]] +name = "getrandom" +version = "0.2.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "190092ea657667030ac6a35e305e62fc4dd69fd98ac98631e5d3a2b1575a12b5" +dependencies = [ + "cfg-if", + "libc", + "wasi", +] + +[[package]] +name = "heck" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "95505c38b4572b2d910cecb0281560f54b440a19336cbbcb27bf6ce6adc6f5a8" + +[[package]] +name = "libc" +version = "0.2.152" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "13e3bf6590cbc649f4d1a3eefc9d5d6eb746f5200ffb04e5e142700b8faa56e7" + +[[package]] +name = "num-bigint" +version = "0.4.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "608e7659b5c3d7cba262d894801b9ec9d00de989e8a82bd4bef91d08da45cdc0" +dependencies = [ + "autocfg", + "num-integer", + "num-traits", +] + +[[package]] +name = "num-integer" +version = "0.1.45" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "225d3389fb3509a24c93f5c29eb6bde2586b98d9f016636dff58d7c6f7569cd9" +dependencies = [ + "autocfg", + "num-traits", +] + +[[package]] +name = "num-traits" +version = "0.2.17" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "39e3200413f237f41ab11ad6d161bc7239c84dcb631773ccd7de3dfe4b5c267c" +dependencies = [ + "autocfg", +] + +[[package]] +name = "opimps" +version = "0.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "857dabe64a7afe2e51ac9962dc3c008e74ae050dd47e21a7e7b1fc69a67a0229" +dependencies = [ + "proc-macro2", + "quote", + "syn", +] + +[[package]] +name = "ppv-lite86" +version = "0.2.17" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5b40af805b3121feab8a3c29f04d8ad262fa8e0561883e7653e024ae4479e6de" + +[[package]] +name = "proc-macro2" +version = "1.0.70" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "39278fbbf5fb4f646ce651690877f89d1c5811a3d4acb27700c1cb3cdb78fd3b" +dependencies = [ + "unicode-ident", +] + +[[package]] +name = "quote" +version = "1.0.33" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5267fca4496028628a95160fc423a33e8b2e6af8a5302579e322e4b520293cae" +dependencies = [ + "proc-macro2", +] + +[[package]] +name = "rand" +version = "0.8.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "34af8d1a0e25924bc5b7c43c079c942339d8f0a8b57c39049bef581b46327404" +dependencies = [ + "libc", + "rand_chacha", + "rand_core", +] + +[[package]] +name = "rand_chacha" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e6c10a63a0fa32252be49d21e7709d4d4baf8d231c2dbce1eaa8141b9b127d88" +dependencies = [ + "ppv-lite86", + "rand_core", +] + +[[package]] +name = "rand_core" +version = "0.6.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ec0be4795e2f6a28069bec0b5ff3e2ac9bafc99e6a9a7dc3547996c5c816922c" +dependencies = [ + "getrandom", +] + +[[package]] +name = "rmath" +version = "0.1.0" +dependencies = [ + "num-bigint", + "num-traits", + "opimps", + "rand", + "sealed", + "trait-group", +] + +[[package]] +name = "sealed" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f4a8caec23b7800fb97971a1c6ae365b6239aaeddfb934d6265f8505e795699d" +dependencies = [ + "heck", + "proc-macro2", + "quote", + "syn", +] + +[[package]] +name = "syn" +version = "2.0.39" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "23e78b90f2fcf45d3e842032ce32e3f2d1545ba6636271dcbf24fa306d87be7a" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "trait-group" +version = "0.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e1b362975c6f0f21a41fbb9ca91fe5dcb7e01e12331360374347476b45f5cb9c" + +[[package]] +name = "unicode-ident" +version = "1.0.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3354b9ac3fae1ff6755cb6db53683adb661634f67557942dea4facebec0fee4b" + +[[package]] +name = "wasi" +version = "0.11.0+wasi-snapshot-preview1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9c8d87e72b64a3b4db28d11ce29237c246188f4f51057d65a7eab63b7987e423" diff --git a/Cargo.toml b/Cargo.toml new file mode 100644 index 0000000..8f966c2 --- /dev/null +++ b/Cargo.toml @@ -0,0 +1,28 @@ +[package] +name = "rmath" +version = "0.1.0" +edition = "2021" +authors = ["Walter Oggioni "] +license = "MIT" +rust-version = "1.60" + +[dependencies] +trait-group = "0.1.0" +sealed = "0.5" +num-traits = "0.2.17" +opimps = "0.2.2" + +[dev-dependencies] +num-bigint = "0.4.4" +rand = "0.8.5" + +[lib] +name = "rmath" +crate-type = ["lib"] +bench = false + +[profile.release] +strip = true +lto = true +debug-assertions = false +codegen-units = 1 diff --git a/examples/benchmark.rs b/examples/benchmark.rs new file mode 100644 index 0000000..d250bd9 --- /dev/null +++ b/examples/benchmark.rs @@ -0,0 +1,42 @@ +use std::env; + +use num_bigint::BigInt; +use num_traits::One; +use num_traits::Zero; +use rmath::HMatrix; +use rmath::Rational; +use rmath::NumericalMatrix; +use rand::rngs::StdRng; +use rand::SeedableRng; +use rand::RngCore; + +fn main( ) { + let size = env::args().skip(1).map(|it| { + it.parse::().unwrap() + }).take(1).collect::>()[0]; + let mut rand = { + let seed = [ + 1,0,1,3, + 2,5,0,0, + 200,1,0,0, + 210,30,0,0, + 78,134,31,0, + 253,11,7,0, + 120,169,89,48, + 200,0,202,0 + ]; + StdRng::from_seed(seed) + }; + let mtx = HMatrix::>::new( + usize::try_from(size).unwrap(),usize::try_from(size).unwrap(), |_| { + Rational::new(BigInt::from(rand.next_u32() % (size * 20)) - size * 10, BigInt::one()) + }); + let b = HMatrix::>::new( + usize::try_from(size).unwrap(),usize::try_from(size).unwrap(), |_| { + Rational::new(BigInt::from(rand.next_u32() % (size * 10)) - (size * 5), BigInt::one()) + }); + let lu = mtx.clone().lu().unwrap(); + let x = lu.solve(&b).unwrap(); + + assert!((mtx * x - b).is_zero()); +} \ No newline at end of file diff --git a/examples/matrix_example.rs b/examples/matrix_example.rs new file mode 100644 index 0000000..d8b3df1 --- /dev/null +++ b/examples/matrix_example.rs @@ -0,0 +1,51 @@ +use rmath::HMatrix; +use rmath::SMatrix; + +fn main() { + let mut acc = 0; + let mut m1 = HMatrix::new(3, 3, |_| { + let res = acc; + acc += 1; + res + }); + println!("{}", m1); + println!("{}", &m1 + &m1); + println!("{}", m1.clone() - m1.clone()); + println!("{}", &m1 * &m1); + + acc = 0; + let mut m2: SMatrix = SMatrix::new(|_| { + let res = acc; + acc += 1; + res + }); + println!("{}", m2); + println!("{}", m2 + m2); + println!("{}", m2 - m2); + println!("{}", m2 * m2); + + for (i, j, v) in m2 { + println!("{} {} {}", i, j, v); + } + for (i, j, v) in &mut m2 { + println!("{} {} {}", i, j, v); + *v = 5; + } + println!("{}", m2); + println!("{}", m1 == m1); + println!("{}", m1 == &m1 + &m1); + println!("{}", m2 == m2); + println!("{}", m2 == m2 + m2); + + for (i, j, v) in &mut m1 { + println!("{} {} {}", i, j, v); + *v = 0; + } + + let a: SMatrix = + SMatrix::new(|(i, j)| i32::try_from(i).unwrap() - i32::try_from(j).unwrap()); + let b: SMatrix = + SMatrix::new(|(i, j)| i32::try_from(i).unwrap() - i32::try_from(j).unwrap()); + println!("{}", a * b); + println!("{}", b * a); +} diff --git a/src/err.rs b/src/err.rs new file mode 100644 index 0000000..5321d7a --- /dev/null +++ b/src/err.rs @@ -0,0 +1,50 @@ +use std::fmt::Debug; + +pub trait RmathError { + fn msg(&self) -> String; +} + +#[derive(Default)] +pub struct SingularMatrixError {} + +impl Debug for SingularMatrixError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + write!(f, "{}", self.msg()) + } +} + +impl RmathError for SingularMatrixError { + fn msg(&self) -> String { + String::from("Matrix is singular") + } +} + +pub struct SizeError { + msg: String, +} + +impl Debug for dyn RmathError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + write!(f, "{}", self.msg()) + } +} + +impl RmathError for SizeError { + fn msg(&self) -> String { + self.msg.clone() + } +} + +impl SizeError { + pub fn new(msg: &str) -> SizeError { + SizeError { + msg: String::from(msg), + } + } +} + +impl Debug for SizeError { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + write!(f, "{}", self.msg()) + } +} diff --git a/src/hludecomposition.rs b/src/hludecomposition.rs new file mode 100644 index 0000000..50a8801 --- /dev/null +++ b/src/hludecomposition.rs @@ -0,0 +1,215 @@ +use num_traits::One; +use num_traits::Zero; +use std::ops::Add; +use std::ops::Div; +use std::ops::DivAssign; +use std::ops::Mul; +use std::ops::MulAssign; +use std::ops::Neg; +use std::ops::SubAssign; + +use crate::err::RmathError; + +use super::hpivot::HPivot; +use super::matrix::LUDecompositionImpl; +use super::matrix::NumericalMatrixImpl; +use super::HMatrix; +use super::LUDecomposition; +use super::Matrix; +use super::NumericalMatrix; +use super::Pivot; +use super::SingularMatrixError; + +use super::misc::EnhancedOption; + +pub struct HLUDecomposition { + matrix: HMatrix, + pivot: HPivot, +} + +impl HLUDecomposition +where TYPE : Zero + SubAssign + Clone, + for<'a> &'a TYPE : Mul<&'a TYPE, Output = TYPE>, + for<'a> TYPE : DivAssign<&'a TYPE> { + pub (crate) fn new(matrix: HMatrix, pivot: Option) -> Self { + let size = matrix.rows(); + HLUDecomposition { + matrix, + pivot: pivot.otherwise(|| HPivot::new(size)), + } + } + + pub fn solve(&self, b : &HMatrix) -> Result, Box> { + let mut x = HMatrix::::new(b.rows(), b.columns(), |_| TYPE::zero()); + for n in 0..b.columns() { + for i in 0..self.matrix.rows() { + x[(i, n)] = b[(self.pivot[i], n)].clone(); + for k in 0..i { + let sub = &self.matrix[(i, k)] * &x[(k, n)]; + x[(i, n)] -= sub; + } + } + for i in (0..self.matrix.rows()).rev() { + for k in (i + 1)..self.matrix.rows() { + let sub = &self.matrix[(i, k)] * &x[(k, n)]; + x[(i, n)] -= sub; + } + if self.matrix[(i, i)].is_zero() { + return Err(Box::::default()); + } else { + x[(i, n)] /= &self.matrix[(i, i)]; + } + } + } + Ok(x) + } +} + +impl LUDecompositionImpl, HPivot> for HLUDecomposition +where + for<'a> TYPE: 'a + One + MulAssign<&'a TYPE> + Neg, +{ + fn get_matrix(&self) -> &HMatrix { + &self.matrix + } + + fn get_pivot(&self) -> &HPivot { + &self.pivot + } +} + +impl LUDecomposition, HPivot> for HLUDecomposition +where + for<'b> TYPE: 'b + + Clone + + Zero + + One + + PartialEq + + PartialOrd + + Add<&'b TYPE, Output = TYPE> + + Neg + + Div + + SubAssign + + SubAssign<&'b TYPE> + + DivAssign + + DivAssign<&'b TYPE> + + MulAssign<&'b TYPE>, + for<'c> &'c TYPE: + Mul + Add + Div + Neg, +{ + fn l(&self) -> HMatrix { + self.matrix.clone().tril_replace(TYPE::one()) + } + + fn u(&self) -> HMatrix { + self.matrix.clone().triu() + } + + fn invert(&self) -> HMatrix { + let mut result = + HMatrix::::new(self.matrix.rows(), self.matrix.columns(), |_| TYPE::zero()); + NumericalMatrixImpl::>::lu_invert( + &self.matrix, + &mut result, + &mut self.pivot.clone(), + ); + result + } + + fn pivot(&self, m: HMatrix) -> HMatrix { + &self.pivot * m + } + + // fn p<'a>(&'a self) -> &'a HPivot { + // &self.pivot + // } +} + + +#[cfg(test)] +mod test { + use rand::rngs::StdRng; + use rand::RngCore; + use rand::SeedableRng; + + use crate::NumericalMatrix; + use crate::Rational; + use crate::HMatrix; + use num_traits::Zero; + + use rand; + + #[test] + fn solve_linear_system() { + let mut rand = { + let seed = [ + 1,0,1,3, + 2,5,0,0, + 200,1,0,0, + 210,30,0,0, + 78,134,31,0, + 253,11,7,0, + 120,169,89,48, + 200,0,202,0 + ]; + StdRng::from_seed(seed) + }; + let mtx = HMatrix::>::new(5, 5, |_| { + Rational::new(i64::try_from(rand.next_u32() % 40).unwrap() - 20, 1) + }); + let b = HMatrix::>::new(5, 5, |_| { + Rational::new(i64::try_from(rand.next_u32() % 20).unwrap() - 10, 1) + }); + let lu = mtx.clone().lu().unwrap(); + let x = lu.solve(&b).unwrap(); + assert!((mtx * x - b).is_zero()); + } + +} + + +#[cfg(test)] +mod test_big_int { + use num_traits::Zero; + use rand::rngs::StdRng; + use rand::RngCore; + use rand::SeedableRng; + + use crate::NumericalMatrix; + use crate::Rational; + use crate::HMatrix; + + use rand; + + #[test] + fn solve_linear_system() { + use num_bigint::BigInt; + use num_traits::One; + + let mut rand = { + let seed = [ + 1,0,1,3, + 2,5,0,0, + 200,1,0,0, + 210,30,0,0, + 78,134,31,0, + 253,11,7,0, + 120,169,89,48, + 200,0,202,0 + ]; + StdRng::from_seed(seed) + }; + let mtx = HMatrix::>::new(10, 10, |_| { + Rational::new(BigInt::from(rand.next_u32() % 200) - 100, BigInt::one()) + }); + let b = HMatrix::>::new(10, 10, |_| { + Rational::new(BigInt::from(rand.next_u32() % 100) - 50, BigInt::one()) + }); + let lu = mtx.clone().lu().unwrap(); + let x = lu.solve(&b).unwrap(); + + assert!((mtx * x - b).is_zero()); + + } + +} \ No newline at end of file diff --git a/src/hmatrix.rs b/src/hmatrix.rs new file mode 100644 index 0000000..8c41c5c --- /dev/null +++ b/src/hmatrix.rs @@ -0,0 +1,1057 @@ +use num_traits::One; +use num_traits::Zero; +use std::fmt::Debug; +use std::fmt::Display; +use std::iter::Enumerate; +use std::ops::Add; +use std::ops::AddAssign; +use std::ops::Div; +use std::ops::DivAssign; +use std::ops::FnMut; +use std::ops::Index; +use std::ops::IndexMut; +use std::ops::Mul; +use std::ops::MulAssign; +use std::ops::Neg; +use std::ops::Sub; +use std::ops::SubAssign; + +use super::matrix::MatrixImpl; +use super::matrix::NumericalMatrixImpl; +use super::misc::debug_fmt; +use super::NumericalMatrix; +use crate::SingularMatrixError; + +use super::hludecomposition::HLUDecomposition; +use super::Matrix; +use super::Pivot; + +use super::hpivot::HPivot; + +pub struct HMatrix { + rows: usize, + columns: usize, + values: Vec, +} + +impl Index<(usize, usize)> for HMatrix { + fn index(&self, index: (usize, usize)) -> &TYPE { + &self.values[self.columns * index.0 + index.1] + } + + type Output = TYPE; +} + +impl IndexMut<(usize, usize)> for HMatrix { + fn index_mut(&mut self, index: (usize, usize)) -> &mut TYPE { + &mut self.values[self.columns * index.0 + index.1] + } +} + +impl Default for HMatrix { + fn default() -> Self { + HMatrix { + rows: 0, + columns: 0, + values: vec![], + } + } +} + +impl Display for HMatrix +where + TYPE: Display, +{ + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + let lim = self.rows() * self.columns() - 1; + writeln!(f, "[")?; + for i in 0..self.rows() { + write!(f, " ")?; + for j in 0..self.columns() { + let index = i * self.columns() + j; + write!(f, "{}", self[(i, j)])?; + if index < lim { + write!(f, ", ")?; + } + } + writeln!(f)?; + } + write!(f, "]") + } +} + +impl Clone for HMatrix +where + TYPE: Clone, +{ + fn clone(&self) -> Self { + HMatrix::new(self.rows(), self.columns(), |position| { + self[position].clone() + }) + } +} + +impl<'a, TYPE> MatrixImpl<'a, TYPE> for HMatrix +where + TYPE: 'a, +{ + type IteratorType = HMatrixIterator; + type RefIteratorType<'b> = HMatrixIteratorRef<'a, TYPE>; + type MutRefIteratorType<'b> = HMatrixIteratorMut<'a, TYPE>; +} + +impl<'a, TYPE> Matrix for HMatrix +where + TYPE: 'a, +{ + fn columns(&self) -> usize { + self.columns + } + + fn rows(&self) -> usize { + self.rows + } + + fn size(&self) -> (usize, usize) { + (self.rows(), self.columns()) + } + + fn transpose(self) -> Result { + let columns = self.columns(); + let rows = self.rows(); + let mut tmp = Vec::::new(); + for (_, _, value) in self { + tmp.push(value); + } + let mut it = tmp.into_iter().rev(); + Ok(Self::new(columns, rows, |_| it.next().unwrap())) + } +} + +impl Add for HMatrix +where + for<'a> TYPE: AddAssign<&'a TYPE>, +{ + fn add(self, rhs: Self) -> Self::Output { + if self.size() != rhs.size() { + panic!( + "Trying to add matrix with size {:?} to matrix with size {:?}, {}", + self.size(), + rhs.size(), + "matrixes need to have the same size in order to perform addition" + ); + } + let mut result = HMatrix { + rows: self.rows(), + columns: self.columns(), + values: self.values, + }; + for (i, j, v) in &mut result { + *v += &rhs[(i, j)] + } + result + } + + type Output = Self; +} + +impl Add for &HMatrix +where + for<'a> &'a TYPE: Add<&'a TYPE, Output = TYPE>, +{ + fn add(self, rhs: Self) -> Self::Output { + if self.size() != rhs.size() { + panic!( + "Trying to add matrix with size {:?} to matrix with size {:?}, {}", + self.size(), + rhs.size(), + "matrixes need to have the same size in order to perform addition" + ); + } + HMatrix::new(self.rows(), self.columns(), |position| { + &self[position] + &rhs[position] + }) + } + + type Output = HMatrix; +} + +impl Sub for HMatrix +where + for<'a> TYPE: SubAssign<&'a TYPE>, +{ + fn sub(self, rhs: Self) -> Self::Output { + if self.size() != rhs.size() { + panic!( + "Trying to subtract matrix with size {:?} to matrix with size {:?}, {}", + self.size(), + rhs.size(), + "matrixes need to have the same size in order to perform subtraction" + ); + } + let mut result = HMatrix { + rows: self.rows(), + columns: self.columns(), + values: self.values, + }; + for (i, j, v) in &mut result { + *v -= &rhs[(i, j)] + } + result + } + + type Output = Self; +} + +impl Sub for &HMatrix +where + TYPE: Sub + Clone, + for<'a> &'a TYPE: Sub<&'a TYPE, Output = TYPE>, +{ + fn sub(self, rhs: Self) -> Self::Output { + if self.size() != rhs.size() { + panic!( + "Trying to subtract matrix with size {:?} to matrix with size {:?}, {}", + self.size(), + rhs.size(), + "matrixes need to have the same size in order to perform subtraction" + ); + } + HMatrix::new(self.rows(), self.columns(), |position| { + &self[position] - &rhs[position] + }) + } + + type Output = HMatrix; +} + +#[opimps::impl_ops(Mul)] +fn mul(self: HMatrix, other: HMatrix) -> HMatrix +where + TYPE: Zero + AddAssign, + for<'a> &'a TYPE: Mul<&'a TYPE, Output = TYPE>, +{ + if self.columns() != other.rows() { + panic!( + "Trying to multiply matrix with size {:?} to matrix with size {:?}, {} {}", + self.size(), + other.size(), + "the number of columns of the first matrix need be equal to the number of rows", + "of the second matrix in order to perform matrix multiplication" + ); + } + let mut result: HMatrix = HMatrix::new(self.rows(), other.columns(), |_| TYPE::zero()); + for i in 0..result.rows() { + for j in 0..result.columns() { + for k in 0..self.columns() { + result[(i, j)] += &self[(i, k)] * &other[(k, j)] + } + } + } + result +} + +impl HMatrix { + pub fn new(rows: usize, columns: usize, mut cb: INITIALIZER) -> Self + where + INITIALIZER: FnMut((usize, usize)) -> TYPE, + { + let mut values = Vec::with_capacity(rows * columns); + for i in 0..rows { + for j in 0..columns { + values.push(cb((i, j))); + } + } + HMatrix { + rows, + columns, + values, + } + } +} + +// impl HMatrix where TYPE : One + Zero { + +// pub fn identity(size : usize) -> HMatrix { +// HMatrix::new(size, size, |(i, j)| if i == j { +// TYPE::one() +// } else { +// TYPE::zero() +// }) +// } +// } + +pub struct HMatrixIterator { + it: Enumerate< as IntoIterator>::IntoIter>, + columns: usize, +} + +impl Iterator for HMatrixIterator { + fn next(&mut self) -> Option { + self.it + .next() + .map(|(index, value)| (index / self.columns, index % self.columns, value)) + } + + type Item = (usize, usize, TYPE); +} + +pub struct HMatrixIteratorRef<'a, TYPE> { + it: Enumerate<<&'a Vec as IntoIterator>::IntoIter>, + columns: usize, +} + +impl<'a, TYPE> Iterator for HMatrixIteratorRef<'a, TYPE> { + fn next(&mut self) -> Option { + self.it + .next() + .map(|(index, value)| (index / self.columns, index % self.columns, value)) + } + + type Item = (usize, usize, &'a TYPE); +} + +pub struct HMatrixIteratorMut<'a, TYPE> { + it: Enumerate<<&'a mut Vec as IntoIterator>::IntoIter>, + columns: usize, +} + +impl<'a, TYPE> Iterator for HMatrixIteratorMut<'a, TYPE> { + fn next(&mut self) -> Option { + self.it + .next() + .map(|(index, value)| (index / self.columns, index % self.columns, value)) + } + + type Item = (usize, usize, &'a mut TYPE); +} + +impl IntoIterator for HMatrix { + type Item = (usize, usize, TYPE); + type IntoIter = HMatrixIterator; + + fn into_iter(self) -> Self::IntoIter { + let columns = self.columns(); + HMatrixIterator { + it: self.values.into_iter().enumerate(), + columns, + } + } +} + +impl<'a, TYPE> IntoIterator for &'a HMatrix { + type Item = (usize, usize, &'a TYPE); + type IntoIter = HMatrixIteratorRef<'a, TYPE>; + + fn into_iter(self) -> Self::IntoIter { + HMatrixIteratorRef { + it: self.values.iter().enumerate(), + columns: self.columns(), + } + } +} + +impl<'a, TYPE> IntoIterator for &'a mut HMatrix { + type Item = (usize, usize, &'a mut TYPE); + type IntoIter = HMatrixIteratorMut<'a, TYPE>; + + fn into_iter(self) -> Self::IntoIter { + let columns = self.columns(); + HMatrixIteratorMut { + it: self.values.iter_mut().enumerate(), + columns, + } + } +} + +impl PartialEq for HMatrix +where + TYPE: PartialEq, +{ + fn eq(&self, other: &Self) -> bool { + self.rows() == other.rows() + && self.columns() == other.columns() + && (0..self.rows()) + .flat_map(|i| (0..self.columns()).map(move |j| (i, j))) + .all(|index| self[index] == other[index]) + } +} + +impl Debug for HMatrix +where + TYPE: Debug, +{ + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + debug_fmt(self, f) + } +} + +impl Zero for HMatrix +where + TYPE: Zero, + for<'a> TYPE: AddAssign<&'a TYPE>, +{ + fn zero() -> Self { + HMatrix::new(1, 1, |_| TYPE::zero()) + } + + fn is_zero(&self) -> bool { + self.into_iter().all(|(_, _, value)| value.is_zero()) + } + + fn set_zero(&mut self) { + for (_, _, value) in self { + *value = TYPE::zero(); + } + } +} + +impl One for HMatrix +where + TYPE: One + Zero + PartialEq + AddAssign, + for<'a> &'a TYPE: Mul<&'a TYPE, Output = TYPE>, +{ + fn one() -> Self { + HMatrix::new(1, 1, |_| TYPE::one()) + } + + fn is_one(&self) -> bool { + self.into_iter().all(|(i, j, value)| { + if i == j { + value.is_one() + } else { + value.is_zero() + } + }) + } + + fn set_one(&mut self) { + for (i, j, value) in self { + if i == j { + *value = TYPE::one(); + } else { + *value = TYPE::zero(); + } + } + } +} + +impl<'a, TYPE> NumericalMatrixImpl<'a, TYPE, HPivot, HLUDecomposition> for HMatrix +where + for<'b> TYPE: Zero + + One + + Clone + + PartialEq + + PartialOrd + + SubAssign + + SubAssign<&'b TYPE> + + Div + + MulAssign<&'b TYPE> + + DivAssign + + DivAssign<&'b TYPE> + + Neg + + Add<&'b TYPE, Output = TYPE> + + 'b, + for<'c> &'c TYPE: Mul<&'c TYPE, Output = TYPE> + + Add<&'c TYPE, Output = TYPE> + + Div<&'c TYPE, Output = TYPE> + + Neg, + TYPE: 'a, + Self: Sized + + Index<(usize, usize), Output = TYPE> + + IndexMut<(usize, usize), Output = TYPE> + + IntoIterator + + 'a, + &'a Self: IntoIterator>, + &'a mut Self: + IntoIterator>, + Self::IteratorType: Iterator, + Self::RefIteratorType<'a>: Iterator, + Self::MutRefIteratorType<'a>: Iterator, +{ + fn identity(size: usize) -> HMatrix { + HMatrix::new( + size, + size, + |(i, j)| if i == j { TYPE::one() } else { TYPE::zero() }, + ) + } +} + +impl NumericalMatrix> for HMatrix +where + for<'b> TYPE: Clone + + 'b + + Zero + + One + + PartialEq + + PartialOrd + + Neg + + Add<&'b TYPE, Output = TYPE> + + Div + + SubAssign + + SubAssign<&'b TYPE> + + DivAssign + + DivAssign<&'b TYPE> + + MulAssign<&'b TYPE>, + for<'c> &'c TYPE: + Mul + Add + Div + Neg, +{ + fn identity(size: usize) -> Self { + NumericalMatrixImpl::>::identity(size) + } + + fn tril(mut self) -> Self { + self.tril_impl(None); + self + } + + fn triu(mut self) -> Self { + self.triu_impl(None); + self + } + + fn tril_replace(mut self, diag_replacement: TYPE) -> Self { + self.tril_impl(Some(diag_replacement)); + self + } + + fn triu_replace(mut self, diag_replacement: TYPE) -> Self { + self.triu_impl(Some(diag_replacement)); + self + } + + fn gauss_jordan_high(mut self) -> Self { + let mut pivot = HPivot::new(self.rows()); + NumericalMatrixImpl::>::gauss_jordan_high( + &mut self, &mut pivot, + ); + self + } + + fn gauss_jordan_low(mut self) -> Self { + let mut pivot = HPivot::new(self.rows()); + NumericalMatrixImpl::>::gauss_jordan_low( + &mut self, &mut pivot, + ); + self + } + + fn det(self) -> TYPE { + let pivot = HPivot::new(self.rows()); + NumericalMatrixImpl::>::det(self, pivot) + } + + fn invert(&mut self) -> Self { + let mut result = + NumericalMatrixImpl::>::identity(self.rows()); + let mut pivot = HPivot::new(self.rows()); + NumericalMatrixImpl::>::invert( + self, + &mut pivot, + &mut result, + ); + result + } + + fn lup(mut self) -> Result, SingularMatrixError> { + let mut pivot = HPivot::new(self.columns()); + NumericalMatrixImpl::>::lup(&mut self, &mut pivot)?; + Ok(HLUDecomposition::::new(self, Some(pivot))) + } + + fn lu(mut self) -> Result, SingularMatrixError> { + NumericalMatrixImpl::>::lu(&mut self)?; + Ok(HLUDecomposition::::new(self, None)) + } +} + +#[cfg(test)] +mod tests { + + use crate::HMatrix; + use crate::LUDecomposition; + use crate::Matrix; + use crate::NumericalMatrix; + use crate::Rational; + + use num_traits::One; + use num_traits::Zero; + + #[test] + fn test_new() { + let mut idx = 0; + let m1 = HMatrix::::new(3, 5, |_| { + let res = idx; + idx += 1; + res + }); + idx = 0; + for (_, _, value) in m1 { + assert_eq!(idx, value); + idx += 1; + } + } + + #[test] + fn test_eq() { + let mut idx = 0; + let m1 = HMatrix::::new(3, 5, |_| { + let res = idx; + idx += 1; + res + }); + idx = 0; + for (_, _, value) in &m1 { + assert_eq!(idx, *value); + idx += 1; + } + let m2 = m1.clone(); + assert!(m1 == m2); + } + + #[test] + fn test_sum() { + let mut idx = 0; + let m1 = HMatrix::::new(3, 5, |_| { + let res = idx; + idx += 1; + res + }); + let m2 = HMatrix::::new(3, 5, |_| { + let res = idx; + idx += 1; + res + }); + let m3 = &m1 + &m2; + for i in 0..m1.rows() { + for j in 0..m1.columns() { + assert_eq!(m3[(i, j)], m1[(i, j)] + m2[(i, j)]); + } + } + } + + #[test] + fn test_sub() { + let mut idx = 0; + let m1 = HMatrix::::new(3, 5, |_| { + let res = idx; + idx += 1; + res + }); + let m2 = HMatrix::::new(3, 5, |_| { + let res = idx; + idx += 1; + res + }); + let m3 = &m1 - &m2; + for i in 0..m1.rows() { + for j in 0..m1.columns() { + assert_eq!(m3[(i, j)], m1[(i, j)] - m2[(i, j)]); + } + } + } + + #[test] + fn test_mul() { + let mut idx = 0; + let m1 = HMatrix::::new(3, 5, |_| { + let res = idx; + idx += 1; + res + }); + idx = 0; + let m2 = HMatrix::::new(5, 2, |_| { + let res = idx; + idx += 1; + res + }); + let m3 = &m1 * &m2; + let expected = { + let mut idx = 0; + let values = [60, 70, 160, 195, 260, 320]; + HMatrix::::new(3, 2, |_| { + let res = values[idx]; + idx += 1; + res + }) + }; + assert!(expected == m3); + } + + #[test] + fn test_one() { + let mut idx = 0; + let m1 = HMatrix::::new(5, 5, |_| { + let res = idx; + idx += 1; + res + }); + let m2 = HMatrix::::identity(5); + let m3 = &m1 * &m2; + assert!(m1 == m3); + } + + #[test] + fn test_inv() { + let arr = [1, 2, 3, 4, 5, 6, 8, 7, 9]; + let mut idx = 0; + let m1 = HMatrix::>::new(3, 3, |_| { + let res = Rational::new(arr[idx], i64::one()); + idx += 1; + res + }); + let m2 = NumericalMatrix::invert(&mut m1.clone()); + let m3 = &m1 * &m2; + assert!(m3.is_one()); + } + + #[test] + fn test_det() { + let arr = [1, 2, 3, 4, 5, 6, 8, 7, 9]; + let mut idx = 0; + let m1 = HMatrix::>::new(3, 3, |_| { + let res = Rational::new(arr[idx], 1); + idx += 1; + res + }); + assert_eq!(Rational::new(-9, 1), NumericalMatrix::det(m1)); + } + + #[test] + fn test_tril() { + let arr = [1, 2, 3, 4, 5, 6, 8, 7, 9]; + let mut idx = 0; + let m1 = HMatrix::::new(3, 3, |_| { + let res = arr[idx]; + idx += 1; + res + }); + let diag_value = 42; + let lt = m1.tril_replace(diag_value); + for (i, j, value) in lt { + if i == j { + assert_eq!(value, diag_value); + } else if i < j { + assert_eq!(value, i64::zero()) + } else { + assert_eq!(value, i64::from(arr[i * 3 + j])) + } + } + } + + #[test] + fn test_triu() { + let arr = [1, 2, 3, 4, 5, 6, 8, 7, 9]; + let mut idx = 0; + let m1 = HMatrix::::new(3, 3, |_| { + let res = arr[idx]; + idx += 1; + res + }); + let diag_value = 42; + let lt = m1.triu_replace(diag_value); + for (i, j, value) in lt { + if i == j { + assert_eq!(value, diag_value); + } else if i > j { + assert_eq!(value, i64::zero()) + } else { + assert_eq!(value, i64::from(arr[i * 3 + j])) + } + } + } + + #[test] + fn test_lu_decomposition() { + let arr = [1, 2, 3, 4, 5, 6, 8, 7, 9]; + let mut idx = 0; + let orig = HMatrix::>::new(3, 3, |_| { + let res = Rational::new(arr[idx], i64::one()); + idx += 1; + res + }); + let m1 = orig.clone(); + let lu = m1.lup().unwrap(); + let l = lu.l(); + let u = lu.u(); + let m2 = lu.pivot(l * u); + assert_eq!(orig, m2); + } + + #[test] + fn test_lu_invert() { + let arr = [1, 2, 3, 4, 5, 6, 8, 7, 9]; + let mut idx = 0; + let orig = HMatrix::>::new(3, 3, |_| { + let res = Rational::new(arr[idx], i64::one()); + idx += 1; + res + }); + let m1 = orig.clone(); + let lu = m1.lup().unwrap(); + let m2 = lu.invert(); + let m3 = orig * m2; + assert!(m3.is_one()); + } + + #[test] + fn test_lu_det() { + let arr = [1, 2, 3, 4, 5, 6, 8, 7, 9]; + let mut idx = 0; + let m1 = HMatrix::>::new(3, 3, |_| { + let res = Rational::new(arr[idx], i32::one()); + idx += 1; + res + }); + let lu = m1.clone().lup().unwrap(); + assert_eq!(Rational::new(-9, i32::one()), lu.det()); + } +} + +#[cfg(test)] +mod big_int_tests { + use num_bigint::BigInt; + + use crate::HMatrix; + use crate::LUDecomposition; + use crate::Matrix; + use crate::NumericalMatrix; + use crate::Rational; + + use num_traits::One; + use num_traits::Zero; + + #[test] + fn test_new() { + let mut idx = 0; + let m1 = HMatrix::::new(3, 5, |_| { + let res = BigInt::from(idx); + idx += 1; + res + }); + idx = 0; + for (_, _, value) in m1 { + assert_eq!(BigInt::from(idx), value); + idx += 1; + } + } + + #[test] + fn test_eq() { + let mut idx = 0; + let m1 = HMatrix::::new(3, 5, |_| { + let res = BigInt::from(idx); + idx += 1; + res + }); + idx = 0; + let m2 = m1.clone(); + assert!(m1 == m2); + for (_, _, value) in &m1 { + assert_eq!(BigInt::from(idx), *value); + idx += 1; + } + } + + #[test] + fn test_sum() { + let mut idx = 0; + let m1 = HMatrix::::new(3, 5, |_| { + let res = BigInt::from(idx); + idx += 1; + res + }); + let m2 = HMatrix::::new(3, 5, |_| { + let res = BigInt::from(idx); + idx += 1; + res + }); + let m3 = &m1 + &m2; + for i in 0..m1.rows() { + for j in 0..m1.columns() { + assert_eq!(m3[(i, j)], &m1[(i, j)] + &m2[(i, j)]); + } + } + } + + #[test] + fn test_sub() { + let mut idx = 0; + let m1 = HMatrix::::new(3, 5, |_| { + let res = BigInt::from(idx); + idx += 1; + res + }); + let m2 = HMatrix::::new(3, 5, |_| { + let res = BigInt::from(idx); + idx += 1; + res + }); + let m3 = &m1 - &m2; + for i in 0..m1.rows() { + for j in 0..m1.columns() { + assert_eq!(m3[(i, j)], &m1[(i, j)] - &m2[(i, j)]); + } + } + } + + #[test] + fn test_mul() { + let mut idx = 0; + let m1 = HMatrix::::new(3, 5, |_| { + let res = BigInt::from(idx); + idx += 1; + res + }); + idx = 0; + let m2 = HMatrix::::new(5, 2, |_| { + let res = BigInt::from(idx); + idx += 1; + res + }); + let m3 = m1 * m2; + let expected = { + let mut idx = 0; + let values = [60, 70, 160, 195, 260, 320]; + HMatrix::::new(3, 2, |_| { + let res = BigInt::from(values[idx]); + idx += 1; + res + }) + }; + assert!(expected == m3); + } + + #[test] + fn test_one() { + let mut idx = 0; + let m1 = HMatrix::::new(5, 5, |_| { + let res = BigInt::from(idx); + idx += 1; + res + }); + let m2 = HMatrix::::identity(5); + let m3 = &m1 * &m2; + assert!(m1 == m3); + } + + #[test] + fn test_inv() { + let arr = [1, 2, 3, 4, 5, 6, 8, 7, 9]; + let mut idx = 0; + let m1 = HMatrix::>::new(3, 3, |_| { + let res = Rational::new(BigInt::from(arr[idx]), BigInt::one()); + idx += 1; + res + }); + let m2 = NumericalMatrix::invert(&mut m1.clone()); + let m3 = m1 * m2; + assert!(m3.is_one()); + } + + #[test] + fn test_det() { + let arr = [1, 2, 3, 4, 5, 6, 8, 7, 9]; + let mut idx = 0; + let m1 = HMatrix::>::new(3, 3, |_| { + let res = Rational::new(BigInt::from(arr[idx]), BigInt::one()); + idx += 1; + res + }); + assert_eq!( + Rational::new(BigInt::from(9), -BigInt::one()), + NumericalMatrix::det(m1) + ); + } + + #[test] + fn test_tril() { + let arr = [1, 2, 3, 4, 5, 6, 8, 7, 9]; + let mut idx = 0; + let m1 = HMatrix::::new(3, 3, |_| { + let res = BigInt::from(arr[idx]); + idx += 1; + res + }); + let diag_value = BigInt::from(42); + let lt = m1.tril_replace(diag_value.clone()); + for (i, j, value) in lt { + if i == j { + assert_eq!(value, diag_value); + } else if i < j { + assert_eq!(value, BigInt::zero()) + } else { + assert_eq!(value, BigInt::from(arr[i * 3 + j])) + } + } + } + + #[test] + fn test_triu() { + let arr = [1, 2, 3, 4, 5, 6, 8, 7, 9]; + let mut idx = 0; + let m1 = HMatrix::::new(3, 3, |_| { + let res = BigInt::from(arr[idx]); + idx += 1; + res + }); + let diag_value = BigInt::from(42); + let lt = m1.triu_replace(diag_value.clone()); + for (i, j, value) in lt { + if i == j { + assert_eq!(value, diag_value); + } else if i > j { + assert_eq!(value, BigInt::zero()) + } else { + assert_eq!(value, BigInt::from(arr[i * 3 + j])) + } + } + } + + #[test] + fn test_lu_decomposition() { + let arr = [1, 2, 3, 4, 5, 6, 8, 7, 9]; + let mut idx = 0; + let orig = HMatrix::>::new(3, 3, |_| { + let res = Rational::new(arr[idx], i64::one()); + idx += 1; + res + }); + let m1 = orig.clone(); + let lu = m1.lup().unwrap(); + let tril = lu.l(); + let triu = lu.u(); + let m2 = lu.pivot(tril * triu); + assert_eq!(orig, m2); + } + + #[test] + fn test_lu_invert() { + let arr = [1, 2, 3, 4, 5, 6, 8, 7, 9]; + let mut idx = 0; + let orig = HMatrix::>::new(3, 3, |_| { + let res = Rational::new(arr[idx], i64::one()); + idx += 1; + res + }); + let m1 = orig.clone(); + let lu = m1.lup().unwrap(); + let m2 = lu.invert(); + let m3 = orig * m2; + assert!(m3.is_one()); + } + + #[test] + fn test_lu_det() { + let arr = [1, 2, 3, 4, 5, 6, 8, 7, 9]; + let mut idx = 0; + let m1 = HMatrix::>::new(3, 3, |_| { + let res = Rational::new(BigInt::from(arr[idx]), BigInt::one()); + idx += 1; + res + }); + let lu = m1.clone().lup().unwrap(); + assert_eq!(Rational::new(BigInt::from(-9), BigInt::one()), lu.det()); + } +} diff --git a/src/hpivot.rs b/src/hpivot.rs new file mode 100644 index 0000000..eb0162b --- /dev/null +++ b/src/hpivot.rs @@ -0,0 +1,87 @@ +use std::ops::Index; +use std::ops::IndexMut; +use std::ops::Mul; + +use super::hmatrix::HMatrix; +use super::matrix::MatrixImpl; +use super::Matrix; +use super::Pivot; + +pub struct HPivot { + data: Vec, + permutations: usize, +} + +impl Index for HPivot { + fn index(&self, index: usize) -> &Self::Output { + &self.data[index] + } + + type Output = usize; +} + +impl IndexMut for HPivot { + fn index_mut(&mut self, index: usize) -> &mut >::Output { + &mut self.data[index] + } +} + +impl Pivot for HPivot { + fn swap(&mut self, i1: usize, i2: usize) { + self.data.swap(i1, i2); + self.permutations += 1; + } + + fn permutations(&self) -> usize { + self.permutations + } + + fn new(size: usize) -> Self { + HPivot { + data: (0..size).collect(), + permutations: 0, + } + } +} + +impl Mul> for HPivot +where + TYPE: Clone, +{ + fn mul(self, matrix: HMatrix) -> Self::Output { + &self * matrix + } + + type Output = HMatrix; +} + +impl Mul> for &HPivot +where + TYPE: Clone, +{ + fn mul(self, matrix: HMatrix) -> Self::Output { + if self.data.len() != matrix.rows() { + panic!("Pivot len must equal matrix rows") + } + let mut result = matrix.clone(); + let mut pclone = self.clone(); + for i in 0..matrix.rows() { + while i != pclone[i] { + let j = pclone[i]; + result.swap_rows_pivot(&mut pclone, i, j); + } + } + result + } + + type Output = HMatrix; +} + +impl Clone for HPivot { + fn clone(&self) -> Self { + HPivot { + data: self.data.clone(), + permutations: self.permutations, + } + } +} diff --git a/src/lib.rs b/src/lib.rs new file mode 100644 index 0000000..a048fff --- /dev/null +++ b/src/lib.rs @@ -0,0 +1,97 @@ +#[macro_use] +extern crate trait_group; + +use num_traits::One; +use num_traits::Zero; +use std::ops::Index; +use std::ops::IndexMut; +use std::ops::MulAssign; +use std::ops::Neg; + +mod rational; +pub use self::rational::Rational; + +mod matrix; +pub use self::matrix::LUDecomposition; + +mod hpivot; +pub use hpivot::HPivot; + +mod hludecomposition; +pub use hludecomposition::HLUDecomposition; + +mod hmatrix; +pub use self::hmatrix::HMatrix; + +mod spivot; +pub use spivot::SPivot; + +mod sludecomposition; +pub use sludecomposition::SLUDecomposition; + +mod smatrix; +pub use self::smatrix::SMatrix; + +mod err; +pub use self::err::SingularMatrixError; +pub use self::err::SizeError; + +mod misc; + +pub trait Pivot: Index + IndexMut + Sized { + fn swap(&mut self, i1: usize, i2: usize); + fn permutations(&self) -> usize; + fn new(size: usize) -> Self; +} + + +pub trait Vector : Index + IndexMut { + fn new TYPE>(size : usize, cb : F) -> Self; + fn length(&self) -> usize; +} + +pub trait NumericalVector : Vector { + fn norm() -> TYPE; +} + +pub trait Matrix +where + Self: Index<(usize, usize), Output = TYPE> + IndexMut<(usize, usize), Output = TYPE> + Sized, +{ + fn rows(&self) -> usize; + fn columns(&self) -> usize; + + fn size(&self) -> (usize, usize) { + (self.rows(), self.columns()) + } + + fn transpose(self) -> Result; +} + +pub trait NumericalMatrix>: + Matrix +where + for<'a> TYPE: Clone + One + Zero + MulAssign<&'a TYPE> + Neg, +{ + fn identity(size: usize) -> Self; + + fn tril(self) -> Self; + + fn triu(self) -> Self; + + fn tril_replace(self, diag_replacement: TYPE) -> Self; + + fn triu_replace(self, diag_replacement: TYPE) -> Self; + + fn gauss_jordan_high(self) -> Self; + + fn gauss_jordan_low(self) -> Self; + + fn det(self) -> TYPE; + + fn invert(&mut self) -> Self; + + fn lup(self) -> Result; + + fn lu(self) -> Result; +} diff --git a/src/matrix.rs b/src/matrix.rs new file mode 100644 index 0000000..e3e7213 --- /dev/null +++ b/src/matrix.rs @@ -0,0 +1,489 @@ +use super::Matrix; +use super::Pivot; +use super::SingularMatrixError; +use super::SizeError; +use num_traits::One; +use num_traits::Zero; +use std::ops::Add; +use std::ops::Div; +use std::ops::DivAssign; +use std::ops::Index; +use std::ops::IndexMut; +use std::ops::Mul; +use std::ops::MulAssign; +use std::ops::Neg; +use std::ops::SubAssign; + +pub trait MatrixImpl<'a, TYPE>: Matrix +where + TYPE: 'a, + Self: Index<(usize, usize), Output = TYPE> + + IndexMut<(usize, usize), Output = TYPE> + + IntoIterator + + 'a, + &'a Self: IntoIterator>, + &'a mut Self: + IntoIterator>, + Self::IteratorType: Iterator, + Self::RefIteratorType<'a>: Iterator, + Self::MutRefIteratorType<'a>: Iterator, +{ + fn swap(&mut self, idx1: (usize, usize), idx2: (usize, usize)) { + unsafe { + // Can't take two mutable loans from one matrix, so instead just cast + // them to their raw pointers to do the swap + let pa: *mut TYPE = &mut self[idx1]; + let pb: *mut TYPE = &mut self[idx2]; + std::ptr::swap(pa, pb); + } + } + + fn transpose_in_place(&mut self) -> Result<(), SizeError> + where + Self: Sized, + { + for i in 0..self.rows() { + for j in 0..self.columns() { + self.swap((i, j), (j, i)); + } + } + Ok(()) + } + + fn swap_rows(&mut self, row1: usize, row2: usize) { + for i in 0..self.columns() { + self.swap((row1, i), (row2, i)) + } + } + + fn swap_rows_pivot(&mut self, pivot: &mut P, row1: usize, row2: usize) { + self.swap_rows(row1, row2); + pivot.swap(row1, row2); + } + + fn swap_rows_mul_pivot( + &mut self, + other: &mut Self, + pivot: &mut P, + row1: usize, + row2: usize, + ) { + self.swap_rows_pivot(pivot, row1, row2); + other.swap_rows(row1, row2); + } + + fn swap_rows_mul(&mut self, other: &mut Self, row1: usize, row2: usize) { + self.swap_rows(row1, row2); + other.swap_rows(row1, row2); + } + + fn swap_columns(&mut self, col1: usize, col2: usize) { + for i in 0..self.rows() { + self.swap((i, col1), (i, col2)) + } + } + + type IteratorType; + type RefIteratorType<'b>; + type MutRefIteratorType<'b>; +} + +pub trait LUDecompositionImpl +where + MATRIX: Matrix, + for<'b> TYPE: One + MulAssign<&'b TYPE> + Neg, +{ + fn get_matrix(&self) -> &MATRIX; + + fn get_pivot(&self) -> &PIVOT; + + fn det(&self) -> TYPE { + let mut result = TYPE::one(); + for i in 0..self.get_matrix().rows() { + result *= &self.get_matrix()[(i, i)]; + } + if self.get_pivot().permutations() % 2 != 0 { + result = -result; + } + result + } +} + +pub trait LUDecomposition, PIVOT: Pivot>: + LUDecompositionImpl +where + for<'b> TYPE: One + MulAssign<&'b TYPE> + Neg, +{ + fn l(&self) -> MATRIX; + + fn u(&self) -> MATRIX; + + fn invert(&self) -> MATRIX; + + fn pivot(&self, m: MATRIX) -> MATRIX; + + fn det(&self) -> TYPE { + LUDecompositionImpl::det(self) + } +} + +pub trait NumericalMatrixImpl<'a, TYPE, PIVOT: Pivot, LU: LUDecomposition>: + MatrixImpl<'a, TYPE> +where + for<'b> TYPE: Zero + + One + + PartialEq + + Neg + + Div + + MulAssign<&'b TYPE>, + for<'c> &'c TYPE: Mul<&'c TYPE, Output = TYPE> + + Add<&'c TYPE, Output = TYPE> + + Div<&'c TYPE, Output = TYPE> + + Neg, + Self: Sized, + TYPE: 'a, + Self: Index<(usize, usize), Output = TYPE> + + IndexMut<(usize, usize), Output = TYPE> + + IntoIterator + + 'a, + &'a Self: IntoIterator>, + &'a mut Self: + IntoIterator>, + Self::IteratorType: Iterator, + Self::RefIteratorType<'a>: Iterator, + Self::MutRefIteratorType<'a>: Iterator, +{ + fn identity(size: usize) -> Self; + + fn add_row(&mut self, source_index: usize, dest_index: usize, factor: &TYPE) + where + TYPE: Mul + Add, + for<'f> &'f TYPE: + Mul + Add + Div + Neg, + { + for i in 0..self.columns() { + let tmp = &(self[(source_index, i)]) * factor; + self[(dest_index, i)] = &(self[(dest_index, i)]) + &tmp; + } + } + + fn add_row_mul( + &mut self, + other: &mut Self, + source_index: usize, + dest_index: usize, + factor: &TYPE, + ) where + TYPE: Add + Mul, + for<'f> &'f TYPE: + Mul + Add + Div + Neg, + { + self.add_row(source_index, dest_index, factor); + other.add_row(source_index, dest_index, factor); + } + + fn gauss_jordan_low(&mut self, pivot: &mut P) + where + TYPE: Zero + + PartialEq + + Add + + Mul + + Neg + + Div, + for<'f> &'f TYPE: + Mul + Add + Div + Neg, + { + for i in 0..self.rows() { + if self[(i, i)] == TYPE::zero() { + for j in (i + 1)..self.columns() { + if self[(j, i)] != TYPE::zero() { + self.swap_rows_pivot(pivot, i, j); + break; + } + } + } + for j in (i + 1)..self.rows() { + if self[(i, i)] != TYPE::zero() { + let neg = -(&self[(j, i)]); + let factor = &neg / &(self[(i, i)]); + self.add_row(i, j, &factor) + } + } + } + } + + fn gauss_jordan_low_mul(&mut self, other: &mut Self, pivot: &mut P) + where + TYPE: Zero + + PartialEq + + Add + + Mul + + Neg + + Div, + for<'f> &'f TYPE: + Mul + Add + Div + Neg, + { + for i in 0..self.rows() { + if self[(i, i)] == TYPE::zero() { + for j in (i + 1)..self.columns() { + if self[(j, i)] != TYPE::zero() { + self.swap_rows_mul_pivot(other, pivot, i, j); + break; + } + } + } + for j in (i + 1)..self.rows() { + if self[(i, i)] != TYPE::zero() { + let neg = -(&self[(j, i)]); + let factor = &neg / &(self[(i, i)]); + self.add_row_mul(other, i, j, &factor) + } + } + } + } + + fn gauss_jordan_high(&mut self, pivot: &mut P) + where + TYPE: Zero + + PartialEq + + Add + + Mul + + Neg + + Div, + for<'f> &'f TYPE: + Mul + Add + Div + Neg, + { + for i in (usize::zero()..self.rows()).rev() { + if self[(i, i)] == TYPE::zero() { + for j in (usize::zero()..i).rev() { + if self[(j, i)] != TYPE::zero() { + self.swap_rows_pivot(pivot, i, j); + break; + } + } + } + for j in (usize::zero()..i).rev() { + if self[(i, i)] != TYPE::zero() { + let neg = -(&self[(j, i)]); + let factor = &neg / &(self[(i, i)]); + self.add_row(i, j, &factor) + } + } + } + } + + fn gauss_jordan_high_mul(&mut self, other: &mut Self, pivot: &mut P) + where + TYPE: Zero + + PartialEq + + Add + + Mul + + Neg + + Div, + for<'f> &'f TYPE: + Mul + Add + Div + Neg, + { + for i in (usize::zero()..self.rows()).rev() { + if self[(i, i)] == TYPE::zero() { + for j in i..usize::zero() { + if self[(j, i)] != TYPE::zero() { + self.swap_rows_mul_pivot(other, pivot, i, j); + break; + } + } + } + for j in (usize::zero()..i).rev() { + if self[(i, i)] != TYPE::zero() { + let neg = -(&self[(j, i)]); + let factor = &neg / &(self[(i, i)]); + self.add_row_mul(other, i, j, &factor) + } + } + } + } + + fn triu_impl(&'a mut self, diag_replacement: Option) + where + TYPE: Zero + Clone, + { + match diag_replacement { + Some(replacement) => { + for (i, j, value) in self { + match &i.cmp(&j) { + std::cmp::Ordering::Equal => *value = replacement.clone(), + std::cmp::Ordering::Greater => *value = TYPE::zero(), + _ => {} + } + } + } + None => { + for (i, j, value) in self { + if i > j { + *value = TYPE::zero() + } + } + } + } + } + + fn tril_impl(&'a mut self, diag_replacement: Option) + where + TYPE: Zero + Clone, + { + match diag_replacement { + Some(replacement) => { + for (i, j, value) in self { + match i.cmp(&j) { + std::cmp::Ordering::Equal => *value = replacement.clone(), + std::cmp::Ordering::Less => *value = TYPE::zero(), + _ => {} + } + } + } + None => { + for (i, j, value) in self { + if i < j { + *value = TYPE::zero() + } + } + } + } + } + + fn invert(&mut self, pivot: &mut PIVOT, result: &mut Self) + where + TYPE: Clone + Zero + One + PartialEq + Neg + Div, + for<'f> &'f TYPE: + Mul + Add + Div + Neg, + { + self.gauss_jordan_low_mul(result, pivot); + self.gauss_jordan_high_mul(result, pivot); + for i in 0..result.rows() { + let f = self[(i, i)].clone(); + for j in 0..result.columns() { + result[(i, j)] = &(result[(i, j)]) / &f; + } + } + } + + fn det(mut self, mut pivot: PIVOT) -> TYPE { + self.gauss_jordan_low(&mut pivot); + (0..self.rows()) + .map(|i| &self[(i, i)]) + .fold(TYPE::one(), |a, b| &a * b) + } + + fn lu_pivot(&mut self, i: usize, pivot: &mut PIVOT) + where + TYPE: PartialOrd + Zero + Clone, + for<'f> &'f TYPE: Neg, + { + let mut max = abs(&self[(i, i)]); + let mut max_index = i; + for j in (i + 1)..self.rows() { + if abs(&self[(j, i)]) > max { + max = abs(&self[(i, j)]); + max_index = j; + } + } + if max_index != i { + self.swap_rows_pivot(pivot, i, max_index) + } + } + + fn lu(&mut self) -> Result<(), SingularMatrixError> + where + TYPE: Zero + Clone, + for<'b> &'b TYPE: Mul<&'b TYPE, Output = TYPE>, + TYPE: SubAssign + DivAssign, + { + for i in 0..self.rows() { + self.lu_row(i)? + } + Ok(()) + } + + fn lup(&mut self, p: &mut PIVOT) -> Result<(), SingularMatrixError> + where + TYPE: Zero + Clone + PartialOrd, + for<'b> &'b TYPE: Mul<&'b TYPE, Output = TYPE> + Neg, + TYPE: SubAssign + DivAssign, + { + for i in 0..self.rows() { + self.lu_pivot(i, p); + self.lu_row(i)? + } + Ok(()) + } + + fn lu_invert(&self, result: &mut Self, p: &mut PIVOT) + where + TYPE: Zero + One + Clone + PartialOrd, + for<'b> &'b TYPE: Mul<&'b TYPE, Output = TYPE>, + for<'c> TYPE: DivAssign<&'c TYPE>, + TYPE: SubAssign, + { + for i in 0..self.rows() { + for j in 0..self.rows() { + if p[j] == i { + result[(j, i)] = TYPE::one(); + } else { + result[(j, i)] = TYPE::zero(); + } + for k in 0..j { + let delta = &self[(j, k)] * &result[(k, i)]; + result[(j, i)] -= delta; + } + } + for j in (0..self.rows()).rev() { + for k in (j + 1)..self.rows() { + let delta = &self[(j, k)] * &result[(k, i)]; + result[(j, i)] -= delta; + } + result[(j, i)] /= &self[(j, j)]; + } + } + } + + fn lu_row(&mut self, i: usize) -> Result<(), SingularMatrixError> + where + TYPE: Zero + Clone, + for<'b> &'b TYPE: Mul<&'b TYPE, Output = TYPE>, + TYPE: SubAssign + DivAssign, + { + if self[(i, i)].is_zero() { + return Err(SingularMatrixError {}); + } + for j in i..self.columns() { + for k in 0..i { + let delta = &self[(i, k)] * &self[(k, j)]; + self[(i, j)] -= delta; + } + } + for j in (i + 1)..self.columns() { + for k in 0..i { + let delta = &self[(j, k)] * &self[(k, i)]; + self[(j, i)] -= delta; + } + let div = self[(i, i)].clone(); + self[(j, i)] /= div; + } + Ok(()) + } + + fn squared_norm2(&'a self) -> TYPE { + self.into_iter().fold(TYPE::zero(), move |acc, (_, _, v)| { acc + v * v }) + } +} + +fn abs(value: &TYPE) -> TYPE +where + TYPE: PartialOrd + Zero + Clone, + for<'a> &'a TYPE: Neg, +{ + if value >= &TYPE::zero() { + value.clone() + } else { + -value + } +} diff --git a/src/misc.rs b/src/misc.rs new file mode 100644 index 0000000..ee13d0b --- /dev/null +++ b/src/misc.rs @@ -0,0 +1,49 @@ +use std::fmt::Debug; +use std::ops::Index; +use std::ops::IndexMut; + +use crate::Matrix; + +pub trait EnhancedOption { + fn otherwise(self, cb: F) -> U + where + F: Fn() -> U; +} + +impl EnhancedOption for Option { + fn otherwise(self, cb: F) -> T + where + F: Fn() -> T, + { + if let Some(value) = self { + value + } else { + cb() + } + } +} + +pub fn debug_fmt<'a, M: Matrix, TYPE: 'a + Debug>( + m: &M, + f: &mut std::fmt::Formatter<'_>, +) -> std::fmt::Result +where + M: Index<(usize, usize), Output = TYPE> + IndexMut<(usize, usize), Output = TYPE>, + M: IntoIterator, +{ + write!(f, "[")?; + for i in 0..m.rows() { + if i > 0 { + write!(f, ", ")?; + } + write!(f, "[")?; + for j in 0..m.columns() { + if j > 0 { + write!(f, ", ")?; + } + write!(f, "{:?}", m[(i, j)])? + } + write!(f, "]")?; + } + write!(f, "]") +} diff --git a/src/rational.rs b/src/rational.rs new file mode 100644 index 0000000..41dad87 --- /dev/null +++ b/src/rational.rs @@ -0,0 +1,791 @@ +use num_traits::sign::Signed; +use num_traits::Num; +use num_traits::NumOps; +use num_traits::One; +use num_traits::Pow; +use num_traits::Signed as SignedOps; +use num_traits::Zero; +use std::cmp::Ordering; +use std::cmp::PartialEq; +use std::cmp::PartialOrd; +use std::fmt::Debug; +use std::fmt::Display; +use std::ops::Add; +use std::ops::AddAssign; +use std::ops::Div; +use std::ops::DivAssign; +use std::ops::Mul; +use std::ops::MulAssign; +use std::ops::Neg; +use std::ops::Rem; +use std::ops::Sub; +use std::ops::SubAssign; + +trait_group! { + pub trait RationalInt + : NumOps + + Debug + + PartialEq + + PartialOrd + + Zero + + One + + Pow + + Neg + + Clone + + Default +} + +fn gcd(a: &T, b: &T) -> T +where + T: Rem + Clone + Zero, +{ + let mut n1 = a.clone(); + let mut n2 = b.clone(); + let mut tmp: T; + while !n2.is_zero() { + tmp = n1; + n1 = n2.clone(); + n2 = tmp % n2; + } + n1 +} + +fn mcm<'a, T>(n1: &'a T, n2: &'a T) -> T +where + T: Div + Rem + Clone + Zero, + &'a T: Mul, +{ + let den = gcd(n1, n2); + n1 * n2 / den +} + +#[derive(Debug)] +pub struct Rational { + num: TYPE, + den: TYPE, +} + +impl Rational +where + TYPE: Zero + One + Clone + Rem + PartialOrd + Neg, + for<'a> &'a TYPE: Div<&'a TYPE, Output = TYPE>, +{ + fn simplify(&mut self) { + if self.num.is_zero() { + self.den = TYPE::one(); + } else if self.den.is_zero() { + } else { + (self.num, self.den) = { + let n1_opt: TYPE; + let n1 = if self.den < TYPE::zero() { + n1_opt = -self.num.clone(); + &n1_opt + } else { + &self.num + }; + let n2_opt: TYPE; + let n2 = if self.den < TYPE::zero() { + n2_opt = -self.den.clone(); + &n2_opt + } else { + &self.den + }; + let num_abs_opt: TYPE; + let num_abs = if self.num < TYPE::zero() { + num_abs_opt = -self.num.clone(); + &num_abs_opt + } else { + &self.num + }; + let gcd = gcd(num_abs, n2); + (n1 / &gcd, n2 / &gcd) + }; + } + } +} + +impl Rational +where + TYPE: RationalInt, +{ + pub fn new(num: TYPE, den: TYPE) -> Rational { + let res = den.partial_cmp(&TYPE::zero()); + match res { + Some(Ordering::Greater) => Rational { num, den }, + Some(Ordering::Less) => Rational { + num: -num, + den: -den, + }, + _ => { + panic!("Rational denominator cannot be zero or NaN") + } + } + } +} + +#[opimps::impl_op(Pow)] +fn pow(self: Rational, exp: u32) -> Rational +where + TYPE: RationalInt, + for<'a> &'a TYPE: Pow, +{ + Rational { + num: self.num.pow(exp), + den: self.den.pow(exp), + } +} + +#[opimps::impl_op(Pow)] +fn pow(self: &Rational, exp: u32) -> Rational +where + TYPE: RationalInt, + for<'a> &'a TYPE: Pow, +{ + Rational { + num: (&self.num).pow(exp), + den: (&self.den).pow(exp), + } +} + +#[opimps::impl_ops(Add)] +fn add(self: Rational, rhs: Rational) -> Rational +where + TYPE: RationalInt, + for<'a> &'a TYPE: NumOps<&'a TYPE, TYPE>, +{ + let den = mcm::(&self.den, &rhs.den); + let a = &self.num * &den; + let b = &rhs.num * &den; + let num = &a / &self.den + &b / &rhs.den; + let mut result: Rational = Rational { num, den }; + result.simplify(); + result +} + +#[opimps::impl_ops(Sub)] +fn sub(self: Rational, rhs: Rational) -> Rational +where + TYPE: RationalInt, + for<'a> &'a TYPE: NumOps<&'a TYPE, TYPE>, +{ + let den = mcm::(&self.den, &rhs.den); + let a = &self.num * &den; + let b = &rhs.num * &den; + let num = &a / &self.den - &b / &rhs.den; + let mut result: Rational = Rational { num, den }; + result.simplify(); + result +} + +#[opimps::impl_ops(Mul)] +fn mul(self: Rational, rhs: Rational) -> Rational +where + TYPE: RationalInt, + for<'a> &'a TYPE: NumOps<&'a TYPE, TYPE>, +{ + let num = &self.num * &rhs.num; + let den = &self.den * &rhs.den; + let mut result = Rational { num, den }; + result.simplify(); + result +} + +#[opimps::impl_ops(Div)] +fn div(self: Rational, rhs: Rational) -> Rational +where + TYPE: RationalInt, + for<'a> &'a TYPE: NumOps<&'a TYPE, TYPE>, +{ + let den = &self.den * &rhs.num; + let num = &self.num * &rhs.den; + let mut result = Rational::new(num, den); + result.simplify(); + result +} + +#[opimps::impl_ops(Rem)] +fn rem(self: Rational, rhs: Rational) -> Rational +where + TYPE: RationalInt, + for<'a> &'a TYPE: NumOps<&'a TYPE, TYPE>, +{ + match self.partial_cmp(&rhs) { + Some(std::cmp::Ordering::Less) => rhs.clone(), + Some(std::cmp::Ordering::Equal) => Rational::::zero(), + Some(std::cmp::Ordering::Greater) => { + let div = self / rhs; + let num = &div.num / &div.den; + div - Rational::new(num, TYPE::one()) + } + None => { + panic!("Cannot compare") + } + } +} + +impl<'a, TYPE> AddAssign<&'a Rational> for Rational +where + TYPE: Zero + + One + + Clone + + Rem + + PartialOrd + + Neg + + Div, + for<'b> &'b TYPE: Mul<&'b TYPE, Output = TYPE> + Div<&'b TYPE, Output = TYPE>, +{ + fn add_assign(&mut self, rhs: &'a Self) { + let den = mcm::(&self.den, &rhs.den); + let a = &self.num * &den; + let b = &rhs.num * &den; + let num = &a / &self.den + &b / &rhs.den; + self.num = num; + self.den = den; + self.simplify(); + } +} + +impl AddAssign for Rational +where + TYPE: RationalInt, + for<'a> &'a TYPE: NumOps<&'a TYPE, TYPE>, +{ + fn add_assign(&mut self, rhs: Self) { + *self += &rhs; + } +} + +impl<'a, TYPE> SubAssign<&'a Rational> for Rational +where + TYPE: Div + + Sub + + Rem + + Neg + + Clone + + Zero + + One + + PartialOrd, + for<'c> &'c TYPE: Mul<&'c TYPE, Output = TYPE> + Div<&'c TYPE, Output = TYPE>, + for<'b> TYPE: SubAssign<&'b TYPE>, +{ + fn sub_assign(&mut self, rhs: &'a Self) { + let den = mcm::(&self.den, &rhs.den); + let a = &self.num * &den; + let b = &rhs.num * &den; + let num = &a / &self.den - &b / &rhs.den; + self.num = num; + self.den = den; + self.simplify(); + } +} + +impl SubAssign for Rational +where + TYPE: Div + + Sub + + Rem + + Neg + + Clone + + Zero + + One + + PartialOrd, + for<'c> &'c TYPE: Mul<&'c TYPE, Output = TYPE> + Div<&'c TYPE, Output = TYPE>, + for<'b> TYPE: SubAssign<&'b TYPE>, +{ + fn sub_assign(&mut self, rhs: Self) { + *self -= &rhs; + } +} + +impl<'a, TYPE> MulAssign<&'a Rational> for Rational +where + for<'b> TYPE: MulAssign<&'b TYPE>, + TYPE: PartialOrd + One + Clone + Zero + Rem + Neg, + for<'c> &'c TYPE: Div<&'c TYPE, Output = TYPE>, +{ + fn mul_assign(&mut self, rhs: &'a Self) { + self.num *= &rhs.num; + self.den *= &rhs.den; + self.simplify(); + } +} + +impl MulAssign for Rational +where + for<'a> TYPE: MulAssign<&'a TYPE>, + TYPE: PartialOrd + One + Clone + Zero + Rem + Neg, + for<'c> &'c TYPE: Div<&'c TYPE, Output = TYPE>, +{ + fn mul_assign(&mut self, rhs: Self) { + *self *= &rhs; + } +} + +impl<'a, TYPE> DivAssign<&'a Rational> for Rational +where + for<'b> TYPE: MulAssign<&'b TYPE>, + TYPE: PartialOrd + One + Clone + Zero + Rem + Neg, + for<'c> &'c TYPE: Div<&'c TYPE, Output = TYPE>, +{ + fn div_assign(&mut self, rhs: &'a Self) { + self.num *= &rhs.den; + self.den *= &rhs.num; + self.simplify() + } +} + +impl DivAssign for Rational +where + for<'a> TYPE: MulAssign<&'a TYPE>, + TYPE: PartialOrd + One + Clone + Zero + Rem + Neg, + for<'c> &'c TYPE: Div<&'c TYPE, Output = TYPE>, +{ + fn div_assign(&mut self, rhs: Self) { + *self /= &rhs; + self.simplify() + } +} + +impl Display for Rational +where + TYPE: RationalInt + Display, +{ + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + write!(f, "{}/{}", self.num, self.den) + } +} + +impl PartialEq for Rational +where + TYPE: RationalInt, + for<'a> &'a TYPE: NumOps<&'a TYPE, TYPE>, +{ + fn eq(&self, other: &Self) -> bool { + self.partial_cmp(other) == Some(Ordering::Equal) + } +} + +impl Eq for Rational +where + TYPE: RationalInt + Eq, + for<'a> &'a TYPE: NumOps<&'a TYPE, TYPE>, +{ +} + +impl PartialOrd for Rational +where + TYPE: RationalInt, + for<'a> &'a TYPE: NumOps<&'a TYPE, TYPE>, +{ + fn partial_cmp(&self, other: &Self) -> Option { + (&self.num * &other.den).partial_cmp(&(&other.num * &self.den)) + } +} + +impl Ord for Rational +where + TYPE: RationalInt + Ord, + for<'a> &'a TYPE: NumOps<&'a TYPE, TYPE>, +{ + fn cmp(&self, other: &Self) -> Ordering { + (&self.num * &other.den).cmp(&(&other.num * &self.den)) + } +} + +impl Clone for Rational +where + TYPE: Clone, +{ + fn clone(&self) -> Self { + Rational { + num: self.num.clone(), + den: self.den.clone(), + } + } +} + +impl Copy for Rational where TYPE: Copy {} + +impl Zero for Rational +where + TYPE: RationalInt, + for<'a> &'a TYPE: NumOps<&'a TYPE, TYPE>, +{ + fn zero() -> Rational { + Rational { + num: TYPE::zero(), + den: TYPE::one(), + } + } + + fn is_zero(&self) -> bool { + self.num.is_zero() && !self.den.is_zero() + } +} + +impl One for Rational +where + TYPE: RationalInt + Display, + for<'a> &'a TYPE: NumOps<&'a TYPE, TYPE>, +{ + fn one() -> Self { + Rational { + num: TYPE::one(), + den: TYPE::one(), + } + } + + fn is_one(&self) -> bool { + self.num == self.den + } +} + +impl Default for Rational +where + TYPE: RationalInt + Display, + for<'a> &'a TYPE: NumOps<&'a TYPE, TYPE>, +{ + fn default() -> Self { + Rational::::zero() + } +} + +#[opimps::impl_uni_op(Neg)] +fn neg(self: Rational) -> Rational +where + TYPE: RationalInt, +{ + Rational::new(-self.num, self.den) +} + +#[opimps::impl_uni_op(Neg)] +fn neg(self: &Rational) -> Rational +where + TYPE: RationalInt, +{ + Rational::new(-self.num.clone(), self.den.clone()) +} + +impl Num for Rational +where + TYPE: Num + RationalInt + Display, + for<'a> &'a TYPE: NumOps<&'a TYPE, TYPE>, +{ + fn from_str_radix(_: &str, _: u32) -> Result { + Result::Err(Rational::::default()) + } + + type FromStrRadixErr = Self; +} + +impl SignedOps for Rational +where + TYPE: Signed + SignedOps + Display + RationalInt, + for<'a> &'a TYPE: NumOps<&'a TYPE, TYPE>, +{ + fn abs(&self) -> Self { + Rational::new(self.num.abs(), self.den.abs()) + } + + fn abs_sub(&self, other: &Self) -> Self { + if self <= other { + Rational::::zero() + } else { + self.clone() - other.clone() + } + } + + fn is_negative(&self) -> bool { + !self.num.is_zero() && self.num.is_negative() ^ self.den.is_negative() + } + + fn is_positive(&self) -> bool { + !self.num.is_zero() && !(self.num.is_positive() ^ self.den.is_positive()) + } + + fn signum(&self) -> Self { + if self.is_positive() { + Rational::::one() + } else if self.is_negative() { + -Rational::::one() + } else { + Rational::::zero() + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + use std::collections::BTreeSet; + + #[test] + fn test_simplify() { + let mut r1 = Rational::new(15, 105); + assert_eq!(Rational::new(1, 7), r1); + r1.simplify(); + assert_eq!(Rational::new(1, 7), r1); + } + + #[test] + fn test_add() { + let r1 = Rational::new(3, 5); + let r2 = Rational::new(4, 10); + let r3 = Rational::new(1, 7); + let one: Rational = Rational::::one(); + assert_eq!(one, r1 + r2); + assert_eq!(r2 + r1, r1 + r2); + assert_eq!(r1 + r2 + r3, r3 + r2 + r1); + } + + #[test] + fn test_add_assign() { + let mut r1 = Rational::new(3, 5); + let r2 = Rational::new(4, 10); + let r3 = Rational::::one(); + r1 += r2; + assert_eq!(r3, r1); + } + + #[test] + fn test_sub() { + let r1 = Rational::new(3, 5); + let r2 = Rational::new(4, 10); + let r3 = Rational::new(1, 7); + assert_eq!(Rational::::zero(), r1 - r1); + assert_eq!(Rational::new(-1, 5), r2 - r1); + assert!(Rational::::zero() > r3 - r1); + assert!(Rational::::zero() > r3 - r2); + assert!(Rational::::zero() > r3 - r1); + } + + #[test] + fn test_sub_assign() { + let mut r1 = Rational::new(3, 5); + let r2 = Rational::new(4, 10); + let r3 = Rational::new(1, 5); + r1 -= r2; + assert_eq!(r3, r1); + } + + #[test] + fn test_mul() { + let r1 = Rational::new(3, 5); + let r2 = Rational::new(-4, 10); + let r3 = Rational::new(1, 7); + assert_eq!(Rational::new(9, 25), r1 * r1); + assert_eq!(Rational::new(-12, 50), r2 * r1); + assert_eq!(Rational::new(2, -35), r2 * r3); + } + + #[test] + fn test_mul_assign() { + let r1 = Rational::new(3, 5); + let r2 = Rational::new(-4, 10); + let r3 = Rational::new(-12, 50); + let mut r4 = r1.clone(); + r4 *= r1; + assert_eq!(Rational::::new(9, 25), r4); + let mut r5 = r1.clone(); + r5 *= r2; + assert_eq!(r3, r5); + } + + #[test] + fn test_div() { + let r1 = Rational::new(3, 5); + let r2 = Rational::new(-4, 10); + let r3 = Rational::new(1, 7); + assert_eq!(Rational::::one(), r1 / r1); + assert_eq!(Rational::new(-2, 3), r2 / r1); + assert_eq!(Rational::new(28, -10), r2 / r3); + } + + #[test] + fn test_div_assign() { + let r1 = Rational::new(3, 5); + let r2 = Rational::new(-4, 10); + let r3 = Rational::new(30, -20); + let mut r4 = r1.clone(); + r4 /= r1; + assert_eq!(Rational::::one(), r4); + let mut r5 = r1.clone(); + r5 /= r2; + assert_eq!(r3, r5); + } + + #[test] + fn test_eq() { + let r1 = Rational::new(3, 5); + let r2 = Rational::new(6, 10); + let r3 = Rational::new(300, 500); + let r4 = Rational::new(300, -500); + assert!(r1 == r2); + assert!(r1 == r3); + assert!(r1 != r4); + } + + #[test] + fn test_ord() { + let r1 = Rational::new(3, 5); + let r2 = Rational::new(7, 10); + let r3 = Rational::new(301, 500); + let r4 = Rational::new(300, -500); + let values: Vec> = BTreeSet::from([r1, r2, r3, r4]) + .iter() + .map(|it| it.clone()) + .collect(); + let expected_order = Vec::from([r4, r1, r3, r2]); + assert_eq!(expected_order, values); + } + + #[test] + fn test_gcd() { + let a = 14; + let b = 21; + assert_eq!(7, gcd(&a, &b)); + } + + #[test] + fn test_mcm() { + let a = 14; + let b = 21; + assert_eq!(42, mcm::(&a, &b)); + } +} + +#[cfg(test)] +mod big_int_tests { + use super::*; + use num_bigint::BigInt; + use std::collections::BTreeSet; + + #[test] + fn test_gcd() { + assert_eq!(6, gcd(&12, &18)); + assert_eq!(6, gcd(&18, &12)); + } + + #[test] + fn test_simplify() { + let mut r1 = Rational::new(BigInt::from(15), BigInt::from(105)); + assert_eq!(Rational::new(BigInt::from(1), BigInt::from(7)), r1); + r1.simplify(); + assert_eq!(Rational::new(BigInt::from(1), BigInt::from(7)), r1); + } + + #[test] + fn test_add() { + let r1 = Rational::new(BigInt::from(3), BigInt::from(5)); + let r2 = Rational::new(BigInt::from(4), BigInt::from(10)); + let r3 = Rational::new(BigInt::from(1), BigInt::from(7)); + assert_eq!(Rational::::one(), &r1 + &r2); + assert_eq!(&r2 + &r1, &r1 + &r2); + assert_eq!(&r1 + &r2 + &r3, &r3 + &r2 + &r1); + } + + #[test] + fn test_add_assign() { + let mut r1 = Rational::new(BigInt::from(3), BigInt::from(5)); + let r2 = Rational::new(BigInt::from(4), BigInt::from(10)); + r1 += r2; + assert_eq!(Rational::::one(), r1); + } + + #[test] + fn test_sub() { + let r1 = Rational::new(BigInt::from(3), BigInt::from(5)); + let r2 = Rational::new(BigInt::from(4), BigInt::from(10)); + let r3 = Rational::new(BigInt::from(1), BigInt::from(7)); + assert_eq!(Rational::::zero(), &r1 - &r1); + assert_eq!(Rational::new(BigInt::from(-1), BigInt::from(5)), &r2 - &r1); + assert!(Rational::::zero() > &r3 - &r1); + assert!(Rational::::zero() > &r3 - &r2); + assert!(Rational::::zero() > r3 - r1); + } + + #[test] + fn test_sub_assign() { + let mut r1 = Rational::new(BigInt::from(3), BigInt::from(5)); + let r2 = Rational::new(BigInt::from(4), BigInt::from(10)); + r1 -= r2; + assert_eq!(Rational::::new(BigInt::one(), BigInt::from(5)), r1); + } + + #[test] + fn test_mul() { + let r1 = Rational::new(BigInt::from(3), BigInt::from(5)); + let r2 = Rational::new(BigInt::from(-4), BigInt::from(10)); + let r3 = Rational::new(BigInt::from(1), BigInt::from(7)); + assert_eq!(Rational::new(BigInt::from(9), BigInt::from(25)), &r1 * &r1); + assert_eq!( + Rational::new(BigInt::from(-12), BigInt::from(50)), + &r2 * &r1 + ); + assert_eq!(Rational::new(BigInt::from(2), BigInt::from(-35)), r2 * r3); + } + + #[test] + fn test_mul_assign() { + let r1 = Rational::new(BigInt::from(3), BigInt::from(5)); + let r2 = Rational::new(BigInt::from(-4), BigInt::from(10)); + let r3 = Rational::new(BigInt::from(-12), BigInt::from(50)); + let mut r4 = r1.clone(); + r4 *= &r1; + assert_eq!( + Rational::::new(BigInt::from(9), BigInt::from(25)), + r4 + ); + let mut r5 = r1.clone(); + r5 *= r2; + assert_eq!(r3, r5); + } + + #[test] + fn test_div() { + let r1 = Rational::new(BigInt::from(3), BigInt::from(5)); + let r2 = Rational::new(BigInt::from(-4), BigInt::from(10)); + let r3 = Rational::new(BigInt::from(1), BigInt::from(7)); + assert_eq!(Rational::::one(), &r1 / &r1); + assert_eq!(Rational::new(BigInt::from(-2), BigInt::from(3)), &r2 / r1); + assert_eq!(Rational::new(BigInt::from(28), BigInt::from(-10)), r2 / r3); + } + + #[test] + fn test_div_assign() { + let r1 = Rational::new(BigInt::from(3), BigInt::from(5)); + let r2 = Rational::new(BigInt::from(-4), BigInt::from(10)); + let r3 = Rational::new(BigInt::from(30), BigInt::from(-20)); + let mut r4 = r1.clone(); + r4 /= &r1; + assert_eq!(Rational::::one(), r4); + let mut r5 = r1.clone(); + r5 /= r2; + assert_eq!(r3, r5); + } + + #[test] + fn test_eq() { + let r1 = Rational::new(BigInt::from(3), BigInt::from(5)); + let r2 = Rational::new(BigInt::from(6), BigInt::from(10)); + let r3 = Rational::new(BigInt::from(300), BigInt::from(500)); + let r4 = Rational::new(BigInt::from(300), BigInt::from(-500)); + assert!(r1 == r2); + assert!(r1 == r3); + assert!(r1 != r4); + } + + #[test] + fn test_ord() { + let r1 = Rational::new(BigInt::from(3), BigInt::from(5)); + let r2 = Rational::new(BigInt::from(7), BigInt::from(10)); + let r3 = Rational::new(BigInt::from(301), BigInt::from(500)); + let r4 = Rational::new(BigInt::from(300), BigInt::from(-500)); + let values: Vec<&Rational> = BTreeSet::from([&r1, &r2, &r3, &r4]) + .iter() + .map(|it| *it) + .collect(); + let expected_order = Vec::from([&r4, &r1, &r3, &r2]); + assert_eq!(expected_order, values); + } +} diff --git a/src/sludecomposition.rs b/src/sludecomposition.rs new file mode 100644 index 0000000..2307a05 --- /dev/null +++ b/src/sludecomposition.rs @@ -0,0 +1,208 @@ +use num_traits::One; +use num_traits::Zero; +use std::ops::Add; +use std::ops::Div; +use std::ops::DivAssign; +use std::ops::Mul; +use std::ops::MulAssign; +use std::ops::Neg; +use std::ops::SubAssign; + +use crate::Matrix; +use crate::SingularMatrixError; + +use super::matrix::LUDecompositionImpl; + +use super::matrix::LUDecomposition; +use super::matrix::NumericalMatrixImpl; +use super::spivot::SPivot; +use super::NumericalMatrix; +use super::Pivot; +use super::SMatrix; + +use super::misc::EnhancedOption; + +pub struct SLUDecomposition { + matrix: SMatrix, + pivot: SPivot, +} + +impl SLUDecomposition +where TYPE : Zero + SubAssign + Clone, + for<'a> &'a TYPE : Mul<&'a TYPE, Output = TYPE>, + for<'a> TYPE : DivAssign<&'a TYPE> { + pub (crate) fn new(matrix: SMatrix, pivot: Option>) -> Self { + SLUDecomposition { + matrix, + pivot: pivot.otherwise(|| SPivot::new(SIZE)), + } + } + + pub fn solve(&self, b : &SMatrix) -> Result, SingularMatrixError> { + let mut x = SMatrix::::new(|_| TYPE::zero()); + for n in 0..b.columns() { + for i in 0..self.matrix.rows() { + x[(i, n)] = b[(self.pivot[i], n)].clone(); + for k in 0..i { + let sub = &self.matrix[(i, k)] * &x[(k, n)]; + x[(i, n)] -= sub; + } + } + for i in (0..self.matrix.rows()).rev() { + for k in (i + 1)..self.matrix.rows() { + let sub = &self.matrix[(i, k)] * &x[(k, n)]; + x[(i, n)] -= sub; + } + if self.matrix[(i, i)].is_zero() { + return Err(SingularMatrixError::default()); + } else { + x[(i, n)] /= &self.matrix[(i, i)]; + } + } + } + Ok(x) + } +} + +impl LUDecomposition, SPivot> + for SLUDecomposition +where + for<'b> TYPE: Clone + + Zero + + One + + PartialEq + + PartialOrd + + Neg + + Div + + SubAssign + + SubAssign<&'b TYPE> + + DivAssign + + DivAssign<&'b TYPE> + + MulAssign<&'b TYPE>, + for<'c> &'c TYPE: + Mul + Add + Div + Neg, +{ + fn l(&self) -> SMatrix { + self.matrix.clone().tril_replace(TYPE::one()) + } + + fn u(&self) -> SMatrix { + self.matrix.clone().triu() + } + + fn invert(&self) -> SMatrix { + let mut result = SMatrix::::new(|_| TYPE::zero()); + self.matrix.lu_invert(&mut result, &mut self.pivot.clone()); + result + } + + fn pivot(&self, m: SMatrix) -> SMatrix { + &self.pivot * m + } + // fn p<'a>(&'a self) -> &'a HPivot { + // &self.pivot + // } +} + +impl LUDecompositionImpl, SPivot> + for SLUDecomposition +where + for<'b> TYPE: One + MulAssign<&'b TYPE> + Neg, +{ + fn get_matrix(&self) -> &SMatrix { + &self.matrix + } + + fn get_pivot(&self) -> &SPivot { + &self.pivot + } +} + + +#[cfg(test)] +mod test { + use rand::rngs::StdRng; + use rand::RngCore; + use rand::SeedableRng; + + use crate::NumericalMatrix; + use crate::Rational; + use crate::SMatrix; + use num_traits::Zero; + + use rand; + + #[test] + fn solve_linear_system() { + let mut rand = { + let seed = [ + 1,0,1,3, + 2,5,0,0, + 200,1,0,0, + 210,30,0,0, + 78,134,31,0, + 253,11,7,0, + 120,169,89,48, + 200,0,202,0 + ]; + StdRng::from_seed(seed) + }; + let mtx = SMatrix::, 5, 5>::new(|_| { + Rational::new(i64::try_from(rand.next_u32() % 40).unwrap() - 20, 1) + }); + let b = SMatrix::, 5, 5>::new(|_| { + Rational::new(i64::try_from(rand.next_u32() % 20).unwrap() - 10, 1) + }); + let lu = mtx.clone().lu().unwrap(); + let x = lu.solve(&b).unwrap(); + assert!((mtx * x - b).is_zero()); + } + +} + + +#[cfg(test)] +mod test_big_int { + use num_traits::Zero; + use rand::rngs::StdRng; + use rand::RngCore; + use rand::SeedableRng; + + use crate::NumericalMatrix; + use crate::Rational; + use crate::SMatrix; + + use rand; + + #[test] + fn solve_linear_system() { + use num_bigint::BigInt; + use num_traits::One; + + let mut rand = { + let seed = [ + 1,0,1,3, + 2,5,0,0, + 200,1,0,0, + 210,30,0,0, + 78,134,31,0, + 253,11,7,0, + 120,169,89,48, + 200,0,202,0 + ]; + StdRng::from_seed(seed) + }; + let mtx = SMatrix::, 10, 10>::new(|_| { + Rational::new(BigInt::from(rand.next_u32() % 200) - 100, BigInt::one()) + }); + let b = SMatrix::, 10, 10>::new(|_| { + Rational::new(BigInt::from(rand.next_u32() % 100) - 50, BigInt::one()) + }); + let lu = mtx.clone().lu().unwrap(); + let x = lu.solve(&b).unwrap(); + + assert!((mtx * x - b).is_zero()); + + } + +} \ No newline at end of file diff --git a/src/smatrix.rs b/src/smatrix.rs new file mode 100644 index 0000000..264f335 --- /dev/null +++ b/src/smatrix.rs @@ -0,0 +1,956 @@ +use num_traits::One; +use num_traits::Zero; + +use std::cmp::PartialEq; +use std::fmt::Debug; +use std::fmt::Display; +use std::iter::Enumerate; +use std::iter::Flatten; +use std::iter::IntoIterator; +use std::ops::Add; +use std::ops::Div; +use std::ops::DivAssign; +use std::ops::FnMut; +use std::ops::Index; +use std::ops::IndexMut; +use std::ops::Mul; +use std::ops::MulAssign; +use std::ops::Neg; +use std::ops::Sub; +use std::ops::SubAssign; + +use super::err::SizeError; + +use super::matrix::MatrixImpl; +use super::matrix::NumericalMatrixImpl; +use super::misc::debug_fmt; +use super::sludecomposition::SLUDecomposition; +use super::Matrix; +use super::NumericalMatrix; +use super::Pivot; +use super::SingularMatrixError; + +use super::spivot::SPivot; + +pub struct SMatrix { + values: [[TYPE; COLUMNS]; ROWS], +} + +impl SMatrix +where + TYPE: Clone, +{ + pub fn transpose(&self) -> SMatrix { + SMatrix::::new(|(i, j)| self[(j, i)].clone()) + } +} + +impl SMatrix { + pub fn transpose_in_place(mut self) + where + Self: Sized, + { + for i in 0..self.rows() { + for j in 0..self.columns() { + self.swap((i, j), (j, i)); + } + } + } +} + +impl Index<(usize, usize)> + for SMatrix +{ + fn index(&self, index: (usize, usize)) -> &TYPE { + &self.values[index.0][index.1] + } + + type Output = TYPE; +} + +impl IndexMut<(usize, usize)> + for SMatrix +{ + fn index_mut(&mut self, index: (usize, usize)) -> &mut TYPE { + &mut self.values[index.0][index.1] + } +} + +impl Default for SMatrix +where + TYPE: Default, +{ + fn default() -> Self { + Self::new(|_| TYPE::default()) + } +} + +impl Clone for SMatrix +where + TYPE: Clone, +{ + fn clone(&self) -> Self { + SMatrix::new(|it| self[it].clone()) + } +} + +impl Copy for SMatrix where + TYPE: Copy +{ +} + +impl<'a, TYPE, const ROWS: usize, const COLUMNS: usize> MatrixImpl<'a, TYPE> + for SMatrix +where + TYPE: 'a, +{ + type IteratorType = SMatrixIterator; + type RefIteratorType<'b> = SMatrixIteratorRef<'a, TYPE, ROWS, COLUMNS>; + type MutRefIteratorType<'b> = SMatrixIteratorMut<'a, TYPE, ROWS, COLUMNS>; +} + +impl Matrix for SMatrix { + fn columns(&self) -> usize { + COLUMNS + } + + fn rows(&self) -> usize { + ROWS + } + + fn transpose(mut self) -> Result { + if ROWS != COLUMNS { + Err(SizeError::new("Matrix must be squared")) + } else { + self.transpose_in_place()?; + Ok(self) + } + } +} + +#[opimps::impl_ops(Add)] +fn add( + self: SMatrix, + rhs: SMatrix, +) -> SMatrix +where + TYPE: Add, + for<'a> &'a TYPE: Add<&'a TYPE, Output = TYPE>, +{ + SMatrix::new(|position| &self[position] + &rhs[position]) +} + +#[opimps::impl_ops(Sub)] +fn sub( + self: SMatrix, + rhs: SMatrix, +) -> SMatrix +where + TYPE: Sub, + for<'a> &'a TYPE: Sub<&'a TYPE, Output = TYPE>, +{ + SMatrix::new(|position| &self[position] - &rhs[position]) +} + +#[opimps::impl_ops(Mul)] +fn mul< + TYPE, + const L_ROWS: usize, + const L_COLUMNS: usize, + const R_ROWS: usize, + const R_COLUMNS: usize, +>( + self: SMatrix, + rhs: SMatrix, +) -> SMatrix +where + TYPE: Mul + Add + Zero, + for<'a> &'a TYPE: Add<&'a TYPE, Output = TYPE> + Mul<&'a TYPE, Output = TYPE>, +{ + let mut result: SMatrix = SMatrix::new(|_| TYPE::zero()); + for i in 0..result.rows() { + for j in 0..result.columns() { + for k in 0..self.columns() { + let prod = &self[(i, k)] * &rhs[(k, j)]; + result[(i, j)] = &result[(i, j)] + ∏ + } + } + } + result +} + +impl Display for SMatrix +where + TYPE: Display, +{ + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + let lim = self.rows() * self.columns() - 1; + writeln!(f, "[")?; + for i in 0..self.rows() { + write!(f, " ")?; + for j in 0..self.columns() { + let index = i * self.columns() + j; + write!(f, "{}", self[(i, j)])?; + if index < lim { + write!(f, ", ")?; + } + } + writeln!(f)?; + } + write!(f, "]") + } +} + +impl SMatrix { + pub fn new(mut cb: INITIALIZER) -> Self + where + INITIALIZER: FnMut((usize, usize)) -> TYPE, + { + let values: [[TYPE; COLUMNS]; ROWS] = + std::array::from_fn(|i| std::array::from_fn(|j| cb((i, j)))); + SMatrix { values } + } +} + +pub struct SMatrixIterator { + it: Enumerate>>, +} + +impl Iterator + for SMatrixIterator +{ + fn next(&mut self) -> Option { + self.it + .next() + .map(|(index, value)| (index / COLUMNS, index % COLUMNS, value)) + } + type Item = (usize, usize, TYPE); +} + +impl IntoIterator for SMatrix { + fn into_iter(self) -> SMatrixIterator { + SMatrixIterator { + it: self.values.into_iter().flatten().enumerate(), + } + } + + type IntoIter = SMatrixIterator; + type Item = (usize, usize, TYPE); +} + +pub struct SMatrixIteratorRef<'a, TYPE, const ROWS: usize, const COLUMNS: usize> { + it: Enumerate>>, +} + +impl<'a, TYPE, const ROWS: usize, const COLUMNS: usize> Iterator + for SMatrixIteratorRef<'a, TYPE, ROWS, COLUMNS> +{ + fn next(&mut self) -> Option { + self.it + .next() + .map(|(index, value)| (index / COLUMNS, index % COLUMNS, value)) + } + type Item = (usize, usize, &'a TYPE); +} + +impl<'a, TYPE, const ROWS: usize, const COLUMNS: usize> IntoIterator + for &'a SMatrix +{ + fn into_iter(self) -> SMatrixIteratorRef<'a, TYPE, ROWS, COLUMNS> { + SMatrixIteratorRef { + it: (self.values).iter().flatten().enumerate(), + } + } + + type IntoIter = SMatrixIteratorRef<'a, TYPE, ROWS, COLUMNS>; + type Item = (usize, usize, &'a TYPE); +} + +pub struct SMatrixIteratorMut<'a, TYPE, const ROWS: usize, const COLUMNS: usize> { + it: Enumerate>>, +} + +impl<'a, TYPE, const ROWS: usize, const COLUMNS: usize> Iterator + for SMatrixIteratorMut<'a, TYPE, ROWS, COLUMNS> +{ + fn next(&mut self) -> Option { + self.it + .next() + .map(|(index, value)| (index / COLUMNS, index % COLUMNS, value)) + } + type Item = (usize, usize, &'a mut TYPE); +} + +impl<'a, TYPE, const ROWS: usize, const COLUMNS: usize> IntoIterator + for &'a mut SMatrix +{ + fn into_iter(self) -> SMatrixIteratorMut<'a, TYPE, ROWS, COLUMNS> { + SMatrixIteratorMut { + it: self.values.iter_mut().flatten().enumerate(), + } + } + + type IntoIter = SMatrixIteratorMut<'a, TYPE, ROWS, COLUMNS>; + type Item = (usize, usize, &'a mut TYPE); +} + +impl PartialEq for SMatrix +where + TYPE: PartialEq, +{ + fn eq(&self, other: &Self) -> bool { + (0..ROWS) + .flat_map(|i| (0..COLUMNS).map(move |j| (i, j))) + .all(|index| self[index] == other[index]) + } +} + +impl Debug for SMatrix +where + TYPE: Debug, +{ + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + debug_fmt(self, f) + } +} + +impl Zero for SMatrix +where + TYPE: PartialEq + One + Zero + Add, + for<'a> &'a TYPE: Add<&'a TYPE, Output = TYPE>, +{ + fn zero() -> Self { + SMatrix::new(|_| TYPE::zero()) + } + + fn is_zero(&self) -> bool { + *self == Self::zero() + } +} + +impl One for SMatrix +where + for<'b> TYPE: Clone + + Zero + + One + + PartialEq + + PartialOrd + + DivAssign + + SubAssign + + Neg + + Div + + SubAssign<&'b TYPE> + + DivAssign<&'b TYPE> + + MulAssign<&'b TYPE>, + for<'a> &'a TYPE: + Mul + Add + Div + Neg, +{ + fn one() -> Self { + , SLUDecomposition>>::identity(SIZE) + } +} + +impl<'a, TYPE, const SIZE: usize> + NumericalMatrixImpl<'a, TYPE, SPivot, SLUDecomposition> + for SMatrix +where + for<'b> TYPE: Clone + + 'a + + Zero + + One + + PartialEq + + PartialOrd + + DivAssign + + SubAssign + + Neg + + Div + + SubAssign<&'b TYPE> + + DivAssign<&'b TYPE> + + MulAssign<&'b TYPE>, + for<'c> &'c TYPE: + Mul + Add + Div + Neg, +{ + fn identity(_: usize) -> Self { + Self::new(|(i, j)| if i == j { TYPE::one() } else { TYPE::zero() }) + } +} + +impl<'a, TYPE, const SIZE: usize> NumericalMatrix, SLUDecomposition> + for SMatrix +where + for<'b> TYPE: Clone + + 'a + + Zero + + One + + PartialEq + + PartialOrd + + DivAssign + + SubAssign + + Neg + + Div + + SubAssign<&'b TYPE> + + DivAssign<&'b TYPE> + + MulAssign<&'b TYPE>, + for<'c> &'c TYPE: + Mul + Add + Div + Neg, +{ + fn identity(size: usize) -> Self { + NumericalMatrixImpl::, SLUDecomposition>::identity(size) + } + + fn tril(mut self) -> Self { + self.tril_impl(None); + self + } + + fn triu(mut self) -> Self { + self.triu_impl(None); + self + } + + fn tril_replace(mut self, diag_replacement: TYPE) -> Self { + self.tril_impl(Some(diag_replacement)); + self + } + + fn triu_replace(mut self, diag_replacement: TYPE) -> Self { + self.triu_impl(Some(diag_replacement)); + self + } + + fn gauss_jordan_high(mut self) -> Self { + let mut pivot = SPivot::::new(self.rows()); + NumericalMatrixImpl::, SLUDecomposition>::gauss_jordan_high( + &mut self, &mut pivot, + ); + self + } + + fn gauss_jordan_low(mut self) -> Self { + let mut pivot = SPivot::::new(self.rows()); + NumericalMatrixImpl::, SLUDecomposition>::gauss_jordan_low( + &mut self, &mut pivot, + ); + self + } + + fn det(self) -> TYPE { + let pivot = SPivot::new(self.rows()); + NumericalMatrixImpl::, SLUDecomposition>::det(self, pivot) + } + + fn invert(&mut self) -> Self { + let mut result = NumericalMatrix::identity(self.rows()); + let mut pivot = SPivot::new(self.rows()); + NumericalMatrixImpl::, SLUDecomposition>::invert( + self, + &mut pivot, + &mut result, + ); + result + } + + fn lup(mut self) -> Result, SingularMatrixError> { + let mut pivot = SPivot::new(self.columns()); + NumericalMatrixImpl::, SLUDecomposition>::lup( + &mut self, &mut pivot, + )?; + Ok(SLUDecomposition::::new(self, Some(pivot))) + } + + fn lu(mut self) -> Result, SingularMatrixError> { + NumericalMatrixImpl::, SLUDecomposition>::lu(&mut self)?; + Ok(SLUDecomposition::::new(self, None)) + } +} + +#[cfg(test)] +mod tests { + + use crate::LUDecomposition; + use crate::Matrix; + use crate::NumericalMatrix; + use crate::Rational; + use crate::SMatrix; + + use num_traits::One; + use num_traits::Zero; + + #[test] + fn test_new() { + let mut idx = 0; + let m1 = SMatrix::::new(|_| { + let res = idx; + idx += 1; + res + }); + idx = 0; + for (_, _, value) in m1 { + assert_eq!(idx, value); + idx += 1; + } + } + + #[test] + fn test_eq() { + let mut idx = 0; + let m1 = SMatrix::::new(|_| { + let res = idx; + idx += 1; + res + }); + idx = 0; + for (_, _, value) in m1 { + assert_eq!(idx, value); + idx += 1; + } + let m2 = m1.clone(); + assert!(m1 == m2); + } + + #[test] + fn test_sum() { + let mut idx = 0; + let m1 = SMatrix::::new(|_| { + let res = idx; + idx += 1; + res + }); + let m2 = SMatrix::::new(|_| { + let res = idx; + idx += 1; + res + }); + let m3 = m1 + m2; + for i in 0..m1.rows() { + for j in 0..m1.columns() { + assert_eq!(m3[(i, j)], m1[(i, j)] + m2[(i, j)]); + } + } + } + + #[test] + fn test_sub() { + let mut idx = 0; + let m1 = SMatrix::::new(|_| { + let res = idx; + idx += 1; + res + }); + let m2 = SMatrix::::new(|_| { + let res = idx; + idx += 1; + res + }); + let m3 = m1 - m2; + for i in 0..m1.rows() { + for j in 0..m1.columns() { + assert_eq!(m3[(i, j)], m1[(i, j)] - m2[(i, j)]); + } + } + } + + #[test] + fn test_mul() { + let mut idx = 0; + let m1 = SMatrix::::new(|_| { + let res = idx; + idx += 1; + res + }); + idx = 0; + let m2 = SMatrix::::new(|_| { + let res = idx; + idx += 1; + res + }); + let m3 = m1 * m2; + let expected = { + let mut idx = 0; + let values = [60, 70, 160, 195, 260, 320]; + SMatrix::::new(|_| { + let res = values[idx]; + idx += 1; + res + }) + }; + assert!(expected == m3); + } + + #[test] + fn test_one() { + let mut idx = 0; + let m1 = SMatrix::::new(|_| { + let res = idx; + idx += 1; + res + }); + let m2 = SMatrix::::one(); + let m3 = m1 * m2; + assert!(m1 == m3); + } + + #[test] + fn test_inv() { + let arr = [1, 2, 3, 4, 5, 6, 8, 7, 9]; + let mut idx = 0; + let m1 = SMatrix::, 3, 3>::new(|_| { + let res = Rational::new(arr[idx], i64::one()); + idx += 1; + res + }); + let m2 = NumericalMatrix::invert(&mut m1.clone()); + let m3 = m1 * m2; + assert!(m3.is_one()); + } + + #[test] + fn test_det() { + let arr = [1, 2, 3, 4, 5, 6, 8, 7, 9]; + let mut idx = 0; + let m1 = SMatrix::, 3, 3>::new(|_| { + let res = Rational::new(arr[idx], 1); + idx += 1; + res + }); + assert_eq!(Rational::new(-9, 1), NumericalMatrix::det(m1)); + } + + #[test] + fn test_tril() { + let arr = [1, 2, 3, 4, 5, 6, 8, 7, 9]; + let mut idx = 0; + let m1 = SMatrix::::new(|_| { + let res = arr[idx]; + idx += 1; + res + }); + let diag_value = 42; + let lt = m1.tril_replace(diag_value); + for (i, j, value) in lt { + if i == j { + assert_eq!(value, diag_value); + } else if i < j { + assert_eq!(value, i64::zero()) + } else { + assert_eq!(value, i64::from(arr[i * 3 + j])) + } + } + } + + #[test] + fn test_triu() { + let arr = [1, 2, 3, 4, 5, 6, 8, 7, 9]; + let mut idx = 0; + let m1 = SMatrix::::new(|_| { + let res = arr[idx]; + idx += 1; + res + }); + let diag_value = 42; + let lt = m1.triu_replace(diag_value); + for (i, j, value) in lt { + if i == j { + assert_eq!(value, diag_value); + } else if i > j { + assert_eq!(value, i64::zero()) + } else { + assert_eq!(value, i64::from(arr[i * 3 + j])) + } + } + } + + #[test] + fn test_lu_decomposition() { + let arr = [1, 2, 3, 4, 5, 6, 8, 7, 9]; + let mut idx = 0; + let orig = SMatrix::, 3, 3>::new(|_| { + let res = Rational::new(arr[idx], i64::one()); + idx += 1; + res + }); + let m1 = orig.clone(); + let lu = m1.lup().unwrap(); + let tril = lu.l(); + let triu = lu.u(); + let m2 = lu.pivot(tril * triu); + assert_eq!(orig, m2); + } + + #[test] + fn test_lu_invert() { + let arr = [1, 2, 3, 4, 5, 6, 8, 7, 9]; + let mut idx = 0; + let orig = SMatrix::, 3, 3>::new(|_| { + let res = Rational::new(arr[idx], i64::one()); + idx += 1; + res + }); + let m1 = orig.clone(); + let lud = m1.lup().unwrap(); + let m2 = lud.invert(); + let m3 = orig * m2; + assert!(m3.is_one()); + } + + #[test] + fn test_lu_det() { + let arr = [1, 2, 3, 4, 5, 6, 8, 7, 9]; + let mut idx = 0; + let m1 = SMatrix::, 3, 3>::new(|_| { + let res = Rational::new(arr[idx], 1); + idx += 1; + res + }); + assert_eq!(Rational::new(-9, 1), m1.lup().unwrap().det()); + } +} + +#[cfg(test)] +mod big_int_tests { + use num_bigint::BigInt; + + use crate::LUDecomposition; + use crate::Matrix; + use crate::NumericalMatrix; + use crate::Rational; + use crate::SMatrix; + + use num_traits::One; + use num_traits::Zero; + + #[test] + fn test_new() { + let mut idx = 0; + let m1 = SMatrix::::new(|_| { + let res = BigInt::from(idx); + idx += 1; + res + }); + idx = 0; + for (_, _, value) in m1 { + assert_eq!(BigInt::from(idx), value); + idx += 1; + } + } + + #[test] + fn test_eq() { + let mut idx = 0; + let m1 = SMatrix::::new(|_| { + let res = BigInt::from(idx); + idx += 1; + res + }); + idx = 0; + let m2 = m1.clone(); + assert!(m1 == m2); + for (_, _, value) in &m1 { + assert_eq!(BigInt::from(idx), *value); + idx += 1; + } + } + + #[test] + fn test_sum() { + let mut idx = 0; + let m1 = SMatrix::::new(|_| { + let res = BigInt::from(idx); + idx += 1; + res + }); + let m2 = SMatrix::::new(|_| { + let res = BigInt::from(idx); + idx += 1; + res + }); + let m3 = &m1 + &m2; + for i in 0..m1.rows() { + for j in 0..m1.columns() { + assert_eq!(m3[(i, j)], &m1[(i, j)] + &m2[(i, j)]); + } + } + } + + #[test] + fn test_sub() { + let mut idx = 0; + let m1 = SMatrix::::new(|_| { + let res = BigInt::from(idx); + idx += 1; + res + }); + let m2 = SMatrix::::new(|_| { + let res = BigInt::from(idx); + idx += 1; + res + }); + let m3 = &m1 - &m2; + for i in 0..m1.rows() { + for j in 0..m1.columns() { + assert_eq!(m3[(i, j)], &m1[(i, j)] - &m2[(i, j)]); + } + } + } + + #[test] + fn test_mul() { + let mut idx = 0; + let m1 = SMatrix::::new(|_| { + let res = BigInt::from(idx); + idx += 1; + res + }); + idx = 0; + let m2 = SMatrix::::new(|_| { + let res = BigInt::from(idx); + idx += 1; + res + }); + let m3 = m1 * m2; + let expected = { + let mut idx = 0; + let values = [60, 70, 160, 195, 260, 320]; + SMatrix::::new(|_| { + let res = BigInt::from(values[idx]); + idx += 1; + res + }) + }; + assert!(expected == m3); + } + + #[test] + fn test_one() { + let mut idx = 0; + let m1 = SMatrix::::new(|_| { + let res = BigInt::from(idx); + idx += 1; + res + }); + let m2 = SMatrix::::one(); + let m3 = &m1 * &m2; + assert!(m1 == m3); + } + + #[test] + fn test_inv() { + let arr = [1, 2, 3, 4, 5, 6, 8, 7, 9]; + let mut idx = 0; + let m1 = SMatrix::, 3, 3>::new(|_| { + let res = Rational::new(BigInt::from(arr[idx]), BigInt::one()); + idx += 1; + res + }); + let m2 = NumericalMatrix::invert(&mut m1.clone()); + let m3 = m1 * m2; + assert!(m3.is_one()); + } + + #[test] + fn test_det() { + let arr = [1, 2, 3, 4, 5, 6, 8, 7, 9]; + let mut idx = 0; + let m1 = SMatrix::, 3, 3>::new(|_| { + let res = Rational::new(BigInt::from(arr[idx]), BigInt::one()); + idx += 1; + res + }); + assert_eq!( + Rational::new(BigInt::from(9), -BigInt::one()), + NumericalMatrix::det(m1) + ); + } + + #[test] + fn test_tril() { + let arr = [1, 2, 3, 4, 5, 6, 8, 7, 9]; + let mut idx = 0; + let m1 = SMatrix::::new(|_| { + let res = BigInt::from(arr[idx]); + idx += 1; + res + }); + let diag_value = BigInt::from(42); + let lt = m1.tril_replace(diag_value.clone()); + for (i, j, value) in lt { + if i == j { + assert_eq!(value, diag_value); + } else if i < j { + assert_eq!(value, BigInt::zero()) + } else { + assert_eq!(value, BigInt::from(arr[i * 3 + j])) + } + } + } + + #[test] + fn test_triu() { + let arr = [1, 2, 3, 4, 5, 6, 8, 7, 9]; + let mut idx = 0; + let m1 = SMatrix::::new(|_| { + let res = BigInt::from(arr[idx]); + idx += 1; + res + }); + let diag_value = BigInt::from(42); + let lt = m1.triu_replace(diag_value.clone()); + for (i, j, value) in lt { + if i == j { + assert_eq!(value, diag_value); + } else if i > j { + assert_eq!(value, BigInt::zero()) + } else { + assert_eq!(value, BigInt::from(arr[i * 3 + j])) + } + } + } + + #[test] + fn test_lu_decomposition() { + let arr = [1, 2, 3, 4, 5, 6, 8, 7, 9]; + let mut idx = 0; + let orig = SMatrix::, 3, 3>::new(|_| { + let res = Rational::new(BigInt::from(arr[idx]), BigInt::one()); + idx += 1; + res + }); + let m1 = orig.clone(); + let lu = m1.lup().unwrap(); + let tril = lu.l(); + let triu = lu.u(); + let m2 = lu.pivot(tril * triu); + assert_eq!(orig, m2); + } + + #[test] + fn test_lu_invert() { + let arr = [1, 2, 3, 4, 5, 6, 8, 7, 9]; + let mut idx = 0; + let orig = SMatrix::, 3, 3>::new(|_| { + let res = Rational::new(BigInt::from(arr[idx]), BigInt::one()); + idx += 1; + res + }); + let m1 = orig.clone(); + let lu = m1.lup().unwrap(); + let m2 = lu.invert(); + let m3 = orig * m2; + assert!(m3.is_one()); + } + + #[test] + fn test_lu_det() { + let arr = [1, 2, 3, 4, 5, 6, 8, 7, 9]; + let mut idx = 0; + let m1 = SMatrix::, 3, 3>::new(|_| { + let res = Rational::new(BigInt::from(arr[idx]), BigInt::one()); + idx += 1; + res + }); + assert_eq!( + Rational::new(BigInt::from(-9), BigInt::one()), + m1.lup().unwrap().det() + ); + } +} diff --git a/src/spivot.rs b/src/spivot.rs new file mode 100644 index 0000000..b4c48ac --- /dev/null +++ b/src/spivot.rs @@ -0,0 +1,94 @@ +use std::ops::Index; +use std::ops::IndexMut; +use std::ops::Mul; + +use super::matrix::MatrixImpl; +use super::smatrix::SMatrix; +use super::Matrix; +use super::Pivot; + +pub struct SPivot { + data: [usize; SIZE], + permutations: usize, +} + +impl Index for SPivot { + fn index(&self, index: usize) -> &Self::Output { + &self.data[index] + } + + type Output = usize; +} + +impl IndexMut for SPivot { + fn index_mut(&mut self, index: usize) -> &mut >::Output { + &mut self.data[index] + } +} + +impl Pivot for SPivot { + fn swap(&mut self, i1: usize, i2: usize) { + self.data.swap(i1, i2); + self.permutations += 1; + } + + fn permutations(&self) -> usize { + self.permutations + } + + fn new(_: usize) -> Self { + SPivot { + data: std::array::from_fn(|i| i), + permutations: 0, + } + } +} + +impl Mul> + for SPivot +where + TYPE: Clone, +{ + fn mul(self, matrix: SMatrix) -> Self::Output { + let mut result = matrix.clone(); + let mut pclone = self.clone(); + for i in 0..ROWS { + while i != pclone[i] { + let j = pclone[i]; + result.swap_rows_pivot(&mut pclone, i, j); + } + } + result + } + + type Output = SMatrix; +} + +impl Mul> + for &SPivot +where + TYPE: Clone, +{ + fn mul(self, matrix: SMatrix) -> Self::Output { + let mut result = matrix.clone(); + let mut pclone = self.clone(); + for i in 0..matrix.rows() { + while i != pclone[i] { + let j = pclone[i]; + result.swap_rows_pivot(&mut pclone, i, j); + } + } + result + } + + type Output = SMatrix; +} + +impl Clone for SPivot { + fn clone(&self) -> Self { + SPivot { + data: self.data, + permutations: self.permutations, + } + } +}