集合幂级数、FMT、FWT学习笔记

发表于 | 更新于 | 分类于 | 博文状态：更新中 | 阅读次数：

受这几天看到的不少高维前缀和题目的影响，我决定系统地学习一下集合幂级数的一套理论了。内容主要来自2015年吕凯风（VFleaKing）国家集训队论文《集合幂级数的性质与应用及其快速算法》（pdf版本会放在附录里），包括集合并卷积、集合对称差卷积、子集卷积、快速莫比乌斯变换、快速莫比乌斯反演、快速沃尔什变换及逆变换等，以及附带进行的一些练习。

(UPD 2021.4.20) 推荐观看：Fourier Analysis of Boolean functions || @ CMU || Lecture 8a of CS Theory Toolkit

[TODO] 把练习第一题搞清楚
[TODO] 把第三道论文题学会

引言及定义

类似数列的生成函数，对于集合，我们引入集合幂级数来解决一些有关集合的动态规划问题。

定义：设 $F$ 是一个域，则称函数 $f:2^U\to F$ 是 F 上的一个形式幂级数，记 $f_S=f(S)$ 为集合幂级数第 $S$ 项的系数。

我们记为：
$$
f=\sum_{S\subseteq 2^U} f_S x^S
$$
显然可以定义集合幂级数的加法：$(f+g)(S) = f(S)+g(S)$ ，减法也类似。这两种运算都是 $O(2^n)$ 的。

如何定义乘法？为了保证乘法对加法的分配律，若 $h=f\cdot g$ ，我们应该有：
$$
\sum_{S\in 2^U} h_S x^S = \big(\sum_{L\in2^U}f_Lx^L\big) \cdot (\sum_{R\in 2^U} g_R x^R) = \sum_{L\in 2^U} \sum_{R\in 2^U} (f_L x^L) \cdot (g_R x^R)
$$
因此我们只需要规定 $(f_L x^L) \cdot (g_R x^R)$ 的运算结果，我们希望它是以某种集合运算乘起来的。我们设一个 $2^U$ 中满足交换律、结合律、空集是单位元的二元运算 $\ast$ ，那么我们就可以定义 $(f_L x^L) \cdot (g_R x^R) = (f_L g_R) x^{L\ast R}$ 。至此集合幂级数形成了一个交换环，并且包含了整个 $F$ 作为子环。

常见定义 $\ast$ 为集合并运算、集合对称差运算以及子集卷积，我们一一解决。

集合并卷积（FWT_or）

原理

即给出两个集合幂级数 $f,g$ ，求 $h = f\ast g$ 满足：
$$
h_S = \sum_{L\subseteq 2^U} \sum_{R\subseteq 2^U} [L\cup R = S] f_L g_R \tag{1}
$$
怎么做呢？暴力做是 $O(4^n)$ 的，一种做法是分治乘法，可是总感觉在理论价值上面会逊色一些。有没有类似 FFT 的一种变换，使得变换后将卷积直接变成点乘呢？莫比乌斯变换做到了这一点。

我们定义集合幂级数 $f$ 的莫比乌斯变换 $\hat{f} $为：
$$
\hat{f_S} = \sum_{T\subseteq S} f_T
$$
反过来，如何求莫比乌斯逆变换呢（常常称为莫比乌斯反演）？有关反演看这里：《炫酷反演魔术》。很容易由容斥原理，我们又有：（这一步建议在纸上推导一下，要不然可能会觉得不直观）
$$
f_S = \sum_{T\subseteq S} (-1)^{|S|-|T|} \hat{f_T}
$$
现在，我们可以解决 (1) 中的问题了。对 (1) 式的左右两边同时做莫比乌斯变换
$$
\begin{aligned}
\hat{h_S} & = \sum_{L\subseteq 2^U} \sum_{R\subseteq 2^U} [L\cup R \subseteq S] f_L g_R\
& = \sum_{L\subseteq 2^U} \sum_{R\subseteq 2^U} [L\subseteq S][ R \subseteq S] f_L g_R\
& = (\sum_{L\subseteq S}f_L) (\sum_{R\subseteq S}g_R)\
& = \hat{f_S}\hat{g_S}
\end{aligned}
$$
因此，我们想要的性质是成立的。只需要把 $f,g$ 做莫比乌斯变换，点乘起来，然后再做莫比乌斯反演即可得到集合并卷积。

