集合幂级数相关

集合幂级数相关算法:快速莫比乌斯变换(FMT),高维前缀和(SOSDP),快速沃尔什变换(FWT),子集卷积(Subset Convolution)

CHANGE LOG

NOI 大纲里没有把位运算卷积如 FMT,FWT,子集卷积等知识点单独列出,但高维前缀和(SOSDP)是应用比较广泛的重要算法。

学习上述算法,首先要理解什么是集合幂级数。

1. 集合幂级数

1.1 定义

集合幂级数最初由吕凯风在他的 2015 年集训队论文《集合幂级数的性质与应用及其快速算法》中提出。这个名字听起来非常高级,但实际上不难理解。如果读者能够有一点点形式幂级数的基础知识会更容易理解。

为方便说明,我们引入一些记号:

  • 定义 \(2 ^ X\) 表示 \(X\) 的所有幂集,即 \(X\) 的所有子集组成的集合。
  • 定义 \(U\) 为全集 \(\{1, 2, \cdots, n\}\),其中 \(n = |U|\),即 \(U\) 的大小。用数字来表示元素是因为我们并不关心集合中的元素具体是什么,数字只是一个易于理解的符号。
  • 若无特殊说明,\(S\) 均为 \(U\) 的子集。

形式化地,域 \(F\) 上的集合幂级数 \(f\)\(2 ^ U \to F\) 的函数,其中对于每个 \(S \subseteq U\),函数 \(f\) 均有对应的取值 \(f_S \in F\)

换句话说,集合幂级数就是从全集的所有子集(全集的幂集)到某个域 \(F\)(OI 中一般研究模某个大质数 \(p\) 意义下的整数域 \(\mathbb Z_p\) 或实数域 \(\mathbb R\))上的 映射。把 \(U\) 的每个子集 \(S\) 带入函数,均有对应的取值 \(f_S\)

通俗地(不严谨),集合幂级数就是把函数的自变量值类型从一个数 \(x\) 换成了一个集合 \(S\),且函数的定义域为全集的所有子集。

  • 对比:形式幂级数由序列导出,而集合幂级数由集合导出。

1.2 表示

有了定义,我们考虑如何表示一个集合幂级数。由于它的定义实在有些抽象,我们需要一个方法清晰地描述集合幂级数 \(f\) 在每个集合 \(S\) 处的取值(总不能说 \(f\)\(\varnothing\) 映射到 \(a\),把 \(\{1\}\) 映射到 \(b\),把 \(\{\cdots\}\) 映射到 \(\cdots\) 吧?这样太麻烦了,而且不方便形式化表述)。

类比形式幂级数,我们可以给每个集合 \(S\) 对应的取值 \(f_S\) 后添加 占位符 \(x ^ S\)。即 \(cx ^ S\) 表示一个将 \(S\) 映射到 \(c\ (f_S = c)\) 的集合幂级数,且对于不等于 \(S\)\(U\) 的子集 \(T\)\(f_T\) 的取值为 \(0\)(因为 \(0x ^ T\) 可以忽略)。因此,\(f = \sum\limits_{S \in 2 ^ U} f_S x ^ S\)

例如,\(f = 3x ^ {\varnothing} + 4x ^ {\{1\}} + 9 ^ {\{1, 2\}}\) 对应一个 \(f_{\varnothing} = 3\)\(f_{\{1\}} = 4\)\(f_{\{2\}} = 0\)\(f_{\{1, 2\}} = 9\) 的集合幂级数 \(f\)

  • 注意,\(cx ^ S\) 是一个整体而非 \(c\times x ^ S\),因为 \(x ^ S\) 本身没有意义,它是占位符。只有它和它前面的系数作为整体时有意义。\(x ^ S\) 表示集合为 \(S\),而 \(c\) 表示 \(f_S\) 对应的值为 \(c\)
  • 下文有些地方会将 \(S\) 看成一个 \(0 \sim 2 ^ n – 1\) 的二进制数 \(i\)。这表示 \(i = \sum\limits_{j \in S} 2 ^ {j – 1}\),或者说 \(S\)\(i\) 在二进制下值为 \(1\) 的位。读者需要在集合 \(S\) 和二进制数 \(i\) 之间建立起快速的对应关系。

1.3 加法与乘法

形式幂级数有加法(对应系数相加)和乘法(加法卷积,相关算法有著名的 FFT 和 NTT),集合幂级数自然也需要定义加法和乘法。

加法比较容易,只要将对应系数相加即可。若集合幂级数 \(h = f + g\),则 \(h_S = f_S + g_S\)。注意,\(f, g\) 对应的全集需相同。

乘法的定义就稍微复杂了一点。考虑集合幂级数 \(h = f \cdot g = \left(\sum\limits_{L \in 2 ^ U} f_L x ^ L \right) \cdot \left(\sum\limits_{R \in 2 ^ U} g_R x ^ R \right)\)。我们希望集合幂级数的乘法对加法有分配律,这样才能推导出更多性质,扩展集合幂级数的应用。因此,\(h = \sum\limits_{L \in 2 ^ U}\sum\limits_{R \in 2 ^ U} f_L x ^ L \cdot g_R x ^ R\)

考虑如何定义 \(f_L x ^ L \cdot g_R x ^ R\)。设其结果为 \((f_L \times g_R) x ^ {L * R}\),其中 \(*\)\(2 ^ U\) 上的二元运算。

  • 为满足集合幂级数的乘法对加法的分配律,\(f_L\)\(g_R\) 应当以域 \(F\) 中的乘法相结合(存在域 \(F\) 不满足乘法对加法的分配律吗?至少笔者没见过,而且这不重要 = . =),即 \(\times\) 是域 \(F\) 中的乘法。一般即模 \(p\) 意义下的整数域 \(\mathbb Z_p\)\(2 \times 3 \equiv 1 \pmod 5\))或有理数域 \(\mathbb R\) 上的乘法(\(2.5 \times 2.5 = 6.25\))。
  • 为满足集合幂级数的乘法交换律,\(*\) 应满足交换律。
  • 为满足集合幂级数的乘法结合律,\(*\) 应满足结合律。
  • 当然,\(*\) 还应该有单位元 \(\varnothing\)。这符合集合之间二元运算关系的基本直觉。

