haar_lib/linalg/mod_p/
determinant.rs

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