如何进行莫比乌斯变换呢？可以使用递推，设 $\hat f^{(i)}S$ 表示只考虑 $S\oplus T\subseteq {1,\cdots,i}$ 的子集 $T$ 时的莫比乌斯变换第 $S$ 项，令 $\hat f_S^{(0)} = f_S$ ，那么对于每个不包含 $i$ 的 $S$ 有：
$$
\hat{f_S^{(i)}} = \hat{f_S}^{(i-1)}\
\hat{f}{S\cup{i}} ^{(i)} = \hat f_{S\cup{i}}^{(i-1)} + \hat f _S ^{(i-1)}
$$
由此递推计算即可。复杂度 $O(n2^n)$

实现

代码非常简短，dmt 变量为 1 时表示正变换，为 -1 时表示逆变换。

1
2
3
4
5
6
vector<int> fmt_or(vector<int> A, int dmt = 1) {
    for(int i = 0; i < n; ++i)
        for(int S = 0; S < m; ++S) if(~S>>i&1)
            A[S|(1<<i)] = mo(A[S|(1<<i)] + dmt*A[S]);
    return A;
}

集合交卷积（FWT_and）

由于 $S\cap T = C_U(\bar S \cup \bar T)$，因此将S、T翻转之后，可以直接用FWT_or导出。可是实际上还有更简洁的实现，如下面代码所示。

复杂度 $O(n2^n)$

1
2
3
4
5
6
vector<int> fmt_and(vector<int> A, int dmt = 1) {
    for(int i = 0; i < n; ++i)
        for(int S = 0; S < m; ++S) if(~S>>i&1)
            A[S] = mo(A[S] + dmt*A[S|(1<<i)]);
    return A;
}

集合对称差卷积（FWT_xor）

原理

这一回把 $\ast$ 定义为集合对称差。即给出两个集合幂级数 $f,g$ 求 $h = f\ast g$ ：
$$
h_S = \sum_{L\subseteq 2^U} \sum_{R\subseteq 2^U} [L\oplus R=S] f_L g_R \tag 2
$$
分治乘法在这里也是可行的，但是对应的变换方法是什么呢？这次我们使用快速沃尔什变换及其逆变换（本质上是高维 FFT）。

快速沃尔什变换：
$$
\hat{f_S} = \sum_{T\subseteq 2^U} f_T (-1) ^{|S\cap T|}
$$
其逆变换为：（我没有找到靠谱证明）
$$
f_S = \frac {1}{2^n}\sum_{T\subseteq 2^U} \hat{f_T}(-1)^{|S\cap T|}
$$
沃尔什变换是具体如何解决（2）式的问题的呢？我们基于下面的事实：
$$
\frac {1}{2^n}\sum_{T\subseteq 2^U} (-1)^{|S\cap T|} = [S=\oslash]
$$
化简过程如下：
$$
\begin{aligned}
h_S & = \sum_{L\subseteq 2^U} \sum_{R\subseteq 2^U} [L\oplus R\oplus S = \oslash] f_L g_R\
& = \sum_{L\subseteq 2^U} \sum_{R\subseteq 2^U} \frac{1}{2^n} \sum_{T\subseteq 2^n} (-1)^{|S\cap (L\oplus R \oplus S)|} f_l g_R\
& = \sum_{L\subseteq 2^U} \sum_{R\subseteq 2^U} \frac{1}{2^n} \sum_{T\subseteq 2^n} (-1)^{|T\cap L|} (-1)^{|T\cap R|}(-1)^{|T\cap S|} f_l g_R\
& = \frac{1}{2^n} \sum_{T\subseteq 2^n} (-1)^{|T\cap S|} \Big(\sum_{L\subseteq 2^U} (-1)^{|T\cap L|} f_L\Big)
\Big(\sum_{R\subseteq 2^U} (-1)^{|T\cap R|} g_R\Big)\
& = \frac{1}{2^n} \sum_{L\subseteq 2^U} (-1)^{|L\cap S|} \hat{f_S} \hat{g_S}
\end{aligned}
$$
即： $\hat{h_S} = \hat{f_S} \cdot \hat{g_S}$

