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\)
- 如果當前位置 \(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)\)
遞迴邊界:若長度為 \(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))\) 即可用等比數列求和計算。