常见的满足上述性质的集合二元运算并不唯一。因此,集合幂级数的乘法类型也并不唯一,常见的有 并卷积对称差卷积 以及 子集卷积。如果将集合的所有子集用二进制数表示,则它们分别对应了位运算卷积中的 或卷积异或卷积子集卷积。相关算法将在接下来介绍。

2. 集合并卷积

2.1 定义

令集合之间的二元运算 \(*\) 为集合并运算即可得到集合并卷积的形式

\[h_S = \sum_{L \in 2 ^ U} \sum_{R \in 2 ^ U} [L \cup R = S] f_L g_R \]

即下标分别为 \(L\)\(R\) 的两个数 \(f_L, g_R\) 相乘后会贡献到下标为 \(L \cup R\) 的值。如果从位运算的角度理解,它就是或卷积,即 \(h_i = \sum\limits_{j | k = i} f_j g_k\)

最暴力地执行卷积,时间复杂度 \(\mathcal{O}(4 ^ n)\),非常劣。

2.2 分治解法

根据 位运算在每一位的独立性,我们考虑分治求解卷积。

\(f = f ^ – + x ^ {\{n\}} \cdot f ^ +\)\(g = g ^ – + x ^ {\{n\}} \cdot g ^ +\)(注意 \(\cdot\) 表示集合并卷积,即两个集合贡献到它们的并)。即我们将第 \(n\) 位单独提出来,\(f ^ -\) 表示 \(f\) 不包含元素 \(n\) 的所有集合对应的集合幂级数,或者说 \(f_0x ^ 0 \sim f_{2 ^ {n – 1} – 1}x ^ {2 ^ {n – 1} – 1}\) 单独提取出来,即 \(f ^ – = \sum\limits_{n \notin S} f_Sx ^ S = \sum\limits_{0\leq i < 2 ^ {n – 1}} f_ix ^ i\)。而 \(f ^ +\) 表示 \(f\) 包含元素 \(n\) 的所有集合对应的集合幂级数,同时将这些集合中的元素 \(n\) 去掉,但系数保留,即 \(f ^ + = \sum\limits_{n\in S} f_Sx ^ {S\backslash\{n\}} = \sum\limits_{2 ^ {n – 1} \leq i < 2 ^ n} f_ix ^ {i – 2 ^ {n – 1}}\)

显然,\(f \cdot g = (f ^ – + x ^ {\{n\}} \cdot f ^ +) \cdot (g ^ – + x ^ {\{n\}} \cdot g ^ +)\)。由于 \(x ^ {\{n\}} \cdot x ^ {\{n\}} = x ^ {\{n\}}\),因此

\[f \cdot g = f ^ – \cdot g ^ – + x ^ {\{n\}} \cdot ((f ^ – + f ^ +) \cdot (g ^ – + g ^ +) – f ^ – \cdot g ^ -) \]

注意到我们将问题分成了两次不包含 \(n\) 的两个集合幂级数的集合并卷积。求解 \(\mathcal{T}(n) = 2\mathcal{T}(\frac n 2) + \mathcal{O}(n)\),得到 \(\mathcal{T}(n) = \mathcal{O}(n \log n)\)。因此,上述分治做法的时间复杂度为 \(\mathcal{O}(n 2 ^ n)\)

vector <int> or_conv(int n, vector <int> &f, vector <int> &g) {
	vector <int> h(1 << n);
	if(n == 0) return h[0] = 1ll * f[0] * g[0] % mod, h;
	vector <int> a(1 << n - 1), b(1 << n - 1), c(1 << n - 1), d(1 << n - 1);
	for(int i = 0; i < 1 << n - 1; i++) {
		a[i] = f[i], b[i] = g[i];
		c[i] = (f[i] + f[i + (1 << n - 1)]) % mod;
		d[i] = (g[i] + g[i + (1 << n - 1)]) % mod; 
	}
	vector <int> l = or_conv(n - 1, a, b), r = or_conv(n - 1, c, d);
	for(int i = 0; i < 1 << n - 1; i++) h[i] = l[i];
	for(int i = 0; i < 1 << n - 1; i++) h[i + (1 << n - 1)] = (r[i] + mod - l[i]) % mod;
	return h;
}

2.3 莫比乌斯变换

上述分治做法已经达到了相当优秀的复杂度,但是由于其涉及到递归和数组复制,导致常数并不理想。对分治做法加以改进,我们得到了不需要递归的快速莫比乌斯变换(FMT)和快速沃尔什变换(FWT,这将在第 4 章提到)。

观察分治做法,可以发现在后半部分,我们相当于先将所有 \(f_i (2 ^ {n – 1} \leq i < 2 ^ n)\) 加上 \(f_{i – 2 ^ {n – 1}}\),相乘后再减去 \(f_{i-2 ^ {n – 1}}\) 相乘之后的结果。前面一部分给我们的感觉是 将所有第 \(n\) 位为 \(0\) 的项前的系数加到其对应的第 \(n\) 位为 \(1\) 的项前的系数上,后面一部分则是 进行了这种变换的两个集合幂级数的对应系数相乘(因为递归到最底层是对应系数相乘)后,将所有第 \(n\) 位为 \(1\) 的项前的系数减去其对应的第 \(n\) 位为 \(0\) 的项前的系数

因此,我们考虑 模拟递归的过程从大到小 枚举每一位 \(j(1\leq j \leq n)\),对所有第 \(j\) 位为 \(1\) 的数 \(i\),执行 \(f_i\gets f_i + f_{i – 2 ^ {j – 1}}\)\(g_i\gets g_i + g_{i – 2 ^ {j – 1}}\)。这样,我们得到了两个新的集合幂级数 \(\hat f\)\(\hat g\)。令 \(\hat h\)\(\hat f, \hat g\) 对应位置系数相乘后得到的集合幂级数,从小到大 枚举每一位 \(j(1 \leq j \leq n)\),对所有第 \(j\) 位为 \(1\) 的数 \(i\),执行 \(\hat h_i \gets \hat h_i – \hat h_{i – 2 ^ {j – 1}}\),得到 \(h\)\(h\) 即为 \(f\cdot g\)

更进一步地,注意到枚举 \(j\) 的顺序实际上并没有影响,因为最终得到的形式均为

\[\begin{aligned} \hat f_S = \sum_{T\subseteq S} f_T \\ \hat f_i = \sum_{j \mid i = i} f_j \end{aligned} \]

  • 注:关于上述操作得到该形式的原因见 2.4 小节。关于该形式的正确性,见本小节最后。