注：在理论上，这种变换还有一个小问题，就是特征为2的 $F$ 上上述算法会失效，因为上述算法涉及乘以 $\frac {1}{2^n}$ 。不过一般情况下不会成问题。

接下来考虑如何计算快速沃尔什变换。依然可以使用递推，设 $\hat f^{(i)}S$ 表示只考虑 $S\oplus T\subseteq {1,\cdots,i}$ 的子集 $T$ 时的沃尔什变换第 $S$ 项，令 $\hat f_S^{(0)} = f_S$ ，那么对于每个不包含 $i$ 的 $S$ 有：
$$
\hat{f_S^{(i)}} = \hat{f_S}^{(i-1)}+\hat{f}{S\cup{i}}^{(i-1)}\
\hat{f}_{S\cup{i}} ^{(i)} = \hat f S ^{(i-1)} - \hat f{S\cup{i}}^{(i-1)}
$$
由此递推计算即可。复杂度 $O(n2^n)$

实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
vector<int> fwt_xor(vector<int> A, int dmt = 1) {
    int l, r;
    for(int i = 0; i < n; ++i)
        for(int S = 0; S < m; ++S) if(~S>>i&1) {
            l = A[S], r = A[S|(1<<i)];
            A[S] = mo(l + r);
            A[S|(1<<i)] = mo(l - r);
        }
    if(dmt==-1) {
        int inv2 = qpow(qpow(2, MOD-2), n);
        for(int S = 0; S < m; ++S)
            A[S] = muln(A[S], inv2);
    }
    return A;
}

子集卷积

这一部分感觉论文讲得有点自闭，推荐 Dance Of Faith 的这篇博客 。

原理

考虑这种形式的卷积如何处理：
$$
h_S = \sum_{T\subseteq S} f_T g_{S-T} \tag 3
$$
实际上，这种卷积等价于定义 $\ast$ 为不相交集合的并，即
$$
L \ast R =
\begin{cases}
\varnothing, & L \cap R \neq \varnothing \
L \cup R, & \text{otherwise}
\end{cases}
$$
我们稍微改写一下 （3） 式：
$$
\begin{aligned}
h_S & = \sum_{L\subseteq 2^U} \sum_{R\subseteq 2^U} [L\cap R = \varnothing] [L\cup R = S] f_L g_R \
& =\sum_{L\subseteq 2^U} \sum_{R\subseteq 2^U} [|L| + |R| = |S|][L\cup R = S] f_L g_R
\end{aligned}
$$
我们只需要处理掉 $|L|+|R| = |S|$ 这一项。实际上我们直接加一维表示集合的大小，暴力处理即可。

具体来讲，初始时，我们先只把 $f_{pc(S),S}$ （pc(S) 表示 popcount ，即 S 二进制表示中有多少位是 1 ）的值赋成原来的 $f_{S}$（$g$ 也使用同样的操作），然后对每一个 $f_i$ 做一遍FMT，点值相乘时这么卷积：$h_{i, S} = \sum\limits_{j = 0}^{i} f_{j,S} * g_{i - j, S}$ 。代码很简洁，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
vector<int> subset_conv(vector<int> A, vector<int> B) {
    vector<int> H(Len);
    vector< vector<int> > siga(Base+1, vector<int>(Len,0)), sigb = siga, sigh = siga;
    for(int S = 0; S < Len; ++S)
        siga[pc[S]][S] = A[S], sigb[pc[S]][S] = B[S];
    for(int i = 0; i <= Base; ++i) {
        fmt_or(siga[i], 1);
        fmt_or(sigb[i], 1);
        for(int S = 0; S < Len; ++S)
            for(int j = 0; j <= i; ++j)
                sigh[i][S] = mo(sigh[i][S] + muln(siga[j][S], sigb[i-j][S]));
        fmt_or(sigh[i], -1);
    }
    for(int S = 0; S < Len; ++S)
        H[S] = sigh[pc[S]][S];
    return H;
}

