用途
Fast Walsh-Hadamard Transform,即FWT,用来解决形如 ci=∑j⊕k=iajbk c i = ∑ j ⊕ k = i a j b k 一类的卷积,其中 ⊕ ⊕ 表示位运算( xor/or/and xor/or/and )。
过程
类比FFT,我们可以构造出一个变换使得 ai a i 和 bj b j 可以直接相乘,再变换回来求得答案。同样的,FWT也分DWT和IDWT两步。
DWT
设 DWT(a)i D W T ( a ) i 表示 a a 经过DWT之后的变换。那么我们可以得到
DWT(a)i=∑j=0n−1ajf(i,j)DWT(a)i∗DWT(b)i=DWT(c)i①②(61)(62)(61)DWT(a)i=∑j=0n−1ajf(i,j)①(62)DWT(a)i∗DWT(b)i=DWT(c)i②
其中 f(i,j) f ( i , j ) 表示变换的系数。把①式代入②中,可以发现 f(i,j)∗f(i,k)=f(i,j⊕k) f ( i , j ) ∗ f ( i , k ) = f ( i , j ⊕ k )
对于不同的位运算,变换是不同的。下面给出结论:
⊕=&:f(i,j)=[i&j=i] ⊕ = & : f ( i , j ) = [ i & j = i ]
⊕=|:f(i,j)=[i&j==j] ⊕ = | : f ( i , j ) = [ i & j == j ]
⊕=^:f(i,j)=(−1)cnt(i&j) ⊕ = ^ : f ( i , j ) = ( − 1 ) c n t ( i & j )
( cnt(i) c n t ( i ) )表示 i i 在二进制下1的个数。
YY一下很对啊,但是不知道它怎么来的啊。
因为位运算每一位的独立性,f(i,j)f(i,j)可以二进制拆分。那么我们就可以开始推式子了:
DWT(a)i=∑j=0n−1ajf(i,j)=∑j=0n2−1ajf(i,j)+∑j=n2n−1ajf(i,j)=f(i0,0)∑j=0n2−1ajf(i1−>l−1,j1−>l−1)+f(i0,1)∑j=n2n−1ajf(i1−>l−1,j1−>l−1)(63)(64)(65) (63) D W T ( a ) i = ∑ j = 0 n − 1 a j f ( i , j ) (64) = ∑ j = 0 n 2 − 1 a j f ( i , j ) + ∑ j = n 2 n − 1 a j f ( i , j ) (65) = f ( i 0 , 0 ) ∑ j = 0 n 2 − 1 a j f ( i 1 − > l − 1 , j 1 − > l − 1 ) + f ( i 0 , 1 ) ∑ j = n 2 n − 1 a j f ( i 1 − > l − 1 , j 1 − > l − 1 )
然后就变成
DWT(A)i=f(0,0)DWT(AL)i+f(0,1)DWT(AR)iDWT(A)i+n2=f(1,0)DWT(AL)i+f(1,1)DWT(AR)i D W T ( A ) i = f ( 0 , 0 ) D W T ( A L ) i + f ( 0 , 1 ) D W T ( A R ) i D W T ( A ) i + n 2 = f ( 1 , 0 ) D W T ( A L ) i + f ( 1 , 1 ) D W T ( A R ) i
递归做就好了。(实际上并不用递归,和FFT一样迭代即可)
IDWT
代回去之后发现解个二元一次方程就好了。
模板
洛谷P4717
// luogu-judger-enable-o2
#include<cctype>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define N 1<<17
#define F inline
using namespace std;
typedef long long LL;
const LL P=;
int n; LL a[][N],b[][N],c[][N],inv[];
F char readc(){
static char buf[],*l=buf,*r=buf;
if (l==r) r=(l=buf)+fread(buf,,,stdin);
return l==r?EOF:*l++;
}
F int _read(){
int x=; char ch=readc();
while (!isdigit(ch)) ch=readc();
while (isdigit(ch)) x=(x<<)+(x<<)+(ch^),ch=readc();
return x;
}
F void writec(LL x){ if (x>) writec(x/); putchar(x%+); }
F void _write(LL x){ writec(x),putchar(' '); }
F void FWT(LL *a,int kd,int f){
for (int k=;k<n;k<<=)
for (int i=;i<n;i+=(k<<))
for (int j=;j<k;j++){
LL x=a[i+j],y=a[i+j+k];
if (kd==) (a[i+j+k]+=f*x)%=P;
else if (kd==) (a[i+j]+=f*y)%=P;
else a[i+j]=(x+y)*inv[-f>>]%P,a[i+j+k]=(x-y)*inv[-f>>]%P;
}
}
int main(){
n=<<_read(),inv[]=,inv[]=P>>|;
for (int i=;i<n;i++)
a[][i]=a[][i]=a[][i]=_read();
for (int i=;i<n;i++)
b[][i]=b[][i]=b[][i]=_read();
for (int i=;i<;i++){
FWT(a[i],i,),FWT(b[i],i,);
for (int j=;j<n;j++)
c[i][j]=a[i][j]*b[i][j]%P;
FWT(c[i],i,-);
}
for (int i=;i<;i++){
for (int j=;j<n;j++)
_write((c[i][j]+P)%P);
puts("");
}
return ;
}