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] 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) + } +}