小结

FWT_OR

1
2
if FMT : f[S|(1<<i)] += f[S]
if IFMT: f[S|(1<<i)] -= f[S]

FWT_AND

1
2
if FMT : f[S] += f[S|(1<<i)]
if IFMT: f[S] -= f[S|(1<<i)]

FWT_XOR：

1
2
3
4
5
l = f[S], r = f[S|(1<<i)]
f[S] = l + r
f[S|(1<<i)] = l - r
if IFMT: 
    f[S] /= qpow(2, n)

子集卷积：

1
2
3
4
5
6
7
8
9
for i in range(0, n + 1):
    F[i] = FWT_OR(f[i])
    G[i] = FWT_OR(g[i])
    for S in range(0, Len):
        for j in range(0, i + 1):
            H[i][S] += F[j][S] * G[i-j][S]
    H[i] = IFWT_OR(H[i])
for S in range(0, Len):
    h[S] = H[bc[S]][S]

完整板子代码如下，由于没有采用类似 FFT 蝴蝶操作的优化，效率不是很高，在LOJ的子集卷积模板题上跑了 4500/5000 ms。Luogu 4717

点击显示/隐藏代码
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
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
#include<bits/stdc++.h>
#define rep(i, n) for(int i = 0, i##_end_ = (n); i < i##_end_; ++i)
using namespace std;
typedef pair<int, int> pii;
typedef long long ll;

// Template starts here
const int MOD = 1000000009;


inline int muln(int x, int y) { return 1LL * x * y % MOD; }
inline int qpow(int x, int y) {
    int ret = 1;
    for(; y; y >>= 1, x = muln(x, x))
        if(y & 1) ret = muln(ret, x);
    return ret;
}
inline int mo(int x) {
    if(x >= MOD) x -= MOD;
    if(x < 0) x += MOD;
    return x;
}

int Len, Base; // Ensure these values are calculated!!!!
vector<int> fmt_or(vector<int> A, int dmt = 1) {
    for(int i = 0; i < Base; ++i)
        for(int S = 0; S < Len; ++S) if(~S>>i&1)
            A[S|(1<<i)] = mo(A[S|(1<<i)] + dmt*A[S]);
    return A;
}
vector<int> fmt_and(vector<int> A, int dmt = 1) {
    for(int i = 0; i < Base; ++i)
        for(int S = 0; S < Len; ++S) if(~S>>i&1)
            A[S] = mo(A[S] + dmt*A[S|(1<<i)]);
    return A;
}
vector<int> fwt_xor(vector<int> A, int dwt = 1) {
    int l, r;
    for(int i = 0; i < Base; ++i)
        for(int S = 0; S < Len; ++S) if(~S>>i&1) {
            l = A[S], r = A[S|(1<<i)];
            A[S] = mo(l + r);
            A[S|(1<<i)] = mo(l - r);
        }
    if(dwt==-1) {
        int inv2 = qpow(Len, MOD-2);
        for(int S = 0; S < Len; ++S)
            A[S] = muln(A[S], inv2);
    }
    return A;
}

vector<int> conv(const vector<int>& A, const vector<int>& B, vector<int>(*fn)(vector<int>,int)) {
    vector<int> aa = fn(A, 1), bb = fn(B, 1);
    for(int i = 0; i < Len; ++i)
        aa[i] = muln(aa[i], bb[i]);
    aa = fn(aa, -1);
    return aa;
}

vector<int> pc;
void get_popcount(int sz) {
    pc.resize(sz, 0);
    for(int i = 1; i < sz; ++i)
        pc[i] = pc[i >> 1] + (i & 1);
}
vector<int> subset_conv(vector<int> A, vector<int> B) {
    get_popcount(Len);
    vector<int> H(Len);
    vector< vector<int> > siga(Base+1, vector<int>(Len,0)), sigb = siga, sigh = siga;
    for(int S = 0; S < Len; ++S)
        siga[pc[S]][S] = A[S], sigb[pc[S]][S] = B[S];
    for(int i = 0; i <= Base; ++i) {
        siga[i] = fmt_or(siga[i], 1);
        sigb[i] = fmt_or(sigb[i], 1);
        for(int S = 0; S < Len; ++S)
            for(int j = 0; j <= i; ++j)
                sigh[i][S] = mo(sigh[i][S] + muln(siga[j][S], sigb[i-j][S]));
        sigh[i] = fmt_or(sigh[i], -1);
    }
    for(int S = 0; S < Len; ++S)
        H[S] = sigh[pc[S]][S];
    return H;
}
// Template ends here

