快速傅立葉變換

lalaouye發表於2024-03-09

FFT

問題：

設 \(A(x)=\sum_{i=0}^n a_ix^i\)，\(B(x)=\sum_{i=0}^m b_ix^i\)。求 \(A(x)\) 和 \(B(x)\) 的卷積。

有一個結論：座標系中 \(n\) 個點確定一個 \(n-1\) 次函式。

可以這樣理解：\(n-1\) 次函式有 \(n\) 個係數，而 \(n\) 個點相當於 \(n\) 個方程。

於是我們可以換一種思路求卷積：將 \(n+m+1\) 個不同的 \(x\) 代入多項式，那麼其卷積的結果的函式 \(C(x)\) 必然經過 （\(x_i,A(x_i)B(x_i)\)），於是再透過插值就可以得到卷積結果了。

但是，這樣子複雜度仍然不變，那該怎麼辦呢？

於是我們考慮用複數解決問題，在一個單位圓上，圓上的點表示 \(z=\cos\theta+i\sin\theta\)（\(0\le \theta< 2\pi\)）。然後我們把圓分成 \(n\) 等份，便有了單位根：\(w_n^k=\cos{\frac{2\pi k}{n}}+i\sin{\frac{2\pi k}{n}}\)，它有很好的性質：

1.\(w_n^a\times w_n^b=w_n^{a+b}\)，其實就是積化和差直接推式子就可以證了，並且也可以透過這個簡單推出 \((w_n^a)^b=w_n^{ab}\)。不過這個式子其實是最最最重要的，也是唯一用到虛數的性質。

2.\(w_n^{k+\frac{n}{2}}=-w_n^k\)，學過三角函式就會發現很顯然。

3.\(w_n^{2k}=w_{\frac{n}{2}}^k\)，顯然。

4.\(w_n^k=w^{k+n}\)，顯然。

然後，考慮如何快速計算多項式的值，透過思考可以得到：

\[\begin{align} A(x)&=a_0+a_1x^1+\cdots+a_{n-1}x^{n-1}+a_nx^n \\ &=(a_0+a_2x^2+\cdots+a_{n-2}x^{n-2})+(a_1+a_3x^2+\cdots+a_{n-1}x^{n-2})x\\ &=A_1(x^2)+A_2(x^2)x\end{align} \]

注意：\(A_1(x)=a_0+a_2x+a_4x^2+\cdots+a_{n-2}x^{\frac{n}{2}-1}\)，\(A_2\) 同理可推，然後在 FFT 中我們規定 \(n=2^{i\in \mathbb {Z}}\)，並且我們在運算中是不會管第 \(n\) 項的，故 \(n\) 必須大於多項式的次數。

但是這並不能最佳化我們的複雜度，於是我認為 FFT 一個很厲害的地方來了，我們將單位根 \(w_n^0,w_n^1,\cdots,w_n^{n-1}\) 分帶入多項式並求值，那麼：

當 \(0\le k < \frac{n}{2}\) 時：

\[A(w_n^k)=A_1(w_{\frac{n}{2}}^k)+w_n^kA_2(w_{\frac{n}{2}}^k) \]

當 \(\frac{n}{2}\le k <n\) 時：

\[A(w_n^k)=A_1(w_{\frac{n}{2}}^k)-w_n^kA_2(w_{\frac{n}{2}}^k) \]

這樣，我們就有了較大的突破，因為這樣子，我們要求的東西就是可轉移的了，它就相當於一個 \(f(n)=f_1(n)+k\times f_2(n)\)，於是，要求它每一項的值就可以透過遞迴實現了，為了更好的理解這裡給出程式碼:

void FFT (int lim, comp F[])
{
	if (lim == 1) return ;
	comp F1[lim >> 1], F2[lim >> 1];
	rep (i, 0, (lim >> 1) - 1) F1[i] = F[i << 1], F2[i] = F[i << 1 | 1];
	FFT (lim >> 1, F1), FFT (lim >> 1, F2);
	comp w = (comp) {1, 0};
	comp wk = (comp) {cos (2.0 * PI / lim), sin (2.0 * PI / lim)};
	for (int i = 0; i < (lim >> 1); ++ i, w = w * wk)
	{
		F[i] = F1[i] + F2[i] * w;
		F[i + (lim >> 1)] = F1[i] - F2[i] * w;
	}
}

