(False)faces | Cerc2009 积和式模 4

题目链接

给定一个两边各有 \(n\) 个点的二分图,判断完美匹配的个数是否是 \(4\) 的倍数。

\(n \le 300\)

完美匹配的个数即积和式。

来自论文的算法:求积和式模 \(4\) 的余数。

积和式定义:考虑所有选 \(n\) 个不同行、不同列的数的方案,对每个方案中 \(n\) 个数的乘积求和。 \[ \text{perm} A = \sum_{p}A_{i,p_i} \]\(A\)\(01\) 矩阵时,有 \[ \text{prem} A = (-1)^n\sum_{x \in \{0,1\}^n}(-1)^{x1+x2+\cdots+x_n}\prod_{i=1}^n(Ax)_i \] 证明可以考虑容斥:钦定某些列没有数被选。

观察到式子中间有一个 \(\prod\) ,由于我们要求这个东西模 \(4\) 的余数,如果它要造成贡献的话,\(Ax\) 就至多有一个位置模 \(2\)\(0\)

考虑 \(Ax\) 每一项模 \(2\) 的余数,由于至多只能有一个 \(0\),因此可以枚举这个东西的取值,它只有 \(n+1\) 种。

对于每种取值,通过高斯消元解出满足条件的所有 \(x\),再将每一组 \(x\) 代入刚刚的式子中求出答案。

问题是,合法的 \(x\) 的个数可能很大,因为需要枚举自由元的取值。

对这个矩阵做一些变换,记为 \(A'\)\[ \begin{bmatrix} A_{1,1}&A_{1,2}&\cdots &A_{1,n}&v_1\\ A_{2,1}&A_{2,2}&\cdots &A_{2,n}&v_2\\ \vdots&&&&\vdots\\ A_{n,1}&A_{n,2}&\cdots &A_{n,n}&v_n\\ 0&0&\cdots&0&1 \end{bmatrix} \] 显然 \(\text{perm} A' = \text{perm} A\)

随机选取 \(v\),则期望 \(O(1)\) 组解:一组 \(x\) 要有贡献,\(x_{n+1}=1\)(注意这一点,不能枚举 \(x_{n+1}=0\),否则复杂度会假),因此翻转 \(v\) 中的任意一位,\(A'x\) 也会翻转对应的一位,所以对于一个确定的 \(A'x\),每一组 \(x\) 成为它的解的概率都是 \(\frac 1{2^n}\)

直接写复杂度是 \(O(n^4)\),使用 bitset 优化,复杂度为 \(O(\frac {n^4}{\omega})\)

笔者想出了一个复杂度更好的算法。

引理:当 \(n\) 阶方阵 \(A\)\(2\) 意义下的秩小于等于 \(n-3\) 时(换句话就是模 \(2\) 意义下高斯消元后有至少 \(3\) 个自由元),有 \(\text{perm} A \equiv 0\pmod 4\)

证明留给读者。

把这种情况判掉就不需要随机化了,因为至多有 \(2\) 个自由元,也就至多有 \(4\) 个解。

之前算法的瓶颈在于枚举 \(Ax\) 后需要 \(O(n^3)\) 求解 \(x\) 并计算贡献,这可以分为四个部分:

  • 高斯消元成上三角矩阵。
  • 判断是否有解(系数全 \(0\) 行的常数项是否为 \(0\))。
  • 枚举每个自由元的取值,并算出每个非自由元。
  • \(x\) 代入容斥式子计算贡献。

注意每次高斯消元的矩阵都是 \(A\),只有常数项不同,因此可以考虑像矩阵求逆一样在 \(A\) 右边放一个单位矩阵进行一次消元,之后就可以 \(O(n^2)\) 地由初始常数项直接计算消元后的常数项。

后三个部分的复杂度都不高于 \(O(n^2)\),对于单个 \(Ax\) 就可以做到 \(O(n^2)\) 的复杂度,总复杂度 \(O(n^3)\)

除高斯消元之外的部分(计算每个 \(Ax\) 的贡献)还可以优化,虽然改变不了复杂度,但实际速度可以快很多倍。

实际上并不需要每次 \(O(n^2)\) 计算消元后的常数项。考虑翻转 \((Ax)_i\),那么常数项中和 \((Ax)_i\) 有关的位也会翻转,我们要求解的 \(Ax\) 都是可以通过全 \(1\) 翻转至多一位得到的,预处理 \(Ax\)\(1\) 的常数项,就可以 \(O(n)\) 得到其他 \(Ax\) 的常数项。

判断是否有解直接判断即可,由于至多有 \(2\) 个自由元,复杂度 \(O(1)\)

求解 \(x\) 也不需要先枚举每个自由元的取值,再算出每个非自由元。考虑消元后的常数项翻转了一些位,如果保持自由元取值不变,\(x\) 也就是翻转了对应的位,在维护常数项的时候顺便维护 \(x\)\(O(n)\) 的。

\(x\) 代入容斥式子计算贡献是 \(O(n^2)\) 的,但是我们只关心贡献模 \(4\),如果 \((Ax)_i\equiv 0\pmod 2\),那么 \(\sum_{x \in \{0,1\}^n}(-1)^{x1+x2+\cdots+x_n}\prod_{i=1}^n(Ax)_i\) 就必然是偶数了,只需要知道它是否是 \(4\) 的倍数,检验一下 \((Ax)_i\) 即可。最后对 \(Ax\)\(1\) 暴力计算即可。