void show(vector<int> v) {
    int first = true;
    for(auto g : v) {
        if(first) first = false;
        else putchar(' ');
        printf("%d", g);
    }
    putchar('\n');
}

int n, m;
vector<int> a, b;

int main() {
    scanf("%d", &n);
    m = Len = 1 << n, Base = n;
    a.resize(m), b.resize(m);
    rep(i, m) scanf("%d", &a[i]);
    rep(i, m) scanf("%d", &b[i]);
    
    show(conv(a, b, fmt_and));
    show(conv(a, b, fmt_or));
    show(conv(a, b, fwt_xor));
    show(subset_conv(a, b));



    return 0;
}

练习

[HAOI2015] 按位或

Description

刚开始你有一个数字0，每一秒钟你会随机选择一个 $[0,2^n-1]$ 的数字，与你手上的数字进行按位或操作。选择数字 $i$ 的概率是 $p_i$。保证 $0\leq p_i \leq1$，$\sum p_i = 1$。问期望多少秒后，你手上的数字变成 $2^n-1$。无穷输出 INF

数据范围：$n\leq 20$

Source: BZOJ4036 vfk论文第一道例题。

Solution

留坑。

Day8 I. 岸边露伴的人生经验

Description

岸边露伴是一个天才漫画家，他经常用自己的替身天堂之门来查看别人的人生经历，为自己的漫画积累素材。最近他学会了将一个人的人生经历编码成一个 $10$ 维的向量，每一维取值为 ${0,1,2}$ 中的一个元素。定义向量$\overrightarrow{V}=(x_1,x_2,\cdots,x_{10})$ 的模长 $|\overrightarrow{V}|$ 为 $\sqrt{x_1^2+x_2^2+\cdots+x_{10}^2}$。令第 $i$ 个人的人生经历对应的向量为$\overrightarrow{V_i}$，则第ii个人和第 $j$ 个人的人生轨迹的差别可以用 $|\overrightarrow{V_i}-\overrightarrow{V_j}|$ 衡量。岸边露伴收集了 $n$ 个人的向量，他想要知道这些人里，人生轨迹差别相同的二元组有多少对，即有多少个四元组 $(i,j,k,l)(1 \le i,j,k,l \le n)$ 满足 $|\overrightarrow{V_i}-\overrightarrow{V_j}|=|\overrightarrow{V_k}-\overrightarrow{V_l}|$

数据范围：$n\leq 10^5$

Source: CCPC-Wannafly Winter Camp Day8 (Div1, onsite) Day 8 Problem I

Solution

考虑每个向量取值只有 ${0,1}$ 的情况，每个人的向量为一个二进制数，装进桶里，直接进行自己卷积自己的异或 FWT ，这时集合幂级数的每一项都代表这个集合被异或出来的 可能方式。再扫一遍每一个状态，按照 bitcount 统计距离，加入另一个数组，平方一下即可。

考虑这道题，多了一个 $2$ ，我们只需要把每一个向量中一个数位拆成两个字符即可，将 1 变为 01 ，2 变为 10 ， 0 变为 00 ，然后把上一种做法求 bitcount 统计距离的方式稍微改一下就可以了。

Code

点击显示/隐藏代码
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
71
72
73
74
75
76
77
78
79
80
81
82
83
/* Generated by powerful Codeforces Tool
 * Author: YangDavid
 * Time: 2019-04-28 19:15:38
*/

#include<bits/stdc++.h>
#define rep(i, n) for(int i = 1, i##_end_ = (n); i <= i##_end_; ++i)
using namespace std;
typedef pair<int, int> pii;
typedef long long ll;