我認為仔細看是可以看懂的，首先我們把係數帶進自己的左右節點的陣列，然後我們發現每一層完成以後自己的陣列變成了每一個 \(F(w_n^k)\) 的答案，然後再轉移即可。

於是每個需要的 (\(x_i,A(x_i)B(x_i)\)) 便求出來了，然後我們需要考慮如何透過已有座標推出多項式係數了。

我們設 \(y_i=A(x_i)B(x_i)\)，\(C(x)=\sum_{i=0}^{n-1} y_ix^i\)，然後我們分別將 \(w_n^0,w_n^{-1},\cdots,w_n^{1-n}\) 帶入進這個多項式，設其結果為 \(z_0,\cdots,z_{1-n}\)，我們開始求 \(z\)：

\[\begin{align} z_k&=\sum_{i=0}^{n-1}y_i(w_n^{-k})^i\notag\\ &=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1} a_j(w_n^i)^j(w_n^{-k})^i\notag\\ &= \sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1} (w_n^{j-k})^i\notag\end{align} \]

此時，當 \(j=k\) 時，內部求和的值顯然為 \(n\)，而當 \(j\neq k\) 時，我們相當於是在求一個等比數列求和，內部求和的值便是 \(\large\frac{(w_n^{j-k})^n-1}{w_n^{j-k}-1}\)，其中可以發現 \((w_n^{j-k})^n\) 就是 \(w_n^0=1\)，於是內部求和的總和就是零，那麼可以得出 \(z_k=na_k\) 即 \(a_k=\frac{z_k}{n}\)。於是我們求出 \(w_n^{-k}\) 就可以求出 \(z_k\) 了，現在考慮怎麼求：

設 \(w_n^k=\cos\theta+isin\theta\)， 而 \(w_n^{-k}=w^{n-k}\)，所以 \(w_n^{-k}=\cos\theta-i\sin\theta\)，這個東西也可以透過 FFT 求，於是我們就學會了 FFT 了！接下來上的程式碼：

void FFT (int lim, comp F[], int op)
{
	if (lim == 1) return ;
	comp F1[lim >> 1], F2[lim >> 1];
	rep (i, 0, (lim >> 1) - 1) F1[i] = F[i << 1], F2[i] = F[i << 1 | 1];
	FFT (lim >> 1, F1, op), FFT (lim >> 1, F2, op);
	comp w = (comp) {1, 0};
	comp wk = (comp) {cos (2.0 * PI / lim), sin (2.0 * PI / lim) * op};
	for (int i = 0; i < (lim >> 1); ++ i, w = w * wk)
	{
		F[i] = F1[i] + F2[i] * w;
		F[i + (lim >> 1)] = F1[i] - F2[i] * w;
	}
}
signed main ()
{
	// freopen ("1.in", "r", stdin);
	// freopen ("1.out", "w", stdout);
	n = rd (), m = rd ();
	int l = 0;
	for (lim = 1; lim <= n + m; lim <<= 1, ++ l) ;
	rep (i, 1, lim) R[i] = (R[i >> 1] >> 1) | (i & 1) << (l - 1);
	rep (i, 0, n) A[i].x = rd ();
	rep (i, 0, m) B[i].x = rd ();
	FFT (lim, A, 1), FFT (lim, B, 1);
	rep (i, 0, lim) A[i] = A[i] * B[i];
	FFT (lim, A, -1);
	rep (i, 0, n + m) printf ("%lld ", (int) (A[i].x / lim + 0.5));
}

這裡的 \(op\) 就是表示帶進去的是 \(w_n^k\) 還是 \(w_n^{-k}\)，然後還有要注意的就是 \(lim\) 其實就是重置成 \(2\) 的次冪的新的 \(n\)，顯然對於答案不會產生任何影響，並且複雜度跟線段樹一樣，\(O(n\log n)\)。

這個寫法效率很慢，有一個叫蝴蝶變換的更優秀的規律可以最佳化它，這裡不多贅述。

相關文章