对于 集合幂级数 \(f\),我们称形如 \(\hat f_i = \sum\limits_{j \mid i = i} f_j\) 的变换为 莫比乌斯变换\(\hat f\) 又记作 \(\mathrm{FMT}(f)\),因为我们通常用快速莫比乌斯变换(Fast Mobius Transform)求解莫比乌斯变换。

考虑 \(\hat h_i \gets \hat h_i – \hat h_{i – 2 ^ {j – 1}}\) 的过程,不难发现对于 \(T\subseteq S\)\(\hat h_T\)\(h_S\) 的贡献的系数为 \((-1) ^ {|S| – |T|}\),这是因为 \(T\) 每减少一个元素,贡献就会乘上 \(-1\)。根据 容斥原理 也可以推导得到该结果。对于集合幂级数 \(f\),我们称形如 \(\hat f_i = \sum\limits_{j\mid i = i} (-1) ^ {\mathrm{popcount}(i) – \mathrm{popcount}(j)}f_j\) 的变换为 莫比乌斯反演(注意与数论的莫比乌斯反演进行区分)。\(\hat f\) 又记作 \(\mathrm{FMI}(f)\)

  • FMT 和 FMI 互为逆变换。

快速莫比乌斯变换和快速莫比乌斯反演的具体算法在上文已经进行了讲解。观察式子,可以总结出一般形式:枚举每一位 \(j\),令第 \(j\) 位为 \(1\)\(i\) 执行 \(f_j\gets f_j + c \times f_{j – 2 ^ {j – 1}}\)。容易发现 FMT 即 \(c = 1\),FMI 即 \(c = -1\)。因此,我们有如下代码(P4717 【模板】快速莫比乌斯/沃尔什变换 (FMT/FWT) 部分代码)。

#include <bits/stdc++.h>
using namespace std;

const int N = 1 << 17;
const int mod = 998244353;

int n, a[N], b[N], A[N], B[N], C[N];
void FMT(int *f, int c) {
	for(int i = 1; i < 1 << n; i <<= 1)
		for(int j = 0; j < 1 << n; j++)
			if(i & j)
				f[j] = (f[j] + 1ll * c * f[j ^ i] + mod) % mod;
}
int main() {
	cin >> n;
	for(int i = 0; i < 1 << n; i++) scanf("%d", &a[i]);
	for(int i = 0; i < 1 << n; i++) scanf("%d", &b[i]);
	memcpy(A, a, N << 2), memcpy(B, b, N << 2);
	FMT(A, 1), FMT(B, 1);
	for(int i = 0; i < 1 << n; i++) C[i] = 1ll * A[i] * B[i] % mod;
	FMT(C, mod - 1);
	for(int i = 0; i < 1 << n; i++) printf("%d ", C[i]);
	return 0;
}

接下来我们证明 FMT 进行集合并卷积的正确性。

据定义,\(h_S = \sum\limits_L \sum\limits_R [L \cup R = S] f_L g_R\),因此 \(\hat h_S = \sum\limits_L \sum\limits_R [L \cup R \subseteq S] f_L g_R\)。由于 \([L \subseteq S][R \subseteq S]\) 等价于 \([L \cup R \subseteq S]\),故 \(\hat h_S = \sum\limits_{L \subseteq S} f_L \sum\limits_{R \subseteq S} g_R = \hat f_L \hat g_R\),证毕。

2.4 高维前缀和

我们将说明 FMT 操作为什么能得到一个集合幂级数 \(f\) 的莫比乌斯变换形式 \(\hat f\)。它和 高维前缀和(SOSDP,Sum over Subsets DP) 有着密不可分的关系。

回忆二维前缀和的做法,我们利用容斥原理得到递推式 \(s_{i, j} = s_{i – 1, j} + s_{i, j – 1} – s_{i – 1, j – 1}\)。但这对于更高维数的前缀和并没有扩展价值,因为对于 \(n\) 维前缀和,利用容斥原理得到的 \(s_i\) 的递推式共有 \(2 ^ n\) 项(包括 \(s_i\) 本身)。此时时间复杂度骤然上升至 \(2 ^ n \times \prod V_i\),其中 \(V_i\) 是第 \(i\) 个维度的大小。

事实上,求解高维前缀和有更一般的做法。即首先枚举每一维,对所有数 仅关于这一维 做前缀和。例如,对于三维前缀和 \(s_{n, m, t}\)。首先对于 \(n\) 这一维枚举所有 \(s_{i, j, k}\),然后令 \(s_{i + 1, j, k} \gets s_{i + 1, j, k} + s_{i, j, k}\),注意 \(i\) 必须升序枚举(显然,做一维前缀和的时候没有人会按下标从大到小枚举),其次再对 \(m\) 这一维和 \(t\) 这一维做类似的操作。容易证明 \(s_{n, m, t} = \sum\limits_{i \in [1, n]}\sum\limits_{j \in [1, m]}\sum\limits_{k \in [1, t]} s_{i, j, k}\)

对于更高的维度同理。这样做高维前缀和的时间复杂度是 \(\sum V_i\times \prod V_i\),非常优秀。

关于 FMT,令第 \(j\) 位为 \(1\)\(i\) 执行 \(f_i\gets f_i + f_{i – 2 ^ {j – 1}}\) 可以看做 \(j\) 这一位执行前缀和,因此变换结束后得到的就是 \(f\) 的高维前缀和。故有 \(\hat f_i = \sum\limits_{j\mid i = i} f_j\) 这一形式。

2.5 FMT 的性质

注意到 FMT 本质上是一个形如 \(\hat f_i = \sum\limits_{0 \leq j < 2 ^ n} c(i, j) f_j\)线性变换,其标准矩阵 \(C\)\(C_{i, j} = [j \mid i = i]\)。故 FMT 具有 线性性。对于 FMI 同理。根据 FMT 和 FMI 互逆,它们的标准矩阵也互逆。

因此 \(\mathrm{FMT}(f + g) = \mathrm{FMT}(f) + \mathrm{FMT}(g)\),且对于任意正整数 \(c\),均有 \(\mathrm{FMT}(cf) = c\mathrm{FMT}(f)\)。更一般地,由 叠加原理,对于若干集合幂级数 \(f_i\)\(\mathrm{FMT}(\sum c_if_i) = \sum c_i\mathrm{FMT}(f_i)\)。这里 \(i\) 表示第 \(i\) 个集合幂级数而非集合幂级数 \(f\)\(x ^ i\) 项系数。

