nalgebra/linalg/
pow.rs

1//! This module provides the matrix exponential (pow) function to square matrices.
2
3use crate::{
4    allocator::Allocator,
5    storage::{Storage, StorageMut},
6    DefaultAllocator, DimMin, Matrix, OMatrix, Scalar,
7};
8use num::{One, Zero};
9use simba::scalar::{ClosedAddAssign, ClosedMulAssign};
10
11impl<T, D, S> Matrix<T, D, D, S>
12where
13    T: Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign,
14    D: DimMin<D, Output = D>,
15    S: StorageMut<T, D, D>,
16    DefaultAllocator: Allocator<D, D> + Allocator<D>,
17{
18    /// Raises this matrix to an integral power `exp` in-place.
19    pub fn pow_mut(&mut self, mut exp: u32) {
20        // A matrix raised to the zeroth power is just the identity.
21        if exp == 0 {
22            self.fill_with_identity();
23        } else if exp > 1 {
24            // We use the buffer to hold the result of multiplier^2, thus avoiding
25            // extra allocations.
26            let mut x = self.clone_owned();
27            let mut workspace = self.clone_owned();
28
29            if exp % 2 == 0 {
30                self.fill_with_identity();
31            } else {
32                // Avoid an useless multiplication by the identity
33                // if the exponent is odd.
34                exp -= 1;
35            }
36
37            // Exponentiation by squares.
38            loop {
39                if exp % 2 == 1 {
40                    self.mul_to(&x, &mut workspace);
41                    self.copy_from(&workspace);
42                }
43
44                exp /= 2;
45
46                if exp == 0 {
47                    break;
48                }
49
50                x.mul_to(&x, &mut workspace);
51                x.copy_from(&workspace);
52            }
53        }
54    }
55}
56
57impl<T, D, S: Storage<T, D, D>> Matrix<T, D, D, S>
58where
59    T: Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign,
60    D: DimMin<D, Output = D>,
61    S: StorageMut<T, D, D>,
62    DefaultAllocator: Allocator<D, D> + Allocator<D>,
63{
64    /// Raise this matrix to an integral power `exp`.
65    #[must_use]
66    pub fn pow(&self, exp: u32) -> OMatrix<T, D, D> {
67        let mut result = self.clone_owned();
68        result.pow_mut(exp);
69        result
70    }
71}