1. 程式人生 > 其它 >FWT & FMT & 子集卷積 & 高維字首和

FWT & FMT & 子集卷積 & 高維字首和

胡亂卷積

1. FWT & FMT

全稱 Fast Walsh Transformation。

它可以在 \(n\log n\) 的時間內解決形如 \(c_k=\sum_{i\oplus j=k}a_ib_j\) 的卷積。其中 \(\oplus\) 可以是與,或和異或運算子。

一些約定:下標從 \(0\) 開始。序列長度用 \(2^n\) 表示,不足的位數補 \(0\),方便處理。

1.1. 或卷積 —— FMT

1.1.1. 引入

\(b_i=\sum_{j|i=i}a_j\)。看起來一臉不可做的樣子。

考慮分治:將長度為 \(2^n\) 的序列 \(a\) 按照第 \(n\) 位為 \(0/1\)

劃分為兩個長 \(2^{n-1}\) 的序列 \(a0,a1\),對於它們分別計算 \(b0,b1\)。假設我們已經得到了 \(b0\)\(b1\),考慮合併它們的答案:

  • 如果當前位置 \(x\) 的第 \(n\) 位為 \(0\),那麼因為 \(0|0=0\)\(1|0\neq 0\),所以 \(b_x=b0_x\)
    • 首先,如果原來(即去掉第 \(n\) 位時) \(y|x\neq x\),那麼無論第 \(n\) 位是什麼,考慮了之後 \(y|x\) 也一定不等於 \(x\)
    • \(y\) 的第 \(n\) 位為 \(0\),那麼因為 \(x\) 的第 \(n\) 位為 \(0\),所以仍然有 \(y|x=x\)
    • \(y\) 的第 \(n\) 位為 \(1\),那麼因為 \(x\) 的第 \(n\) 位為 \(0\),所以 \(y|x\neq x\)
    • 好比 \((2+0)|(3+0)=(3+0)\),而 \((2+4)|(3+0)\neq (3+0)\)
  • 如果當前位置 \(x\) 的第 \(n\) 位為 \(1\),那麼因為 \(0|1=1\)\(1|1=1\),所以 \(b_x=b0_{x-2^{n-1}}+b1_{x-2^{n-1}}\)
    • 原因與上方類似。
    • 好比 \((2+0)|(3+4)=(3+4)\)\((2+4)|(3+4)=(3+4)\)

綜上,\(b=\mathrm{merge}(b0,b0+b1)\)

\(\mathrm{merge}\) 表示連線兩個序列。

遞迴邊界:若長度為 \(1\),則直接 \(b_i=a_i\)

實際上這就是 FWT 的正變換。實現程式碼如下:

inline void OR(int *a){
    for(int d=2,k=1;d<=n;d<<=1,k<<=1)
        for(int i=0;i<n;i+=o)
            for(int j=0;j<k;j++)
                f[i+j+k]+=f[i+j];
}

1.1.2. 逆變換

\(b=\mathrm{merge}(b0,b0+b1)\) 顯然有 \(a=\mathrm{merge}(a0,a1-a0)\)

將正變換和逆變換合併,程式碼如下:

inline void OR(int *a,int x){
    for(int d=2,k=1;d<=n;d<<=1,k<<=1)
        for(int i=0;i<n;i+=o)
            for(int j=0;j<k;j++)
                f[i+j+k]+=f[i+j]*x;
}

1.1.3. 正確性證明

原理:如果 \(i|k=k\)\(j|k=k\),那麼顯然 \((i|j)|k=k\)

