FFT學得還是有點模糊，原理那些基本還是算有所理解了吧，不過自己推這個推不動。
看的資料主要有這兩個：
http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform
https://www.zybuluo.com/397915842/note/37965
這兒簡單做做筆記。
多項式點值表示
首先$FFT$可以用來快速計算兩個多項式的乘積。
一個$n$次多項式（最高次為$n$），可以用系數表示法表示和點值表示法表示。
- 系數表示法：$A(x)=\sum_{i=0}^{n}c_ix^i$，可以知道用系數表示法進行多項式乘法時間復雜度是$\Theta (n^2)$。
- 點值表示法：用$n+1$個點$(x_k,y_k)，k=0\dots n$，其中$x_k$互不相同，$y_k=A(x_k)$。
這$n+1$個點是可以唯一表示一個$n$次多項式的，因為：$$\begin{bmatrix} 1 & x_0 & x_0^2 & ... & x_0^n \\ 1 & x_1 & x_1^2 & ... & x_1^n \\\vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n & x_n^2 & ... & x_n^n \\\end{bmatrix}\begin{bmatrix} c_0 \\ c_1 \\ \vdots \\c_n \\\end{bmatrix} = \begin{bmatrix} y_0 \\ y_1 \\ \vdots \\y_n \\\end{bmatrix}$$
前面那個方陣是范德蒙德矩陣，由於$x_k$互不相同，其行列式可以知道是非0，所以該方陣存在逆矩陣，方程有解且有唯一解。
現在考慮兩個n次多項式$A(x)$和$B(x)$，其乘積的多項式結果是$C(x)$，那么$C(x)$的次數為$2n$。
不妨給$A$和$B$取$2n+1$個點，這樣求得兩個多項式的點值表示為$(x_k,A(x_k))、(x_k,B(y_k))，k=0\dots 2n$，而$C(x)$的點值表示在$\Theta (n)$就能直接求出，即為$(x_k,A(x_k)B(y_k))$。
這樣進行多項式乘法運算，原本用系數表示多項式需要$\Theta (n^2)$的時間復雜度，而點值則只需要$\Theta (n)$。
不過將一個多項式轉化成點值的形式，需要枚舉$2n$個點，然后對於每一個點$\Theta (n)$求得多項式的值，這樣時間復雜度是$\Theta (n^2)$的；另外，求得結果后也是一個點值的表示，需要還原回系數形式，直接解上面矩陣表示的方程，這要立方的時間復雜度。
這時就是$FFT$做的事了。不妨把前面轉化為點值形式看作$DFT$，后面逆轉化回系數形式看作$IDFT$，而這兩個過程利用$FFT$這個算法就能快速進行並得出結果。
DFT
對於多項式$A(x)$，比如我們要取$n$個點$(x_k,y_k)$，$x_k$選取的是$n$次單位根（滿足$z^n=1$的復數解$z$），即：$$e^{\frac{2\pi ki}{n}}, k = 0\cdots n - 1, i=\sqrt{-1}$$
另外，計算這個可以用歐拉公式：$$e^{\theta i}=\cos\theta + i\sin\theta$$
記$\omega_n=e^{\frac{2\pi i}{n}}$，那么點值表示即為$$(\omega_n^0, A(\omega_n^0)), (\omega_n^1, A(\omega_n^1)),\cdots ,(\omega_n^{n-1}, A(\omega_n^{n-1}))$$
現在如何求$A(\omega_n^0), A(\omega_n^1), \cdots, A(\omega_n^{n-1})$？
看下面：$$A(x) = c_0 + c_1x + c_2x^2 + c_3x^3 + \cdots + c_nx^{n-1}\\A^{[0]}(x) = c_0 + c_2x + c_4x^2 + \cdots + c_{n - 2}x^{n/2 - 1} \\A^{[1]}(x) = c_1 + c_3x + c_5x^2 + \cdots + c_{n - 1}x^{n/2 - 1}$$
然后，可以發現$$A(x) = A^{[0]}(x^2) + xA^{[1]}(x^2)$$
而這邊是將$\omega_n^k，k=0\cdots n-1$代入$x$，其中$x^2$也就是$(\omega_n^k)^2$可以化簡：$$(\omega_n^k)^2=\left(e^{\frac{2\pi ki}{n}}\right)^2=e^{\frac{2\pi ki}{n/2}}=\omega_{\frac{n}{2}}^k$$而$\omega_{\frac{n}{2}}^k$復數根有$n/2$個。
那么可以看到，把$A(x)$分成兩個子問題，而這兩個子問題選取的點也縮減了一半。這樣時間復雜度就是$\Theta (nlogn)$了。於是可以用遞歸去實現這個分治的算法。另外一開始$n$拓展成$2$的次冪，多出來的補$0$即可。
不過，事實上可以改用迭代實現，而不用遞歸。畫出遞歸的每一次划分。。$此處沒圖$。。
然后對比最初和最后一層各個下標的二進制，可以發現相當於是每一個下標的二進制反轉了。
那么可以一開始直接反轉，從最后一層開始，一層層一個個迭代上去去求得結果。
另外，有個現成的$Rader$雷德算法可以直接求得反轉結果。
IDFT
上面就在$\Theta (nlogn)$下將系數表示轉化成點值表示；而現在需要逆過來，把點值表示還原成所需要的結果，系數表示。$$ \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} \begin{bmatrix} c_0 \\ c_1 \\ \vdots \\ c_{n-1} \end{bmatrix} = \begin{bmatrix} A(\omega_n^0) \\ A(\omega_n^1) \\ \vdots \\ A(\omega_n^{n-1}) \end{bmatrix}$$
其實就是求$c_k$。
而那個范德蒙德矩陣的逆也有結論的，就是：$$\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^{-(n-1)})^0 & (\omega_n^{-(n-1)})^1 & \cdots & (\omega_n^{-(n-1)})^{n-1} \end{bmatrix}$$
證明過程。。兩個矩陣相乘。。我就不搬運了= =。。
這樣的話等式兩邊同時乘上方陣的逆：$$\begin{bmatrix} c_0 \\ c_1 \\ \vdots \\ c_{n-1} \end{bmatrix} = \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^{-(n-1)})^0 & (\omega_n^{-(n-1)})^1 & \cdots & (\omega_n^{-(n-1)})^{n-1} \end{bmatrix} \begin{bmatrix} A(\omega_n^0) \\ A(\omega_n^1) \\ \vdots \\ A(\omega_n^{n-1}) \end{bmatrix}$$
那么可以發現，這可以用$DFT$一樣的過程去求了，不同的是選的點是$\omega_n^{-k}$，最后結果除以$n$。
總之
總的過程大概就是：
- 進行一些預處理，把次數$\times 2$，並拓展成$2$的次冪，高位多出來的設置成$0$；
- 用$FFT$進行$DFT$的過程，在$\Theta (nlogn)$分別將兩個多項式$A(x)$和$B(x)$轉化成點值的形式；
- 利用點值表示在$\Theta (n)$直接就能求得$A(x)\times B(x)$的點值表示；
- 用$FFT$進行$IDFT$的過程，在$\Theta (nlogn)$將乘積結果的點值表示還原成系數表示。
代碼的實現抄的，稍微改了點，寫了一遍：
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
#define MAXN 262144
const double PI=acos(-1.0);
struct Complex{
double real,imag;
Complex(){}
Complex(double _real,double _imag):real(_real),imag(_imag){}
Complex operator-(const Complex &cp) const{
return Complex(real-cp.real,imag-cp.imag);
}
Complex operator+(const Complex &cp) const{
return Complex(real+cp.real,imag+cp.imag);
}
Complex operator*(const Complex &cp) const{
return Complex(real*cp.real-imag*cp.imag,real*cp.imag+imag*cp.real);
}
void setValue(double _real=0,double _imag=0){
real=_real; imag=_imag;
}
};
int len;
Complex wn[MAXN],wn_anti[MAXN];
void FFT(Complex y[],int op){
// Rader, 位逆序置換
for(int i=1,j=len>>1,k; i<len-1; ++i){
if(i<j) swap(y[i],y[j]);
k=len>>1;
while(j>=k){
j-=k;
k>>=1;
}
if(j<k) j+=k;
}
// h=1，Wn^0=1
for(int h=2; h<=len; h<<=1){
// Wn = e^(2*PI*i/n)，如果是插值Wn = e^(-2*PI*i/n)，i虛數單位
Complex Wn = (op==1 ? wn[h] : wn_anti[h]);
for(int i=0; i<len; i+=h){
Complex W(1,0);
for(int j=i; j<i+(h>>1); ++j){
Complex u=y[j],t=W*y[j+(h>>1)];
y[j]=u+t;
y[j+(h>>1)]=u-t;
W=W*Wn;
}
}
}
if(op==-1){
for(int i=0; i<len; ++i) y[i].real/=len;
}
}
void Convolution(Complex A[],Complex B[],int n){
// 初始化
for(len=1; len<(n<<1); len<<=1);
for(int i=n; i<len; ++i){
A[i].setValue();
B[i].setValue();
}
// e^(θi) = cosθ+isinθ -> Wn = cos(2*PI/n)+isin(2*PI/n)
for(int i=0; i<=len; ++i){ // 預處理
wn[i].setValue(cos(2.0*PI/i),sin(2.0*PI/i));
wn_anti[i].setValue(wn[i].real,-wn[i].imag);
}
FFT(A,1); FFT(B,1);
for(int i=0; i<len; ++i){
A[i]=A[i]*B[i];
}
FFT(A,-1);
}