@演算法 - [email protected] 位運算卷積 —— FWT
目錄
- @0 - 參考資料@
- @1 - 異或卷積概念及性質@
- @2 - 快速沃爾什正變換（異或）@
- @3 - 快速沃爾什逆變換（異或）@
- @4 - 與卷積、或卷積@
- @5 - 參考程式碼實現@
- @6 - 例題與應用@
@0 - 參考資料@
@1 - 異或卷積概念及性質@
現在已知兩個 n 維的向量 \(A=(a_0,a_1,...,a_{n-1}),B=(b_0,b_1,...,b_{n-1})\)。定義異或卷積式：
\[A\oplus B=(\sum_{i\oplus j=0}a_i*b_j,\sum_{i\oplus j=1}a_i*b_j,\dots,\sum_{i\oplus j=n-1}a_i*b_j)\]
定義加法 \(A\pm B=(a_0\pm b_0,a_1\pm b_1,...,a_{n-1}\pm b_{n-1})\)。
定義乘積 \(A*B=(a_0*b_0,a_1*b_1,...,a_{n-1}*b_{n-1})\)。
定義\(A_0=(a_0,a_1,...a_{n/2-1}),A_1=(a_{n/2},a_{n/2+1},...,a_{n-1})\)，則有 \(A=(A_0,A_1)\)。可以把 \(A_0,A_1\) 理解為向量 \(A\) 的前半部分與後半部分，也可以理解為 \(A\) 中二進位制最高位為 0 的部分與二進位制最高位為 1 的部分。
異或卷積有著良好的遞迴性質：
\[A\oplus B=(A_0\oplus B_0+A_1\oplus B_1,A_1\oplus B_0+A_0\oplus B_1)\]
怎麼理解呢？根據我剛剛說的 \(A_0,A_1\) 的定義，討論最高位為 0 或 1 就可以得到如上的式子。
這些運算有以下幾個性質。
（1）\(A\oplus B=B\oplus A,A+B=B+A\)（交換律）
（2）\(A\oplus B\oplus C=A\oplus (B\oplus C),A+B+C=A+(B+C)\)（結合律）
（3）\(A\oplus (B+C)=A\oplus B+A\oplus C\)（分配律）
證明可以根據定義，拆 \(\sum\) 得到。
@2 - 快速沃爾什正變換（異或）@
回想我們的加法卷積 \(C_k=\sum_{i+j=k}a_i*b_j\)
我們根據 FFT 的優化思路，來得到異或卷積的高效演算法 FWT。
在搞 FFT 的時候，我們是先把原多項式 \(A(x)\) 由係數表示轉為點值表示，再直接逐個相乘，最後由點值表示轉回成係數表示。
將係數表示的多項式抽象成一個向量 \(A=(a_0,a_1,...,a_{n-1})\)
將點值表示的多項式也抽象成一個向量 \(A'=(a'_0,a'_1,...,a'_{n-1})\)，其中\(a'_i=A(w_n^i)\)。
則 FFT 的優化思路就可以抽象出來：將原向量\(A\)轉化為新向量\(A'\)，使原向量的卷積\(A\oplus B\)對應著新向量的乘積\(A'*B'\)，最後把乘積轉化為卷積。
我們定義一個函式 \(FWT(A) = A'\)，一個定義域和值域都為向量的函式，使得 \(FWT(A\oplus B)=FWT(A)*FWT(B)\)。再定義它的反函式 \(IFWT(A') = A\)。
則我們的 FFT 相當於選用了一個特殊的 \(FWT\) 函式，可以讓我們快速求解\(FWT(A),IFWT(A')\)。
對於異或卷積而言，我們構造性地選取 \(FWT\)。
定義異或卷積的 \(FWT\) 函式為：
\[FWT(A)=\begin{cases} (FWT(A_0)+FWT(A_1),FWT(A_0)-FWT(A_1)) & (n>0)\\ A & (n=0)\end{cases}\]
【敲黑板，劃重點，這個必考】
我們來驗證一下這玩意具有上面那個性質。
性質（1）：\(FWT(A+B)=FWT(A)+FWT(B)\)
證明：\(FWT(A+B)\) 中的每一項都是線性組合，故滿足分配律。
性質（2）：\(FWT(A\oplus B)=FWT(A)*FWT(B)\)
證明：數學歸納法。【集中注意力！！！】
當 n = 1 時，顯然成立。
假設對於 k <= n-1 結論成立，則：
\[\begin{aligned} FWT(A\oplus B)&=FWT((A_0\oplus B_0+A_1\oplus B_1,A_1\oplus B_0+A_0\oplus B_1))&\\ &=(FWT(A_0\oplus B_0+A_1\oplus B_1)+FWT(A_1\oplus B_0+A_0\oplus B_1),FWT(A_0\oplus B_0+A_1\oplus B_1)-FWT(A_1\oplus B_0+A_0\oplus B_1))&\\ &=(FWT(A_0\oplus B_0+A_1\oplus B_1+A_1\oplus B_0+A_0\oplus B_1),FWT(A_0\oplus B_0+A_1\oplus B_1-A_1\oplus B_0-A_0\oplus B_1))&\\ &=(FWT((A_0+A_1)\oplus (B_0+B_1)),FWT((A_0-A_1)\oplus (B_0-B_1))&\\ &=(FWT(A_0+A_1)*FWT(B_0+B_1),FWT(A_0-A_1)*FWT(B_0-B_1))&\\ &=(FWT(A_0+A_1),FWT(A_0-A_1))*(FWT(B_0+B_1),FWT(B_0-B_1))&\\ &=FWT(A)*FWT(B) \end{aligned}\]
稍微分析我們每個等號都用了哪些性質。
第一個：用了異或卷積本身的遞迴性。
第二個：用了 FWT 的定義。
第三個：用了性質 \(FWT(A+B)=FWT(A)+FWT(B)\)。
第四個：用了因式分解。
第五個：用了我們數學歸納法的歸納前提。
第六個：用了一個恆等變換（類因式分解），可以驗證等式左右兩邊是相等的。
第七個：逆用了 FWT 的定義。
然後得證我們的結論。
……
只能說很壯觀。當初想出來這玩意兒的人真的是不容易。
所以我選擇直接背定義。
根據 FWT 的定義，我們可以進行迭代求解。時間複雜度與 FFT 相同，為O(nlog n)。
@3 - 快速沃爾什逆變換（異或）@
好了，正變換搞定了，我們可以寫 FWT 了……
稍等。好像還差些什麼。
……
好的我們來看看逆變換怎麼弄。
由定義可得：\(FWT(A)=(FWT(A_0)+FWT(A_1),FWT(A_0)-FWT(A_1))=A'\)。
因此有\(\begin{cases} A'_0=FWT(A_0)+FWT(A_1)&\\ A'_1=FWT(A_0)-FWT(A_1) \end{cases}\)
解一下這個二元方程組就可得：\(\begin{cases} FWT(A_0)=\frac{A'_0+A'_1}{2}&\\ FWT(A_1)=\frac{A'_0-A'_1}{2} \end{cases}\)
兩邊同時進行逆變換可以得到：\(\begin{cases} A_0=IFWT(\frac{A'_0+A'_1}{2})&\\ A_1=IFWT(\frac{A'_0-A'_1}{2}) \end{cases}\)
所以：\(IFWT(A')=A=(A_0,A_1)=(IFWT(\frac{A'_0+A'_1}{2}),IFWT(\frac{A'_0-A'_1}{2}))\)
簡化一下，加上邊界： \[IFWT(A)=\begin{cases} (IFWT(\frac{A'_0+A'_1}{2}),IFWT(\frac{A'_0-A'_1}{2})) & (n>0)\\ A & (n=0)\end{cases}\]
因此，我們也可以在 O(nlog n) 的時間複雜度內實現逆變換。
@4 - 與卷積、或卷積@
要知道二進位制的運算可不止異或。
定義與卷積，或卷積如下：
\[A|B=(\sum_{i|j=0}a_i*b_j,\sum_{i|j=1}a_i*b_j,\dots,\sum_{i|j=n-1}a_i*b_j)\\ A\&B=(\sum_{i\&j=0}a_i*b_j,\sum_{i\&j=1}a_i*b_j,\dots,\sum_{i\&j=n-1}a_i*b_j)\]
我們直接給出與卷積，或卷積的 FWT 與 IFWT 的表示式。
如果想要證明可以自行參照異或的證明思路進行證明 （或者直接背） 。
或卷積：
\[FWT(A)=\begin{cases} (FWT(A_0),FWT(A_1)+FWT(A_0)) & (n>0)\\ A & (n=0)\end{cases}\\ IFWT(A)=\begin{cases} (IFWT(A_0),IFWT(A_1)-IFWT(A_0)) & (n>0)\\ A & (n=0)\end{cases} \]
與卷積：
\[FWT(A)=\begin{cases} (FWT(A_0)+FWT(A_1),FWT(A_1)) & (n>0)\\ A & (n=0)\end{cases}\\ IFWT(A)=\begin{cases} (IFWT(A_0)-IFWT(A_1),IFWT(A_1)) & (n>0)\\ A & (n=0)\end{cases} \]
@5 - 參考程式碼實現@
本程式碼為 洛谷P4717 的 AC 程式碼。
#include<cstdio>
const int MOD = 998244353;
const int INV = (998244353 + 1)/2;
const int MAXN = (1<<17);
void fwt_or(int *a, int n, int type) {
for(int s=2;s<=n;s<<=1) {
int t = (s>>1);
for(int i=0;i<n;i+=s) {
for(int j=0;j<t;j++) {
int x = a[i+j], y = a[i+j+t];
a[i+j] = (type == -1) ? x : x;
a[i+j+t] = (type == -1) ? (y + MOD - x)%MOD : (y + x)%MOD;
}
}
}
}
void fwt_and(int *a, int n, int type) {
for(int s=2;s<=n;s<<=1) {
int t = (s>>1);
for(int i=0;i<n;i+=s) {
for(int j=0;j<t;j++) {
int x = a[i+j], y = a[i+j+t];
a[i+j] = (type == -1) ? (x + MOD - y)%MOD : (x + y)%MOD;
a[i+j+t] = (type == -1) ? y : y;
}
}
}
}
void fwt_xor(int *a, int n, int type) {
for(int s=2;s<=n;s<<=1) {
int t = (s>>1);
for(int i=0;i<n;i+=s) {
for(int j=0;j<t;j++) {
int x = a[i+j], y = a[i+j+t];
a[i+j] = (type == -1) ? 1LL*(x + y)%MOD*INV%MOD : (x + y)%MOD;
a[i+j+t] = (type == -1) ? 1LL*(x + MOD - y)%MOD*INV%MOD : (x + MOD - y)%MOD;
}
}
}
}
int A[MAXN + 5], B[MAXN + 5], C[MAXN + 5];
int main() {
int N; scanf("%d", &N); N = (1<<N);
for(int i=0;i<N;i++)
scanf("%d", &A[i]);
for(int i=0;i<N;i++)
scanf("%d", &B[i]);
fwt_or(A, N, 1); fwt_or(B, N, 1);
for(int i=0;i<N;i++)
C[i] = 1LL*A[i]*B[i]%MOD;
fwt_or(A, N, -1); fwt_or(B, N, -1); fwt_or(C, N, -1);
for(int i=0;i<N;i++)
printf("%d ", C[i]);
puts("");
fwt_and(A, N, 1); fwt_and(B, N, 1);
for(int i=0;i<N;i++)
C[i] = 1LL*A[i]*B[i]%MOD;
fwt_and(A, N, -1); fwt_and(B, N, -1); fwt_and(C, N, -1);
for(int i=0;i<N;i++)
printf("%d ", C[i]);
puts("");
fwt_xor(A, N, 1); fwt_xor(B, N, 1);
for(int i=0;i<N;i++)
C[i] = 1LL*A[i]*B[i]%MOD;
fwt_xor(A, N, -1); fwt_xor(B, N, -1); fwt_xor(C, N, -1);
for(int i=0;i<N;i++)
printf("%d ", C[i]);
puts("");
}@6 - 例題與應用@
和 FFT 一樣，留坑，待填。