用于解决一类位运算快速卷积问题:
\(C_{k}=\sum_{i\oplus j}A_i\times B_j\),其中 \(\oplus\) 为某种位运算,可以取按位或、按位与或按位异或。
P4717 【模板】快速莫比乌斯/沃尔什变换 (FMT/FWT)
概论
模仿(笔记)多项式基础 FFT 快速傅里叶变换 NTT 快速数论变换,我们尝试把需要操作的值都放在一个多项式的系数上,然后用一些快速的方法把它们变成点值的形式(其中我们会需要每个 \(x^k\),这是我们需要构造的),\(\Theta(n)\) 相乘,然后再用逆变换把这个东西变回系数的形式。
FWT
根据概论,我们不妨令点值 \(\text{FWT}(A)[i]=\sum_{j=0}^{n-1} A_j\times c(i,j)\),其中 \(c(i,j)\) 是 \(A_j\) 对 \(\text{FWT}(A)[i]\) 的贡献系数,类比 FFT 可以理解为幂级数中的 \(x^i\)。这个东西和 \(x^i\) 同时又有本质不同,我们不需要依赖 \(n\) 个本质不同的单位根分别带入求到点值再求逆变换,而是在位运算过程中可以分治处理。
根据点值性质,\(\text{FWT}(A)[i]\cdot\text{FWT}(B)[i]=\text{FWT}(C)[i]\),拆开可以得到:
根据 \(C_{k}=\sum_{i\oplus j}A_i\times B_j\):
性质 1,这是一条关键性质。
同时,由于位运算每位的独立性,\(c\) 还具有性质 2: \(c(i,j)=c(i_0,j_0)c(i_1,j_1)c(i_2,j_2)\dots\)。为什么是累乘是因为我们选用了乘法作为卷积的内置运算符号,这和上面那条性质是相符的。
接下来我们考虑求出一些合法的 \(c\) 后如何快速计算 \(\text{FWT}\)。
暴力是 \(O(n^2)\) 的,我们预期达到 \(O(n\log n)\),那直接从低到高按位考虑合并答案。由于位运算每一位的独立性,我们不需要像 FFT 那样奇偶拆半(单位根加速的需要),而是可以直接利用最高位为 \(0/1\) 左右拆半。这样考虑的好处是去掉了蝴蝶变换这一优化,且可以利用二进制位直接处理。
由于性质 2,让 \(i',j'\) 分别表示 \(i,j\) 去掉最高位后的数,我们可以如此按位拆开每个分治过程:
于是每次计算出 \(\sum_{j=0}^{n/2-1}c(i',j')A_j\) 这样的子问题即可。
转移矩阵构造
根据性质 2,我们只需要构造出 \(c([0,1],[0,1])\) 就相当于构造出了所有 \(c\)。我们还要解决的是 FWT 如何求逆?有如下公式:
类似莫比乌斯反演和二项式反演,我没法给出很严谨的推导过程,但是不妨代入下式到上式理解一下。此外,一个矩阵 \(A\) 的逆定义为 \(A\times B=B\times A=E\) 的矩阵 \(B\),这也代表了 \(A\) 不能存在某一行或某一列全为 \(0\),否则没有逆。
下文 \(\oplus\) 统一表示按位异或,\(\cup\) 表示按位或,\(\cap\) 表示按位与。
可以记一下结论,构造过程参考位运算卷积(FWT) & 集合幂级数。
或矩阵
我们构造矩阵需要满足 \(c(i,j)c(i,k)=c(i,j\cup k)\)。
比较自然地,我们发现这个东西实际上就是高位为 \(0\) 只取 \(0\) 的贡献,高位为 \(1\) 取 \(0\) 和 \(1\) 的贡献,也就是这个东西实际上是在做子集卷积((笔记)Sum over Subsets SOS 子集 DP 高维前缀和)。
主包没学过矩阵相关,只能从定义入手推一下逆:
大致可以列几个式子:
因此:
即逆矩阵:
感性理解构造其实可以拆开分治过程:
有点像前缀和和差分的关系。
与矩阵
构造:
逆:
异或矩阵
逆:
Code
点击查看代码
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL MOD=998244353,inv2=(MOD+1)/2;
const int N=17;
const LL cOR[2][2]={{1,0},{1,1}},cAND[2][2]={{1,1},{0,1}},cXOR[2][2]={{1,1},{1,MOD-1}},cnOR[2][2]={{1,0},{MOD-1,1}},cnAND[2][2]={{1,MOD-1},{0,1}},cnXOR[2][2]={{inv2,inv2},{inv2,MOD-inv2}};
int n;
LL A[1<<N],B[1<<N],FWTA[1<<N],FWTB[1<<N];
LL C[1<<N],FWTC[1<<N];
void FWT(LL *R,LL *FWTR,const LL c[2][2],int len){if(len==1){FWTR[0]=R[0];return ;}int half=len>>1;FWT(R,FWTR,c,half);FWT(R+half,FWTR+half,c,half);for(int i=0;i<half;i++){LL pre=FWTR[i],suf=FWTR[i+half];FWTR[i]=(c[0][0]*pre%MOD+c[0][1]*suf%MOD)%MOD;FWTR[i+half]=(c[1][0]*pre%MOD+c[1][1]*suf%MOD)%MOD;}
}
void solve(const LL c[2][2],const LL cn[2][2]){FWT(A,FWTA,c,n);FWT(B,FWTB,c,n);for(int i=0;i<n;i++)FWTC[i]=FWTA[i]*FWTB[i]%MOD;FWT(FWTC,C,cn,n);
}
int main(){ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);cin>>n;n=1<<n;for(int i=0;i<n;i++)cin>>A[i];for(int i=0;i<n;i++)cin>>B[i];solve(cOR,cnOR);for(int i=0;i<n;i++)cout<<C[i]<<' ';cout<<'\n';solve(cAND,cnAND);for(int i=0;i<n;i++)cout<<C[i]<<' ';cout<<'\n';solve(cXOR,cnXOR);for(int i=0;i<n;i++)cout<<C[i]<<' ';return 0;
}
Q&A
- 性质 2 成立的必要性?
它出现的必要性是我们需要保证 \(c(i_t,j_t)c(i_t,k_t)=c(i_t,j_t\oplus k_t)\iff c(i,j)c(i,k)=c(i,j\oplus k)\)。
- 异或卷积是否依赖模数?
是,否则将使用 double
。同时与或卷积可以不带取模,但是需要数据足够小。