目錄

• 快速傅立葉變換FFT
  • 用途
  • 前置知識
    • 係數表示法
    • 點值表示法
    • 單位負根
      • 定義
      • 性質
  • 快速傅立葉變換FFT
  • 快速傅立葉逆變換
    • 作用
    • 方法與推導

快速傅立葉變換FFT

用途

\(\operatorname{FFT}\)演算法支援在\(O(n log n)\)時間內計算兩個\(n\)度的多項式的乘法。也可以用來加速大整數乘法運算。

---

前置知識

係數表示法

用一個多項式的各項係數來表達這個多項式（升冪順序）：

\[\large f(x) = a_0 + a_1x + a_2x^2 + \ldots + a_nx^n \iff f(x) = \{a_0,a_1, \ldots ,a_n \} \]

則它的計算公式為：

\[\large f(x) = \sum _{i = 0}^{n} a_i * x^i \]

時間複雜度：\(O(n^2)\)

點值表示法

把一個多項式看成一個座標系中的函式影象，從影象上選取\(n +1\)個點，利用這\(n + 1\)個點來唯一地表示這個函式。

為什麼\(n + 1\)個點可以唯一地表示\(f(n)\)

\[由小學數學可知，兩點確定一條直線（一次函式），三點確定一條拋物線（二次函式）... 以此類推，n + 1個點可以確定一個n次函式，就是變數x的最高次冪為n的函式。 \]

知道了正確性，那麼我們可以設：

\[\large f(x_0) = y_0 = a_0 + a_1x_0 + a_2x_0^2 + \ldots +a_nx_0^n \]\[\large f(x_1) = y_1 = a_0+ a_1x_1 + a_2x_1^2 + \ldots + a_nx_1^n \]\[\large f(x_2)= y_2= a_0+a_1x_2+ a_2x_2^2+\ldots+ a_nx_2^n \]\[\ldots \]\[\large f(x_n) = y_n = a_0 + a_1x_n + a_2x_n^2+\dots + a_nx_n^n \]

那麼用點值表示法表示\(f(x)\)

如下：

\[\large f(x) = a_0 + a_1x + a_2x^2 + a_3x^3 + \ldots +a_nx^n \iff f(x) = \{(x_0,y_0),(x_1,y_1),(x_2,y_2),(x_3,y_3),\ldots,(x_n,y_n)\} \]

它的計算公式為：

\[\large f(i) = \sum _{j = 0}^{n - 1} a_j * x_i^j \]

時間複雜度：\(O(n^2)\)

單位負根

定義

\(x^n = 1\)在複數意義下的解釋\(n\)次負根。這樣的解有\(n\)個。

設：\(\large \omega_n = e^{\frac{2 \pi i}{n}}\)

則：\(x^n\)的解集表示為：

\[\large \{\omega_n | k = 0,1,\ldots,n - 1\} \]

稱\(\omega_n\)是\(n\)次單位負根，其他負根均可以用單位負根的冪表示。

由於我數學太菜了所以不會證明。

根據尤拉公式可得：

\[\large \omega_n = e^{\frac{2 \pi i}{n}} = cos(\frac{2\pi}{n}) + i * sin(\frac{2\pi}{n}) \]

柿子很好記，可是我不知道怎麼證，可是記住就行了。

性質

單位負根。對於任意正整數\(n\)和整數\(k\)：

\[\large \omega ^n_n = 1 \]\[\large \omega^k_n = \omega^{2k}_{2n} \]\[\large \omega^{k + n}_{2n} = -\omega^k_{2n} \]

---

快速傅立葉變換FFT

\(\operatorname{FFT}\)的基本思想是分治。

先看\(\operatorname{DFT}\)，它分治地來求當\(x = \omega^k_n\)的時候\(f(x)\)的值，分治思想體現在將多項式分為奇次項和偶次項處理。

令\(f(x)\)為一個\(n\)次的多項式。

令\(G(x)\)表示偶次項係數建立的新函式：

\[\large G(x) = a_0 + a_2x + a_4x^2 +\ldots+ a_{n -2}x^{\frac{n - 2}{2}} \]

令\(H(x)\)表示奇次項係數建立的新函式：