const int maxn = 202000, maxs = 1 << 20, BB = 20, MOD = 998244353;
int n, dif[10200];
vector<int> a(maxs, 0), val(maxs, 0);

int muln(int x, int y) { return 1LL * x * y % MOD; }
int mo(int x) {
	if(x >= MOD) x -= MOD;
	if(x < 0) x += MOD;
	return x;
}
int qpow(int x, int y) {
	int ret = 1;
	for(; y; y >>= 1, x = muln(x, x))
		if(y & 1) ret = muln(ret, x);
	return ret;
}


void fwt(int dwt) {
	int l, r;
	for(int i = 0; i < BB; ++i) {
		int mask = (maxs-1) ^ (1<<i);
		for(int S = mask; S >= 0; S = (S-1)&mask) {
			l = a[S], r = a[S | (1<<i)];
			a[S] = mo(l + r);
			a[S|(1<<i)] = mo(l - r);
			if(S == 0) break;
		}
	}
	if(dwt == -1) {
		int fac = qpow(qpow(2, MOD - 2), BB);
		for(int i = 0; i < maxs; ++i)
			a[i] = muln(a[i], fac);
	}
}

int main() {
    scanf("%d", &n);
    for(int i = 0; i < n; ++i) {
    	int msk = 0, x;
    	rep(j, 10) {
    		scanf("%d", &x);
    		msk <<= 2;
    		if(x == 1) msk += 1;
    		else if(x == 2) msk += 2;
    	}
		a[msk]++;
    }
    fwt(1);
    for(int i = 0; i < maxs; ++i)
    	a[i] = muln(a[i], a[i]);
    fwt(-1);

    val[0] = 0, val[1] = 1, val[2] = 4, val[3] = 1;
    for(int S = 4; S < maxs; S += 4) {
    	int t = val[S >> 2];
    	val[S] = t;
    	val[S+1] = t + 1;
    	val[S+2] = t + 4;
    	val[S+3] = t + 1;
    }
    for(int S = 0; S < maxs; ++S) {
    	dif[val[S]] += a[S];
    }
    ll res = 0;
    for(int i = 0; i <= 100; ++i) {
    	res = mo(res + 1LL * dif[i] * dif[i] % MOD);
    }
    printf("%lld\n", res);
    
    return 0;
}

CF1119 H. 难题，精妙的 FWT

Description

给定 $a,b,c$ ，以及 $n$ 个 $k$ 位 bitmask $A_i,B_i,C_i$ ，求下面 $n$ 个集合幂级数的异或 FWT，即：
$$
\prod_{1\leq i\leq n} (ax^{A_i}+bx^{B_i}+cx^{C_i})
$$
数据范围：$n\leq 10^5,k\leq 17,0\leq A_i,B_i,C_i < 2^k,$

Source: Codeforces Global Round 2 Problem H

Solution

按照正常的 FWT 思路，就是把这 $n$ 个集合幂级数分别求出 FWT，然后再点乘起来，然后再做 IFWT。复杂度是 $O(nk2^k)$ 的，无法通过此题。

但是这道题需要注意的一点是每一个集合幂级数都只有三项，并且 $a,b,c$ 是固定的。回顾 FWT 的公式：
$$
\hat{f}(S) = \sum_{T\subseteq 2^U} f(T) (-1)^{|S\cap T|}
$$
因此每个集合幂级数 FWT 之后也仅仅会有 $\pm a \pm b \pm c$ 这么 8 种项。为了使问题更加简单，我们进行这样的变换： $B_i:=B_i \oplus A_i,C_i := C_i \oplus A_i, A_i := 0$ 。变换之后我们得出的结果的第 $S$ 项就等价于答案的第 $S\oplus xorsum$ 项，其中 $xorsum = \oplus_{i=1}^n A_i$。可以发现，这样做之后只剩下 $a\pm b \pm c$ 这四项了。

