快速沃尔什变换
快速沃尔什变换是求这样的式子:
对序列 (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 ag{1}
]
[operatorname{UFWT}(operatorname{FWT}(C)) = C ag{2}
]
这样, 我们就可以把序列通过 (operatorname{FWT}) 转换, 把位运算转换为点值乘法, 再通过 (operatorname{FWT}) 转换回来.
构造
或运算
举或运算为例.
考虑倍增地构造 (FWT(F)). 记长为 (n) 的序列 (A) 的低 (frac n2) 位为 (A_0), 高 (frac n2) 位为 (A_1).
假设已经得到了 (FWT(A_0)) 和 (FWT(A_1)), 我们构造
[egin{cases}
FWT(A)_0 = FWT(A_0) \
FWT(A)_1 = FWT(A_0)+FWT(A_1) \
end{cases}
]
可以证明这样构造出来的式子满足 ((1)) 式.
然后我们可以根据 ((2)) 解出 (UFWT(F)) (在式子两边套上 (UFWT()) 即可):
[egin{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
]
[egin{cases}
FWT(A)_0 = FWT(A_0)+FWT(A_1) \
FWT(A)_1 = FWT(A_1) \
end{cases}
]
[egin{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$ 的个数.
[egin{cases}
FWT(A)_0 = FWT(A_0)+FWT(A_1) \
FWT(A)_1 = FWT(A_0)-FWT(A_1) \
end{cases}
]
[egin{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<<'
';
mul(a,b,c,n,fwtand);
rep(i,0,n-1)cout<<(c[i]+nmod)%nmod<<' ';
cout<<'
';
mul(a,b,c,n,fwtxor);
rep(i,0,n-1)cout<<(c[i]+nmod)%nmod<<' ';
cout<<'
';
return 0;
}