\[\large H(x) = a_1 + a_3x + a_5x^2 +\ldots+a_{n - 1}x^{\lfloor \frac{n - 1}{2}\rfloor} \]

那麼原來的\(f(x)\)用新的兩個函式表示為：

\[\large f(x) = G(x^2) + x * H(x ^ 2) \]

利用單位負根的性質推柿子：

\[\large \operatorname{DFT}(f(\omega^k_n)) = \large\operatorname{DFT}(G((\omega^k_n)^2)) + \operatorname{DFT}(H((\omega^k_n)^2)) \]\[=\large\operatorname{DFT}(G(\omega^{2k}_n)) + \omega^k_n * \operatorname{DFT}(H(\omega^{2k}_n)) \]\[=\large\operatorname{DFT}(G(\omega^k_{\frac{n}{2}})) + \omega^k_n * \operatorname{DFT}(H(\omega^k_{\frac{n}{2}})) \]

同理：

\[\large \operatorname{DFT}(f(\omega^{k + \frac{n}{2}}_n)) = \large\operatorname{DFT}(G((\omega^k_n)^2)) + \operatorname{DFT}(H((\omega^k_n)^2)) \]\[=\large\operatorname{DFT}(G(\omega^{2k}_n)) + \omega^k_n * \operatorname{DFT}(H(\omega^{2k}_n)) \]\[=\large\operatorname{DFT}(G(\omega^k_{\frac{n}{2}})) + \omega^k_n * \operatorname{DFT}(H(\omega^k_{\frac{n}{2}})) \]

