FWT(快速沃尔什变换) 学习笔记
前言
因为最近模拟赛考到了,学习了一下 \(FWT\) 的相关知识。感觉还是非常有用的。
参考文章 点这里 。 焦神还是太强了!!!
位运算卷积
一般的卷积形如: \(h_{i} = \sum\limits_{0 \leq j \leq i} f_j g_{i - j}\)。我们将它称作 \((+, \times)\) 卷积。
闵可夫斯基和可以维护特殊的卷积: \(h_i = \max\limits_{0 \leq j \leq i}(f_j + g_{i - j})\)。我们将它称作 \((\max, +)\) 卷积。
由此可以发现 \((*, *)\) 卷积的一个符号是 贡献方式,第二个符号是两个数组的 运算方式。而下标变换一般默认是 求和。
那么 \(FWT\) 维护的是什么形式的卷积呢?
对于下标为 \(0 \sim 2^n - 1\) 的数组 \(f,g, h\)而言,
形如: \(h_{i} = \sum_{j \oplus k = i} f_jg_k\),其中 \(\oplus\) 为某种位运算。
那么贡献方式是 \(+\),运算方式是 \(\times\),只不过下标变换不是求和而是位运算了。可以称作 \((+, \times )\) 位运算卷积。
\(FWT\) 的基本思路:
- 找到一种变换 \(FWT(A)\) 满足 \(FWT(h) = FWT(f) \cdot FWT(g)\)。这里的点乘是每位分别乘的意思。
- 在 \(O(n \log n)\) 的复杂度求出 \(FWT(f)\) 和 \(FWT(g)\)。
- 在 \(O(n)\) 的复杂度点乘求出 \(FWT(h)\)。
- 在 \(O(n \log n)\) 的复杂度将 \(FWT(h)\) 还原为 \(h\)。
我们希望 \(FWT(A)\) 是一个线性变换,也就是满足:
\(FWT(A + B) = FWT(A) + FWT(B)\),\(kFWT(A) = FWT(kA)\)。
假设我们想让 \(A, B\) 数组位运算卷积得到 \(C\):
不妨设 \(FWT(A)_i = \sum\limits_{j}f(i, j)A_j\),其中 \(f(i, j)\) 表示 \(j\) 对 \(i\) 的贡献系数。
那么有:
又因为:
所以
那么对比式子,得到: \(f(i, j)f(i, k) = f(i, j \oplus k)\)。
注意我们现在并没有用到 位运算的性质,因此上述关系对于多项式同样是适用的。
那么现在我们的任务就是构造一个 \(2^n \times 2^n\) 的系数矩阵 \(f\),满足 \(f(i, j)f(i, k) = f(i, j \oplus k)\) 的性质,那么向量 \(A \times f = FWT(A)\)。由于我们还要还原,因此矩阵 \(f\) 还需要有逆矩阵。并且还要能快速求出 \(A \times f\) 的结果。
用上位运算的性质:每位运算独立。尝试给 \(f\) 加上更强的限制:任意两个数的 \(f(i, j)\) 都等于它们每一位上的数的 \(f\) 乘起来。
这样有什么好处呢?如果满足这个性质那么我们可以快速的求出 \(A \times f\):
考虑一个分治的过程:
\(FWT(A)_i = \sum\limits_{j = 0}^{2^{n - 1}-1}f(i, j)A_j + \sum\limits_{j = 2^{n - 1}}^{2^n - 1}f(i, j)A_j\)
\(FWT(A)_i = f(c, 0)\sum\limits_{j = 0}^{2^{n - 1} - 1}f(i',j')A_j + f(c, 1)\sum\limits_{j = 2^{n - 1}}^{2^n - 1}f(i', j')A_j\)
其中 \(c\) 为 \(A_i\) 最高位上的数字,\(i',j'\) 为 \(i,j\) 去掉最高位后的数。
这样可以转化成一个子问题,可以用线段树的结构去理解。
这样分治的时间复杂度是 \(O(n \log_2 n)\) 的。
如何还原?
称还原的过程为 \(IFWT\)。
由于 \(A \times f = FWT(A)\),因此 \(FWT(A) \times f^{-1} = A\)。
那么如何构造 \(f^{-1}\) 呢?
具体证明我也不懂:对于 \(f\) 在一位上的 \(2 \times 2\) 矩阵,求出它的逆就是 \(f^{-1}\) 在每一位上的矩阵。也就是 \(f^{-1}\) 同样满足数值等于每一位上分别求 \(f^{-1}\) 再乘起来。
那么只需要将这个 \(2 \times 2\) 的矩阵求出逆矩阵后再跑一遍 \(FWT\) 就能还原数组了。
还有一个需要注意的点:
由于每一位的贡献是独立的,那么我们可以以任意顺序求出每一位的贡献。
可以理解成一个数 \(x\) 对一个数 \(y\) 的贡献系数只取决于每一位上两个数字的 \(f\) 的乘积。如果将每一位看作一个维度,每一位上的转移看作在相应维度上走了一步,那么无论先走哪个维度,最后 \(x\) 是一定能走到 \(y\) 的,并且这个贡献系数在走的过程中累乘起来也是对的。
基础位运算卷积
或卷积
要满足 \(f(i, j)f(i, k)=f(i, j|k)\)。
注意到如果 \(i|j = i\),\(i|k = i\),那么 \(i|(j|k) = i\)。
因此令 \(f(i, j) = [i|j = i]\),也就是 \(j\) 是 \(i\) 的子集。
可以拆成每位的乘积,构造 \(2 \times 2\) 的矩阵为:
相当于分治左区间对右区间对应位置的贡献是 \(1\)。
它的逆为:
相当于分治左区间对右区间对应位置的贡献是 \(-1\)。
可以合并成一个函数,代码如下:
inline void OR(LL *f, int n, LL opt) { // or 卷积
for(int o = 2, k = 1; o <= n; o <<= 1, k <<= 1)
for(int i = 0; i < n; i += o)
for(int j = 0; j < k; j ++ )
f[i + j + k] = (f[i + j + k] + f[i + j] * opt % mod) % mod;
}
与卷积
要满足 \(f(i, j)f(i, k) = f(i, j\&k)\)
注意到如果 \(i \& j = i\),\(i \& k = i\),那么 \(i \& (j \& k) = i\)
因此令 \(f(i, j) = [i \& j = i]\),也就是 \(j\) 是 \(i\) 的超集。
每位上 \(2 \times 2\) 的矩阵:
相当于右边对左边的系数是 \(1\)。
逆矩阵:
相当于右边对左边的系数是 \(-1\)。
同样可以合并成一个函数:
inline void AND(LL *f, int n, LL opt) {
for(int o = 2, k = 1; o <= n; o <<= 1, k <<= 1)
for(int i = 0; i < n; i += o)
for(int j = 0; j < k; j ++ )
f[i + j] = (f[i + j] + f[i + j + k] * opt % mod) % mod;
}
异或卷积
异或比较麻烦。
还是注意到 \((-1)^{|i \& j|}(-1)^{|i \& k|} = (-1)^{|i \& (j \oplus k)|}\)(注意力惊人 )。
这是因为如果 \(i, j\) 在某一位上的数相同,在这一位上 \(i \& j\) 与 \(i \& k\) 一定相同,因此乘起来一定是 \(1\),可以直接异或消掉。
如果不同,那么这一位是 \((-1)\) 等价于 \(i\) 在这一位上是 \(1\),左右仍然相等。
由此也能得到一个十分重要的式子:
构造矩阵:
相当于
逆矩阵比较特殊:
符号每边,系数变成了 \(\frac{1}{2}\)。
也可以合并成一个函数:
inline void XOR(LL *f, int n, LL opt) {
for(int o = 2, k = 1; o <= n; o <<= 1, k <<= 1)
for(int i = 0; i < n; i += o)
for(int j = 0; j < k; j ++ )
f[i + j] = (f[i + j] + f[i + j + k]) % mod,
f[i + j + k] = (f[i + j] - 2 * f[i + j + k] % mod + mod) % mod,
f[i + j] = f[i + j] * opt % mod, f[i + j + k] = f[i + j + k] * opt % mod;
}
最后的封装版:
namespace FWT {
inline void OR(LL *f, int n, LL opt) { // or 卷积
for(int o = 2, k = 1; o <= n; o <<= 1, k <<= 1)
for(int i = 0; i < n; i += o)
for(int j = 0; j < k; j ++ )
f[i + j + k] = (f[i + j + k] + f[i + j] * opt % mod) % mod;
}
inline void AND(LL *f, int n, LL opt) { // and 卷积
for(int o = 2, k = 1; o <= n; o <<= 1, k <<= 1)
for(int i = 0; i < n; i += o)
for(int j = 0; j < k; j ++ )
f[i + j] = (f[i + j] + f[i + j + k] * opt % mod) % mod;
}
inline void XOR(LL *f, int n, LL opt) { // xor 卷积
for(int o = 2, k = 1; o <= n; o <<= 1, k <<= 1)
for(int i = 0; i < n; i += o)
for(int j = 0; j < k; j ++ )
f[i + j] = (f[i + j] + f[i + j + k]) % mod,
f[i + j + k] = (f[i + j] - 2 * f[i + j + k] % mod + mod) % mod,
f[i + j] = f[i + j] * opt % mod, f[i + j + k] = f[i + j + k] * opt % mod;
}
}
子集卷积
给你两个长度为 \(2^n\) 的数组 \(A,B\),你需要求出数组 \(C\),满足:
相当于 \(j, k\) 是 \(i\) 的两个不交子集并且它们的并等于 \(i\)。
\(0 \leq n \leq 20\)。
首先可以枚举子集暴力计算子集卷积,时间复杂度 \(O(3^n)\),但这不是我们想要的。
如果只有 \(j | k = i\) 这一个条件,那就是正常的或卷积。注意到满足 \(j | k = i\) 条件时, \(j \& k = 0\) 等价于 \(|j| + |k| = |i|\)。因此我们可以多记一维表示大小。
设 \(A_{i, j} = A_j[|j| = i]\),\(B_{i, j} = B_j[|j| = i]\),\(C_{i, j} = \sum\limits_{L|R = j,|L|+|R| =i}A_LB_R\)。
那么 \(C_{|i|, i}\) 就是答案。
考虑 \(C_{i, j}\) 怎么计算:
设 \(C'_{i, j} = \sum\limits_{p | j = j}C_{i, p}\), \(A'_{i, j} = \sum\limits_{p | j = j}A_{i, p}\), \(B'_{i, j} = \sum\limits_{p | j = j}B_{i, p}\)
那么有 \(C'_{i, j} = \sum\limits_{k \leq i}A'_{k, j}B'_{i -k, j}\)。
原因是这样:
考虑数组的实际意义:那么 \(A'_{k, j}\) 就是所有大小为 \(k\) 且是 \(j\) 的子集的 \(A_p\) 之和, \(B'_{i -k, j}\) 就是所有大小为 \(i - k\) 且是 \(j\) 的子集的 \(B_p\) 之和。那么它们相乘就得到了 所有是 \(j\) 子集的 \(C_{i, p} 之和\),这与 \(C'_{i, p}\) 的含义相同。
求 \(A'_{i, j}\) 和 \(B'_{i, j}\) 就是做 \(n\) 次 \(FWTor\),时间复杂度 \(O(n^2 2^n)\)。
求 \(C'_{i, j}\) 就是枚举 \(k\) 然后点乘,时间复杂度 \(O(n^2 2^n)\)。
最后还原 \(C_{i, j}\) 是做 \(n\) 次 \(IFWT\),复杂度也是 \(O(n^2 2^n)\)。
总时间复杂度 \(O(n^2 2^n)\),空间复杂度 \(O(n2^n)\)。
CODE:
// 子集卷积模板
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL mod = 1e9 + 9;
inline int sz(int x) {return __builtin_popcount(x);}
int n;
LL A[21][1 << 20], B[21][1 << 20], C[21][1 << 20];
namespace FWT {
inline void OR(LL *f, int n, LL opt) {
for(int o = 2, k = 1; o <= n; o <<= 1, k <<= 1)
for(int i = 0; i < n; i += o )
for(int j = 0; j < k; j ++ )
f[i + j + k] = (f[i + j + k] + f[i + j] * opt % mod) % mod;
}
}
int main() {
ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
cin >> n;
for(int i = 0; i < (1 << n); i ++ ) cin >> A[sz(i)][i];
for(int i = 0; i < (1 << n); i ++ ) cin >> B[sz(i)][i];
for(int i = 0; i <= n; i ++ ) FWT::OR(A[i], (1 << n), 1);
for(int i = 0; i <= n; i ++ ) FWT::OR(B[i], (1 << n), 1);
for(int i = 0; i <= n; i ++ )
for(int j = 0; j <= i; j ++ )
for(int k = 0; k < (1 << n); k ++ )
C[i][k] = (C[i][k] + A[j][k] * B[i - j][k] % mod) % mod;
for(int i = 0; i <= n; i ++ ) FWT::OR(C[i], (1 << n), mod - 1);
for(int i = 0; i < (1 << n); i ++ ) cout << C[sz(i)][i] << ' ';
cout << endl;
return 0;
}
\(K\) 进制卷积
只有构造 \(f\) 是不同的,正常考虑一位的 \(f\) 怎么构造就行。
\(K\) 进制 \(min\) 卷积
含义为:每位取 \(min\)
可以看作 与卷积 的拓展,因为 与卷积实际上就是二进制下的 \(min\) 卷积。
有一种合法构造为: \(f(a, b) = [a \leq b]\)。以 \(K = 4\) 为例,可得到一下矩阵:
逆矩阵为:
\(K\) 进制 \(max\) 卷积
含义为: 每位取 \(max\)
同样可以看作 或卷积 的拓展。
有一种合法构造为:\(f(a, b) = [a \geq b]\),以 \(K = 4\) 为例:
逆矩阵为:
\(K\) 进制异或卷积(不进位加法)
含义:每位相加后对 \(K\) 取模。
也就是需要满足 \(f(a, b)f(a, c) = f(a, (b + c) \% K)\)。
有一种合法构造为:\(f(a, b) = w_{K}^{b}\)。但是这样会导致矩阵每一行都一样,因此我们调整为 \(f(a, b) = w_{K}^{ab}\)。这里我们认为矩阵的行和列都从 \(0\) 开始编号。
那么可以得到以下矩阵:
发现矩阵很熟悉?没错,就是 \(FFT\) 中的 范德蒙德矩阵,它的逆比较特殊,为:
但是这样有一个问题: 如果答案要取模,那么单位根在模意义可能不存在逆元!!
如果存在逆元,那么变成 \(NTT\) 那样即可。
否则可能需要用到 扩域。
这个科技我也不懂,放一篇教学和一道题目:教学,题目。有生之年再学吧。
练习题
BZOJ4589 Hard Nim
CF662C Binary Table
WC2018 州区划分
THUPC 2019 找树 / findtree
题意:
定义一个作用在 \(w\) 位二进制数上的新位运算 \(\oplus\):\(a \oplus b\) 的每一位等于 \(a_i \oplus_i b_i\) 。其中 \(a_i, b_i\) 表示 \(a,b\) 在第 \(i\) 位上的数字。 \(\oplus_i\) 表示运算 \(\oplus\) 在第 \(i\) 位上的运算符,其中 \(oplus_i \in \{或,与,异或\}\)。 运算 \(\oplus\) 最开始读入一个长度位 \(w\) 的字符串得到。
给出一张 \(n\) 个点 \(m\) 条边的无向图,每一条边的边权是一个 \(w\) 位的二进制数。请找出原图的一棵生成树。设 \(v_1,v_2,...,v_{n - 1}\) 为生成树的边权,你需要最大化 \(v1 \oplus v_2 ... \oplus v_{n - 1}\)。
\(1 \leq n \leq 70\),\(1 \leq m \leq 5000\), \(1 \leq w \leq 12\)。
分析:
神题。感觉每一步都比较难想。但是如果见过类似的那么就比较套路。
首先这题不是一个最优化问题,考虑对每种权值的生成树计数,最后的答案就是最大的权值满足这种权值的生成树数量大于 \(0\)。
由于要对生成树计数,考虑 矩阵树定理。 矩阵树定理不会的可以参考这个博客
我们知道矩阵树定理求出的是 \(\sum\limits_{T}\prod\limits_{e \in T}w_e\),如果 \(w_e\) 等于 \(1\) 时求出的就是生成树的数量。那么怎么对每种权值的生成树都求出答案呢?
首先需要知道矩阵树定理的适用条件:
设边权的数集为 \(W\)(这里的 \(W\) 可以是整数集,实数集,复数集。也可以是 多项式集 甚至自己定义的元素集)。设 \(W\) 中的加法运算为 \(+\),乘法运算为 \(\times\)。(这里的 \(+,\times\) 的运算法则都可以自己定义) 那么只要满足 \((W, +, \times)\) 构成一个 环,矩阵树定理就成立。
什么是环? 通俗的说就是元素之间的 \(+\) 运算满足 交换律,结合律, \(\times\) 运算满足 结合律, \(\times\) 对 \(+\) 可分配。
那么接着来看上面的问题:套路的,考虑将边权变成一个集合幂级数(多项式)。
对于一条边权为 \(v\) 的边,将它的边权修改成 \(x^v\)。那么如果此时的度数矩阵 \(D\) 和邻接矩阵 \(E\) 的每一项都变成了多项式,构造得到 \(Lapace\) 矩阵(\(L = D - E\)) 的每一项也变成了多项式。定义多项式的 \(\times\) 运算为 多项式卷积,\(+\) 运算为逐位累加。那么这时候求出的行列式也是一个多项式,并且 \(x^v\) 前面的系数就是 边权和 为 \(v\) 的生成树的数量。
多项式矩阵不能 高斯消元,怎么求行列式? 实际上,由于上面的 \(\times\) 和 \(+\) 是正常的多项式运算,可以通过 拉插 还原系数的方式将矩阵变成整数矩阵。还有一种方法下面会说。
回到本题,一棵生成树的权值不是 边权和,而是边权 \(\oplus\) 和,怎么办?
也就是 \(x^{v_1}\) 和 \(x^{v_2}\) 的系数乘积应该贡献给 \(x^{v_1 \oplus v_2}\),那么就是 \(h_{x} = \sum\limits_{a \oplus b = x}f_ag_b\) 的形式。因此需要将多项式卷积变成 位运算卷积。
我们修改 \(\times\) 的运算法则为两个多项式做 位运算卷积。那么求出的行列式 \(x^v\) 前的系数就是边权经过位运算后为 \(v\) 的生成树数量。
考虑怎么求行列式:
首先先考虑 \(n!\) 求行列式的方法,枚举排列,发现每次相当于是若干个多项式做位运算卷积,然后把每一位的系数对应加起来。
那么 \(FWT\) 之后就变成了每位上点乘,最后 \(IFWT\) 后将每一位系数加起来。
但是考虑到 \(FWT\) 为线性变换,满足 \(IFWT(A + B) = IFWT(A) + IFWT(B)\)。也就说,可以将所有排列得到的多项式对应位加起来之后再 \(IFWT\) 回去。
发现这时候每一位是独立的!!!可以每次将某一位上的系数拿出来单独求一遍行列式,此时就可以高斯消元了!
刚才说的正常多项式卷积也可以用这个思路求: \(FFT\) 之后就变成了点值相乘,每一位单独拿出来求行列式即可。
还有一个问题:每一位上的运算不同,怎么 \(FWT\)?
由于求 \(FWT\) 时每一位的贡献是独立的,因此每一位用对应运算的系数矩阵即可。不懂可以看看代码。
由于方案数太多,因此需要对大质数取模。 总时间复杂度就是 \(O(n^32^w)\) 的,虽然很大,但我没咋卡常就过了。
// so hard
// 第一步:想到对每种权值的生成树计数, 权值最大的并且数量不为 0 的就是答案
// 第二步:对生成树计数,考虑矩阵树定理。套路的,我们要对每种权值的生成树求出数量,因此将边权转化为集合幂级数(多项式) v -> x^v
// 第三步:修改多项式相乘和相加的定义: 由于矩阵树定理求的是 Σ_{T} prod_{e \in T} w_e, 一棵树的权值应该是 x^{v1 \oplus v2 ...},因此将 × 定义为位运算卷积,将 + 定义为逐位相加
// 第四步:考虑每一个元素为多项式的行列式怎么计算结果:这个应该是比较套路的。先考虑 n! 计算行列式的方法,发现是若干多项式进行位运算卷积,那么将每个多项式FWT之后就变成了按位点乘,然后在按位加起来,最后在还原(线性变换)。那么这时候每一位上就是独立的,因此单独拿出来高斯消元
// 补充:至于位运算卷积变成了广义FWT,由于每一位贡献独立因此正常做就行
#include<bits/stdc++.h>
using namespace std;
const int N = 75;
const int M = 5010;
const int mod = 1e9 + 7;
int n, m, w;
int bin[N];
char str[N];
int a[M], b[M], v[M];
int c[3][2] = {{1, mod - 1}, {1, mod - 1}, {1, 500000004}};
inline int Del(int x, int y) {return x - y < 0 ? x - y + mod : x - y;}
inline int Add(int x, int y) {return x + y >= mod ? x + y - mod : x + y;}
struct F { // 多项式
int f[1 << 12];
friend F operator - (F a, F b) {
F c; for(int i = 0; i < (1 << w); i ++ ) c.f[i] = Del(a.f[i], b.f[i]);
return c;
}
};
F deg[N][N], E[N][N]; // 度数矩阵和邻接矩阵
F ans;
namespace FWT {
inline void transform(int *f, int n, int opt) {
for(int o = 2, k = 1, l = 0; o <= n; o <<= 1, k <<= 1, l ++ )
for(int i = 0; i < n; i += o )
for(int j = 0; j < k; j ++ )
if(str[l] == '|') f[i + j + k] = Add(f[i + j + k], 1LL * f[i + j] * c[0][opt] % mod);
else if(str[l] == '&') f[i + j] = Add(f[i + j], 1LL * f[i + j + k] * c[1][opt] % mod);
else
f[i + j] = Add(f[i + j], f[i + j + k]),
f[i + j + k] = Add(f[i + j], 1LL * (mod - 2LL) * f[i + j + k] % mod),
f[i + j] = 1LL * f[i + j] * c[2][opt] % mod,
f[i + j + k] = 1LL * f[i + j + k] * c[2][opt] % mod;
}
}
namespace Det { // 计算行列式
int a[N][N];
inline int sign(int x) {return (x & 1) ? mod - 1 : 1;}
inline int Pow(int x, int y) {
int res = 1, k = x;
while(y) {
if(y & 1) res = 1LL * res * k % mod;
y >>= 1;
k = 1LL * k * k % mod;
}
return res;
}
int calc(int n) { // 高斯消元
int c = 0; int res = 1;
for(int i = 1; i <= n; i ++ ) {
int id = i;
for(int j = i; j <= n; j ++ )
if(a[j][i] != 0) id = j;
if(id != i) c ++;
for(int j = 1; j <= n; j ++ ) swap(a[i][j], a[id][j]);
int inv = Pow(a[i][i], mod - 2LL); res = 1LL * res * a[i][i] % mod;
for(int j = 1; j <= n; j ++ ) a[i][j] = 1LL * a[i][j] * inv % mod;
for(int j = i + 1; j <= n; j ++ ) {
int tmp = a[j][i];
for(int k = 1; k <= n; k ++ )
a[j][k] = Del(a[j][k], 1LL * a[i][k] * tmp % mod);
}
}
return 1LL * res * sign(c) % mod;
}
}
int Find(int x) {return x == bin[x] ? x : bin[x] = Find(bin[x]);}
inline void Merge(int x, int y) {
int f1 = Find(x), f2 = Find(y);
if(f1 == f2) return ;
bin[f1] = f2;
}
int main() {
scanf("%d%d", &n, &m);
scanf("%s", str); w = strlen(str);
for(int i = 1; i <= n; i ++ ) bin[i] = i;
for(int i = 1; i <= m; i ++ ) {
scanf("%d%d%d", &a[i], &b[i], &v[i]);
deg[a[i]][a[i]].f[v[i]] ++; deg[b[i]][b[i]].f[v[i]] ++;
E[a[i]][b[i]].f[v[i]] ++; E[b[i]][a[i]].f[v[i]] ++;
Merge(a[i], b[i]);
}
int ct = 0;
for(int i = 1; i <= n; i ++ ) ct += (Find(i) == i);
if(ct > 1) {puts("-1"); return 0;}
for(int i = 1; i <= n; i ++ )
for(int j = 1; j <= n; j ++ )
deg[i][j] = deg[i][j] - E[i][j];
for(int i = 2; i <= n; i ++ ) // 去掉一行一列
for(int j = 1; j <= n; j ++ )
deg[i - 1][j] = deg[i][j];
for(int i = 2; i <= n; i ++ )
for(int j = 1; j < n; j ++ )
deg[j][i - 1] = deg[j][i];
n --;
for(int i = 1; i <= n; i ++ )
for(int j = 1; j <= n; j ++ )
FWT::transform(deg[i][j].f, 1 << w, 0); // 0 是正变化,1 是变换回来
for(int i = 0; i < (1 << w); i ++ ) { // 点乘, 每一位是独立的,拿出来高斯消元
for(int j = 1; j <= n; j ++ )
for(int k = 1; k <= n; k ++ )
Det::a[j][k] = deg[j][k].f[i];
ans.f[i] = Det::calc(n);
}
FWT::transform(ans.f, 1 << w, 1);
int ret = 0;
for(int i = 0; i < (1 << w); i ++ )
if(ans.f[i] > 0) ret = max(ret, i);
cout << ret << endl;
return 0;
}
ARC132F Takahashi The Strongest
题意:\(A, B, C\) 三个人在玩剪刀石头布游戏,一共玩了 \(K\) 局。设 \(P\) 表示布,\(R\) 表示石头, \(S\) 表示剪刀。
给你 \(n\) 个长度为 \(K\) 只有 \(\{P,R,S\}\)字符串表示 \(A\) 的策略集合 和 \(m\) 个长度为 \(K\) 的字符串表示 \(B\) 的策略集合。
一个策略表示每局出什么手势。\(A,B\) 会从他们的策略集合中挑一种作为自己的策略。
对 \(C\) 的 \(3^K\) 种策略,问 \(nm\) 种组合中有多少种使得 \(C\) 在这种策略下至少有一局完全赢过 \(A,B\)。
完全赢过的意思是 \(B,C\) 手势相同且 \(A\) 能胜 \(B\)。
\(1 \leq K \leq 12\), \(1 \leq n,m \leq 3^ K\)。
分析:
首先将 \(P,R,S\) 分别看作 \(0,1,2\)。定义 \(to_i\) 为能赢过手势 \(i\) 的手势编号。
设 \(f(S) = to_{S_1} + to_{S_2}\dots + to_{S_k}\),相当于就是每一局都能胜过策略 \(S\) 的策略。
那么 \(A,B\) 的一个策略组合 \(a,b\) 只有在 \(f(a) \cap f(b) \ne \emptyset\) 时才会对某个状态右贡献。
我们希望先将 \(A,B\) 的每种组合贡献到它们的交上。
定义新运算 \(\oplus\) 为
可以将 \(3\) 看作这一位上交集为空的标识。
那么首先要求的就是 \(h_{i} = \sum\limits_{j \oplus k = i} f_i g_j\)。 \(f_i, g_i\) 分别表示 \(A,B\) 策略集合中 \(i, j\) 策略的数量。
如果能求出 \(h_i\),那么只需要将每个 \(h_i\) 贡献到与它有交的状态 \(ans_j\) 上即可。
发现 \(h_i\) 是可以位运算卷积求出的。我们套用理论,想要构造一个 \(4 \times 4\) 的系数矩阵满足 \(f(i, j)f(i, k) = f(i, j \oplus k)\)。
那么可以定义 \(f(i, j) = [i = j 或 i = 3]\),不难发现是满足条件的。
也可以理解为规定 \(3 \subset 0, 3 \subset 1, 3 \subset 2\),\(1,2, 3\) 之间无子集关系。那么相当于一个求超集贡献。
然后可以 \(FWT\) 之后点乘,然后 \(IFWT\) 还原。
接下来考虑怎么求 \(h\) 对 \(ans\) 的贡献:
我们想让\(h_j\) 贡献到与 \(i,j\) 有交的 \(ans_i\) 上:考虑二项式反演的形式容斥。
假设交集大小为 \(n\),那么有 \(\sum\limits_{i = 0}^{n}(-1)^i\binom{n}{i} = (1 - 1)^n = [n = 0]\)。
那么我们想让交集的每个子集都以 \((-1)^{sz}\) 为系数贡献。考虑先对 \(h\) 求一遍超集贡献,得到 \(h'\),然后令 \(h‘_{S} \gets h’_{s} \times (-1)^{|S|}\)。
接着只要对 \(h'\) 做一遍子集贡献就可以得到与每个状态不交的方案数,拿 \(nm\) 减去这个数就是答案,
注意到 \(h'\) 就是点乘之后的结果,因此点乘之后直接乘容斥系数然后做一遍求子集贡献的 \(FWT\) 即可。
时间复杂度 \(O(4^k \times k)\)。
CODE:
// P -> 0, R -> 1, S -> 2。 不匹配 -> 3
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 540000;
char str[20];
int val[200], to[5];
int K, n, m, a[N], b[N], sz[1 << 24], mi[13];
LL f[1 << 24], g[1 << 24];
inline int sign(int x) {return (x & 1) ? -1 : 1;}
inline int turn(int x) { // 把一个三进制数转化成四进制数
int ret = 0, p = 0;
while(x) {
ret = ret + (x % 3) * (1 << p * 2);
x /= 3; p ++;
}
return ret;
}
namespace FWT { // 四进制下求子集贡献和超集贡献
inline void transform(LL *f, int n, int opt) {
if(opt == 0) { // 求超集贡献
for(int o = 4, k = 1; o <= n; o <<= 2, k <<= 2 )
for(int i = 0; i < n; i += o )
for(int j = 0; j < k; j ++ )
f[i + j + 3 * k] = (f[i + j] + f[i + j + k] + f[i + j + 2 * k]);
}
else {
for(int o = 4, k = 1; o <= n; o <<= 2, k <<= 2 )
for(int i = 0; i < n; i += o )
for(int j = 0; j < k; j ++ )
f[i + j] += f[i + j + 3 * k],
f[i + j + k] += f[i + j + 3 * k],
f[i + j + 2 * k] += f[i + j + 3 * k];
}
}
}
int main() {
val['P'] = 0; val['R'] = 1; val['S'] = 2;
to[0] = 2, to[1] = 0, to[2] = 1;
scanf("%d%d%d", &K, &n, &m);
mi[0] = 1; for(int i = 1; i <= K; i ++ ) mi[i] = mi[i - 1] * 3;
for(int i = 0; i < (1 << 2 * K); i ++ ) {
int j = i;
for(int k = 0; k < K; k ++ ) {
if(j % 4 != 3) sz[i] ++;
j /= 4;
}
}
for(int i = 1; i <= n; i ++ ) {
int x = 0;
scanf("%s", str + 1); // K 位
for(int j = 1; j <= K; j ++ ) x = x * 4 + to[val[str[j]]];
f[x] ++;
}
for(int i = 1; i <= m; i ++ ) {
int x = 0;
scanf("%s", str + 1);
for(int j = 1; j <= K; j ++ ) x = x * 4 + to[val[str[j]]];
g[x] ++;
}
FWT::transform(f, 1 << 2 * K, 0);
FWT::transform(g, 1 << 2 * K, 0);
for(int i = 0; i < (1 << 2 * K); i ++ ) f[i] = f[i] * g[i] * sign(sz[i]);
FWT::transform(f, 1 << 2 * K, 1);
for(int i = 0; i < mi[K]; i ++ ) printf("%lld\n", 1LL * n * m - f[turn(i)]);
return 0;
}
P7930 [COCI 2021/2022 #1] Set
题意:
定义有序多元组 \(a\) 的第 \(i\) 项为 \(a_i\)。
给定 \(n\) 个有序 \(m\) 元组 \(b_i\),要从这些 \(m\) 元组中选出 \(3\) 个,设这三个 \(m\) 元组的下标为 \(i, j, k\)。需要满足:
- \(i < j < k\)
- \(\forall 1 \leq z \leq m\),\(b_i, b_j, b_k\) 在第 \(z\) 位 相同 或者 两两不同。
\(1 \leq m \leq 12, 1 \leq n \leq 3^m\),\(b_i\) 互不相同,\(1 \leq b_{i, j} \leq 3\)。
分析:
考虑合法三元组的等价条件:发现需要它们在每一位上的和 \(\% \ 3\) 等于 \(0\)。
将每一位上的数字都减一,等价条件不变。我们将每个 \(m\) 元组看作一个 \(m\) 位的三进制数。
那么相当于是一个 \(3\) 进制下异或卷积的形式。写出式子:
其中 \(\oplus\) 符号表示 \(K\) 进制异或卷积,\(g_{a}\) 表示是否存在一个数为 \(a\) 的元组。
那么答案就是 \(\frac{f_{0} - n}{6}\)。 \(-n\) 是因为会有三个数选到同一个位置的情况,\(/6\) 是因为不考虑顺序。
那么只需要将 \(g\) \(FWT\) 一次,然后每一个位置变为原来的三次方,然后在 \(IFWT\) 回去就能得到 \(f\)。这个不用取模直接构造单位根矩阵即可。
复杂度 \(O(m \times3^m)\)。
CODE:
// 合法的三元组每一位上的和 mod 3 = 0
// K 进制 XOR 卷积
#include<bits/stdc++.h>
#define pi acos(-1)
using namespace std;
typedef long long LL;
const int N = 550000;
struct comp { // 定义复数类
double x, y;
friend comp operator + (comp a, comp b) {return (comp) {a.x + b.x, a.y + b.y};}
friend comp operator * (comp a, comp b) {return (comp) {a.x * b.x - a.y * b.y, a.x * b.y + b.x * a.y};}
};
comp f[N];
int n, m, mi[20];
char b[N];
inline void print(comp a) {
printf("%.3lf + %.3lf i\n", a.x, a.y);
}
namespace FWT {
inline void XOR(comp *f, int n, int op) {
for(int o = 3, k = 1; o <= n; o *= 3, k *= 3 ) {
for(int i = 0; i < n; i += o )
for(int j = 0; j < k; j ++ ) {
comp tmp[3] = {f[i + j], f[i + j + k], f[i + j + k * 2]};
for(int t = i + j, l = 0; l < 3; t += k, l ++ ) {
comp wk = (comp) {cos(2.0 * pi * l / 3.0), op * sin(2.0 * pi * l / 3.0)};
comp w = (comp) {1, 0};
f[t] = (comp) {0, 0};
for(int h = 0; h < 3; h ++, w = w * wk) {
f[t] = f[t] + (w * tmp[h]);
}
if(op == -1) f[t] = f[t] * (comp) {(double)1.0 / 3.0, 0};
}
}
}
}
}
int main() {
scanf("%d%d", &n, &m);
mi[0] = 1; for(int i = 1; i <= m; i ++ ) mi[i] = mi[i - 1] * 3;
for(int i = 1; i <= n; i ++ ) {
scanf("%s", b + 1);
int x = 0;
for(int j = 1; j <= m; j ++ ) x = x * 3 + (b[j] - '1');
f[x].x = 1.0;
}
FWT::XOR(f, mi[m], 1); // FWT
for(int i = 0; i < mi[m]; i ++ ) f[i] = f[i] * f[i] * f[i]; // 卷两次
FWT::XOR(f, mi[m], -1); // IFWT
LL res = (LL)(round)(f[0].x);
res -= n; res /= 6;
cout << res << endl;
return 0;
}
【UNR #2】黎明前的巧克力
非常有意思的一道题。
题意:
给你一个大小为 \(n\) 的可重数集 \(\{a_i\}\),可认为每个元素都有编号。你需要从中选出两个不交子集 \(A,B\),满足 \(A\) 中所有元素 \(a_i\) 的异或和等于 \(B\) 中所有 \(a_i\) 的异或和,\(A,B\) 至少有一个不为空。求方案数对 \(998244353\) 取模后的结果。
\(1 \leq n \leq 10^6\), \(1 \leq a_i \leq 10^6\)。
分析:
设全集为 \(U\)。
直接考虑两个集合的选法是比较难的。但是我们发现 \(U / (A \cup B)\) 中元素的异或和一定与 \(U\) 中所有数的异或和相等,这是因为 \(A,B\) 异或和相同,异或后为 \(0\)。
设原来所有数的异或和为 \(v\),那么我们可以先从 \(U\) 中选出一个异或和为 \(v\) 的子集 \(S\),然后将 \(2^{|U - S|}\) 表示剩下的元素任意分到 \(A\) 或 \(B\) 的方案数贡献到答案中。容易得知其正确性。
进一步的,我们给每个元素赋一个 \(\frac{1}{2}\) 的权值。那么只需要计算所有异或和为 \(v\) 的子集 \(S\) 内部元素权值乘积的和,最后乘上 \(2^n\) 就是答案了。
考虑一个值域上的 \(dp\):
设 \(dp_{i, j}\) 表示从小到大考虑到了值为 \(i\) 的数,当前所有异或和为 \(j\) 的集合 \(S\) 的元素权值乘积之和。那么有转移:
- \(dp_{i, j} \times \sum\limits_{k = 0}^{2k \leq cnt_{i + 1}} \binom{cnt_{i + 1}}{2k}(\frac{1}{2})^{2k} \to dp_{i + 1, j}\)
- \(dp_{i, j} \times\ \sum\limits_{k = 0}^{2k + 1 \leq cnt_{i + 1}} \binom{cnt_{i + 1}}{2k + 1}(\frac{1}{2})^{2k + 1} \to dp_{i + 1, j \oplus (i + 1)}\)
设 \(g(x) = \sum\limits_{k = 0}^{2k \leq x} \binom{x}{2k}(\frac{1}{2})^{2k}\),\(h(x) = \sum\limits_{k = 0}^{2k + 1 \leq x} \binom{x}{2k + 1}(\frac{1}{2})^{2k + 1}\)。
那么转移可以写成:
- \(dp_{i, j} \times g(cnt_{i + 1}) \to dp_{i + 1, j}\)
- \(dp_{i, j} \times h(cnt_{i + 1}) \to dp_{i + 1, j \oplus (i + 1)}\)
这个其实就是 位运算卷积 的形式了:
看作有一个 \(0\) 位置上是 \(g(cnt_{i + 1})\), \(i + 1\) 位置上是 \(h(cnt_{i + 1})\),其余位置都是 \(0\) 的数组 \(F_{i + 1}\)。
转移可以写作: \(dp_{i} * F_{i + 1} \to dp_{i + 1}\)。其中 \(*\) 表示位运算卷积。
但是显然不可能每次做一遍 \(FWT\),这样复杂度是 \(O(n \times V)\) 的,\(V\) 代表最小的大于 \(max(a_i)\) 的 \(2\) 的整次幂。
我们来观察\(FWT(F_{i})\):
根据 异或卷积 的矩阵,可以知道 \(FWT(F_i)\) 的每个位置上都是 \(g(cnt_i) + h(cnt_i)\) 或 \(g(cnt_i) - h(cnt_i)\) 的形式。
然后惊奇的发现 \(g(cnt_i) + h(cnt_i) = (\frac{3}{2})^{cnt_i}\), \(g(cnt_i) - h(cnt_i) = (\frac{1}{2})^{cnt_i}\)。
具体的来说,对于一个 \(i\) 而言,如果 \(|i \& j|\) 为奇数,那么 \(FWT(F_i)_j = (\frac{1}{2})^{cnt_i}\),为偶数则是 \((\frac{3}{2})^{cnt_i}\)。
那么问题就很简单了:只要我们能快速求出 \(dp\) 到最后的 \(FWT\) 数组,然后 \(IFWT\) 一遍还原回去就行了。至于求末尾状态的 \(FWT\) 数组,根据刚才的结论,只要求出每个位置上 \(\frac{1}{2}\) 的幂次 \(y\) 和 \(\frac{3}{2}\) 的幂次 \(x\) 即可。可以先将初始的 \(cnt\) 数组 \(FWT\) 一遍,求出来每个位置上 \(x - y\) 的值,然后 \(x + y\) 一定等于 \(n\),因此可以直接解的 \(x, y\)。
时间复杂度 \(O(V \log_2 V)\)。
CODE:
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 1e6 + 10;
const int mod = 998244353;
const int inv_2 = 499122177;
inline int add(int x, int y) {return x + y >= mod ? x + y - mod : x + y;}
inline int del(int x, int y) {return x - y < 0 ? x - y + mod : x - y;}
inline int mul(int x, int y) {return 1LL * x * y % mod;}
int n, a[N], f[1 << 20], V, mi_2[N], mi_3_2[N], mi_1_2[N];
namespace FWT {
inline void XOR(int *f, int n, int op) {
for(int o = 2, k = 1; o <= n; o <<= 1, k <<= 1 )
for(int i = 0; i < n; i += o )
for(int j = 0; j < k; j ++ )
f[i + j] = add(f[i + j], f[i + j + k]),
f[i + j + k] = del(f[i + j], add(f[i + j + k], f[i + j + k])),
f[i + j] = mul(f[i + j], op), f[i + j + k] = mul(f[i + j + k], op);
}
}
int main() {
scanf("%d", &n);
mi_2[0] = mi_3_2[0] = mi_1_2[0] = 1;
for(int i = 1; i <= n; i ++ ) {
mi_2[i] = mul(mi_2[i - 1], 2);
mi_3_2[i] = mul(mi_3_2[i - 1], add(inv_2, 1));
mi_1_2[i] = mul(mi_1_2[i - 1], inv_2);
}
for(int i = 1; i <= n; i ++ ) {
scanf("%d", &a[i]);
f[a[i]] ++;
V ^= a[i];
}
FWT::XOR(f, 1 << 20, 1);
for(int i = 0; i < (1 << 20); i ++ ) { // 虽然是在指数上 FWT, 但是计算出的 x, y 仍然是真实值,不用变换模数
int x = add(f[i], n) / 2, y = n - x;
f[i] = mul(mi_3_2[x], mi_1_2[y]);
}
FWT::XOR(f, 1 << 20, inv_2);
int res = del(mul(f[V], mi_2[n]), 1);
cout << res << endl;
return 0;
}
【集训队互测2015】胡策的统计
题意:
对于一张无向图,定义它的联通值为连通块数的阶乘。
给你一张 \(n\) 个点,\(m\) 条边的简单无向图 \(G\),你需要求出它所有生成子图的联通值之和。
生成子图:原图删掉若干条边(可以是 \(0\) 条)后得到的图。
\(1 \leq n \leq 20,1\leq m \leq \frac{n \times (n - 1)}{2}\)。
分析:
设 \(g_S\) 表示只考虑 \(S\) 中的点形成的诱导子图,这些点连通的方案数。 \(E_S\) 表示 \(s\) 的诱导子图中的边数。那么可以容斥求 \(g_S\):
其中 \(mn(S)\) 表示 \(S\) 中编号最小的点,这里相当于枚举了编号最小的点所在的连通块,所以不重不漏。
令 \(T_S = 2^{E_S}\)。那么转移就变成了:
如果没有 \(mn(S') = mn(S)\) 的条件,那么 \(g_S\) 的转移就是一个常见的 子集卷积。
实际上有这个条件只需要在 \(FWT\) 的过程特判 \(lowbit(i) = lowbit(j)\) 才让 \(i\) 往 \(j\) 转移即可。
那么求 \(g\) 可以 半在线子集卷积 做到 \(O(n^2 \times 2^n)\)。
有了 \(g_S\) 考虑怎么求答案:
设 \(f_{S}\) 表示 \(S\) 中的点形成的诱导子图的所有生成子图的联通值之和。转移显然是枚举编号最小的点所在连通块。
但是我们需要乘上排列的系数,因此可能还需要加入一维。这样时间复杂度是 \(O(n3^n)\),即使使用子集卷积也是 \(O(n^32^n)\) 的。
但是我们注意到 排列 就是 不钦定枚举的连通块包含编号最小的点,一种方案会算重复的次数!!
因此直接平凡的子集卷积就是对的:
半在线子集卷积可做到 \(O(n^2 \times 2^n)\)。总复杂度 \(O(n^2 \times 2^n)\),需要卡常。
CODE:
// 两次子集卷积。 能过吗?
// 话说求 f 的时候去掉钦定直接就把系数给乘上了也很妙啊
#include<bits/stdc++.h>
using namespace std;
const int N = 21;
const int mod = 998244353;
inline int add(int x, int y) {return x + y >= mod ? x + y - mod : x + y;}
inline int del(int x, int y) {return x - y < 0 ? x - y + mod : x - y;}
inline int mul(int x, int y) {return 1LL * x * y % mod;}
int n, m, mi[N * N], eg[N][N];
int E[1 << 20], G[1 << 20], g[21][1 << 20], f[21][1 << 20]; // G[s] 表示 s 内部的点联通的方案数, 这个可以容斥
inline int lowbit(int x) {return (x & -x);}
inline int sz(int x) {return __builtin_popcount(x);}
namespace FWT {
inline void OR0(int *f, int n, int op) { // o = 1 代表有限制, o = 0 代表没限制
for(int o = 2, k = 1; o <= n; o <<= 1, k <<= 1 )
for(int i = 0; i < n; i += o )
for(int j = 0; j < k; j ++ )
if(op == 1) f[i + j + k] = add(f[i + j + k], f[i + j]);
else f[i + j + k] = del(f[i + j + k], f[i + j]);
}
inline void OR1(int *f, int n, int op) { // o = 1 代表有限制, o = 0 代表没限制
for(int o = 2, k = 1; o <= n; o <<= 1, k <<= 1 )
for(int i = 0; i < n; i += o )
for(int j = 0; j < k; j ++ )
if(lowbit(i + j) == lowbit(i + j + k)) {
if(op == 1) f[i + j + k] = add(f[i + j + k], f[i + j]);
else f[i + j + k] = del(f[i + j + k], f[i + j]);
}
}
}
int main() {
scanf("%d%d", &n, &m);
for(int i = 1; i <= m; i ++ ) {
int u, v; scanf("%d%d", &u, &v);
eg[u][v] = eg[v][u] = 1;
}
mi[0] = 1; for(int i = 1; i <= m; i ++ ) mi[i] = mul(mi[i - 1], 2);
for(int i = 0; i < (1 << n); i ++ ) {
if(!i) continue;
int p;
for(int j = 0; j < n; j ++ )
if((i >> j & 1)) {p = j; break;}
int s = (i ^ (1 << p)), ct = E[s];
for(int j = 0; j < n; j ++ )
if((s >> j & 1) && eg[p + 1][j + 1]) ct ++;
E[i] = ct;
}
for(int i = 0; i < (1 << n); i ++ ) f[sz(i)][i] = mi[E[i]];
for(int i = 1; i <= n; i ++ ) FWT::OR0(f[i], 1 << n, 1); // 2^n * n^2
for(int i = 0; i <= n; i ++ ) { // 按照 sz 大小依次确定 g
if(i > 1) {
for(int s = 0; s < (1 << n); s ++ ) // 2^n * n^2
for(int j = 1; j <= min(i - 1, sz(s)); j ++ )
g[i][s] = add(g[i][s], mul(g[j][s], f[i - j][s]));
}
if(i > 1) FWT::OR1(g[i], 1 << n, -1); // 2^n * n^2
for(int j = 0; j < (1 << n); j ++ ) {
if(sz(j) == i) G[j] = g[i][j] = del(mi[E[j]], g[i][j]);
else g[i][j] = 0;
}
if(i > 0) FWT::OR1(g[i], 1 << n, 1); // 正着 FWT 一遍 2^n * n^2
}
for(int i = 0; i <= n; i ++ ) {
memset(g[i], 0, sizeof g[i]);
memset(f[i], 0, sizeof f[i]);
}
for(int i = 0; i < (1 << n); i ++ ) g[sz(i)][i] = G[i];
for(int i = 1; i <= n; i ++ ) FWT::OR0(g[i], 1 << n, 1); // 2^n * n^2
for(int i = 0; i <= n; i ++ ) { // 依次求 f
if(i == 0) f[0][0] = 1;
else if(i == 1) {for(int j = 0; j < n; j ++ ) f[1][1 << j] = 1;}
else {
for(int s = 0; s < (1 << n); s ++ ) // 2^n * n^2
for(int j = 1; j <= min(i, sz(s)); j ++ )
f[i][s] = add(f[i][s], mul(g[j][s], f[i - j][s]));
}
if(i > 1) FWT::OR0(f[i], 1 << n, -1); // 2^n * n^2
for(int j = 0; j < (1 << n); j ++ ) {
if(sz(j) == i) G[j] = f[i][j];
else f[i][j] = 0;
}
FWT::OR0(f[i], 1 << n, 1); // 2^n * n^2
}
printf("%lld\n", G[(1 << n) - 1]);
return 0;
}

浙公网安备 33010602011771号