FFT,FWT演算法比對——HDU6966. I love sequences(FWT)
卷積定義:現有\(A,B\)兩個序列,\(A*B\)為其卷積
\[(A*B)_k=\sum\limits_{i+j=k}A_iB_j \]對其進行拓展,可以將\(i+j=k\)中的加號換成其它的二元運算子,例如與,或,異或等等。
FFT——快速傅立葉變換
快速傅立葉變換就是針對\(i+j=k\)型別的卷積的。
設\(T\)為傅立葉變換,我們要求這個變換滿足
\[T(A)_kT(B)_k=T(A*B)_k \]這樣計算\(A*B\)就可以通過計算\(T(A),T(B)\)得到\(T(A*B)\),再對其進行\(T^{-1}\)變換得到。
\(T\)是一種線性變換,所以我們可以將\(T(A)_k\)
要使得\(T(A)_kT(B)_k=T(A*B)_k\)成立,就要
\[\begin{aligned} \left(\sum\limits_{i=0}^{n-1}A_if(k,i)\right)\left(\sum\limits_{j=0}^{n-1}B_jf(k,j)\right)&=\sum\limits_{m=0}^{n-1}\sum\limits_{i+j=m}A_iB_jf(k,m)\\ \sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}A_iB_jf(k,i)f(k,j)&=\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}A_iB_jf(k,i+j)\\ \end{aligned} \]即需要
顯然我們可以令\(f(k,i)=a^i\),理論上\(a\)的取值可以任意。
在FFT中,選擇了\(n\)次單位根的冪次\(\omega_n^m(0\leq m<n)\)作為\(a\)。
這樣保證了有\(n\)組相互之間線性無關的線性變換,這保證了逆變換還能變得回來。
FFT採用分治演算法進行,根據奇偶,將原序列分成兩部分,\(A_{(\dots 0)_2}\)以及\(A_{(\dots 1)_2}\)
\[\begin{aligned} T(A)_k&=\sum\limits_{i=0}^{n-1}A_i\omega_n^{ki}\\ &=\sum\limits_{i=0}^{\frac{n}{2}-1}A_{2i}\omega_n^{2ki}+\sum\limits_{i=0}^{\frac{n}{2}-1}A_{2i+1}\omega_n^{2ki+k}\\ &=\sum\limits_{i=0}^{\frac{n}{2}-1}A_{2i}\omega_{\frac{n}{2}}^{ki}+\omega_n^k\sum\limits_{i=0}^{\frac{n}{2}-1}A_{2i+1}\omega_\frac{n}{2}^{ki}\\ &=T(A_{(\dots 0)_2})_k+\omega_n^kT(A_{(\dots 1)_2})_k&(k<\frac{n}{2})\\ T(A)_k&=\sum\limits_{i=0}^{n-1}A_i\omega_n^{ki}\\ &=\sum\limits_{i=0}^{\frac{n}{2}-1}A_{2i}\omega_n^{2ki}+\sum\limits_{i=0}^{\frac{n}{2}-1}A_{2i+1}\omega_n^{2ki+k}\\ &=\sum\limits_{i=0}^{\frac{n}{2}-1}A_{2i}\omega_{\frac{n}{2}}^{ki}+\omega_n^k\sum\limits_{i=0}^{\frac{n}{2}-1}A_{2i+1}\omega_\frac{n}{2}^{ki}\\ &=\sum\limits_{i=0}^{\frac{n}{2}-1}A_{2i}\omega_{\frac{n}{2}}^{(k-\frac{n}{2})i}+\omega_n^k\sum\limits_{i=0}^{\frac{n}{2}-1}A_{2i+1}\omega_\frac{n}{2}^{(k-\frac{n}{2})i}\\ &=T(A_{(\dots 0)_2})_{k-\frac{n}{2}}+\omega_n^kT(A_{(\dots 1)_2})_{k-\frac{n}{2}}\\ &=T(A_{(\dots 0)_2})_{k-\frac{n}{2}}-\omega_n^{k-\frac{n}{2}}T(A_{(\dots 1)_2})_{k-\frac{n}{2}}&(k\geq\frac{n}{2})\\ \end{aligned} \]這樣就可以進行FFT
最終形成了一個線性變換,
不難得出,最終的變換矩陣為
\[A=\begin{bmatrix} (\omega_n^0)^0&(\omega_n^0)^1&\cdots&(\omega_n^0)^{n-1}\\ (\omega_n^1)^0&(\omega_n^1)^1&\cdots&(\omega_n^1)^{n-1}\\ \vdots&\vdots&\ddots&\vdots\\ (\omega_n^{n-1})^0&(\omega_n^{n-1})^1&\cdots&(\omega_n^{n-1})^{n-1} \end{bmatrix}\\ \]其逆矩陣
\[A^{-1}=\frac{1}{n}\begin{bmatrix} (\omega_n^0)^0&(\omega_n^0)^1&\cdots&(\omega_n^0)^{n-1}\\ (\omega_n^{-1})^0&(\omega_n^{-1})^1&\cdots&(\omega_n^{-1})^{n-1}\\ \vdots&\vdots&\ddots&\vdots\\ (\omega_n^{1-n})^0&(\omega_n^{1-n})^1&\cdots&(\omega_n^{1-n})^{n-1} \end{bmatrix}\\ \]FWT——快速沃爾什變換
快速沃爾什變換是針對\(i\oplus j=k\),這裡的\(\oplus\)不一定指異或,可以是任意的二元按位運算子(不一定是二進位制,但下面以二進位制為例)
設\(T\)為沃爾什變換,我們要求這個變換滿足
\[T(A)_kT(B)_k=T(A*B)_k \]這樣計算\(A*B\)就可以通過計算\(T(A),T(B)\)得到\(T(A*B)\),再對其進行\(T^{-1}\)變換得到。
\(T\)是一種線性變換,所以我們可以將\(T(A)_k\)表示成如下的式子
\[T(A)_k=\sum\limits_{i=0}^{n-1}A_if(k,i) \]要使得\(T(A)_kT(B)_k=T(A*B)_k\)成立,就要
\[\begin{aligned} \left(\sum\limits_{i=0}^{n-1}A_if(k,i)\right)\left(\sum\limits_{j=0}^{n-1}B_jf(k,j)\right)&=\sum\limits_{m=0}^{n-1}\sum\limits_{i\oplus j=m}A_iB_jf(k,m)\\ \sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}A_iB_jf(k,i)f(k,j)&=\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}A_iB_jf(k,i\oplus j)\\ \end{aligned} \]即需要
\[f(k,i)f(k,j)=f(k,i\oplus j) \]這裡以二進位制異或為例,認為\(\oplus\)為二進位制異或。
由於是按位運算,我們不妨逐位考慮,令\(f(k,i)=\prod\limits_x f(k,i_x)\)
那麼\(f(k,i)f(k,j)=f(k,i\oplus j)\)成立當且僅當下面\(3\)個式子成立
\[f(k,0)f(k,0)=f(k,0)\\ f(k,1)f(k,1)=f(k,0)\\ f(k,1)f(k,0)=f(k,1) \]可以解得這\(3\)組解(並不是任取一組就好的,具體原因下面說)
\[f(k,0)=0,f(k,1)=0\\ f(k,0)=1,f(k,1)=-1\\ f(k,0)=1,f(k,1)=1 \]同樣地,可以採用分治的辦法,根據最高位可以是\(0\),\(1\)可以將序列分為\(A_{(0\dots)_2},A_{(1\dots)_2}\)
\[\text{$i_0$為$i$的最高位,$i_1$為$i$去掉最高位}\\ \begin{aligned} T(A)_k&=\sum\limits_{i=0}^{n-1}A_if(k,i)\\ &=\sum\limits_{i=0}^{\frac{n}{2}-1}A_if(k,i)+\sum\limits_{i=\frac{n}{2}}^{n-1}A_if(k,i)\\ &=f(k_0,0)\sum\limits_{i=0}^{\frac{n}{2}-1}A_if(k_1,i_1)+f(k_0,1)\sum\limits_{i=0}^{\frac{n}{2}-1}A_if(k_1,i_1)\\ &=f(k_0,0)T(A_{(0\dots)_2})_{k_1}+f(k_0,1)T(A_{(1\dots)_2})_{k_1} \end{aligned} \]注意,式子形式雖然一樣,但是對於\(k<\frac{n}{2}\)和\(k\geq\frac{n}{2}\)兩種情況,如果\(f(k_0,0),f(k_0,1)\)取的值分別相同,那麼會造成前一半和後一半是一樣的,顯然變換無法保持可逆。
所以我們要在上面解出的若干個解中,選出\(2\)個線性無關的解,這裡別無選擇,只能選擇
\[f(k,0)=1,f(k,1)=-1\\ f(k,0)=1,f(k,1)=1 \]讓其中之一(隨便哪一個)作為\(k<\frac{n}{2}\)時的\(f(k_0,0),f(k_0,1)\),另外一組作為\(k\geq\frac{n}{2}\)時的\(f(k_0,0),f(k_0,1)\)
最終線性變換不太好寫,可以寫每一步的矩陣
\[A=\begin{bmatrix} 1&-1\\ 1&1 \end{bmatrix} \]逆矩陣
\[A^{-1}=\begin{bmatrix} \frac{1}{2}&\frac{1}{2}\\ -\frac{1}{2}&\frac{1}{2} \end{bmatrix} \]HDU6966. I love sequences
題面:
問題描述:
給定三個序列\(a=[a_1,a_2,\dots,a_n],b=[b_1,b_2,\dots,b_n],c=[c_1,c_2,\dots,c_n]\)每一個含有\(n\)個非負整數。你需要計算
\[\sum\limits_{p=1}^{n}\sum\limits_{k=1}^{+\infty}d_{p,k}c_p^k \]其中序列\(d=[d_{p,1},d_{p,2},\dots,d_{p,n}]\)按以下的規則計算:
對於每一個整數\(p(1\leq p\leq n)\),每一個整數\(k(1\leq k\leq n)\),有
\[d_{p,k}=\sum\limits_{k=i\oplus j}a_ib_j\\ \text{$i,j$的範圍為$1\leq i\leq \frac{n}{p},1\leq j\leq \frac{n}{p}$} \]其中\(\oplus\)代表一個三進位制二元按位運算子,滿足\(k_x=\gcd(i_x,j_x)\)。
比如\(12\oplus 7=(110)_3\oplus (021)_3=(111)_3=13\)
輸入:
這個問題只有一個樣例。
第一行包含一個整數\(n(n\leq2\times10^5)\),表示序列\(a,b,c\)的長度
第二行包含\(n\)個整數,表示\(a_1,a_2,\dots,a_n\)
第三行包含\(n\)個整數,表示\(b_1,b_2,\dots,b_n\)
第四行包含\(n\)個整數,表示\(c_1,c_2,\dots,c_n\)
輸出:
輸出一個整數,為答案模\(10^9+7\)
樣例輸入:
4
2 3 4 5
1 2 3 4
9 8 7 6
樣例輸出:
1643486
分析:
位運算卷積考慮FWT,考慮\(f(k,i)f(k,j)=f(k,i\oplus j)\)成立的充要條件
\[f(k,0)f(k,0)=f(k,0)\\ f(k,1)f(k,1)=f(k,1)\\ f(k,2)f(k,2)=f(k,2)\\ f(k,0)f(k,1)=f(k,1)\\ f(k,0)f(k,2)=f(k,2)\\ f(k,1)f(k,2)=f(k,1)\\ \]解得
\[f(k,0)=0,f(k,1)=0,f(k,2)=0\\ f(k,0)=1,f(k,1)=0,f(k,2)=0\\ f(k,0)=1,f(k,1)=0,f(k,2)=1\\ f(k,0)=1,f(k,1)=1,f(k,2)=1\\ \]同樣地,使用分治
\[\text{$i_0$為$i$的最高位,$i_1$為$i$去掉最高位}\\ \begin{aligned} T(A)_k&=\sum\limits_{i=0}^{n-1}A_if(k,i)\\ &=\sum\limits_{i=0}^{\frac{n}{3}-1}A_if(k,i)+\sum\limits_{i=\frac{n}{3}}^{\frac{2n}{3}-1}A_if(k,i)+\sum\limits_{i=\frac{2n}{3}}^{n-1}A_if(k,i)\\ &=f(k_0,0)\sum\limits_{i=0}^{\frac{n}{3}-1}A_if(k_1,i_1)+f(k_0,1)\sum\limits_{i=\frac{n}{3}}^{\frac{2n}{3}-1}A_if(k_1,i_1)+f(k_0,2)\sum\limits_{i=\frac{2n}{3}}^{n-1}A_if(k_1,i_1)\\ &=f(k_0,0)T(A_{(0\dots)_3})_{k_1}+f(k_0,1)T(A_{(1\dots)_3})_{k_1}+f(k_0,2)T(A_{(2\dots)_3})_{k_1} \end{aligned} \]\(f(k_0,x)\)需要根據\(k<\frac{n}{3},\frac{n}{3}\leq k<\frac{2n}{3},\frac{2n}{3}\leq k<n\)三種情況,選擇三個線性無關的取值。
所以只能取
\[f(k,0)=1,f(k,1)=0,f(k,2)=0\\ f(k,0)=1,f(k,1)=0,f(k,2)=1\\ f(k,0)=1,f(k,1)=1,f(k,2)=1\\ \]形成矩陣
\[A=\begin{bmatrix} 1&0&0\\ 1&0&1\\ 1&1&1 \end{bmatrix} \]其逆矩陣
\[A^{-1}=\begin{bmatrix} 1&0&0\\ 0&-1&1\\ -1&1&0 \end{bmatrix} \]程式碼:
#include <cstdio>
typedef long long Lint;
const int maxn = 540000;
const int mod = 1e9 + 7;
int a[maxn], b[maxn], c[maxn];
int aa[maxn], bb[maxn], cc[maxn];
void FWT(int* a, int n, int flag) {
for (int len = 1; len < n; len *= 3) {
for (int i = 0; i < n; i += len * 3) {
for (int j = 0; j < len; j++) {
int x = a[i + j], y = a[i + j + len], z = a[i + j + 2 * len];
if (flag == 1) {
a[i + j + len] = (x + z) % mod;
a[i + j + 2 * len] = ((x + y) % mod + z) % mod;
} else {
a[i + j + len] = (mod - y + z) % mod;
a[i + j + 2 * len] = (mod - x + y) % mod;
}
}
}
}
}
int cal(int n, int c) {
int tot = 1;
while (tot <= n)
tot *= 3;
for (int i = 1; i <= n; i++) {
aa[i] = a[i];
bb[i] = b[i];
}
for (int i = n + 1; i < tot; i++) {
aa[i] = bb[i] = 0;
}
FWT(aa, tot, 1), FWT(bb, tot, 1);
for (int i = 0; i < tot; i++) {
cc[i] = (Lint)aa[i] * bb[i] % mod;
}
FWT(cc, tot, -1);
int res = 0;
for (int k = 1, ck = c; k < tot; k++, ck = (Lint)ck * c % mod) {
res = (res + (Lint)cc[k] * ck % mod) % mod;
}
return res;
}
int main() {
int n;
scanf("%d", &n);
for (int i = 1; i <= n; i++) {
scanf("%d", a + i);
}
for (int i = 1; i <= n; i++) {
scanf("%d", b + i);
}
for (int i = 1; i <= n; i++) {
scanf("%d", c + i);
}
int ans = 0;
for (int p = 1; p <= n; p++) {
ans = (ans + cal(n / p, c[p])) % mod;
}
printf("%d\n", ans);
return 0;
}