haar_lib/linalg/mod_2/
lineq.rs

1//! 連立一次方程式$A \boldsymbol{x} = \boldsymbol{b} \pmod 2$を解く。
2//!
3//! # Problems
4//! - <https://judge.yosupo.jp/problem/system_of_linear_equations_mod_2>
5use crate::{ds::bitset::Bitset, linalg::mod_2::gaussian_elim::*};
6
7/// 連立一次方程式$A \boldsymbol{x} = \boldsymbol{b} \pmod 2$を解く。
8///
9/// ここで、$A$は$n \times m$の行列、$\boldsymbol{x}$は$m$行の縦ベクトル、$\boldsymbol{b}$は$n$行の縦ベクトル。
10///
11/// 連立方程式が解をもたないとき、`None`を返す。
12/// そうでなければ、`Some((sol, bases))`を返す。
13///
14/// `sol`は$m$行のベクトル、`bases`は`dim`個の$m$行のベクトルで、
15/// 連立方程式の解は、`bases`の要素の線型結合と`sol`の和で表される。
16pub fn lineq(mut a: Vec<Bitset>, b: Vec<bool>) -> Option<(Bitset, Vec<Bitset>)> {
17    let n = a.len();
18    assert_eq!(b.len(), n);
19
20    let Some(m) = a.first().map(|a| a.len()) else {
21        panic!("行数は0以上でなければならない。")
22    };
23    assert!(a.iter().all(|r| r.len() == m));
24
25    for (r, bi) in a.iter_mut().zip(b.iter()) {
26        r.push(*bi);
27    }
28
29    let (rank, mut a) = gaussian_elim(a);
30
31    let dim = m - rank;
32
33    let b: Vec<_> = a.iter_mut().map(|r| r.pop().unwrap()).collect();
34
35    if rank > 0 && a[rank - 1].count_ones() == 0 {
36        return None;
37    }
38
39    let mut index_zero = vec![];
40    let mut index_one = vec![];
41    {
42        let mut k = 0;
43        for ai in a.iter().take(rank) {
44            for j in (0..).take(m).skip(k) {
45                if ai.test(j) {
46                    index_one.push(j);
47                    k = j + 1;
48                    break;
49                }
50                index_zero.push(j);
51            }
52        }
53        for j in k..m {
54            index_zero.push(j);
55        }
56    }
57
58    assert_eq!(index_zero.len(), dim);
59    assert_eq!(index_one.len(), rank);
60
61    let mut sol = Bitset::new(m);
62    for (i, x) in b.into_iter().take(rank).enumerate() {
63        sol.set(index_one[i], x);
64    }
65
66    let mut bases = vec![Bitset::new(m); dim];
67    for i in 0..rank {
68        for (j, &k) in index_zero.iter().enumerate() {
69            bases[j].set(index_one[i], a[i].test(k));
70        }
71    }
72
73    for i in 0..dim {
74        bases[i].set(index_zero[i], true);
75    }
76
77    Some((sol, bases))
78}