因此，我們求出\(\operatorname{DFT}(G(\omega^k_{\frac{n}{2}}))\)和\(DFT(H(\omega^k_{\frac{n}{2}}))\)後，就可以同時求出\(\operatorname{DFT}(f(\omega^k_n)\)和\(\operatorname{DFT}(f(\omega^{k + \frac{n}{2}}_n))\)，於是對\(G\)和\(H\)分別遞迴\(\operatorname {DFT}\)即可。

注意：考慮到分治\(\operatorname {DFT}\)能處理的多項式長度只能是\(2^m(m \in N^*)\)，否則在分治的時候左式和右式長度不等，右式就取不到係數了。所以要在第一次\(\operatorname{DFT}\)之前就把序列向上補成長度為\(2^m(m \in N^*)\)，最高項次數為\(2^m - 1\)的多項式（高次係數補0）。

在代入值的時候，因為要代入\(n\)個不同值，所以我們代入\(\omega^0_n,\omega^1_n,\omega^2_n\ldots,\omega^{n - 1}_n(n = 2^m(m\in N^*))\)一共\(2^m\)個不同值。

---

快速傅立葉逆變換

作用

把點值表示法轉化為係數表示法。

方法與推導

考慮原本的多項式：

\[\large f(x) = a_0 + a_1x + a_2x^2+\ldots+a_{n - 1}x^{n - 1} = \sum^{n - 1}_{i = 0}a_ix^i \]

考慮構造法。

我們已知\(y_i = f(\omega^i_n),i\in\{0,1,\ldots,n - 1\}\)，要求\(\{a_0,a_1,\ldots,a_{n - 1}\}\)

構造多項式：

\[\large A(x) = \sum^{n - 1}_{i = 0}y_ix^i \]

也就是把\(\{y_0,y_1,\ldots,y_{n - 1}\}\)當做多項式\(A\)的係數表示法。

設\(b_i = \omega^{-i}_n\)，則多項式\(A\)在\(x = b_0,b_1,\ldots,b_{n-1}\)處的點值表示法為：\(\{A(b_0),A(b_1),\ldots,A(b_{n - 1})\}\)

對\(A(x)\)的定義式做一下變換，可以將\(A(b_k)\)表示為：

\[\large A(b_k) = \sum^{n-1}_{i=0}f(\omega^i_n)\omega^{-ik}_n \]\[\large =\sum^{n-1}_{i=0}\omega^{-ik}_{n}\sum^{n-1}_{j=0}a_j(\omega^i_n)^j \]\[\large =\sum^{n-1}_{i=0}\sum^{n-1}_{j=0}a_j\omega^{i(j-k)}_{n} \]\[\large =\sum^{n-1}_{j=0}a_j\sum^{n-1}_{i=0}(\omega^{j-k}_{n})^i \]

記\(S( \omega^a_n)=\sum^{n-1}{i=0}(\omega^a_n)^i\)

當\(a=0(\bmod n)\)時，\(S(\omega^a_n) = n\)

當\(a\ne0(\bmod n)\)時，錯位相減：

\[\large S(\omega^a_n)=\sum^{n-i}_{i=0}(\omega^a_n)^i \]\[\large\omega^a_nS(\omega^a_n)=\sum^n_{i=1}(\omega^a_n)^i \]\[\large S(\omega^a_n)=\dfrac{(\omega^a_n)^n-(\omega^a_n)^0}{\omega^a_n-1}=0 \]

也就是：

\[\large S(\omega^a_n) = \begin{cases}n&a=0\\0&a\ne0\end{cases} \]

代回原式得：

\[\large A(b_k)=\sum^{n-1}_{j=0}a_jS(\omega^{j-k}_n)=a_k\cdot n \]

那麼多項式\(A\)的點值表示法為：

\[\begin{aligned}\large \{(b_0,A(b_0),(b_1,A(b_1)),\ldots,(b_{n-1},A(b_{n-1})))\}\\=\large \{(b_0,a_0\cdot n),(b_1,a_1\cdot n),\ldots,(b_{n-1},a_{n-1}\cdot n)\}\end{aligned} \]

總結：我們取單位根為其倒數，對\(\{y_0,y_1,\ldots,y_{n-1}\}\)跑一遍\(\operatorname{FFT}\)，然後除以\(n\)即可得到\(f(x)\)的係數表示。

---

程式碼實現：luogu P3803[模板]多項式乘法(FFT)

#include <bits/stdc++.h>
using namespace std;
//#define ll long long
#define ri register int
const int maxn = 5e6 + 5;
const double pi = acos(-1.0);
int n,m,l,r[maxn],lim = 1;
inline int read() {
	int s = 0,w = 1; char ch = getchar();
	while (ch < '0' || ch > '9') {if (ch == '-') w = -1; ch = getchar();}
	while (ch >= '0' && ch <= '9') s = s * 10 + ch - '0',ch = getchar();
	return s * w;
}
struct cpx {
	double x,y;
	cpx(double a = 0,double b = 0) {x = a; y = b;}
	cpx operator + (const cpx &a) {return cpx(x + a.x,y + a.y);}
	cpx operator - (const cpx &a) {return cpx(x - a.x,y - a.y);}
	cpx operator * (const cpx &a) {return cpx(x * a.x - y * a.y,x * a.y + y * a.x);}
}a[maxn],b[maxn];
inline void fft(cpx *f,int type) {
	for (ri i = 0;i < lim;i++) if (i < r[i]) swap(f[i],f[r[i]]);
	for (ri mid = 1;mid < lim;mid <<= 1) {
		cpx W(cos(pi / mid),type * sin(pi / mid));
		int len = (mid << 1);
		for (ri j = 0;j < lim;j += len) {
			cpx pw(1,0);
			for (ri k = 0;k < mid;k++,pw = pw * W) {
				cpx x = f[j + k],y = pw * f[j + k + mid];
				f[j + k] = x + y;
				f[j + k + mid] = x - y;
			}
		}
	}
}
int main() {
	n = read(); m = read();
	for (ri i = 0;i <= n;i++) a[i].x = read();
	for (ri i = 0;i <= m;i++) b[i].x = read();
	while (lim <= n + m) {lim <<= 1; l++;}
	for (ri i = 0;i < lim;i++) r[i] = (r[i >> 1] >> 1) | (i & 1) << (l - 1);
	fft(a,1); fft(b,1);
//	for (int i = 0;i <= lim;i++) cout << a[i].x << " " << a[i].y << endl;
	for (ri i = 0;i <= lim;i++) a[i] = a[i] * b[i];
	fft(a,-1);
	for (ri i = 0;i <= n + m;i++) printf("%d ",(int)(a[i].x / lim + 0.5));
	return 0;
}

---

本文大部分內容來自：oi-wiki