haar_lib/linalg/mod_p/
determinant.rs

1//! 行列式 (mod 素数)
2use crate::num::{ff::FFElem, one_zero::*};
3
4/// 素数mod p上での行列式を求める。
5///
6/// **Time complexity** $O(n^3)$
7pub fn determinant<T>(mut a: Vec<Vec<T>>) -> T
8where
9    T: FFElem + Copy + Zero + One,
10{
11    let n = a.len();
12
13    assert!(a.iter().all(|r| r.len() == n));
14
15    let mut s = 0;
16    for i in 0..n {
17        if a[i][i].value() == 0 {
18            if let Some(j) = (i + 1..n).find(|&j| a[j][i].value() != 0) {
19                a.swap(i, j);
20                s ^= 1;
21            } else {
22                return T::zero();
23            }
24        }
25
26        let d = a[i][i].inv();
27        let ai = a.swap_remove(i);
28
29        for aj in a.iter_mut().skip(i) {
30            let t = aj[i] * d;
31            for (x, y) in aj.iter_mut().zip(ai.iter()) {
32                *x -= *y * t;
33            }
34        }
35
36        a.push(ai);
37        a.swap(i, n - 1);
38    }
39
40    let mut ret = T::one();
41    for (i, a) in a.into_iter().enumerate() {
42        ret *= a[i];
43    }
44
45    if s == 1 {
46        ret = -ret;
47    }
48
49    ret
50}
51
52#[cfg(test)]
53mod tests {
54    use super::*;
55    use crate::num::const_modint::*;
56
57    fn convert<U, T>(a: Vec<Vec<T>>) -> Vec<Vec<U>>
58    where
59        U: From<T>,
60    {
61        a.into_iter()
62            .map(|b| b.into_iter().map(From::from).collect())
63            .collect()
64    }
65
66    #[test]
67    fn test() {
68        const P: u32 = 998244353;
69
70        let a = vec![vec![3, 1, 4], vec![1, 5, 9], vec![2, 6, 5]];
71        let a = convert::<ConstModInt<P>, _>(a);
72        assert_eq!(determinant(a).value(), 998244263);
73
74        let a = vec![vec![1, 2, 3], vec![4, 5, 6], vec![7, 8, 9]];
75        let a = convert::<ConstModInt<P>, _>(a);
76        assert_eq!(determinant(a).value(), 0);
77
78        let a = vec![vec![0, 1], vec![1, 0]];
79        let a = convert::<ConstModInt<P>, _>(a);
80        assert_eq!(determinant(a).value(), 998244352);
81    }
82}