1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

//! 行列累乗

#![allow(clippy::needless_range_loop)]

use crate::utils::num_traits::{One, Zero};

#[derive(Debug, Clone)]
pub struct Matrix<const N: usize, T>(pub [[T; N]; N]);

impl<const N: usize, T> Matrix<N, T>
where
    T: One + Zero + Copy,
{
    /// 単位行列を返す
    pub fn id() -> Self {
        let mut res = [[T::zero(); N]; N];
        for i in 0..N {
            res[i][i] = T::one();
        }
        Self(res)
    }

    /// ベクトル`x`と行列`A`について、`Ax`を返す
    pub fn apply(&self, v: &[T; N]) -> [T; N] {
        let mut res = [T::zero(); N];
        for i in 0..N {
            for j in 0..N {
                res[i] = res[i] + self.0[i][j] * v[j];
            }
        }
        res
    }

    /// 行列の累乗を返す（繰り返し2乗法）
    pub fn pow(&self, mut p: usize) -> Self {
        let mut res = Self::id();
        let mut tmp = self.clone();
        while p > 0 {
            if p & 1 == 1 {
                res = tmp.dot(&res);
            }
            tmp = tmp.dot(&tmp);
            p >>= 1;
        }
        res
    }

    /// 行列のドット積
    pub fn dot(&self, other: &Self) -> Self {
        let mut res = [[T::zero(); N]; N];
        for i in 0..N {
            for j in 0..N {
                for k in 0..N {
                    res[i][j] = res[i][j] + self.0[i][k] * other.0[k][j];
                }
            }
        }
        Self(res)
    }
}