FMT 还有一个重要的性质,那就是 \(\prod f_i = \mathrm{FMI}(\prod \mathrm{FMT}(f_i))\)。等式左边的求积符号表示 集合并卷积,右边的求积符号表示 对应项系数相乘。同样的,\(i\) 表示第 \(i\) 个集合幂级数而非集合幂级数 \(f\)\(x ^ i\) 项系数。

这使得我们在计算若干集合幂级数的集合并卷积时,可以先全部求莫比乌斯变换,对应位置相乘后再莫比乌斯反演回来,而不用每计算两个集合幂级数的卷积就要莫比乌斯反演一遍。显然,一个集合幂级数 FMI 之后再 FMT,所得结果不变(互逆性)。

注意,这和若干形式幂级数求加法卷积通常使用的分治 + FFT(而非分治 FFT!)不同。这是因为两个 \(n\) 次多项式相乘会得到 \(2n\) 次多项式(除非对 \(x ^ {n + 1} – 1\) 取模,这相当于循环卷积),次数快速增长,所以求 DFT 相乘之后要赶紧 IDFT 回去。但两个 \(n\) 维集合幂级数相乘,结果仍然是 \(n\) 维。

3. 集合对称差卷积

3.1 定义

定义二元集合运算 对称差 表示 \(A\oplus B = \{x \mid (x \in A) \oplus (x \in B)\}\),即只属于其中一个集合的所有元素构成的集合。它相当于逻辑运算中的 异或 运算符。

令集合之间的二元运算 \(*\) 为集合对称差运算即可得到集合对称差卷积的形式

\[h_S = \sum_{L \in 2 ^ U} \sum_{R \in 2 ^ U} [L \oplus R = S] f_L g_R \]

即下标分别为 \(L\)\(R\) 的两个数 \(f_L, g_R\) 相乘后会贡献到下标为 \(L \oplus R\) 的值。如果从位运算的角度理解,它就是异或卷积,即 \(h_i = \sum\limits_{j \oplus k = i} f_j g_k\)

最暴力地执行卷积,时间复杂度 \(\mathcal{O}(4 ^ n)\),非常劣。

3.2 分治解法

类似地,利用位运算在每一维独立的性质,考虑分治。

\(f = f ^ – + x ^ {\{n\}} \cdot f ^ +\)\(g = g ^ – + x ^ {\{n\}} \cdot g ^ +\),类似集合并卷积的推导,我们有

\[\begin{aligned} fg & = (f ^ – + x ^ {\{n\}}f ^ +)(g ^ – + x ^ {\{n\}}g ^ +) \\ & = (f ^ – g ^ – + f ^ + g ^ +) + x ^ {\{n\}}(f ^ – g ^ + + f ^ + g ^ -) \end{aligned} \]

如果直接拆成四个子卷积,时间复杂度仍然是 \(\mathcal{O}(4 ^ n)\)。考虑计算 \(p = (f ^ – + f ^ +)(g ^ – + g ^ +)\) 以及 \(q = (f ^ – – f ^ +)(g ^ – – g ^ +)\)(当然,\((f ^ + – f ^ -)(g ^ + – g ^ -)\) 也可以,可以说明两者是等价的)。那么 \(fg = \dfrac {p + q} 2 + x ^ {\{n\}} \dfrac {p – q} 2\)。时间复杂度 \(\mathcal{O}(n 2 ^ n)\)

3.3 沃尔什变换

以下内容均来自于吕凯风的论文。

我们进行一步非常神仙的转化:\([S = \varnothing] = \dfrac 1 {2 ^ n}\sum\limits_{T \subseteq 2 ^ U} (-1) ^ {|S \cap T|}\)。若 \(S \neq \varnothing\),则选出 \(S\) 中任意一个元素 \(x\),包含 \(x\)\(T\) 和不包含 \(x\)\(T\) 一一对应且贡献相反。若 \(S = \varnothing\) 显然成立。

注意到 \(L \oplus R = S\) 等价于 \(L \oplus R \oplus S = \varnothing\),且 \((L \oplus R) \cap S = (L \cap S) \oplus (R \cap S)\),所以

\[\begin{aligned} h_S & = \sum_{L \in 2 ^ U} \sum_{R \in 2 ^ U} [L \oplus R = S] f_L g_R \\ & = \dfrac 1 {2 ^ n} \sum_{L \in 2 ^ U} \sum_{R \in 2 ^ U} \sum_{T \in 2 ^ U} (-1) ^ {|(S \oplus L \oplus R) \cap T|} f_L g_R \\ & = \dfrac 1 {2 ^ n} \sum_{T \in 2 ^ U} (-1) ^ {|S \cap T|} \left(\sum_{L \in 2 ^ U} (-1) ^ {|L \cap T|} f_L \right) \left(\sum_{R \in 2 ^ U} (-1) ^ {|R \cap T|} g_R \right) \\ \end{aligned} \]

因此,考虑令 \(\hat f_S = \sum\limits_{T \in 2 ^ U} f_T(-1) ^ {|S \cap T|}\),它被称为 \(f\)沃尔什变换\(\hat f\) 也记作 \(\mathrm{FWT}(f)\),因为我们用快速沃尔什变换(Fast Walsh-Hadamard Transform)求解沃尔什变换。

同理,根据 \([S = \varnothing] = \dfrac 1 {2 ^ n}\sum\limits_{T \subseteq 2 ^ U} (-1) ^ {|S \cap T|}\),我们定义 \(\hat f\)沃尔什逆变换\(f_S = \dfrac 1 {2 ^ n} \sum\limits_{T \subseteq 2 ^ U} \hat f_T(-1) ^ {|S\cap T|}\)\(f\) 也记作 \(\mathrm{IFWT}(\hat f)\),因为我们用快速沃尔什逆变换(Inverse Fast Walsh-Hadamard Transform)求解沃尔什你变换。

