[模板] 快速沃尔什变换
快速沃尔什变换
快速沃尔什变换是求这样的式子:
对序列 \(A\), \(B\), 求序列 \(C\), 使得
\[C_{i}=\sum_{j \oplus k} A_{j} B_{k}
\]
其中 \(\oplus\) 是任意位运算.
算法
类似 FFT, 构造从序列映射到序列的函数 \(FWT(F)\) 和 \(\operatorname{UFWT}()\), 满足
\[\operatorname{FWT}(C)_i = \operatorname{FWT}(A)_i \cdot \operatorname{FWT}(B)_i \tag{1}
\]
\[\operatorname{UFWT}(\operatorname{FWT}(C)) = C \tag{2}
\]
这样, 我们就可以把序列通过 \(\operatorname{FWT}\) 转换, 把位运算转换为点值乘法, 再通过 \(\operatorname{FWT}\) 转换回来.
构造
或运算
举或运算为例.
考虑倍增地构造 \(FWT(F)\). 记长为 \(n\) 的序列 \(A\) 的低 \(\frac n2\) 位为 \(A_0\), 高 \(\frac n2\) 位为 \(A_1\).
假设已经得到了 \(FWT(A_0)\) 和 \(FWT(A_1)\), 我们构造
\[\begin{cases}
FWT(A)_0 = FWT(A_0) \\
FWT(A)_1 = FWT(A_0)+FWT(A_1) \\
\end{cases}
\]
可以证明这样构造出来的式子满足 \((1)\) 式.
然后我们可以根据 \((2)\) 解出 \(UFWT(F)\) (在式子两边套上 \(UFWT()\) 即可):
\[\begin{cases}
UFWT(A)_0 = UFWT(A_0) \\
UFWT(A)_1 = UFWT(A_1)-UFWT(A_0) \\
\end{cases}
\]
我们还可以给出 \(FWT()\) 函数的显式:
\[FWT(A)_i=\sum_{i=i \lor j} A_j
\]
与运算
类似的, 可以对与运算构造:
\[FWT(A)_i=\sum_{i=i \land j} A_j
\]
\[\begin{cases}
FWT(A)_0 = FWT(A_0)+FWT(A_1) \\
FWT(A)_1 = FWT(A_1) \\
\end{cases}
\]
\[\begin{cases}
UFWT(A)_0 = UFWT(A_0)-UFWT(A_1) \\
UFWT(A)_1 = UFWT(A_1) \\
\end{cases}
\]
异或运算
\[FWT(A)_i=\sum_{j} (-1)^{ popcount(i \land j)} A_j
\]
其中 $popcount(n)$ 是 $n$ 的二进制表示中 $1$ 的个数.
\[\begin{cases}
FWT(A)_0 = FWT(A_0)+FWT(A_1) \\
FWT(A)_1 = FWT(A_0)-FWT(A_1) \\
\end{cases}
\]
\[\begin{cases}
UFWT(A)_0 = \frac{UFWT(A_0)+UFWT(A_1)}2 \\
UFWT(A)_1 = \frac{UFWT(A_0)-UFWT(A_1)}2 \\
\end{cases}
\]
其他运算
容易发现其他运算可以通过 \(\lnot\) 转换成以上三种运算.
因此, 只要对 \(A\), \(B\), 或者 \(C\) 下标取反 (\(C_i' = C_{\lnot i}\)) 即可.
例如, 对于 \(nor\) 运算, 只要计算 \(or\) 运算, 然后对 \(C\) 下标取反即可.
事实上, 由于任何位运算都可以用 \(\lnot\) 和 \(\lor\) 表示出来, 仅仅用 \(\lor\) 的 FWT 就可以计算所有位运算卷积. (然而难写)
代码
#include<cstdio>
#include<iostream>
#include<cmath>
#include<cstring>
#include<algorithm>
#include<queue>
#include<map>
#include<set>
using namespace std;
#define rep(i,l,r) for(register int i=(l);i<=(r);++i)
#define repdo(i,l,r) for(register int i=(l);i>=(r);--i)
#define il inline
typedef long long ll;
typedef double db;
//---------------------------------------
const int nsz=3e5+50;
const int nmod=998244353,inv2=499122177;
int n,l,a[nsz],b[nsz],c[nsz];
int c1[nsz],c2[nsz];
void cp(int *a,int *b,int n){memcpy(a,b,n*sizeof(int));}
void fwtor(int *a,int n,int fl){
for(int i=1;i<n;i<<=1){
for(int j=0,p=i<<1;j<n;j+=p){
for(int k=0;k<i;++k){
int x=a[j+k],y=a[j+k+i];
a[j+k+i]=(y+fl*x)%nmod;
}
}
}
}
void fwtand(int *a,int n,int fl){
for(int i=1;i<n;i<<=1){
for(int j=0,p=i<<1;j<n;j+=p){
for(int k=0;k<i;++k){
int x=a[j+k],y=a[j+k+i];
a[j+k]=(x+fl*y)%nmod;
}
}
}
}
void fwtxor(int *a,int n,int fl){
int val=(fl==1?1:inv2);
for(int i=1;i<n;i<<=1){
for(int j=0,p=i<<1;j<n;j+=p){
for(int k=0;k<i;++k){
int x=a[j+k],y=a[j+k+i];
a[j+k]=(x+y)*(ll)val%nmod;
a[j+k+i]=(x-y)*(ll)val%nmod;
}
}
}
}
void mul(int *a,int *b,int *c,int n,void (*fwt)(int*,int,int)){
//align n to 2^k
cp(c1,a,n),cp(c2,b,n);
fwt(c1,n,1),fwt(c2,n,1);
rep(i,0,n-1)c1[i]=(ll)c1[i]*c2[i]%nmod;
fwt(c1,n,-1);
cp(c,c1,n);
}
int main(){
ios::sync_with_stdio(0),cin.tie(0);
cin>>l;
n=(1<<l);
rep(i,0,n-1)cin>>a[i];
rep(i,0,n-1)cin>>b[i];
mul(a,b,c,n,fwtor);
rep(i,0,n-1)cout<<(c[i]+nmod)%nmod<<' ';
cout<<'\n';
mul(a,b,c,n,fwtand);
rep(i,0,n-1)cout<<(c[i]+nmod)%nmod<<' ';
cout<<'\n';
mul(a,b,c,n,fwtxor);
rep(i,0,n-1)cout<<(c[i]+nmod)%nmod<<' ';
cout<<'\n';
return 0;
}
