sim_lib_discrete_algebra/
power.rs

1//! Matrix power by square-and-multiply. `A^k[i][j]` over the counting semiring
2//! counts walks of length `k`; over min-plus it bounds `k`-edge shortest paths.
3
4use crate::error::AlgebraError;
5use crate::matrix::{AlgebraLimits, Matrix};
6use crate::semiring::Semiring;
7
8impl<S: Semiring> Matrix<S> {
9    /// Raise a square matrix to the `k`-th power over its semiring.
10    ///
11    /// `k == 0` returns the identity. Returns [`AlgebraError::ShapeMismatch`]
12    /// for non-square input and [`AlgebraError::LimitExceeded`] when the
13    /// dimension exceeds `limits.max_dim`.
14    pub fn power(&self, k: usize, limits: AlgebraLimits) -> Result<Self, AlgebraError> {
15        self.validate()?;
16        if !self.is_square() {
17            return Err(AlgebraError::ShapeMismatch(format!(
18                "power requires a square matrix, got {}x{}",
19                self.rows, self.cols
20            )));
21        }
22        if self.rows > limits.max_dim {
23            return Err(AlgebraError::LimitExceeded(format!(
24                "power: dimension {} exceeds max_dim {}",
25                self.rows, limits.max_dim
26            )));
27        }
28        let n = self.rows;
29        let mut result = Matrix::try_identity_with_limits(n, limits)?;
30        let mut base = self.clone();
31        let mut e = k;
32        while e > 0 {
33            if e & 1 == 1 {
34                result = result.matmul(&base)?;
35            }
36            e >>= 1;
37            if e > 0 {
38                base = base.matmul(&base)?;
39            }
40        }
41        Ok(result)
42    }
43}
44
45#[cfg(test)]
46mod tests {
47    use super::*;
48    use crate::Counting;
49
50    #[test]
51    fn power_zero_is_identity() {
52        let a = Matrix::from_rows(vec![
53            vec![Counting::from_u64(1), Counting::from_u64(1)],
54            vec![Counting::from_u64(0), Counting::from_u64(1)],
55        ])
56        .unwrap();
57        assert_eq!(
58            a.power(0, AlgebraLimits::default()).unwrap(),
59            Matrix::identity(2)
60        );
61    }
62
63    #[test]
64    fn power_counts_walks() {
65        // Triangle cycle 0->1->2->0. A^2[i][j] = number of length-2 walks.
66        let mut a = Matrix::new(3, 3);
67        a.set(0, 1, Counting::from_u64(1)).unwrap();
68        a.set(1, 2, Counting::from_u64(1)).unwrap();
69        a.set(2, 0, Counting::from_u64(1)).unwrap();
70        let a2 = a.power(2, AlgebraLimits::default()).unwrap();
71        // 0->1->2 is the only length-2 walk from 0, ending at 2.
72        assert_eq!(a2.get(0, 2).unwrap(), &Counting::from_u64(1));
73        assert_eq!(a2.get(0, 0).unwrap(), &Counting::from_u64(0));
74        // A^3 returns to the start: exactly one closed length-3 walk per node.
75        let a3 = a.power(3, AlgebraLimits::default()).unwrap();
76        assert_eq!(a3.get(0, 0).unwrap(), &Counting::from_u64(1));
77    }
78
79    #[test]
80    fn power_non_square_fails() {
81        let a: Matrix<Counting> = Matrix::new(2, 3);
82        assert!(matches!(
83            a.power(2, AlgebraLimits::default()),
84            Err(AlgebraError::ShapeMismatch(_))
85        ));
86    }
87}
sim_lib_discrete_algebra/power.rs

sim_lib_discrete_algebra/
power.rs