\[\begin{aligned} f_S & = \dfrac 1 {2 ^ n} \sum\limits_{T \subseteq 2 ^ U} \hat f_T(-1) ^ {|S\cap T|} \\ & = \dfrac 1 {2 ^ n} \sum\limits_{T \subseteq 2 ^ U} \sum\limits_{X \subseteq 2 ^ U} (-1) ^ {|S\cap T|} f_X (-1) ^ {|T \cap X|} \\ & = \sum\limits_{X \subseteq 2 ^ U} f_X \dfrac 1 {2 ^ n} \sum\limits_{T \subseteq 2 ^ U} (-1) ^ {|(S\oplus X)\cap T|} \\ & = \sum\limits_{X \subseteq 2 ^ U} f_X [S\oplus X = \varnothing] \\ & = f_S \end{aligned} \]

因此沃尔什变换和沃尔什逆变换互逆。


考虑如何计算沃尔什(逆)变换。直接枚举子集,时间复杂度 \(\mathcal{O}(3 ^ n)\),比一开始的 \(4 ^ n\) 更优,但仍然不够快。我们尝试一位一位地构造。

设当前考虑第 \(j\) 位,则对于两个除了第 \(j\) 位以外均相同的位置 \(p, q(p + 2 ^ {j – 1} = q)\),由于 \(p\) 这一位是 \(0\),由于 \(0\)\(0 / 1\) 的与均为 \(0\),因此贡献了 \(0\)\(-1\),故 \(f_p \gets f_p + f_q\)。同理,因为 \(1 \land 0 = 0\)\(1\land 1 = 1\),可以推出 \(f_q \gets f_p – f_q\)。因此对每一位的每两对这样的 \(p, q\) 做如上变换,即可得到 \(f\) 的沃尔什变换。这相当于依次把每一位 \(j\)(每一个元素)放到 \(f_i = \sum\limits_{0\leq k < n} (-1) ^ {\mathrm{popcount}(i\ \mathrm{and}\ k)} f_k\) 的考虑范围内。

  • 方便起见,我们将 \(\mathrm{popcount}(i\ \mathrm{and}\ k)\)\(i\)\(k\) 的按位与的 \(1\) 的个数的奇偶性记作 \(i\odot k\)

没听懂?没关系,我们换一种理解方式。在计算 \(f\) 的沃尔什变换时,设 $f = f ^ – + x ^ {{n}} \cdot f ^ + $,先计算 \(f ^ -\)\(f ^ +\) 的沃尔什变换,再考虑将它们合并。对于 \(f_i\),若 \(i\) 的第 \(n\) 位为 \(0\),由于 \(0 \land 0 = 0\)\(0\land 1 = 0\),所以合并

\[\begin{aligned} \hat f_i & = \left(\sum\limits_{0\leq j < 2 ^ {n – 1}} f_j \times (-1) ^ {i\odot j} \right) + \left(\sum\limits_{2 ^ {n – 1}\leq j < 2 ^ n} f_j \times (-1) ^ {i \odot j} \right) \\ & = \left(\sum\limits_{0\leq j < 2 ^ {n – 1}} f_j \times (-1) ^ {i\odot j} \right) + \left(\sum\limits_{2 ^ {n – 1}\leq j < 2 ^ n} f_j \times (-1) ^ {i\odot (j – 2 ^ {n – 1})} \right) \\ & = \hat f ^ -_j + \hat f ^ +_j \end{aligned} \]

同理,当 \(i\) 的第 \(n\) 位为 \(1\) 时,由于 \(1 \land 0 = 0\)\(1\land 1 = 1\),故

\[\begin{aligned} \hat f_i & = \left(\sum\limits_{0\leq j < 2 ^ {n – 1}} f_j \times (-1) ^ {i\odot j} \right) + \left(\sum\limits_{2 ^ {n – 1}\leq j < 2 ^ n} f_j \times (-1) ^ {i \odot j} \right) \\ & = \left(\sum\limits_{0\leq j < 2 ^ {n – 1}} f_j \times (-1) ^ {(i – 2 ^ {j – 1})\odot j} \right) {\color{red}-} \left(\sum\limits_{2 ^ {n – 1}\leq j < 2 ^ n} f_j \times (-1) ^ {(i – 2 ^ {n – 1})\odot (j – 2 ^ {n – 1})} \right) \\ & = \hat f ^ -_j – \hat f ^ +_j \end{aligned} \]

因此,我们可以考虑化递归为递推,直接从小到大枚举每一位执行 \(p \gets p + q\)\(q\gets p – q\) 即可。实际上按任意顺序枚举位均可,因为各位独立且等价。代码实现类似 DFT 的写法。

对于沃尔什逆变换,只需在最后对每一项乘以 \(2 ^ n\) 的逆元即可。需要保证 \(2\) 的逆元存在。

  • 上述快速沃尔什变换的算法和递归算法本质相同,注意我们令 \(p = f ^ – + f ^ +\)\(q = f ^ – – f ^ +\) 相当于沃尔什变换,而 \(f \cdot g = \dfrac{p + q} 2 + x ^ {\{n\}} \dfrac {p – q} 2\) 相当于逆沃尔什变换,并将 \(2 ^ n\) 的逆元分摊在了所有 \(n\) 层当中。

P4717 【模板】快速莫比乌斯/沃尔什变换 (FMT/FWT) 部分代码。

#include <bits/stdc++.h>
using namespace std;

const int N = 1 << 17;
const int mod = 998244353;

int n, a[N], b[N], A[N], B[N], C[N];
void FWT(int *f, int n, int coef) {
	for(int k = 1; k < 1 << n; k <<= 1)
		for(int i = 0; i < 1 << n; i += k << 1)
			for(int j = 0, x, y; j < k; j++)
				x = f[i | j], y = f[i | j | k],
				f[i | j] = (x + y) % mod, f[i | j | k] = (x - y + mod) % mod;
	for(int i = 0; i < 1 << n; i++) f[i] = 1ll * f[i] * coef % mod;
}
int main() {
	cin >> n;
	for(int i = 0; i < 1 << n; i++) scanf("%d", &a[i]);
	for(int i = 0; i < 1 << n; i++) scanf("%d", &b[i]);
	memcpy(A, a, N << 2), memcpy(B, b, N << 2);
	FWT(A, n, 1), FWT(B, n, 1);
	for(int i = 0; i < 1 << n; i++) C[i] = 1ll * A[i] * B[i] % mod;
	int inv = mod + 1 >> 1;
	for(int i = 1; i < n; i++) inv = 1ll * inv * (mod + 1 >> 1) % mod;
	FWT(C, n, inv);
	for(int i = 0; i < 1 << n; i++) printf("%d ", C[i]);
	cout << endl;
	return 0;
}

