haar_lib/linalg/mod_p/
pfaffian.rs

1//! $\mathbb{Z} / p \mathbb{Z}$($p$は素数)上の行列のパフィアン
2//!
3//! # References
4//! - <https://en.wikipedia.org/wiki/Pfaffian>
5//! # Problems
6//! - <https://judge.yosupo.jp/problem/pfaffian_of_matrix>
7
8use crate::misc::swap::swap_vv;
9use crate::num::ff::*;
10
11fn swap<T>(a: &mut [Vec<T>], i: usize, j: usize)
12where
13    T: FFElem,
14{
15    assert!(i < j);
16
17    for k in 0..i {
18        swap_vv(a, i, k, j, k);
19    }
20
21    for k in i + 1..j {
22        swap_vv(a, k, i, j, k);
23        a[k][i] = -a[k][i];
24        a[j][k] = -a[j][k];
25    }
26
27    for k in j + 1..a.len() {
28        swap_vv(a, k, i, k, j);
29    }
30
31    a[j][i] = -a[j][i];
32}
33
34fn add<T>(a: &mut [Vec<T>], s: T, i: usize, j: usize)
35where
36    T: FFElem,
37{
38    assert!(i < j);
39
40    for k in 0..i {
41        unsafe {
42            let t = *a.get_unchecked(i).get_unchecked(k);
43            *a.get_unchecked_mut(j).get_unchecked_mut(k) += t * s;
44        }
45    }
46
47    for k in i + 1..j {
48        unsafe {
49            let t = *a.get_unchecked(k).get_unchecked(i);
50            *a.get_unchecked_mut(j).get_unchecked_mut(k) -= t * s;
51        }
52    }
53
54    for k in j + 1..a.len() {
55        unsafe {
56            let t = *a.get_unchecked(k).get_unchecked(i);
57            *a.get_unchecked_mut(k).get_unchecked_mut(j) += t * s;
58        }
59    }
60}
61
62/// $\mathbb{Z} / p \mathbb{Z}$($p$は素数)上の行列のパフィアンを求める。
63///
64/// 入力の行列は、$n \times n$($n$は偶数)の[歪対称行列](https://en.wikipedia.org/wiki/Skew-symmetric_matrix)である。
65///
66/// **Time complexity** $O(n^3)$
67pub fn pfaffian<F>(mut a: Vec<Vec<F::Element>>, modulo: &F) -> F::Element
68where
69    F: FF,
70    F::Element: FFElem + std::fmt::Debug,
71{
72    let n = a.len();
73
74    assert_eq!(n % 2, 0);
75    assert!(a.iter().all(|r| r.len() == n));
76
77    for i in 0..n {
78        for j in 0..n {
79            assert_eq!(a[i][j], -a[j][i]);
80        }
81    }
82
83    for (i, r) in a.iter_mut().enumerate() {
84        r.truncate(i + 1);
85    }
86
87    let mut minus = false;
88    for i in (0..n).step_by(2) {
89        if a[i + 1][i].value() == 0 {
90            if let Some(j) = (i + 2..n).find(|&j| a[j][i].value() != 0) {
91                swap(&mut a, i + 1, j);
92                minus = !minus;
93            } else {
94                return modulo.zero();
95            }
96        }
97
98        assert_ne!(a[i + 1][i], modulo.zero());
99        let t = a[i + 1][i].inv();
100        for j in i + 2..n {
101            let c = -a[j][i] * t;
102            add(&mut a, c, i + 1, j);
103
104            let c = a[j][i + 1] * t;
105            add(&mut a, c, i, j);
106        }
107    }
108
109    let mut ret = modulo.one();
110
111    for i in (0..n).step_by(2) {
112        ret *= -a[i + 1][i];
113    }
114    if minus {
115        ret = -ret;
116    }
117
118    ret
119}
120
121#[cfg(test)]
122mod tests {
123    use super::*;
124    use crate::num::modint::*;
125
126    fn check(a: Vec<Vec<i64>>, m: u32, ans: u32) {
127        let m = ModIntBuilder::new(m);
128        let a = a
129            .into_iter()
130            .map(|b| b.into_iter().map(|x| m.from_i64(x)).collect())
131            .collect();
132        assert_eq!(pfaffian(a, &m).value(), ans);
133    }
134
135    #[test]
136    fn test() {
137        check(
138            vec![
139                vec![0, 1, 2, 3],
140                vec![-1, 0, 4, 5],
141                vec![-2, -4, 0, 6],
142                vec![-3, -5, -6, 0],
143            ],
144            998244353,
145            8,
146        );
147        check(vec![vec![0, 1], vec![-1, 0]], 998244353, 1);
148    }
149}
haar_lib/linalg/mod_p/pfaffian.rs

haar_lib/linalg/mod_p/
pfaffian.rs