我们的一个重要观察是，将 FWT 之后的这 $n$ 个集合幂级数乘起来的结果中，对于一个固定的位置 $S$ ，这一项的系数一定是 $(a+b+c)^x (a+b-c)^y (a-b+c)^z (a-b-c)^w$ 的形式。如果我们能够把 $x,y,z,w$ 解出来，整个问题就解决了。首先，我们有一个最朴素的关系式：
$$
x+y+z+w = n
$$
之后有两种说法，一种是官方题解所说的考虑所有 $n$ 个 FWT 之后的幂级数的第 $S$ 项之和，然后再多考虑一个条件即可得到四个方程；另一种则非常简洁、对称、优美，可是我对这种方法还没有简洁优美的证明。这种方法是这样的：

  • 将 $A_i\oplus B_i$ 的值统计入一个数组 $f$，对 $f$ FWT之后，对于其第 $S$ 项有：$x+y-z-w=f[S]$。
  • 将 $A_i\oplus C_i$ 的值统计入一个数组 $g$，对 $g$ FWT之后，对于其第 $S$ 项有：$x-y+z-w=g[S]$。
  • 将 $B_i\oplus C_i$ 的值统计入一个数组 $h$，对 $h$ FWT之后，对于其第 $S$ 项有：$x-y-z+w=h[S]$。

然后消元法即可。

Code

点击显示/隐藏代码
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
/* Generated by powerful Codeforces Tool
 * Author: YangDavid
 * Time: 2019-05-02 16:57:53
*/

#include<bits/stdc++.h>
#define rep(i, n) for(int i = 1, i##_end_ = (n); i <= i##_end_; ++i)
using namespace std;
typedef pair<int, int> pii;
typedef long long ll;

const int maxn = 102000, MS = (1 << 17) + 2, MOD = 998244353;
int A[maxn], B[maxn], C[maxn], n, k, full;
int AB[MS], AC[MS], BC[MS], ans[MS], out[MS];

int mo(int x) { if(x >= MOD) x -= MOD; if(x < 0) x += MOD; return x; }
int muln(int x, int y) { return 1LL * x * y % MOD; }
int qpow(int x, int y) {
    int ret = 1;
    for(; y; y >>= 1, x = muln(x, x))
        if(y & 1) ret = muln(ret, x);
    return ret;
}

const int inv2 = qpow(2, MOD - 2), inv4 = qpow(4, MOD - 2);

void fwt(int* arr, int dwt = 1) {
    int l, r;
    for(int i = 0; i < k; ++i)
        for(int S = 0; S < full; ++S) if(~S>>i&1) {
            l = arr[S], r = arr[S|(1<<i)];
            arr[S] = mo(l + r);
            arr[S|(1<<i)] = mo(l - r);
        }

    int iv = qpow(full, MOD - 2);
    if(dwt == -1) for(int S = 0; S < full; ++S)
        arr[S] = muln(arr[S], iv);
}

int gg[4], a, b, c, xorsum;
int gett(ll s, ll t, ll u, ll v) {
    int x = muln(mo(mo(s+t)+mo(u+v)), inv4);
    int y = mo(muln(mo(s+t), inv2) - x);
    int z = mo(muln(mo(s+u), inv2) - x);
    int w = mo(muln(mo(s+v), inv2) - x);
    return muln( muln(qpow(gg[0],x),qpow(gg[1],y)) , muln(qpow(gg[2],z),qpow(gg[3],w)) );
}

int main() {
    scanf("%d%d%d%d%d", &n, &k, &a, &b, &c); full = (1 << k);
    gg[0] = mo(mo(a+b)+c), gg[1] = mo(mo(a+b)-c), gg[2] = mo(mo(a-b)+c), gg[3] = mo(mo(a-b)-c);
    for(int i = 0; i < n; ++i) {
        scanf("%d%d%d", &A[i], &B[i], &C[i]);
        AB[A[i]^B[i]]++, AC[A[i]^C[i]]++, BC[B[i]^C[i]]++;
        xorsum ^= A[i];
    }
    fwt(AB), fwt(AC), fwt(BC);
    for(int S = 0; S < full; ++S) {
        ans[S] = gett(n, AB[S], AC[S], BC[S]);
    }
    fwt(ans, -1);
    for(int S = 0; S < full; ++S)
        out[S] = ans[S ^ xorsum];
    for(int S = 0; S < full; ++S)
        printf("%d%c", out[S], " \n"[S==full-1]);

    return 0;
}