4. 快速沃尔什变换

本章将从另一个角度详细讲解快速沃尔什变换的原理和推导过程,需要一点点线性代数的知识(相关信息见 线性代数学习笔记 第一章)。内容部分来自于 command_block 的博客 位运算卷积(FWT) & 集合幂级数

4.1 核心思想

依据 DFT 的思想,我们 化卷积为乘积,考虑对集合幂级数进行 线性变换 \(T\),设其标准矩阵为 \(C\),即令 \(\hat f_i = \sum\limits_{j = 0} ^ {2 ^ n – 1} C_{i, j}f_j\)。若要求 \(\hat f\)\(\hat g\) 各位相乘后等于 \(\hat f\),不难发现需要有 \(C(i, j) C(i, k) = C(i, j\oplus k)\)。这是因为

\[\begin{aligned} \hat h_i & = \sum\limits_{j = 0} ^ {2 ^ n – 1} C_{i, j} h_j \\ & = \sum\limits_{j = 0} ^ {2 ^ n – 1} \sum\limits_{k = 0} ^ {2 ^ n – 1}\sum\limits_{l = 0} ^ {2 ^ n – 1} [k \oplus l = j]C_{i, j} f_kg_l \\ & = \sum\limits_{k = 0} ^ {2 ^ n – 1}\sum\limits_{l = 0} ^ {2 ^ n – 1} C_{i, k\oplus l} f_k g_l \end{aligned} \]

\[\begin{aligned} \hat h_i & = \hat f_i \hat g_i \\ & = \left(\sum\limits_{k = 0} ^ {2 ^ n – 1} C_{i, k} f_k\right) \left(\sum\limits_{l = 0} ^ {2 ^ n – 1} C_{i, l} g_l \right)\\ & = \sum\limits_{k = 0} ^ {2 ^ n – 1}\sum\limits_{l = 0} ^ {2 ^ n – 1} C_{i, k} C_{i, l} f_k g_l\\ \end{aligned} \]

我们利用位运算的一个非常重要的性质,即 各位独立性,将 \(i, j\) 二进制分解,从而将 \(C(i, j)\) 表示成若干个 \(c(i_p, j_p)\) 相乘,其中 \(i_p, j_p\in [0, 1]\)。即设 \(i\) 的各位分别为 \(i_1, i_2, \cdots, i_p,\cdots,i_n\)\(j\) 同理,则 \(C(i, j) = \prod\limits_{p = 1} ^ n c(i_p, j_p)\)。根据各位独立性,\(C(i, j) C(i, k) = C(i, j\oplus k)\) 等价于 \(c(i_p, j_p) c(i_p, k_p) = c(i_p, j_p \oplus k_p)\),我们可以从这一结论出发,反过来推导 \(c\) 的具体值。

  • 注意,\(\oplus\) 可以是 任何位运算符,即按位与,按位或和按位异或。
  • 注意,矩阵 \(c\)存在逆元,否则没有逆变换。

考虑得到矩阵 \(c\) 之后如何快速实现变换。继续使用按位讨论的思路,设 $f = f ^ – + x ^ {{n}} \cdot f ^ + $,当 \(i\) 的第 \(n\) 位为 \(0\) 时,我们有

\[\begin{aligned} \hat f_i & = \left(c_{0, 0} \times \sum\limits_{0\leq j < 2 ^ {n – 1}}C_{i, j} f_j\right) + \left(c_{0, 1} \times \sum\limits_{2 ^ {n – 1} \leq j < 2 ^ n} C_{i, j – 2 ^ {n – 1}}f_j \right) \\ & = c_{0, 0} \hat f ^ -_j + c_{0, 1} \hat f ^ +_j \end{aligned} \]

同理,当 \(i\) 的第 \(n\) 位为 \(1\) 时,我们有 \(\hat f_i = c_{1, 0} \hat f ^ -_j + c_{1, 1} \hat f ^ +_j\)。因此,我们只需递归计算 \(\hat f ^ -\)\(\hat f ^ +\) 即可 \(\mathcal{O}(2 ^ n)\) 合并,时间复杂度 \(\mathcal{O}(n 2 ^ n)\)。可以通过从低到高按位枚举化为递推。

4.2 或卷积

\(\oplus\) 是按位或时,我们有如下构造:

  • \(0 \lor 0 = 0\),因此 \(c_{0, 0}c_{0, 0} = c_{0, 0}\),以及 \(c_{1, 0}c_{1, 0} = c_{1, 0}\)。这说明 \(c_{0, 0}\)\(c_{1, 0}\) 均等于 \(0\)\(1\),但不同时等于 \(0\),否则 \(c\) 无逆。
  • \(0 \lor 1 = 1\),因此 \(c_{0, 0}c_{0, 1} = c_{0, 1}\),以及 \(c_{1, 0}c_{1, 1} = c_{1, 1}\)。这说明 \(c_{0, 0} = 1\)\(c_{0, 1} = 0\),并且 \(c_{1, 0} = 1\)\(c_{1, 1} = 0\)
  • \(1\lor 0 = 1\),这和上面一条限制等价。
  • \(1 \lor 1 = 1\),因此 \(c_{0, 1}c_{0, 1} = c_{0, 1}\),以及 \(c_{1, 1}c_{1, 1} = c_{1, 1}\)。这说明 \(c_{0, 1}\)\(c_{1, 1}\) 均等于 \(0\)\(1\),但不同时等于 \(0\),否则 \(c\) 无逆。

根据上述限制,我们容易枚举得到所有合法的矩阵 \(c\)\(\begin{bmatrix} 1 & 1 \\ 1 & 0\end{bmatrix}\) 以及 \(\begin{bmatrix} 1 & 0 \\ 1 & 1\end{bmatrix}\)。其中第二个矩阵相当于子集求和,因为当且仅当 \(i_p = 0\)\(j_p = 1\)\(i_p \lor j_p \neq i_p\),此时 \(c(i_p, j_p) = c_{0, 1} = 0\),即不属于 \(i\) 的子集的集合对位置 \(i\) 的贡献为 \(0\),其它情况均为 \(1\)

