haar_lib/linalg/mod_m/
determinant.rs

1//! $\mathbb{Z} / m \mathbb{Z}$上の行列式
2use crate::num::zz::{ZZElem, ZZ};
3
4/// $\mathbb{Z} / m \mathbb{Z}$上で行列式を求める。
5pub fn determinant<R>(mut a: Vec<Vec<R::Element>>, ring: &R) -> R::Element
6where
7    R: ZZ + Copy,
8    R::Element: ZZElem,
9{
10    let n = a.len();
11
12    assert!(a.iter().all(|r| r.len() == n));
13
14    let mut minus = false;
15
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                minus = !minus;
21            } else {
22                return ring.zero();
23            }
24        }
25
26        let mut ai = a.swap_remove(i);
27
28        for aj in a.iter_mut().skip(i) {
29            let t = aj[i];
30            if t.value() == 0 {
31                continue;
32            }
33
34            loop {
35                if aj[i].value() == 0 {
36                    break;
37                }
38
39                if ai[i].value() > aj[i].value() {
40                    std::mem::swap(&mut ai, aj);
41                    minus = !minus;
42                }
43
44                let t = ring.from_u64((aj[i].value() / ai[i].value()) as u64);
45
46                for (x, y) in aj.iter_mut().skip(i).zip(ai.iter().skip(i)) {
47                    *x -= *y * t;
48                }
49            }
50        }
51
52        a.push(ai);
53        a.swap(i, n - 1);
54    }
55
56    let mut ret = ring.one();
57    for (i, ai) in a.into_iter().enumerate() {
58        ret *= ai[i];
59    }
60
61    if minus {
62        ret = -ret;
63    }
64
65    ret
66}
haar_lib/linalg/mod_m/determinant.rs

haar_lib/linalg/mod_m/
determinant.rs
