淺析快速沃爾什變換FWT(or,and,xor)
技術標籤:多項式全家桶
一通百通,先說最簡單的或運算吧
求 c i = ∑ j ∣ k = i a j b k c_i=\sum\limits_{j|k=i}a_jb_k ci=j∣k=i∑ajbk
先構造一個 F W T FWT FWT規則
也就是 F W T ( a ) = ∑ j ∣ i = i a j FWT(a)=\sum\limits_{j|i=i}a_j FWT(a)=j∣i=i∑aj
也就是
F
W
T
(
b
)
=
∑
k
∣
i
=
k
b
k
FWT(b)=\sum\limits_{k|i=k}b_k
FWT(b)=k∣i=k∑bk
那麼 F W T ( a ) ∗ F W T ( b ) FWT(a)*FWT(b) FWT(a)∗FWT(b)
= ∑ j ∣ i = i a j ∑ k ∣ i = i b k =\sum\limits_{j|i=i}a_j\sum\limits_{k|i=i}b_k =j∣i=i∑ajk∣i=i∑bk
= ∑ ( j ∣ k ) ∣ i a j b k =\sum\limits_{(j|k)|i}a_jb_k =(j∣k)∣i∑ajbk
考慮到陣列
c
c
c是
∑
j
∣
k
=
i
a
j
b
k
\sum\limits_{j|k=i}a_jb_k
j∣k=i∑
也就是 F W T ( a ) ∗ F W T ( b ) = F W T ( c ) FWT(a)*FWT(b)=FWT(c) FWT(a)∗FWT(b)=FWT(c)
那麼思路出來了,先把 a , b a,b a,b陣列做一遍 F W T FWT FWT得到 F W T ( c ) FWT(c) FWT(c)
再使用逆變換把 F W T ( c ) FWT(c) FWT(c)變回 c c c得到答案。
Ⅰ . a − − − − F W T ( a ) \color{Red}Ⅰ.a----FWT(a) Ⅰ.a−−−−FWT(a)
F
W
T
(
a
)
i
=
∑
j
∣
i
=
i
a
j
FWT(a)_i=\sum\limits_{j|i=i}a_j
現在假定陣列 a a a的長度是 2 n 2^n 2n
如果 O ( 2 n ∗ 2 n ) O(2^n*2^n) O(2n∗2n)列舉太慢了(這樣變換有個錘子用!!!)
考慮分治法,把陣列 a a a先按照下標二進位制最高位是 0 0 0還是 1 1 1分成兩個序列
最高位是 0 0 0的作為序列 a 0 a_0 a0,最高位是 1 1 1的作為序列 a 1 a_1 a1
那麼有 F W T ( a ) = [ F W T ( a 0 ) , F W T ( a 0 ) + F W T ( a 1 ) ] FWT(a)=[FWT(a_0),FWT(a_0)+FWT(a_1)] FWT(a)=[FWT(a0),FWT(a0)+FWT(a1)]
這裡務必清楚, F W T ( a ) FWT(a) FWT(a)是一個長 2 n 2^n 2n的序列,所以等式右邊的兩部分分別長 2 n − 1 2^{n-1} 2n−1
首先 F W T ( a ) FWT(a) FWT(a)的前 2 n − 1 2^{n-1} 2n−1項就是 F W T ( a 0 ) FWT(a_0) FWT(a0)
因為當下標 i i i小於 2 n − 1 2^{n-1} 2n−1時,和任何一項最高位有 1 1 1的或起來都不會等於 i i i了
通俗點說,由於 a 1 a_1 a1中的所有下標最高位都是 1 1 1,所以當 i i i的最高位沒有 1 1 1時,不滿足 j ∣ i = i j|i=i j∣i=i
然後 F W T ( a ) FWT(a) FWT(a)的後 2 n − 1 2^{n-1} 2n−1項就是 F W T ( a 0 ) + F W T ( a 1 ) FWT(a_0)+FWT(a_1) FWT(a0)+FWT(a1)
比如對於 F W T ( a 0 ) FWT(a_0) FWT(a0)來說,因為滿足 j ∣ i = i j|i=i j∣i=i
那麼也會滿足 ( j + 2 n − 1 ) ∣ ( i + 2 n − 1 ) = i + 2 n − 1 (j+2^{n-1})|(i+2^{n-1})=i+2^{n-1} (j+2n−1)∣(i+2n−1)=i+2n−1
所以後 2 n − 1 2^{n-1} 2n−1項係數涵蓋 F W T ( a 0 ) FWT(a_0) FWT(a0),也涵蓋了 F W T ( a 1 ) FWT(a_1) FWT(a1)
所以求 F W T ( a ) FWT(a) FWT(a),就先求出左半段 F W T ( a 0 ) FWT(a_0) FWT(a0),然後累加到右半段上去
void OR(int f[],int type)
{
for(int mid=2,k=1;mid<=mx;mid<<=1,k<<=1)
for(int i=0;i<mx;i+=mid )//長度每mid作為一個大區間
for(int j=0;j<k;j++)//[i,i+k)到[i+k,i+2k)分別作為左右區間,左邊加到右邊去
f[i+j+k] = ( f[i+j+k]+f[i+j]*type )%mod;
}
F W T ( a ) − − − − a \color{Red}FWT(a)----a FWT(a)−−−−a
反演回去怎麼搞 ? ? ?? ??
首先我們知道當陣列長度為 1 1 1時, a = F W T ( a ) a=FWT(a) a=FWT(a)
否則, U F W T ( a ) = [ U F W T ( a 0 ) , U F W T ( a 1 ) − U F W T ( a 0 ) ] UFWT(a)=[UFWT(a_0),UFWT(a_1)-UFWT(a_0)] UFWT(a)=[UFWT(a0),UFWT(a1)−UFWT(a0)]
感性理解,就是之前多加的部分,現在減回去
再提一下與運算 & \& &
求 c i = ∑ j & k = i a j b k c_i=\sum\limits_{j\&k=i}a_jb_k ci=j&k=i∑ajbk
構造 F W T ( a ) ∗ F W T ( b ) FWT(a)*FWT(b) FWT(a)∗FWT(b)
= ∑ j & i = i a j ∑ k & i = i b k =\sum\limits_{j\&i=i}a_j\sum\limits_{k\&i=i}b_k =j&i=i∑ajk&i=i∑bk
= ∑ ( j & k ) & i = i a j b k =\sum\limits_{(j\&k)\&i=i}a_jb_k =(j&k)&i=i∑ajbk
= F W T ( c ) =FWT(c) =FWT(c)
你會發現沒什麼區別,因為 ∣ | ∣和 & \& &都滿足結合律
然後快速求 a − − − > F W T ( a ) a--->FWT(a) a−−−>FWT(a)也是按照最高位分成 a 0 a_0 a0和 a 1 a_1 a1
讓我們來 y y yy yy一下,當 i > = 2 n − 1 i>=2^{n-1} i>=2n−1時滿足(此時 i i i的最高位是 1 1 1)
j & i = i j\&i=i j&i=i說明 j j j的最高位必須是 1 1 1,所以後 2 n − 1 2^{n-1} 2n−1項是 F W T ( a 1 ) FWT(a_1) FWT(a1)
當 i < 2 n − 1 i<2^{n-1} i<2n−1最高位沒有 1 1 1了,滿足 j ∣ i = i j|i=i j∣i=i的 j j j最高位可以是 0 0 0也可以是 1 1 1
至於低位怎麼樣滿足條件,丟給 F W T ( a 0 ) FWT(a_0) FWT(a0)和 F W T ( a 1 ) FWT(a_1) FWT(a1)即可
F W T ( a ) = [ F W T ( a 0 ) + F W T ( a 1 ) , F W T ( a 1 ) ] FWT(a)=[FWT(a_0)+FWT(a_1),FWT(a_1)] FWT(a)=[FWT(a0)+FWT(a1),FWT(a1)]
所以程式也和或運算差不多
void AND(int f[],int type)
{
for(int mid=2,k=1;mid<=mx;mid<<=1,k<<=1)
for(int i=0;i<mx;i+=mid )
for(int j=0;j<k;j++)
f[i+j] = ( f[i+j]+f[i+j+k]*type )%mod;
}
今天就到這裡吧,異或下次再來寫。下面是三種運算的程式碼。
#include <bits/stdc++.h>
using namespace std;
#define int long long
const int N = (1<<17),mod = 998244353;
int n,mx,a[N],b[N],ta[N],tb[N];
void OR(int f[],int type)
{
for(int mid=2,k=1;mid<=mx;mid<<=1,k<<=1)
for(int i=0;i<mx;i+=mid )
for(int j=0;j<k;j++)
f[i+j+k] = ( f[i+j+k]+f[i+j]*type )%mod;
}
void AND(int f[],int type)
{
for(int mid=2,k=1;mid<=mx;mid<<=1,k<<=1)
for(int i=0;i<mx;i+=mid )
for(int j=0;j<k;j++)
f[i+j] = ( f[i+j]+f[i+j+k]*type )%mod;
}
void XOR(int f[],int type)
{
for(int mid=2,k=1;mid<=mx;mid<<=1,k<<=1)
for(int i=0;i<mx;i+=mid)
for(int j=0;j<k;j++)
{
int t = i+j, v = i+j+k;
f[t] = ( f[t]+f[v] )%mod;
f[v] = ( f[t]-f[v]-f[v] )%mod;
f[t]=f[t]*type%mod, f[v]=f[v]*type%mod;
}
}
void init(){ for(int i=0;i<mx;i++) ta[i]=a[i],tb[i]=b[i]; }
void print(int a[]){ for(int i=0;i<mx;i++) printf("%lld%c",(a[i]%mod+mod)%mod,(i==mx-1)?'\n':' '); }
void mul(int a[],int b[] ){ for(int i=0;i<mx;i++) a[i] = a[i]*b[i]%mod; }
signed main()
{
scanf("%lld",&n);
mx = (1<<n);
for(int i=0;i<mx;i++) scanf("%lld",&a[i] ),a[i]%=mod;
for(int i=0;i<mx;i++) scanf("%lld",&b[i] ),b[i]%=mod;
init(); OR(ta,1); OR(tb,1); mul( ta,tb ); OR(ta,-1); print(ta);
init(); AND(ta,1); AND(tb,1); mul( ta,tb ); AND(ta,-1); print(ta);
init(); XOR(ta,1); XOR(tb,1); mul( ta,tb ); XOR(ta,(mod+1)>>1); print(ta);
}