上述两个矩阵对应的逆矩阵分别为 \(\begin{bmatrix} 0 & 1 \\ -1 & 1\end{bmatrix}\)\(\begin{bmatrix} 1 & 0 \\ -1 & 1\end{bmatrix}\),后者相当于求子集和的逆变换(互逆的标准矩阵显然对应着互逆的线性变换)。

因此,通过 FWT 实现的或卷积和通过 FMT 实现的或卷积 本质相同

对于与卷积同理,读者可自行推导,这里给出结论:常用的矩阵是 \(c = \begin{bmatrix} 1 & 1 \\ 0 & 1\end{bmatrix}\),这相当于对 \(i\) 的超集求和。

4.3 异或卷积

\(\oplus\) 是按位异或时,我们有如下构造:

  • 对于任意 \(i, j\),均有 \(c_{i, 0}c_{i, j} = c_{i, j}\)。所以 \(c_{0, 0} = c_{1, 0} = 1\)
  • \(0 \oplus 1 = 1\),因此 \(c_{0, 0}c_{0, 1} = c_{0, 1}\),以及 \(c_{1, 0}c_{1, 1} = c_{1, 1}\)。由于上面一条限制,该限制自然满足。
  • \(1 \oplus 1 = 0\),因此 \(c_{0, 1}c_{0, 1} = c_{0, 0}\),以及 \(c_{1, 1}c_{1, 1} = c_{1, 0}\)。这说明 \(c_{0, 1}\)\(c_{1, 1}\) 等于 \(1\)\(-1\),但不相等,否则 \(c\) 无逆。

综上,矩阵 \(\begin{bmatrix} 1 & -1 \\ 1 & 1 \end{bmatrix}\)\(\begin{bmatrix} 1 & 1 \\ 1 & -1 \end{bmatrix}\) 合法。其中第二个矩阵所对应的的变换和 3.3 小节所讲解的快速沃尔什变换相等,因为 \(\begin{bmatrix} 1 & 1 \\ 1 & -1 \end{bmatrix} \begin{bmatrix} p\\ q \end{bmatrix} = \begin{bmatrix} p + q \\ p – q \end{bmatrix}\)

对第二个矩阵求逆(可以直接手算),得到 \(\begin{bmatrix} 0.5 & 0.5 \\ 0.5 & -0.5 \end{bmatrix}\),恰对应着快速沃尔什逆变换。根据上文的分析,我们可以将 \(\dfrac 1 2\) 提出,放到最后乘以 \(\dfrac 1 {2 ^ n}\),以减少取模次数,减小常数。

5. 子集卷积

5.1 定义

将集合之间的二元运算定义为 不相交集合并,即 \(x ^ L \cdot x ^ R = \begin{cases} 0 & (L \cap R \neq \varnothing) \\ x ^ {L \cup R} & (L \cap R = \varnothing)\end{cases}\),我们得到了 子集卷积 的形式

\[h_S = \sum_{L \in 2 ^ U} \sum_{R \in 2 ^ U} [L \cup R = S][L \cap R = \varnothing] f_L g_R \]

它和集合并卷积的唯一区别在于两个集合不能有交。暴力枚举子集,时间复杂度 \(\mathcal{O}(3 ^ n)\)

5.2 快速算法

注意到当 \(L \cup R = S\) 时,\(L \cap R = \varnothing\)充要条件\(|L| + |R| = |S|\)。因此,我们设 \(f_{i, S} = [|S| = i]f_S\)\(g\) 同理,令 \(H_i = \sum\limits_{j + k = i} f_j \cdot g_k\)。这里 \(H_i\) 表示一个集合幂级数而非集合幂级数 \(H\) 的第 \(i\) 项,\(\cdot\) 表示集合并卷积。

根据性质我们有 \(h_S = H_{|S|, S}\),因为只有所有大小之和为 \(|S|\)\(L, R\) 才会贡献到 \(H_{|S|}\) 上,而取 \(h_{|S|}\) 位置 \(S\) 上的值保证了贡献到该位置上的集合 \(L, R\) 的并为 \(S\),因此符合限制。

将等式两边同时取 FMT,得到 \(\hat H_i = \sum\limits_{j + k = i} \hat f_j \cdot \hat g_k\),这里 \(\cdot\) 表示对应位置相乘。因此,我们首先 \(\mathcal{O}(n ^ 2 2 ^ n)\) 计算出每个 \(f_j\)\(g_k\) 的莫比乌斯变换(共有 \(n\) 个集合幂级数,每做一次 FMT 需要 \(\mathcal{O}(n 2 ^ n)\) 的时间),然后 \(\mathcal{O}(n ^ 2 2 ^ n)\) 按式子计算出 \(\hat H_i\)(共要做 \(n ^ 2\) 次对应位置相乘,每次需要 \(\mathcal{O}(2 ^ n)\) 的时间),最后再将每个 \(\hat H_i\) 逆变换回来即可。三个部分时间复杂度均为 \(\mathcal{O}(n ^ 2 2 ^ n)\),因此这也是总时间复杂度。

P6097【模板】子集卷积 代码如下。

#include <bits/stdc++.h>
using namespace std;

inline int read() {
	int x = 0; char s = getchar();
	while(!isdigit(s)) s = getchar();
	while(isdigit(s)) x = x * 10 + s - '0', s = getchar();
	return x;
}
inline void print(int x) {
	if(x >= 10) print(x / 10);
	putchar(x % 10 + '0');
}

const int N = 20, mod = 1e9 + 9;
long long f[N + 1][1 << N], g[N + 1][1 << N], h[1 << N], ans[1 << N];
int n;
template <class T> void FMT(T *f, int op) {
	for(int j = 1; j < 1 << n; j <<= 1)
		for(int i = 0; i < 1 << n; i++)
			if(i & j)
				f[i] += op ? f[i ^ j] : -f[i ^ j];
	for(int i = 0; i < 1 << n; i++) f[i] = (f[i] % mod + mod) % mod;
}
int main() {
	cin >> n;
	for(int i = 0; i < 1 << n; i++) f[__builtin_popcount(i)][i] = read();
	for(int i = 0; i < 1 << n; i++) g[__builtin_popcount(i)][i] = read();
	for(int i = 0; i <= n; i++) FMT(f[i], 1), FMT(g[i], 1);
	ans[0] = f[0][0] * g[0][0] % mod;
	for(int i = 1; i <= n; i++) {
		memset(h, 0, sizeof(h));
		for(int j = i; ~j; j--) {
			for(int p = 0; p < 1 << n; p++) h[p] += f[j][p] * g[i - j][p];
			if(j % 9 == 0) for(int p = 0; p < 1 << n; p++) h[p] %= mod;
		}
		FMT(h, 0);
		for(int p = 0; p < 1 << n; p++) if(__builtin_popcount(p) == i) ans[p] = h[p];
	}
	for(int p = 0; p < 1 << n; p++) print(ans[p]), putchar(' ');
	return 0;
}