因此 \(a’_kb’_k=\sum_{i|k=k}\sum_{j|k=k}a_ib_j=\sum_{(i|j)|k=k}a_ib_j=c'_k\)

1.2. 與卷積 —— FMT

與或卷積類似,\(b=\mathrm{merge}(b0+b1,b1)\)\(a=\mathrm{merge}(a0-a1,a1)\)

inline void AND(int *a,int x){
    for(int d=2,k=1;d<=n;d<<=1,k<<=1)
        for(int i=0;i<n;i+=o)
            for(int j=0;j<k;j++)
                f[i+j]+=f[i+j+k]*x;
}

1.3. 異或卷積 —— FWT

異或看得我一臉異或。這是人能想出來的麼。

\(i\odot j=\mathrm{popcount}(i\& j)\bmod 2\)\(b_i=\sum_{j\odot i=0}a_j-\sum_{j\odot i=1}a_j\),在分治計算的時候:

  • 如果當前位置 \(x\) 的第 \(n\) 位為 \(0\),那麼因為 \(0\odot 0=0\)\(1\odot 0=0\),所以 \(b_x=b0_x+b1_x\)
    • 因為 \(x\) 的第 \(n\) 位為 \(0\),而 \(\odot\) 的定義是內部是 \(i\& j\),所以如果去掉第 \(n\) 位時 \(y\odot x=k\),那麼考慮第 \(n\) 位時仍然有 \(y\odot x=k\)
  • 如果當前位置 \(x\) 的第 \(n\) 位為 \(1\),那麼因為 \(0\odot 1=0\)\(1\odot 1=1\),所以 \(b_x=b0_{x-2^{n-1}}-b1_{x-2^{n-1}}\)
    • \(y\) 的第 \(n\) 位為 \(1\),那麼 \(x\odot y\) 相對於原來的 \(x\odot y\) 發生了改變,由 \(0\) 變成 \(1\) 或由 \(1\) 變成 \(0\),所以根據 \(b\) 的定義,\(b1_x\) 也變成了原來的相反數。

如果

綜上,\(b=\mathrm{merge}(b0+b1,b0-b1)\)。同理,逆變換 \(a=\mathrm{merge}(\frac{a0+a1}{2},\frac{a0-a1}{2})\)

void XOR(int *a,int x){
    for(int d=2,k=1;d<=n;d<<=1,k<<=1)
        for(int i=0;i<n;i+=o)
            for(int j=0;j<k;j++)
                f[i+j]+=f[i+j+k],
                f[i+j+k]=f[i+j]-f[i+j+k]-f[i+j+k],
                f[i+j]*=x,f[i+j+k]*=x;
}

程式碼倒是挺短的。

1.4. 參考文章

1.5. 例題

I. P4717 【模板】快速莫比烏斯/沃爾什變換 (FMT/FWT)

板子題。

#include <bits/stdc++.h>
using namespace std;

#define ll long long
#define mcpy(x,y) memcpy(x,y,sizeof(x))

const int mod=998244353;
const int N=1<<17;

ll n,lim,a[N],b[N],c[N],x[N],y[N];
void OR(ll *a,ll x){
	for(int d=2,k=1;d<=lim;d<<=1,k<<=1)
		for(int i=0;i<lim;i+=d)
			for(int j=0;j<k;j++)
				a[i+j+k]=(a[i+j+k]+a[i+j]*x)%mod; 
}
void AND(ll *a,ll x){
	for(int d=2,k=1;d<=lim;d<<=1,k<<=1)
		for(int i=0;i<lim;i+=d)
			for(int j=0;j<k;j++)
				a[i+j]=(a[i+j]+a[i+j+k]*x)%mod;
}
void XOR(ll *a,ll x){
	for(int d=2,k=1;d<=lim;d<<=1,k<<=1)
		for(int i=0;i<lim;i+=d)
			for(int j=0;j<k;j++)
				a[i+j]=(a[i+j]+a[i+j+k])%mod,
				a[i+j+k]=(a[i+j]-a[i+j+k]-a[i+j+k]+mod+mod)%mod,
				a[i+j]=a[i+j]*x%mod,a[i+j+k]=a[i+j+k]*x%mod;
}
int main(){
	cin>>n,lim=1<<n;
	for(int i=0;i<lim;i++)cin>>a[i];
	for(int i=0;i<lim;i++)cin>>b[i];
	mcpy(x,a),mcpy(y,b),OR(x,1),OR(y,1);
	for(int i=0;i<lim;i++)c[i]=x[i]*y[i]%mod;
	OR(c,mod-1); for(int i=0;i<lim;i++)cout<<c[i]<<" "; cout<<endl;
	mcpy(x,a),mcpy(y,b),AND(x,1),AND(y,1);
	for(int i=0;i<lim;i++)c[i]=x[i]*y[i]%mod;
	AND(c,mod-1); for(int i=0;i<lim;i++)cout<<c[i]<<" "; cout<<endl;
	mcpy(x,a),mcpy(y,b),XOR(x,1),XOR(y,1);
	for(int i=0;i<lim;i++)c[i]=x[i]*y[i]%mod;
	XOR(c,mod+1>>1); for(int i=0;i<lim;i++)cout<<c[i]<<" "; cout<<endl;
	return 0;
}

II. CF1530F Bingo

非正解。

首先補集轉化,求沒有任何一行 / 列 / 對角線的事件全部發生的概率。不妨把 “沒有發生” 視為 \(1\),“發生” 視為 \(0\)​,將所有​列以及對角線的狀態壓縮成 \(n+2\) 位二進位制數,對於每一行,列舉它 \(2^n-1\) 個狀態(\(0\) 不行,因為至少有一個事件沒有發生)所形成的 \(2^{n+2}\) 位進位制數並求出發生的概率,不難發現最終答案就是每一行做完 OR 卷積之後 \(2^{n+2}-1\) 位上的值。

時間複雜度 \(\mathcal{O}(2^nn^2)\)\(4\) 倍常數,將 long long 換成 int 勉強能過。

const int N=21;

int n,c,a[N][N],f[1<<N+2],g[1<<N+2],q[1<<N],p[N][N],iv[N][N];
int add(int &x,int y){x=(x+y>=mod?x+y-mod:x+y);}
void FMT(int *a){
	for(register int k=1,d=2;d<=c;k<<=1,d<<=1)
		for(register int i=0;i<c;i+=d)
			for(register int j=0;j<k;j++)
				a[i+j+k]=(a[i+j+k]+a[i+j])%mod;
}
void REV(int *a){
	for(int k=1,d=2;d<=c;k<<=1,d<<=1)
		for(int i=0;i<c;i+=d)
			for(int j=0;j<k;j++)
				a[i+j+k]=(a[i+j+k]-a[i+j]+mod)%mod;
}
int main(){
	cin>>n,c=1<<n+2;
	for(int i=0;i<n;i++)for(int j=0;j<n;j++){
		cin>>a[i][j],a[i][j]=a[i][j]*inv(10000)%mod;
		iv[i][j]=inv(a[i][j]),p[i][j]=(mod+1-a[i][j])%mod;
	}
	for(int i=0;i<c;i++)g[i]=1;
	
	for(int i=0;i<n;i++){
		q[0]=1;
		mem(f,0);
		for(int j=0;j<n;j++)q[0]=q[0]*a[i][j]%mod;
		for(int j=1;j<1<<n;j++){
			for(int k=0;k<n;k++)
				if((j>>k)&1){
					q[j]=q[j-(1<<k)]*p[i][k]%mod*iv[i][k]%mod;
					break;
				}
			int nw=j;
			if((nw>>i)&1)nw+=1<<n;
			if((nw>>n-i-1)&1)nw+=1<<n+1;
			f[nw]=q[j];
		} FMT(f);
		for(int j=0;j<c;j++)g[j]=g[j]*f[j]%mod;
	}
	REV(g);
	cout<<(1-g[c-1]+mod)%mod<<endl;
	return 0;
}

III. ABC212H Nim Counting

不難發現答案為 \(\sum_{i=1}^n\mathrm{IFWT}(\mathrm{FWT}^i(A))\)​。變成 \(\mathrm{IFWT}(\sum_{i=1}^n\mathrm{FWT}^i(A))\) 即可用等比數列求和計算。