FWT 重学
不会忘记的 FWT。
去年好像学了 FWT，但是过了不久就忘的差不多了，甚至想用的时候只能去复制之前写的板子，于是决定重学。
由于广为流传的推法我记不住，cmd 博客中写的矩阵推法也感觉不是很详细（或许是因为我太菜了吧），所以决定自己写一篇学习笔记。
FWT 用于将 \(O(n^2)\) 的位运算卷积优化到 \(O(n\log n)\)。
核心思想是将序列 \(A\) 对应到序列 \(F(A)\)，使得 \(F(A)\cdot F(B) = F(A*B)\)，做到将卷积变为点乘。
由于这个变换是线性变换，设 \(F(A)_i=\sum C_{i,j}A_j\)。
\(C\) 只需满足 \(C_{i,j}C_{i,k}=C_{i,j\oplus k}\)，其中 \(\oplus\) 为你要做的卷积中的位运算。
证明
\(\because F(A)\cdot F(B)=F(A*B)\)
\(\therefore\sum_j C_{i,j}A_j \sum_k C_{i,k}B_k = \sum_t C_{i,t}\sum_{x\oplus y=t}A_xB_y\)
\(\therefore\sum_{j,k}C_{i,j}C_{i,k}A_jB_k=\sum_{x,y}C_{i,x\oplus y}A_xB_y\)
\(\therefore C_{i,j}C_{i,k}=C_{i,j\oplus k}\)
在满足这个条件的前提下，\(C\) 矩阵当然越简单越好。
由于位运算中每位相互独立，所以实际上只需要给出 \(C_{0,0},C_{0,1},C_{1,0},C_{1,1}\) 的值，对于其他位置考虑下标的二进制分解，将该位置的值定为二进制下每一位的 \(C\) 值之积。
这样做有什么好处呢？其实是为了快速计算 \(F(A)\)。
根据上面的定义，我们在做矩阵 \(C\) 乘向量 \(A\) 时，可以将矩阵竖着切成两半，向量横着切成两半，那么前一半和后一半在二进制下有且仅有最高位不同。
考虑分治下去做，合并时前一半乘上系数 \(C_{i_0,0}\)，后一半乘上系数 \(C_{i_0,1}\)，其中 \(i_0\) 为 \(i\) 的最高位。
由于递归比较慢，我们将其改成自底向上完成分治过程。时间复杂度 \(O(n\log n)\)。
实现
inline void FWT(int*a,const int C[2][2]){
for(int len=1;len<n;len<<=1)for(int p=0;p<n;p+=len<<1)fr(i,p,p+len-1){
ll l=a[i],r=a[i+len];
a[i]=(C[0][0]*l%P+C[0][1]*r%P)%P,a[i+len]=(C[1][0]*l%P+C[1][1]*r%P)%P;
}
}
对于给定的位运算 \(\oplus\)，可以根据 \(C\) 矩阵的性质手推出 \(C\)。
下面给出几个常用的 \(C\) 矩阵（每个括号中左边是原矩阵，右边是它的逆）：
-
或卷积：\(\left(\begin{bmatrix}1&0\\1&1\end{bmatrix},\begin{bmatrix}1&0\\-1&1\end{bmatrix}\right)\) 或 \(\left(\begin{bmatrix}1&1\\1&0\end{bmatrix},\begin{bmatrix}0&1\\1&-1\end{bmatrix}\right)\)
-
与卷积：\(\left(\begin{bmatrix}1&1\\0&1\end{bmatrix},\begin{bmatrix}1&-1\\0&1\end{bmatrix}\right)\) 或 \(\left(\begin{bmatrix}0&1\\1&1\end{bmatrix},\begin{bmatrix}-1&1\\1&0\end{bmatrix}\right)\)
-
异或卷积：\(\left(\begin{bmatrix}1&1\\1&-1\end{bmatrix},\begin{bmatrix}0.5&0.5\\0.5&-0.5\end{bmatrix}\right)\) 或 \(\left(\begin{bmatrix}1&-1\\1&1\end{bmatrix},\begin{bmatrix}0.5&0.5\\-0.5&0.5\end{bmatrix}\right)\)
模板题完整代码
#include<bits/stdc++.h>
#define fr(i,l,r) for(int i(l),_##i(r);i<=_##i;i++)
template<class T>inline T rd(T&a){
T x=0;bool f=1;char c=getchar_unlocked();
for(;c<48|c>57;c=getchar_unlocked())f&=c!=45;
for(;c>47&c<58;c=getchar_unlocked())x=10*x+c-48;
return a=f?x:-x;
}template<class T,class...V>inline void rd(T&x,V&...v){rd(x),rd(v...);}
using namespace std;using ll=long long;
const int N=1<<17|3,P=998244353,inv2=P+1>>1;
int n,a[N],b[N],A[N],B[N];
const int OR[][2][2]={{{1,0},{1,1}},{{1,0},{P-1,1}}},
AND[][2][2]={{{1,1},{0,1}},{{1,P-1},{0,1}}},
XOR[][2][2]={{{1,1},{1,P-1}},{{inv2,inv2},{inv2,P-inv2}}};
inline void FWT(int*a,const int C[2][2]){
for(int len=1;len<n;len<<=1)for(int p=0;p<n;p+=len<<1)fr(i,p,p+len-1){
ll l=a[i],r=a[i+len];
a[i]=(C[0][0]*l%P+C[0][1]*r%P)%P,a[i+len]=(C[1][0]*l%P+C[1][1]*r%P)%P;
}
}
inline void conv(int*a,int*b,const int mat[][2][2]){
FWT(a,mat[0]),FWT(b,mat[0]);
fr(i,0,n-1)a[i]=1ll*a[i]*b[i]%P;
FWT(a,mat[1]);
}
int main(){
n=1<<rd(n);
fr(i,0,n-1)rd(A[i]);fr(i,0,n-1)rd(B[i]);
memcpy(a,A,n<<2),memcpy(b,B,n<<2),conv(a,b,OR);
fr(i,0,n-1)cout<<a[i]<<' ';puts("");
memcpy(a,A,n<<2),memcpy(b,B,n<<2),conv(a,b,AND);
fr(i,0,n-1)cout<<a[i]<<' ';puts("");
memcpy(a,A,n<<2),memcpy(b,B,n<<2),conv(a,b,XOR);
fr(i,0,n-1)cout<<a[i]<<' ';puts("");
return 0;
}
参考资料：
浙公网安备 33010602011771号