前言
首先,我们来看看多项式乘法。
f ( x ) = a 0 + a 1 x + a 2 x 2 + ⋯ + a n x n f(x)=a_0+a_1x+a_2x^2+\dots+a_nx^n f(x)=a0+a1x+a2x2+⋯+anxn
g ( x ) = b 0 + b 1 x + b 2 x 2 + ⋯ + b n x n g(x)=b_0+b_1x+b_2x^2+\dots+b_nx^n g(x)=b0+b1x+b2x2+⋯+bnxn
而 h ( x ) = f ( x ) × g ( x ) h(x)=f(x)\times g(x) h(x)=f(x)×g(x),如果
h ( x ) = c 0 + c 1 x + c 2 x 2 + ⋯ + c n x n h(x)=c_0+c_1x+c_2x^2+\dots+c_nx^n h(x)=c0+c1x+c2x2+⋯+cnxn
那么显然 c k = ∑ i = 0 k a i × b k − i c_k=\sum\limits_{i=0}^ka_i\times b_{k-i} ck=i=0∑kai×bk−i
我们可以用 F F T FFT FFT(快速傅里叶变换)来 O ( n log n ) O(n\log n) O(nlogn)求出多项式乘法。
为了方便,我们可以用 A A A来表示多项式 f ( x ) f(x) f(x)。
A = ( a 0 , a 1 , a 2 , … , a n ) A=(a_0,a_1,a_2,\dots,a_n) A=(a0,a1,a2,…,an)
同样地,多项式 g ( x ) g(x) g(x)和 h ( x ) h(x) h(x)同样可以用 B B B, C C C来表示。
那么,多项式乘法可以表示为
C k = ∑ i + j = k A i × B j C_k=\sum\limits_{i+j=k}A_i\times B_j Ck=i+j=k∑Ai×Bj
那如果式子变成这样呢?
C k = ∑ i ⊕ j = k A i × B j C_k=\sum\limits_{i\oplus j=k}A_i\times B_j Ck=i⊕j=k∑Ai×Bj
其中 ⊕ \oplus ⊕可以为 & , ∣ , xor \&,|,\text{xor} &,∣,xor。
这个时候, F F T FFT FFT就不管用了,我们就需要使用 F W T FWT FWT(快速沃尔什变换)。
约定
A + B = ( A 0 + B 0 , A 1 + B 1 , A 2 + B 2 , … ) A+B=(A_0+B_0,A_1+B_1,A_2+B_2,\dots) A+B=(A0+B0,A1+B1,A2+B2,…)
A − B = ( A 0 + B 0 , A 1 + B 1 , A 2 + B 2 , … ) A-B=(A_0+B_0,A_1+B_1,A_2+B_2,\dots) A−B=(A0+B0,A1+B1,A2+B2,…)
A × B = ( A 0 × B 0 , A 1 × B 1 , A 2 × B 2 , … ) A\times B=(A_0\times B_0,A_1\times B_1,A_2\times B_2,\dots) A×B=(A0×B0,A1×B1,A2×B2,…)
A ⊕ B = ( ∑ i ⊕ j = 0 A i × B j , ∑ i ⊕ j = 1 A i × B j , ∑ i ⊕ j = 2 A i × B j , … ) A\oplus B=(\sum\limits_{i\oplus j=0}A_i\times B_j,\sum\limits_{i\oplus j=1}A_i\times B_j,\sum\limits_{i\oplus j=2}A_i\times B_j,\dots) A⊕B=(i⊕j=0∑Ai×Bj,i⊕j=1∑Ai×Bj,i⊕j=2∑Ai×Bj,…)
注意 ⊕ \oplus ⊕在此处并不一定表示异或,还可以是 & , ∣ \&,| &,∣。
快速沃尔什变换
对于求上面这类运算的卷积,我们可以用 F F T FFT FFT的思路,将两个多项式先正变换,再按位相乘,最后逆变换回来。这就是 F W T FWT FWT的思路。
首先,令 A A A, B B B的最高次项的次数为 m m m,我们取一个整数 n n n,使 2 n ≥ m 2^n\geq m 2n≥m且 n n n最小。那么让 A A A, B B B的长度为 2 n 2^n 2n。这个操作和 F F T FFT FFT中的是一样的。
对于每一个卷积,我们需要找到一个变换 T T T,使得
- T ( A ⊕ B ) = T ( A ) × T ( B ) T(A\oplus B)=T(A)\times T(B) T(A⊕B)=T(A)×T(B)
- I T ( T ( A ) × T ( B ) ) = A ⊕ B IT(T(A)\times T(B))=A\oplus B IT(T(A)×T(B))=A⊕B
其中 I T IT IT表示 T T T的逆变换。
注意 T ( A ) T(A) T(A)和 I T ( B ) IT(B) IT(B)也表示多项式。
或卷积
C k = ∑ i ∣ j = k A i × B j C_k=\sum\limits_{i|j=k}A_i\times B_j Ck=i∣j=k∑Ai×Bj,也可记作 C = A ∣ B C=A|B C=A∣B。
构造 T ( A ) i = ∑ j ⊆ i A j T(A)_i=\sum\limits_{j\subseteq i}A_j T(A)i=j⊆i∑Aj,其中我们将二进制位的关系看作集合关系,即 j ⊆ i j\subseteq i j⊆i表示 ( j ∣ i ) = i (j|i)=i (j∣i)=i,下同。
则有
T ( A ∣ B ) p = ∑ k ⊆ p ∑ i ∣ j = k A i × B j = ∑ i ⊆ p A i × ∑ j ⊆ p B j = T ( A ) p × T ( B ) p T(A|B)_p=\sum\limits_{k\subseteq p}\sum\limits_{i|j=k}A_i\times B_j=\sum\limits_{i\subseteq p}A_i\times \sum\limits_{j\subseteq p}B_j=T(A)_p\times T(B)_p T(A∣B)p=k⊆p∑i∣j=k∑Ai×Bj=i⊆p∑Ai×j⊆p∑Bj=T(A)p×T(B)p
其中 ( ( i ∣ j ) ⊆ p ) ⇔ ( i ⊆ p ) ∧ ( j ⊆ p ) ((i|j)\subseteq p)\Leftrightarrow(i\subseteq p)\wedge(j\subseteq p) ((i∣j)⊆p)⇔(i⊆p)∧(j⊆p)
设当前的多项式 A A A有 2 n 2^n 2n项,用 A 0 A_0 A0表示前 2 n − 1 2^{n-1} 2n−1项, A 1 A_1 A1表示后 2 n − 1 2^{n-1} 2n−1项。
那么有递推式
T ( A ) = ( T ( A 0 ) , T ( A 0 ) + T ( A 1 ) ) T(A)=(T(A_0),T(A_0)+T(A_1)) T(A)=(T(A0),T(A0)+T(A1))
中间的逗号表示将 T ( A 0 ) T(A_0) T(A0)和 T ( A 1 ) T(A_1) T(A1)拼接在一起。
为什么是这样呢?看上面的 T T T的式子 T ( A ) i = ∑ j ⊆ i A j T(A)_i=\sum\limits_{j\subseteq i}A_j T(A)i=j⊆i∑Aj,第 n n n位为 1 1 1的部分对第 n n n位为 0 0 0的部分没有影响,而第 n n n位为 0 0 0的部分可以作为第 n n n位为 1 1 1的部分的子集。
那么, T T T的逆变换 I T IT IT如下
I T ( A ′ ) = ( I T ( A 0 ′ ) , I T ( A 1 ′ ) − I T ( A 0 ′ ) ) IT(A')=(IT(A_0'),IT(A_1')-IT(A_0')) IT(A′)=(IT(A0′),IT(A1′)−IT(A0′))
其中 A , A 0 ′ , A 1 ′ A,A_0',A_1' A,A0′,A1′表示经过变换 T T T得到的多项式。
正变换和逆变换的实现可以合并在一起。
code
void fwt_or(long long *w,int fl){for(int s=2;s<=1<<n;s<<=1){int mid=s>>1;for(int v=0;v<1<<n;v+=s){for(int i=0;i<mid;i++){w[v+mid+i]=(w[v+mid+i]+fl*w[v+i]+mod)%mod;}}}
}
与卷积
C k = ∑ i & j = k A i × B j C_k=\sum\limits_{i\& j=k}A_i\times B_j Ck=i&j=k∑Ai×Bj,也可记作 C = A & B C=A\& B C=A&B。
构造 T ( A ) i = ∑ i ⊆ j A j T(A)_i=\sum\limits_{i\subseteq j}A_j T(A)i=i⊆j∑Aj。
则有
T ( A & B ) p = ∑ p ⊆ k ∑ i & j = k A i × B j = ∑ p ⊆ i A i × ∑ p ⊆ j B j = T ( A ) p × T ( B ) p T(A\& B)_p=\sum\limits_{p\subseteq k}\sum\limits_{i\& j=k}A_i\times B_j=\sum\limits_{p\subseteq i}A_i\times \sum\limits_{p\subseteq j}B_j=T(A)_p\times T(B)_p T(A&B)p=p⊆k∑i&j=k∑Ai×Bj=p⊆i∑Ai×p⊆j∑Bj=T(A)p×T(B)p
其中 ( p ⊆ ( i & j ) ) ⇔ ( p ⊆ i ) ∧ ( p ⊆ j ) (p\subseteq (i\& j))\Leftrightarrow(p\subseteq i)\wedge(p\subseteq j) (p⊆(i&j))⇔(p⊆i)∧(p⊆j)
设当前的多项式 A A A有 2 n 2^n 2n项,用 A 0 A_0 A0表示前 2 n − 1 2^{n-1} 2n−1项, A 1 A_1 A1表示后 2 n − 1 2^{n-1} 2n−1项。
那么有递推式
T ( A ) = ( T ( A 0 ) , , + T ( A 1 ) , T ( A 1 ) ) T(A)=(T(A_0),,+T(A_1),T(A_1)) T(A)=(T(A0),,+T(A1),T(A1))
为什么是这样呢?也看上面的 T T T的式子 T ( A ) i = ∑ i ⊆ j A j T(A)_i=\sum\limits_{i\subseteq j}A_j T(A)i=i⊆j∑Aj,第 n n n位为 0 0 0的部分对第 n n n位为 1 1 1的部分没有影响,而第 n n n位为 1 1 1的部分可以作为第 n n n位为 0 0 0的部分的超集。
那么, T T T的逆变换 I T IT IT如下
I T ( A ′ ) = ( I T ( A 0 ′ ) − I T ( A 1 ′ ) , I T ( A 1 ′ ) ) IT(A')=(IT(A_0')-IT(A_1'),IT(A_1')) IT(A′)=(IT(A0′)−IT(A1′),IT(A1′))
同样地,正变换和逆变换的实现可以合并在一起。
code
void fwt_and(long long *w,int fl){for(int s=2;s<=1<<n;s<<=1){int mid=s>>1;for(int v=0;v<1<<n;v+=s){for(int i=0;i<mid;i++){w[v+i]=(w[v+i]+fl*w[v+mid+i]+mod)%mod;}}}
}
异或卷积
C k = ∑ i ⊕ j = k A i × B j C_k=\sum\limits_{i\oplus j=k}A_i\times B_j Ck=i⊕j=k∑Ai×Bj,也可记作 C = A ⊕ B C=A\oplus B C=A⊕B。
在这里, ⊕ \oplus ⊕就表示异或了。
首先,两个二进制数的异或,不会改变这两个数中二进制为 1 1 1的位的数量和的奇偶性。
设 d i d_i di表示 i i i的二进制为 1 1 1的位的个数的奇偶性(值为 0 0 0或 1 1 1),则有
d i ⊕ j = d i ⊕ d j d_{i\oplus j}=d_i\oplus d_j di⊕j=di⊕dj
显然有
d i & k ⊕ d j & k = d ( i & k ) ⊕ ( j & k ) = d ( i ⊕ j ) & k d_{i\& k}\oplus d_{j\& k}=d_{(i\& k)\oplus (j\& k)}=d_{(i\oplus j)\& k} di&k⊕dj&k=d(i&k)⊕(j&k)=d(i⊕j)&k
于是,我们可以构造
T ( A ) i = ∑ ( − 1 ) j & i A j T(A)_i=\sum (-1)^{j\& i}A_j T(A)i=∑(−1)j&iAj
则有
T ( A ⊕ B ) p = ∑ ( − 1 ) k & p ∑ i ⊕ j = k A i × B j = ∑ ( − 1 ) d i & p ⊕ d j & p A i × B j T(A\oplus B)_p=\sum(-1)^{k\& p}\sum\limits_{i\oplus j=k}A_i\times B_j=\sum\limits(-1)^{d_{i\& p}\oplus d_{j\& p}}A_i\times B_j T(A⊕B)p=∑(−1)k&pi⊕j=k∑Ai×Bj=∑(−1)di&p⊕dj&pAi×Bj
= ∑ ( − 1 ) d i & p A i × ∑ ( − 1 ) d j & p B j = T ( A ) × T ( B ) =\sum (-1)^{d_{i\& p}}A_i\times \sum(-1)^{d_{j\& p}}B_j=T(A)\times T(B) =∑(−1)di&pAi×∑(−1)dj&pBj=T(A)×T(B)
设当前的多项式 A A A有 2 n 2^n 2n项,用 A 0 A_0 A0表示前 2 n − 1 2^{n-1} 2n−1项, A 1 A_1 A1表示后 2 n − 1 2^{n-1} 2n−1项。
那么有递推式
T ( A ) = ( T ( A 0 ) + T ( A 1 ) , T ( A 0 ) − T ( A 1 ) ) T(A)=(T(A_0)+T(A_1),T(A_0)-T(A_1)) T(A)=(T(A0)+T(A1),T(A0)−T(A1))
为什么呢?
对于 0 ≤ i < 2 n − 1 0\leq i < 2^{n-1} 0≤i<2n−1
T ( A ) i = ∑ ( − 1 ) d j & i A j = T ( A 0 ) i + ∑ j = 0 2 n − 1 − 1 ( − 1 ) d j & i ( A 1 ) j T(A)_i=\sum (-1)^{d_{j\& i}}A_j=T(A_0)_i+\sum\limits_{j=0}^{2^{n-1}-1}(-1)^{d_{j\& i}}(A_1)_j T(A)i=∑(−1)dj&iAj=T(A0)i+j=0∑2n−1−1(−1)dj&i(A1)j
而
T ( A 1 ) i = ∑ j = 0 2 n − 1 − 1 ( − 1 ) d j & i ( A 1 ) j T(A_1)_i=\sum\limits_{j=0}^{2^{n-1}-1} (-1)^{d_{j\& i}}(A_1)_j T(A1)i=j=0∑2n−1−1(−1)dj&i(A1)j
所以
T ( A ) = T ( A 0 ) + T ( A 1 ) T(A)=T(A_0)+T(A_1) T(A)=T(A0)+T(A1)
注意因为 0 ≤ i < 2 n − 1 0\leq i < 2^{n-1} 0≤i<2n−1,所以 ( j + 2 n − 1 ) & i = j & i (j+2^{n-1})\& i=j\& i (j+2n−1)&i=j&i
对于 2 n − 1 ≤ i < 2 n 2^{n-1}\leq i<2^n 2n−1≤i<2n
T ( A ) i = ∑ ( − 1 ) d j & i A j = ∑ j = 0 2 n − 1 − 1 ( − 1 ) d j & ( i − 2 n − 1 ) ( A 0 ) j + ∑ j = 0 2 n − 1 − 1 ( − 1 ) d j & ( i − 2 − 1 ) + 1 ( A 1 ) j T(A)_i=\sum (-1)^{d_{j\& i}}A_j=\sum\limits_{j=0}^{2^{n-1}-1}(-1)^{d_{j\& (i-2^{n-1})}}(A_0)_j+\sum\limits_{j=0}^{2^{n-1}-1}(-1)^{d_{j\&(i-2^{-1})}+1}(A_1)_j T(A)i=∑(−1)dj&iAj=j=0∑2n−1−1(−1)dj&(i−2n−1)(A0)j+j=0∑2n−1−1(−1)dj&(i−2−1)+1(A1)j
所以
T ( A ) = T ( A 0 ) + ( − 1 ) × T ( A 1 ) = T ( A 0 ) − T ( A 1 ) T(A)=T(A_0)+(-1)\times T(A_1)=T(A_0)-T(A_1) T(A)=T(A0)+(−1)×T(A1)=T(A0)−T(A1)
T T T的逆变换 I T IT IT如下
I T ( A ′ ) = ( I T ( A 0 ′ ) + I T ( A 1 ′ ) 2 , I T ( A 0 ′ ) − I T ( A 1 ′ ) 2 ) IT(A')=(\dfrac{IT(A_0')+IT(A_1')}{2},\dfrac{IT(A_0')-IT(A_1')}{2}) IT(A′)=(2IT(A0′)+IT(A1′),2IT(A0′)−IT(A1′))
将正变换和逆变换的实现合并在一起。
code
void fwt_xor(long long *w,int fl){for(int s=2;s<=1<<n;s<<=1){int mid=s>>1;for(int v=0;v<1<<n;v+=s){for(int i=0;i<mid;i++){long long t1=w[v+i],t2=w[v+mid+i];w[v+i]=(t1+t2)%mod;w[v+mid+i]=(t1-t2+mod)%mod;if(fl==-1){w[v+i]=w[v+i]*ny2%mod;w[v+mid+i]=w[v+mid+i]*ny2%mod;}}}}
}
例题
P4717 【模板】快速莫比乌斯/沃尔什变换 (FMT/FWT)
或卷积、与卷积、异或卷积的模板。
code
#include<bits/stdc++.h>
using namespace std;
int n;
long long a[1<<18],b[1<<18],c[1<<18],v1[1<<18],v2[1<<18];
long long mod=998244353,ny2=499122177;
void pt(){for(int i=0;i<1<<n;i++){a[i]=v1[i];b[i]=v2[i];}
}
void gt(){for(int i=0;i<1<<n;i++) c[i]=a[i]*b[i]%mod;
}
void out(){for(int i=0;i<1<<n;i++) printf("%lld ",c[i]);printf("\n");
}
void fwt1(long long *w,int fl){for(int s=2;s<=1<<n;s<<=1){int mid=s>>1;for(int v=0;v<1<<n;v+=s){for(int i=0;i<mid;i++){w[v+mid+i]=(w[v+mid+i]+fl*w[v+i]+mod)%mod;}}}
}
void fwt2(long long *w,int fl){for(int s=2;s<=1<<n;s<<=1){int mid=s>>1;for(int v=0;v<1<<n;v+=s){for(int i=0;i<mid;i++){w[v+i]=(w[v+i]+fl*w[v+mid+i]+mod)%mod;}}}
}
void fwt3(long long *w,int fl){for(int s=2;s<=1<<n;s<<=1){int mid=s>>1;for(int v=0;v<1<<n;v+=s){for(int i=0;i<mid;i++){long long t1=w[v+i],t2=w[v+mid+i];w[v+i]=(t1+t2)%mod;w[v+mid+i]=(t1-t2+mod)%mod;if(fl==-1){w[v+i]=w[v+i]*ny2%mod;w[v+mid+i]=w[v+mid+i]*ny2%mod;}}}}
}
int main()
{scanf("%d",&n);for(int i=0;i<1<<n;i++) scanf("%lld",&v1[i]);for(int i=0;i<1<<n;i++) scanf("%lld",&v2[i]);pt();fwt1(a,1);fwt1(b,1);gt();fwt1(c,-1);out();pt();fwt2(a,1);fwt2(b,1);gt();fwt2(c,-1);out();pt();fwt3(a,1);fwt3(b,1);gt();fwt3(c,-1);out();return 0;
}
子集卷积
C k = ∑ i ⊆ k A i × B k − i C_k=\sum\limits_{i\subseteq k}A_i\times B_{k-i} Ck=i⊆k∑Ai×Bk−i,也可以写成 C k = ∑ i ∣ j = k , i & j = 0 A i × B j C_k=\sum\limits_{i|j=k,i\& j=0}A_i\times B_j Ck=i∣j=k,i&j=0∑Ai×Bj.
相对于或卷积来说,它多了一个限制,即必须是两个二进制位没有交集的下标进行或运算。
可以按二进制位为 1 1 1的数量来分类,设 d i d_i di表示 i i i中二进制位为 1 1 1的位数。那么
C k = ∑ i ∣ j = k , d i + d j = d k A i × B j C_k=\sum\limits_{i|j=k,d_i+d_j=d_k}A_i\times B_j Ck=i∣j=k,di+dj=dk∑Ai×Bj
多项式 V i V_i Vi表示二进制位为 1 1 1的个数为 i i i的 A A A值,二进制位为 1 1 1的个数不等于 i i i的数在 V i V_i Vi中值为零。
对于每一个 V i V_i Vi,做一次 F W T FWT FWT。若 k k k中二进制位为 1 1 1的个数为 i i i,则 C k = ∑ j = 0 i ( V j ) k × ( V i − j ) k C_k=\sum\limits_{j=0}^i(V_j)_k\times (V_{i-j})_k Ck=j=0∑i(Vj)k×(Vi−j)k。
例题
P6097 【模板】子集卷积
子集卷积的模板。
code
#include<bits/stdc++.h>
using namespace std;
int n;
long long a[1<<21],b[1<<21],c[1<<21],ct[1<<21],ans[1<<21],ta[22][1<<21],tb[22][1<<21];
const long long mod=1e9+9;
void fwt(long long *w,long long fl){for(int s=2;s<=(1<<n);s<<=1){int mid=s>>1;for(int v=0;v<(1<<n);v+=s){for(int i=0;i<mid;i++){w[v+mid+i]=(w[v+mid+i]+fl*w[v+i]+mod)%mod;}}}
}
int main()
{scanf("%d",&n);for(int i=1;i<(1<<n);i++){ct[i]=ct[i-(i&(-i))]+1;}for(int i=0;i<(1<<n);i++){scanf("%lld",&a[i]);ta[ct[i]][i]=a[i];}for(int i=0;i<(1<<n);i++){scanf("%lld",&b[i]);tb[ct[i]][i]=b[i];}for(int i=0;i<=n;i++){fwt(ta[i],1);fwt(tb[i],1);}for(int i=0;i<=n;i++){for(int j=0;j<=i;j++){for(int k=0;k<(1<<n);k++){c[k]=(c[k]+ta[j][k]*tb[i-j][k]%mod)%mod;}}fwt(c,-1);for(int j=0;j<(1<<n);j++){if(ct[j]==i) ans[j]=c[j];c[j]=0;}}for(int i=0;i<(1<<n);i++){printf("%lld ",ans[i]);}return 0;
}
总结
F W T FWT FWT主要是要找到变换 T T T,然后用一种快速的方法变换,将两个多项式相乘后再用逆变换变回来,即可得出答案。
一般情况下,时间复杂度为 O ( n log n ) O(n\log n) O(nlogn)。