积和式模 2 的次幂

题目链接

给定一个两边各有 $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})$ 代码：

code

#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})$ 代码：

code

#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());
}