题目链接
给定一个两边各有 \(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());
}
|