1. 程式人生 > 其它 >淺析快速沃爾什變換FWT(or,and,xor)

淺析快速沃爾什變換FWT(or,and,xor)

技術標籤:多項式全家桶

洛谷FWT模板

一通百通,先說最簡單的或運算

c i = ∑ j ∣ k = i a j b k c_i=\sum\limits_{j|k=i}a_jb_k ci=jk=iajbk

先構造一個 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)=ji=iaj

也就是 F W T ( b ) = ∑ k ∣ i = k b k FWT(b)=\sum\limits_{k|i=k}b_k FWT(b)=ki=kbk

那麼 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 =ji=iajki=ibk

= ∑ ( j ∣ k ) ∣ i a j b k =\sum\limits_{(j|k)|i}a_jb_k =(jk)iajbk

考慮到陣列 c c c ∑ j ∣ k = i a j b k \sum\limits_{j|k=i}a_jb_k jk=i

ajbk

也就是 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) .aFWT(a)

F W T ( a ) i = ∑ j ∣ i = i a j FWT(a)_i=\sum\limits_{j|i=i}a_j

FWT(a)i=ji=iaj

現在假定陣列 a a a的長度是 2 n 2^n 2n

如果 O ( 2 n ∗ 2 n ) O(2^n*2^n) O(2n2n)列舉太慢了(這樣變換有個錘子用!!!)

考慮分治法,把陣列 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} 2n1

首先 F W T ( a ) FWT(a) FWT(a)的前 2 n − 1 2^{n-1} 2n1項就是 F W T ( a 0 ) FWT(a_0) FWT(a0)

因為當下標 i i i小於 2 n − 1 2^{n-1} 2n1時,和任何一項最高位有 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 ji=i

然後 F W T ( a ) FWT(a) FWT(a)的後 2 n − 1 2^{n-1} 2n1項就是 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 ji=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+2n1)(i+2n1)=i+2n1

所以後 2 n − 1 2^{n-1} 2n1項係數涵蓋 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=iajbk

構造 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=iajk&i=ibk

= ∑ ( j & k ) & i = i a j b k =\sum\limits_{(j\&k)\&i=i}a_jb_k =(j&k)&i=iajbk

= 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>=2n1時滿足(此時 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} 2n1項是 F W T ( a 1 ) FWT(a_1) FWT(a1)

i < 2 n − 1 i<2^{n-1} i<2n1最高位沒有 1 1 1了,滿足 j ∣ i = i j|i=i ji=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); 
}