CF662C 超级好题

Description

你有一个 $n$ 行 $m$ 列的 01 矩阵，你可以进行翻转行与翻转列两种操作（翻转即 0 变 1，1 变 0）任意多次，问最后得到的矩阵的 1 的个数最少是多少。

数据范围：$n\leq 20,m\leq 100000$

Source: CROC 2016 - Final Round Problem C. Binary Table

Solution

注意 $n\leq 20$ 的条件。先考虑 $O(m2^n )$ 暴力怎么做。记第 $j$ 列的数的 bitmask 为 $A_j$ ，我们枚举行的翻转状态的 bitmask为 $S$ ，然后对于每一种状态 $S$，我们查看每一列，第 $j$ 列被变为了 $A_j \oplus S$。考虑这一列是否翻转，显见答案增加了 $\min{pop(A_j \oplus S), n-pos(A_J\oplus S)}$，其中 $pop(S)$ 是指 $S$ 二进制表示中 1 的个数。我们记$f(S) = \min{pop(S),n-pop(S)}$ ，那么对于行的枚举状态 $S$ 我们得到的答案就是这个：
$$
\sum_{i=1}^m f(S\oplus A_i)
$$
由于 $A_j$ 的值域也是 $2^n$，我们记 $A$ 中 $x$ 这个值出现的次数为 $c(x)$，那么上式可以改写为：
$$
\sum_{T \subseteq 2^n} f(T \oplus S) \cdot c(T)
$$
是不是有点 FWT 的形式了？进一步地：
$$
\begin{aligned}
ans(S) &= \sum_{T \oplus U = S} f(T) c(U)
\end{aligned}
$$
这不正是 FWT 吗？于是我们求出 $f$ 与 $c$ 的异或 FWT，然后扫一边 ans 数组找最小值即可。

Code

点击显示/隐藏代码
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
#include<bits/stdc++.h>
using namespace std;
typedef pair<int, int> pii;
typedef long long ll;

const int maxs = (1 << 20) + 2333, MOD = 998244353;
int c[maxs], g[maxs], n, m, full , pop[maxs];
char s[22][101010];

inline int mo(int x) { if(x >= MOD) return x-MOD; if(x < 0) return x+MOD; return x; }
inline int muln(int x, int y) { return (ll)x * y % MOD; }
inline int qpow(int x, int y) {
    int ret = 1;
    for(; y; y >>= 1, x = muln(x, x))
        if(y & 1) ret = muln(ret, x);
    return ret;
}

int Base, Len;
int* fwt(int* A, int dwt = 1) {
    int l, r;
    for(int i = 0; i < Base; ++i)
        for(int S = 0; S < Len; ++S) if(~S>>i&1) {
            l = A[S], r = A[S|(1<<i)];
            A[S] = mo(l + r);
            A[S|(1<<i)] = mo(l - r);
        }
    if(dwt==-1) {
        int inv2 = qpow(Len, MOD-2);
        for(int S = 0; S < Len; ++S)
            A[S] = muln(A[S], inv2);
    }
    return A;
}

int main() {
    scanf("%d%d", &n, &m);
    full = (1 << n) - 1;
    for(int i = 0; i < n; ++i)
        scanf("%s", s[i]);
    for(int i = 0; i < m; ++i) {
        int res = 0;
        for(int j = 0; j < n; ++j)
            res = 2 * res + s[j][i] - '0';
        c[res]++;
    }
    for(int S = 1; S <= full; ++S) {
        pop[S] = pop[S>>1] + (S&1);
        g[S] = min(pop[S], n - pop[S]);
    }
    Len = 1<<n, Base = n;
    fwt(c, 1);
    fwt(g, 1);
    for(int S = 0; S <= full; ++S)
        c[S] = muln(c[S], g[S]);
    fwt(c, -1);

    int ans = 0x3f3f3f3f;
    for(int S = 0; S <= full; ++S)
        ans = min(ans, c[S]);
    printf("%d\n", ans);


    return 0;
}

[HDU5909]

[WC2018]州区划分