From ab85766b5a1a58e3a6a04de4e915315b70c00459 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Violeta=20Hern=C3=A1ndez?= Date: Mon, 5 Apr 2021 12:17:49 -0500 Subject: [PATCH 1/5] Added pow function. I'll try adding some unchecked and in place variants soon. --- src/linalg/mod.rs | 2 ++ src/linalg/pow.rs | 52 +++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 54 insertions(+) create mode 100644 src/linalg/pow.rs diff --git a/src/linalg/mod.rs b/src/linalg/mod.rs index 075d91c2..0c272494 100644 --- a/src/linalg/mod.rs +++ b/src/linalg/mod.rs @@ -19,6 +19,7 @@ pub mod householder; mod inverse; mod lu; mod permutation_sequence; +mod pow; mod qr; mod schur; mod solve; @@ -41,6 +42,7 @@ pub use self::full_piv_lu::*; pub use self::hessenberg::*; pub use self::lu::*; pub use self::permutation_sequence::*; +pub use self::pow::*; pub use self::qr::*; pub use self::schur::*; pub use self::svd::*; diff --git a/src/linalg/pow.rs b/src/linalg/pow.rs new file mode 100644 index 00000000..8769b144 --- /dev/null +++ b/src/linalg/pow.rs @@ -0,0 +1,52 @@ +//! This module provides the matrix exponential (pow) function to square matrices. + +use std::ops::DivAssign; + +use crate::{allocator::Allocator, DefaultAllocator, DimMin, DimMinimum, MatrixN}; +use num::PrimInt; +use simba::scalar::ComplexField; + +impl MatrixN +where + D: DimMin, + DefaultAllocator: Allocator + + Allocator<(usize, usize), DimMinimum> + + Allocator + + Allocator + + Allocator, +{ + /// Raises a matrix to an integer power using exponentiation by squares. + /// Returns `None` only when the matrix is non-invertible and raised to a + /// negative power. + pub fn pow(&self, mut e: T) -> Option { + let zero = T::zero(); + + if e == zero { + let mut i = self.clone(); + i.fill_with_identity(); + return Some(i); + } + + let mut acc; + if e < zero { + acc = self.clone().try_inverse()?; + } else { + acc = self.clone(); + } + + let one = T::one(); + let two = T::from(2u8).unwrap(); + let mut multiplier = acc.clone(); + + while e != zero { + if e % two == one { + acc *= &multiplier; + } + + e /= two; + multiplier *= multiplier.clone(); + } + + Some(acc) + } +} From 06b657ad49adc35d8ed9feb0713e1bf293e85762 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Violeta=20Hern=C3=A1ndez?= Date: Mon, 5 Apr 2021 12:32:12 -0500 Subject: [PATCH 2/5] Added pow_mut. Actually, I think this will do. --- src/linalg/pow.rs | 46 ++++++++++++++++++++++++++-------------------- 1 file changed, 26 insertions(+), 20 deletions(-) diff --git a/src/linalg/pow.rs b/src/linalg/pow.rs index 8769b144..27cab313 100644 --- a/src/linalg/pow.rs +++ b/src/linalg/pow.rs @@ -2,51 +2,57 @@ use std::ops::DivAssign; -use crate::{allocator::Allocator, DefaultAllocator, DimMin, DimMinimum, MatrixN}; +use crate::{allocator::Allocator, DefaultAllocator, DimMin, MatrixN}; use num::PrimInt; use simba::scalar::ComplexField; impl MatrixN where D: DimMin, - DefaultAllocator: Allocator - + Allocator<(usize, usize), DimMinimum> - + Allocator - + Allocator - + Allocator, + DefaultAllocator: Allocator, { - /// Raises a matrix to an integer power using exponentiation by squares. - /// Returns `None` only when the matrix is non-invertible and raised to a - /// negative power. - pub fn pow(&self, mut e: T) -> Option { + /// Attempts to raise this matrix to an integer power in-place. Returns + /// `false` and leaves `self` untouched if the power is negative and the + /// matrix is non-invertible. + pub fn pow_mut(&mut self, mut e: T) -> bool { let zero = T::zero(); if e == zero { - let mut i = self.clone(); - i.fill_with_identity(); - return Some(i); + self.fill_with_identity(); + return true; } - let mut acc; if e < zero { - acc = self.clone().try_inverse()?; - } else { - acc = self.clone(); + if !self.try_inverse_mut() { + return false; + } } let one = T::one(); let two = T::from(2u8).unwrap(); - let mut multiplier = acc.clone(); + let mut multiplier = self.clone(); while e != zero { if e % two == one { - acc *= &multiplier; + *self *= &multiplier; } e /= two; multiplier *= multiplier.clone(); } - Some(acc) + true + } + + /// Raise this matrix to an integer power. Returns `None` only if the power + /// is negative and the matrix is non-invertible. + pub fn pow(&self, e: T) -> Option { + let mut clone = self.clone(); + + if clone.pow_mut(e) { + Some(clone) + } else { + None + } } } From 15a63cb8927f729ee4f9ea4f356cdb6a118b6f61 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Violeta=20Hern=C3=A1ndez?= Date: Fri, 9 Apr 2021 23:43:59 -0500 Subject: [PATCH 3/5] Memory improvements, extra comments. The result of `multiplier ^ 2` is now written into a single buffer. --- src/linalg/pow.rs | 49 +++++++++++++++++++++++++++++++++++++---------- 1 file changed, 39 insertions(+), 10 deletions(-) diff --git a/src/linalg/pow.rs b/src/linalg/pow.rs index 27cab313..5ed3ae7f 100644 --- a/src/linalg/pow.rs +++ b/src/linalg/pow.rs @@ -10,18 +10,38 @@ impl MatrixN where D: DimMin, DefaultAllocator: Allocator, + DefaultAllocator: Allocator, { - /// Attempts to raise this matrix to an integer power in-place. Returns - /// `false` and leaves `self` untouched if the power is negative and the - /// matrix is non-invertible. + /// Computes the square of this matrix and writes it into a given buffer. + fn square_buf(&mut self, buf: &mut Self) { + // We unroll the first iteration to avoid new_uninitialized. + let mut aux_col = self.column(0).clone_owned(); + aux_col = &*self * aux_col; + buf.column_mut(0).copy_from(&aux_col); + + // We multiply the matrix by its i-th column, + for i in 1..self.ncols() { + aux_col.copy_from(&self.column(i)); + aux_col = &*self * aux_col; + self.column_mut(i).copy_from(&aux_col); + } + } + + /// Attempts to raise this matrix to an integral power `e` in-place. If this + /// matrix is non-invertible and `e` is negative, it leaves this matrix + /// untouched and returns `false`. Otherwise, it returns `true` and + /// overwrites this matrix with the result. pub fn pow_mut(&mut self, mut e: T) -> bool { let zero = T::zero(); + // A matrix raised to the zeroth power is just the identity. if e == zero { self.fill_with_identity(); return true; } + // If e is negative, we compute the inverse matrix, then raise it to the + // power of -e. if e < zero { if !self.try_inverse_mut() { return false; @@ -30,22 +50,31 @@ where let one = T::one(); let two = T::from(2u8).unwrap(); - let mut multiplier = self.clone(); - while e != zero { + // We use the buffer to hold the result of multiplier ^ 2, thus avoiding + // extra allocations. + let mut multiplier = self.clone(); + let mut buf = self.clone(); + + // Exponentiation by squares. + loop { if e % two == one { *self *= &multiplier; } e /= two; - multiplier *= multiplier.clone(); - } + multiplier.square_buf(&mut buf); + multiplier.copy_from(&buf); - true + if e == zero { + return true; + } + } } - /// Raise this matrix to an integer power. Returns `None` only if the power - /// is negative and the matrix is non-invertible. + /// Attempts to raise this matrix to an integral power `e`. If this matrix + /// is non-invertible and `e` is negative, it returns `None`. Otherwise, it + /// returns the result as a new matrix. Uses exponentiation by squares. pub fn pow(&self, e: T) -> Option { let mut clone = self.clone(); From 81f2fc38d7f7260b9dd1d84b168aaa247cbd69ad Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Violeta=20Hern=C3=A1ndez?= Date: Fri, 9 Apr 2021 23:59:22 -0500 Subject: [PATCH 4/5] Use mul_to instead of square_buf Didn't realize that this was something that was already implemented. --- src/linalg/pow.rs | 17 +---------------- 1 file changed, 1 insertion(+), 16 deletions(-) diff --git a/src/linalg/pow.rs b/src/linalg/pow.rs index 5ed3ae7f..4b73784b 100644 --- a/src/linalg/pow.rs +++ b/src/linalg/pow.rs @@ -12,21 +12,6 @@ where DefaultAllocator: Allocator, DefaultAllocator: Allocator, { - /// Computes the square of this matrix and writes it into a given buffer. - fn square_buf(&mut self, buf: &mut Self) { - // We unroll the first iteration to avoid new_uninitialized. - let mut aux_col = self.column(0).clone_owned(); - aux_col = &*self * aux_col; - buf.column_mut(0).copy_from(&aux_col); - - // We multiply the matrix by its i-th column, - for i in 1..self.ncols() { - aux_col.copy_from(&self.column(i)); - aux_col = &*self * aux_col; - self.column_mut(i).copy_from(&aux_col); - } - } - /// Attempts to raise this matrix to an integral power `e` in-place. If this /// matrix is non-invertible and `e` is negative, it leaves this matrix /// untouched and returns `false`. Otherwise, it returns `true` and @@ -63,7 +48,7 @@ where } e /= two; - multiplier.square_buf(&mut buf); + multiplier.mul_to(&multiplier, &mut buf); multiplier.copy_from(&buf); if e == zero { From 341091f64788da1a06b23f013cff91df6a27ee58 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Violeta=20Hern=C3=A1ndez?= Date: Sat, 10 Apr 2021 00:12:26 -0500 Subject: [PATCH 5/5] `pow_mut` now returns `Result`. --- src/linalg/pow.rs | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/linalg/pow.rs b/src/linalg/pow.rs index 4b73784b..edcf55dd 100644 --- a/src/linalg/pow.rs +++ b/src/linalg/pow.rs @@ -14,22 +14,23 @@ where { /// Attempts to raise this matrix to an integral power `e` in-place. If this /// matrix is non-invertible and `e` is negative, it leaves this matrix - /// untouched and returns `false`. Otherwise, it returns `true` and + /// untouched and returns `Err(())`. Otherwise, it returns `Ok(())` and /// overwrites this matrix with the result. - pub fn pow_mut(&mut self, mut e: T) -> bool { + #[must_use] + pub fn pow_mut(&mut self, mut e: T) -> Result<(), ()> { let zero = T::zero(); // A matrix raised to the zeroth power is just the identity. if e == zero { self.fill_with_identity(); - return true; + return Ok(()); } // If e is negative, we compute the inverse matrix, then raise it to the // power of -e. if e < zero { if !self.try_inverse_mut() { - return false; + return Err(()); } } @@ -52,7 +53,7 @@ where multiplier.copy_from(&buf); if e == zero { - return true; + return Ok(()); } } } @@ -63,10 +64,9 @@ where pub fn pow(&self, e: T) -> Option { let mut clone = self.clone(); - if clone.pow_mut(e) { - Some(clone) - } else { - None + match clone.pow_mut(e) { + Ok(()) => Some(clone), + Err(()) => None, } } }