FFT（快速傅立葉變換）
- 概念引入
- 點值表示
對於一個$n - 1$次多項式$A(x)$，可以通過確定$n$個點與值（即$x$和$y$）來表示這唯一的$A(x)$
- 複數
對於一元二次方程
$$x^2 + 1 = 0$$
在實數範圍內無解，那麼我們將實數範圍擴充，就得到了複數，再令$i$為該方程的解，即
$$i^2 = - 1$$
那麼就定義$z = a + bi$的數為複數，則有
當$b = 0$時，$z$為實數
當$b \neq 0$時，$z$為虛數
當$a = 0 \ \&\& \ b \neq 0$時，$z$為純虛數
其中，複數滿足加法交換律、結合律、乘法分配率等
複數相乘：$z_1 = a_1 + b_1i, \ z_2 = a_2 + b_2i$，則$z_1 × z_2 = (a_1 + b_1i)(a_2 + b_2i) = (a_1a_2 - b_1b_2) + (a_1b_2 + b_1a_2)i$，它們相乘還是一個複數，在複平面上理解，就是模長相乘，幅角相加
共軛複數：當兩個複數實部相同，虛部為相反數時，兩個複數被稱為共軛複數
對於一個複數$z = a + bi$，還可以在複數平面上用向量表示出來，即有$\overrightarrow{OZ}$一一對應每個$z = a + bi$，那麼就有複數的模等於$\overline{Z}$，即$|z| = \sqrt{a^2 + b^2}$
複數直接比較大小沒有意義，除非它是一個實數
- $FFT$作用
那麼，有了點值表示，對於多項式$A(x)$和$B(x)$相乘，就不需要$O(n^2)$，而只需要$O(n)$，因為$C(x_i) = A(x_i) * B(x_i)$，然後列舉$x_i$即可
於是現在的複雜度癥結就在於將多項式轉化成點值表示的$O(n^2)$了，於是就有$FFT$，可以實現在$O(n \ logn)$的時間內轉化
- 離散傅立葉變換
於是傅立葉規定，點值中的$x$為模長為$1$的複數
至於為什麼要取複數而不是實數，因為它有很神奇的性質
那麼對於這$n$個複數的取法，取的是複平面上半徑為$1$的一個圓，將其$n$等分後得到的點，令第$i$個點表示為$\omega_n^k(0 \le k \le n - 1)$，那麼這個點在複平面中的表示即為$(cos\frac{k}{n}2\pi, \ sin\frac{k}{n}2\pi)$，那麼根據複數乘法在複平面上的意義為模長相乘，幅長相加，即$\omega_n^k$相當於$(\omega_n^1)^k$，那麼稱$\omega_n^1$為單位根
我們就把這$\omega_n^0, \ \omega_n^1, \ ..., \ \omega_n^{n - 1}$作為點值表示的$x$，稱作離散傅立葉變換的結果
先給出關於傅立葉逆變換的結論：
- 將多項式$A(x)$的離散傅立葉變換結果作為係數代入多項式$B(x)$，再將每個離散傅立葉變換結果的倒數，即$\omega_n^0, \ \omega_n^{- 1}, \ ..., \ \omega_n^{- (n - 1)}$作為$x$代入$B(x)$得到點值表示，那麼這些表示除以$n$就得到了$A(x)$的原係數 [至於證明：不存在的]
這就是傅立葉變換神奇的性質
- $FFT$
有了傅立葉變換，雖然多項式與點值表示相互轉化已經很輕鬆了，但是複雜度仍然不理想，就有了快速傅立葉變換
結論：
- $\omega_{2n}^{2k} = \omega_n^k$ [證明：代進原公式顯然]
- $\omega_n^{k + \frac{n}{2}} = - \omega_n^k$ [證明：關於複平面原點中心對稱]
對於多項式$A(x) = \sum\limits_{i = 0}^{n - 1} a_ix^i$，將它的奇偶項拆開，並將$x$轉為$x^2$得到（此處先假定$n$為偶數）
$$N(x) = a_0 + a_2x + a_4x^2 + ... + a_{n - 2}x^{\frac{n}{2} - 1}$$
$$M(x) = a_1 + a_3x + a_5x^2 + ... + a_{n - 1}x^{\frac{n}{2} - 1}$$
當然此時代入的是$N(x^2), \ M(x^2)$，則有
$$A(x) = N(x^2) + xM(x^2)$$
此時就需要把$\omega_n^k$代入，分情況討論：
若$k < \frac{n}{2}$，則有
$$A(\omega_n^k) = N(\omega_n^{2k}) + \omega_n^kM(\omega_n^{2k})$$
$$\ \ \ \ \ \ \ \ \ \ \ = N(\omega_{\frac{n}{2}}^k) + \omega_n^kM(\omega_{\frac{n}{2}}^k)$$
反之
$$A(\omega_n^k) = N(\omega_n^{2k + n}) + \omega_n^{k + \frac{n}{2}}M(\omega_n^{2k + n})$$
$$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ = N(\omega_n^{2k} × \omega_n^n) + \omega_n^{k + \frac{n}{2}}M(\omega_n^{2k} × \omega_n^n)$$（這一步是通過複數乘法複平面上的意義化簡的）
$$= N(\omega_{\frac{n}{2}}^k) - \omega_n^kM(\omega_{\frac{n}{2}}^k)$$
於是，我們只要求得$\{\omega_{\frac{n}{2}}^0, \ \omega_{\frac{n}{2}}^1, \ ..., \ \omega_{\frac{n}{2}}^{\frac{n}{2} - 1}\}$，就可以得到$A(x)$的所有關於所有離散傅立葉變換結果的點值表示了，可用分治實現，複雜度$O(n \ logn)$
- 程式碼（分治FFT）
1 #include <iostream> 2 #include <cstdio> 3 #include <cstring> 4 #include <algorithm> 5 #include <cmath> 6 7 using namespace std; 8 9 const int MAXN = 2e06 + 10; 10 11 const double PI = acos (- 1.0); 12 13 int N, M; 14 15 struct mcomplex { 16 double x, y; 17 18 mcomplex () {} 19 mcomplex (double fx, double fy) : 20 x (fx), y (fy) {} 21 22 mcomplex operator + (const mcomplex& p) const { 23 return mcomplex (x + p.x, y + p.y); 24 } 25 mcomplex operator - (const mcomplex& p) const { 26 return mcomplex (x - p.x, y - p.y); 27 } 28 mcomplex operator * (const mcomplex& p) const { 29 return mcomplex (x * p.x - y * p.y, x * p.y + y * p.x); 30 } 31 } ; 32 mcomplex omega (int n, int k, int inv) { 33 return mcomplex (cos (2.0 * PI * (double) k / (double) n), 1.0 * inv * sin (2.0 * PI * (double) k / (double) n)); // 當一個複數的模長為1時，它的共軛複數即為它的倒數 34 } 35 mcomplex A[MAXN], B[MAXN]; 36 37 void FFT (mcomplex* a, int n, int inv) { // inv表示是否取倒數 38 if (n == 1) 39 return ; 40 mcomplex a1[n >> 1], a2[n >> 1]; 41 int m = n >> 1; 42 for (int i = 0; i <= n; i += 2) { 43 a1[i >> 1] = a[i]; 44 a2[i >> 1] = a[i + 1]; 45 } 46 FFT (a1, m, inv); 47 FFT (a2, m, inv); 48 for (int i = 0; i < m; i ++) { 49 mcomplex x = omega (n, i, inv); 50 a[i] = a1[i] + x * a2[i]; 51 a[i + m] = a1[i] - x * a2[i]; 52 } 53 } 54 55 int getnum () { 56 int num = 0; 57 char ch = getchar (); 58 59 while (! isdigit (ch)) 60 ch = getchar (); 61 while (isdigit (ch)) 62 num = (num << 3) + (num << 1) + ch - '0', ch = getchar (); 63 64 return num; 65 } 66 67 int main () { 68 N = getnum (), M = getnum (); 69 for (int i = 0; i <= N; i ++) 70 A[i].x = (double) getnum (); 71 for (int i = 0; i <= M; i ++) 72 B[i].x = (double) getnum (); 73 74 int n; 75 for (n = 1; n <= N + M; n <<= 1); 76 FFT (A, n, 1); 77 FFT (B, n, 1); 78 for (int i = 0; i <= n; i ++) 79 A[i] = A[i] * B[i]; 80 FFT (A, n, - 1); 81 for (int i = 0; i <= N + M; i ++) { 82 if (i) 83 putchar (' '); 84 printf ("%d", (int) (A[i].x / n + 0.5)); 85 } 86 puts (""); 87 88 return 0; 89 } 90 91 /* 92 1 2 93 1 2 94 1 2 1 95 */ 96 97 /* 98 5 5 99 1 7 4 0 9 4 100 8 8 2 4 5 5 101 */分治FFT
- 參考資料