快速莫比乌斯 / 沃尔什变换 (FMT / FWT) / 子集卷积
luogu P4717 【模板】快速莫比乌斯 / 沃尔什变换 (FMT / FWT)
给定长 \(2^n\) 的两个序列 \(a,b\)，设：
\[c_k=\sum_{i\otimes j=k}a_ib_j \]分别求出 \(\otimes=\operatorname{and},\operatorname{or},\operatorname{xor}\) 时的 \(c\)。
引入
以下记序列长度 \(m=2^n\)，\(n\) 为维数。
FMT/FWT 在 OI 中，实际上是在做 \(\mathcal O\left(n2^n\right)\)的位运算卷积。
但如果从集合的角度来理解，也是可以理解的。
FMT 常常用于或卷积、与卷积，而 FWT 用于异或卷积、同或卷积。
过程有点类似于 FFT/NTT，大概就是将原序列变换后直接乘，乘完了再做逆变换，从而避免 \(\mathcal O\left(m^2\right)\) 计算。
FMT
或卷积
考虑计算：
考虑 \(i\operatorname{or}k=k,j\operatorname{or}k=k\implies(i\operatorname{or}j)\operatorname{or}k=k\)。
于是对于序列 \(a\)，构造 FMT 变换：
于是有：
于是计算 \(\operatorname{FMT}(a)_k\cdot\operatorname{FMT}(b)_k=\operatorname{FMT}(c)_k\) 是 \(\mathcal O(2^n)\) 的，再做一遍逆变换即可，即高维差分。
计算 \(\operatorname{FMT}(a)_i\) 可以用高维前缀和做到 \(\mathcal O\left(n2^n\right)\)。
也可以考虑集合表示，\(i\subseteq k,j\subseteq k\implies i\cup j\implies k\)。
于是可以理解这就是求子集的和，做高维前缀和即可。对于逆变换，就是高维差分。
有趣的事情是，可以发现这里印证了学长的那句「莫比乌斯变换本质是高维差分」。
因为每次只对维度 \(i\) 考虑，所以从小到大遍历和从大到小不重要。
代码实现上，\(\textit{op}=1\) 为正变换，\(\textit{op}=-1\) 为逆变换。
void FMTOR(int a[],int n,int op){
for(int i=0;i<n;i++){
for(int j=0;j<(1<<n);j++){
if(j&(1<<i)){
a[j]=(a[j]+op*a[j^(1<<i)])%P;
}
}
}
}
与卷积
考虑计算：
类似的考虑 \(i\operatorname{and} k=k,j\operatorname{and}k=k\implies (i\operatorname{and}j)\operatorname{and}k=k\)，构造：
证明是几乎一模一样：
同样可以写集合：
可以发现这就是求超集的和，也可以类似高维前缀和（其实是高维后缀和）。
void FMTAND(int a[],int n,int op){
for(int i=0;i<n;i++){
for(int j=0;j<(1<<n);j++){
if(!(j&(1<<i))){
a[j]=(a[j]+op*a[j^(1<<i)])%P;
}
}
}
}
FWT
异或卷积
考虑求：
对于序列 \(a\)，定义：
于是有：
考虑到：
不妨假定 \(i,j,k\in\set{0,1}\)，把每一位分开考虑。
- \(k=0\) 时，原式等价于 \(0+0\equiv0\pmod2\) 显然成立。
- \(k=1\) 时，原式等价于 \(i+j\equiv((i+j)\bmod2)\pmod2\)，因为 \(\operatorname{xor}\) 也是不进位加法。这显然成立。
于是有：
考虑集合语言描述：
于是这么构造是正确的，考虑如何实现 FWT 变换。
扫每一维 \(i\)，设 \(j\operatorname{and}2^i=0,j\operatorname{or}2^i\)。
那么对于正变换有：
考虑 \(j\) 的 \(i\) 位为 \(0\)，答案在 \(\operatorname{popcount}(k\operatorname{and} j)=0\) 的时候贡献为 \(1\)。
否则指数增加了 \(1\)，贡献为负数。
对于逆变换有：
逆变换的时候要除以 \(2\) 剔除重复贡献。
void FWTXOR(int a[],int n,int op){
for(int i=0;i<n;i++){
for(int j=0;j<(1<<n);j++){
if(!(j&(1<<i))){
int x=a[j],y=a[j^(1<<i)];
if(op==1){
a[j]=(x+y)%P;
a[j^(1<<i)]=(x-y)%P;
}else{
a[j]=1ll*(x+y)*inv2%P;
a[j^(1<<i)]=1ll*(x-y)*inv2%P;
}
}
}
}
}
同或卷积
和异或卷积类似。把变换的时候 \(j,j\operatorname{or}2^i\) 交换一下即可。
void FWTXNOR(int a[],int n,int op){
for(int i=0;i<n;i++){
for(int j=0;j<(1<<n);j++){
if(j&(1<<i)){
int x=a[j],y=a[j^(1<<i)];
if(op==1){
a[j]=(x+y)%P;
a[j^(1<<i)]=(x-y)%P;
}else{
a[j]=1ll*(x+y)*inv2%P;
a[j^(1<<i)]=1ll*(x-y)*inv2%P;
}
}
}
}
}
子集卷积
给定长 \(2^n\) 的序列 \(a,b\)，求：
\[c_k=\sum_{i\operatorname{and}j=0\\i\operatorname{or}j=k}a_ib_j \]
将其视为集合：
设或卷积运算 \(*\)，有 \(c=a*b\)。
子集卷积并没有像 FMT/FWT 那样的 \(\mathcal O\left(n2^n\right)\) 的简单方法。
\(i\cup j=k\) 的限制可以用 FMT 或卷积维护，但是处理不了 \(i\cap j=\varnothing\) 的限制。
考虑运算，若 \(i\cap j=\varnothing,i\cup j=k\)，可以得到 \(\vert i\vert+\vert j\vert=\vert k\vert\)。并且 \(\vert i\vert+\vert j\vert=\vert k\vert\) 当且仅当 \(i\cap j=\varnothing\)。
考虑强行升维，设 \(a_{i,s}\) 表示大小为 \(i\) 的集合 \(s\) 的值，\(b_{i,s},c_{i,s}\) 同理，答案即 \(c_{\vert s\vert,s}\)。
因此有：
具体而言，对 \(a_{i,s},b_{i,s}\) 进行 FMT 变换后，有：
于是可以 \(\mathcal O\left(n^22^n\right)\) 计算。
参考代码
//#include<bits/stdc++.h>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<iomanip>
#include<cstdio>
#include<string>
#include<vector>
#include<cmath>
#include<ctime>
#include<deque>
#include<queue>
#include<stack>
#include<list>
using namespace std;
constexpr const int N=20,P=1e9+9;
int n,a[N+1][1<<N|1],b[N+1][1<<N|1],c[N+1][1<<N|1];
void FMTOR(int a[],int n,int op){
for(int i=0;i<n;i++){
for(int j=0;j<(1<<n);j++){
if(j&(1<<i)){
a[j]=(a[j]+op*a[j^(1<<i)])%P;
}
}
}
}
int main(){
/*freopen("test.in","r",stdin);
freopen("test.out","w",stdout);*/
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
cin>>n;
for(int i=0;i<(1<<n);i++){
cin>>a[__builtin_popcount(i)][i];
}
for(int i=0;i<(1<<n);i++){
cin>>b[__builtin_popcount(i)][i];
}
for(int i=0;i<=n;i++){
FMTOR(a[i],n,1);
FMTOR(b[i],n,1);
}
for(int i=0;i<=n;i++){
for(int k=0;k<=i;k++){
for(int j=0;j<(1<<n);j++){
c[i][j]=(c[i][j]+1ll*a[k][j]*b[i-k][j])%P;
}
}
}
for(int i=0;i<=n;i++){
FMTOR(c[i],n,-1);
}
for(int i=0;i<(1<<n);i++){
cout<<(c[__builtin_popcount(i)][i]+P)%P<<' ';
}
cout.flush();
/*fclose(stdin);
fclose(stdout);*/
return 0;
}
浙公网安备 33010602011771号