6. ln 和 exp

该部分内容过于高深,笔者能力有限(笔者对形式幂级数的 ln 和 exp 都不甚了解),因此挖坑待填。

7. 例题

I. CF1530F Bingo

非正解。考虑补集转化,求没有任何一行 / 列(将对角线看成两个特殊的列)的事件全部发生的概率。设 \(f_{i, S}\) 表示前 \(i\) 行已经有集合 \(S\) 中的列和对角线存在一个事件没有发生。转移则考虑枚举当前行的 \(2 ^ n\) 个状态,设没有发生事件的列和对角线为 \(S\),需满足 \(S \neq \varnothing\) 否则这一行的事件就全部发生了,不符合限制。求出 \(g_S\) 表示没有发生事件的列和对角线恰好为 \(S\) 的可能性(尽管 \(S\)\(2 ^ {n + 2}\) 种可能,但 \(g\) 只有 \(2 ^ n – 1\) 个位置上有值,因为它依赖于第 \(i\) 行的 \(2 ^ n – 1\) 个合法状态),转移 \(f_i \cdot g \to f_{i + 1}\) 是 OR 卷积的形式。最终答案即 \(1 – f_{n, 2 ^ {n + 2} – 1}\),因为所有列和对角线都需要存在一个事件没有发生。

时间复杂度 \(\mathcal{O}(2^nn^2)\)\(4\) 倍常数。把 long long 换成 int 加上 7s 时限,勉强能过。代码

II. ABC212H Nim Counting

\(f_i\) 表示 \(i\)\(A\) 中出现的次数。根据 NIM 游戏的结论,先手获胜当且仅当初始石子异或和个数异或和不为 \(0\)。故答案为 \(\sum\limits_{i = 1} ^ n \mathrm{IFWT} (\mathrm{FWT} ^ i(f))\)\(x ^ 1 \sim x ^ {2 ^ {16} – 1}\) 处的取值。根据 FWT 的性质,变成 \(\mathrm{IFWT} \left(\sum\limits_{i = 1} ^ n\mathrm{FWT} ^ i (f)\right)\) 即可用等比数列求和计算,注意特判 \(f_i = 1\) 的情况等比数列求和不管用。时间复杂度线性对数。代码

*III. 某联考题 /kk

\(k\)\(A_{i,j}\gets A_{i,j}\oplus A_{i,j+1}\oplus A_{i+1,j}\oplus A_{i+1,j+1}\ (0\leq i<n,0\leq j<m)\) 的变换,求 \(A_{0,0}\)\(q\) 组询问。

\(nm \leq 5\times 10 ^ 6\)\(a_{i, j}, k < 2 ^ {31}\)\(q \leq 5 \times 10 ^ 7\)。时间限制 2.5s,空间限制 1GB。

\((i, j)\) 对于 \((0, 0)\) 的贡献次数是 \(\dbinom k i \dbinom k j\),因为 \(k\) 次操作里有 \(i\) 次选行 \(+1\)(剩下选 \(0\)),\(j\) 次选列 \(+1\)。根据 Lucas 定理,当且仅当 \(k \& i = i\)\(k \& j = j\)\(k \& (i | j) = i | j\) 时上式为奇数。高维前缀异或和即可。时间复杂度 \(\mathcal{O}(\max(n, m) \log\max(n, m) + q)\)

IV. AT4168 [ARC100C] Or Plus Max

ARC 也会出裸题?

记录 \(f_S, g_S\) 表示 \(a_i\ (i\in S)\) 的最大值和次大值,用高维前缀和做。\(K = i\) 时的答案 \(ans_i\)\(\max(ans_{i – 1}, f_i + g_i)\)代码

*V. P4221 [WC2018]州区划分

首先判断一个点集 \(S\) 是否合法:若不连通显然合法,否则根据欧拉定理可知,若原图以 \(S\) 为点集的点导出子图中存在度数为奇数的点则合法,反之不合法。

\(q_S\) 表示 \(\sum_{i \in S} w_i\)\(f_S\) 表示点集为 \(S\) 时的满意度之和,则转移方程可写为

\[f_S=\sum_{T\subset S}f_T\times\left(\dfrac{q_{S\setminus T}}{q_S}\right)^p=\dfrac{1}{q_S^p}\sum_{T\subset S}f_Tq_{S\setminus T}^p \]

一个显然的子集卷积形式。

这个 \(\dfrac{1}{q_S^p}\) 不太好处理,我们可以这样:子集卷积需要额外记录一维 \(i\) 表示集合大小,而 \(i\) 只会从比它小的 \(j\) 转移过来,因此我们可以按 \(i\) 从小到大 DP,每 DP 完一层就将该层的 \(f_S\ (|S| = i)\) 做一遍 FMT。这样就不会出现系数的问题了。

时间复杂度 \(\mathcal{O}(2^nn^2)\)代码

VI. CF914G Sum the Fibonacci

三合一题。

注意到 \(s_a | s_b = i\)\(s_a\&s_b = 0\) 是子集卷积,\(s_d \oplus s_e = i\) 是异或卷积,外层是与卷积。于是暴力卷起来即可。时间复杂度 \(\mathcal{O}(n + s\log ^ 2 s)\)代码

VII. P5387 [Cnoi2019]人形演舞

打表求出每个数的 SG 值,发现是 \(i – 2 ^ {\lfloor \log_2 i\rfloor} + 1\)。根据 nim 游戏 SG 值非零先手必胜的结论,直接 FWT 即可,时间复杂度 \(m(\log p + \log m)\)

VIII. CC21JAN ORAND

将存在性问题转化为计数问题,直接 FWT 即可。一直 FWT 到答案不再更新为止,复杂度为 \(\mathcal{O}(TV \log ^ 2 V)\),因为最多进行 \(\log V\) 次位运算卷积(考虑每次卷积必然更新至少一位)。实际远远达不到上界,跑得很快。代码

8. 参考文章