使用 bitset 优化,复杂度 \(O(\frac {n^3}{\omega})\),除了读入和高斯消元外,其他部分复杂度不高于 \(O(\frac {n^2}{\omega})\)

当计算积和式模 \(8\) 时,其他部分复杂度达到 \(O(\frac {n^3}{\omega})\),总复杂度仍然是 \(O(\frac {n^3}{\omega})\)

\(O(\frac {n^4}{\omega})\) 代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
#include <bits/stdc++.h>
#define rep(i, l, r) for(int i = (l); i <= (r); i++)
#define per(i, r, l) for(int i = (r); i >= (l); i--)
#define mem(a, b) memset(a, b, sizeof a)
#define For(i, l, r) for(int i = (l), i##e = (r); i < i##e; i++)
#define pb push_back
#pragma GCC target("popcnt")

using namespace std;

const int N = 305;

int T, n;
char s[N], ans;
bitset<320> A[N], bas[N], x;

void dfs(int i) {
if(i < 0) {
char prod = 1;
rep(j, 0, n) prod *= (A[j] & x).count();
ans += x.count() & 1 ? -prod : prod;
} else if(bas[i][i]) dfs(i - 1);
else rep(j, 0, 1) {
x[i] = j, dfs(i - 1);
rep(k, 0, i) x[k] = x[k] ^ bas[k][i];
}
}
void calc() {
rep(i, 0, n) bas[i].reset();
rep(i, 0, n) {
auto nw = A[i];
rep(j, 0, n) if(nw[j]) {
if(bas[j][j]) nw ^= bas[j];
else { bas[j] = nw; break; }
}
if(nw[n + 1] && nw.count() == 1) return;
}
per(i, n, 0) if(bas[i][i]) rep(j, i + 1, n) if(bas[j][j] & bas[i][j]) bas[i] ^= bas[j];
rep(i, 0, n) x[i] = bas[i][n + 1];
dfs(n);
}
void solve() {
cin >> n, ans = 0, x.reset();
rep(i, 0, n) A[i].reset();
For(i, 0, n) {
cin >> s, A[i][n] = rand() & 1, A[i][n + 1] = 1;
For(j, 0, n) A[i][j] = s[j] - 48;
}
A[n][n] = A[n][n + 1] = 1;
rep(i, 0, n) A[i][n + 1] = i == n, calc(), A[i][n + 1] = 1;
cout << (ans & 3 ? "NO\n" : "YES\n");
}
int main() {
cin.tie(0)->sync_with_stdio(0);
for(cin >> T; T--; solve());
}

\(O(\frac {n^3}{\omega})\) 代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
#include <bits/stdc++.h>
#define rep(i, l, r) for(int i = (l); i <= (r); i++)
#define per(i, r, l) for(int i = (r); i >= (l); i--)
#define mem(a, b) memset(a, b, sizeof a)
#define For(i, l, r) for(int i = (l), i##e = (r); i < i##e; i++)
#define pb push_back
#pragma GCC target("popcnt")

using namespace std;

const int N = 305;

int T, n, idx;
char ans, s[N], x[N];
bitset<640> a[N], bas[N], fre[3], sol[5], X;

void dfs(int i) {
if(!i) {
sol[++idx].reset();
rep(i, 0, n) sol[idx][i] = x[i];
} else if(bas[i][i]) dfs(i - 1);
else {
x[i] = 0, dfs(i - 1);
rep(k, 0, i) if(bas[k][i]) x[k] ^= 1;
x[i] = 1, dfs(i - 1);
rep(k, 0, i) if(bas[k][i]) x[k] ^= 1;
}
}
void solve() {
cin >> n;
rep(i, 1, n) a[i].reset(), bas[i].reset();
rep(i, 1, n) {
cin >> s, a[i][n + i] = a[i][n + n + 2] = 1;
rep(j, 1, n) a[i][j] = s[j - 1] - 48;
}
int p = 0;
rep(i, 1, n) {
auto nw = a[i];
rep(j, 1, n) if(nw[j]) {
if(bas[j][j]) nw ^= bas[j];
else { bas[j] = nw; goto L; }
}
if(p == 2) return cout << "YES\n", void();
fre[++p] = nw; L:;
}
if(p == 0) return cout << "NO\n", void();
per(i, n, 1) if(bas[i][i]) rep(j, i + 1, n) if(bas[j][j] & bas[i][j]) bas[i] ^= bas[j];
rep(i, 1, n) x[i] = bas[i][n + n + 2];
idx = 0, X.reset(), ans = 0, dfs(n);
bitset<640> y;
rep(k, n + 1, n + n + 1) {
bool ban = 0;
rep(i, 1, p) ban |= fre[i][k] != fre[i][n + n + 2];
if(ban) continue;
rep(i, 1, n) y[i] = bas[i][k];
rep(i, 1, idx) {
X = sol[i] ^ y;
if(k > n + n) {
char prod = 1;
rep(j, 1, n) prod *= (a[j] & X).count();
ans += X.count() & 1 ? -prod : prod;
} else ans += (a[k - n] & X).count();
}
}
cout << (ans & 3 ? "NO\n" : "YES\n");
}
int main() {
cin.tie(0)->sync_with_stdio(0);
for(cin >> T; T--; solve());
}