FFT&DFT簡單總結
前言
相信大家都知道大（chou）名（ming）鼎（zhao）鼎（zhu）的FFT（fake_fake_true）(fast_fast_tle)，並且都有過被它各種玄學操作虐待的經歷（大佬請繞路），那么希望這篇詳細的FFT簡介能夠幫到你。
注：本文中的多項式的次數默認為2的整數次冪，實際中如果不是，需在后面補0。
What it is？
快速傅里葉變換（FFT）是一種能在$O(nlogn)$時間內將一個多項式轉換為它的點值表示的算法。、
什么？你不知道點值表示？好吧我也不知道
划重點
設$A(x)$為一個$n-1$次多項式，我們將$n$個不同的$x$帶入，得到$n$個$y$，這$n$對$(x,y)$就可以確定1個$n-1$次多項式。並不會證明
引入
FFT一般用來加速多項式乘法，也就是處理一個A(x)*B(x)的乘積C(x)。
考慮暴力算法，即$O(n^{2})$枚舉每一項並相乘。這樣做嘛，太慢了！！！
既然我們剛剛學了點值表示，不妨用點值表示試一試？
一番嘗試后我們發現，對於點值表示中的每個$x_{i}$，都有$C(x_{i})=A(x_{i})*B(x_{i})$，而$A(x_{i})$和$B(x_{i})$就是點值表示中的$y_{i}$！如此一來，只要$O(n)$枚舉$x_{i}$就可以算出$C(x)$的點值表示，從而得到答案了！
我們成功的將這個過程提速到了$O(n)$！耶！完結撒花！！！
你以為到此就結束了嗎？並沒有！
我們仔細一想，發現暴力求解$A(x)$和$B(x)$的點值表示居然要$O(n^{2})$！（卒，死因：自閉）
所以，我們還有一個問題：加速求解多項式的點值表示。
FFT：這個我會！！
DFT
誒誒，不是FFT嗎？？沒錯，這是FFT的前身。
在講DFT之前，先講講復數。
復數
復數，可以理解為一個點的坐標，分兩部分，一個是實部，一個是虛部。$a+bi$中，$a$是實部，$b$是虛部($i$是$\sqrt{-1}$)。這里將實部作為橫坐標，虛部作為縱坐標。
但與點不同的是，它還是個數，可以代入，可以運算。
這里特別說一點，復數的乘法：
在幾何意義上，定義為模長相乘，幅角相加。模長是向量的模長，幅角是從x軸正方向逆時針旋轉到與當前向量重合經過的角。
在代數上，
$(a+bi)*(c+di)$
$=ac+bci+adi+bdi^{2}$
$=ac+bci+adi-bd$
$=(ac-bd)+(bc+ad)i$
正文
FFT（DFT）點值表示中用到的$n$個$x$不是隨便找的，而是在單位圓（圓心為原點，半徑為1的圓）上$n$等分后的$n$個等分點所表示的虛數。
從(1,0)開始，將$n$個點從0編號，設第k號節點代表的虛數為$w^{k}_{n}$，即它的實部為$cos \frac{k}{n} 2\pi$，虛部為$sin \frac{k}{n} 2\pi$。根據復數相乘時模長相乘，幅角相加，$w^{k}_{n}$的模長為1，幅角為$\frac{k}{n} 2\pi$，所以$w^{k}_{n}$為$w^{1}_{n}$的$k$次方。我們將這n個特殊的$x$代入，得到了一種特殊的點值表示。
【公式恐懼症患者請速速離開現場】
特殊點性質：
$1. w^{dk}_{dn}=w^{k}_{n}$
(代表的虛數一樣)
$2. w^{k+\frac{n}{2}}_{n}=-w^{k}_{n}$
(對應點關於原點對稱)
不明情況的吃瓜群眾：所以為啥要代入特殊點？
當然是因為代入特殊點后的點值表示有特殊性質啦。
將多項式$A(x)$代入特殊點后的點值表示作為多項式$B(x)$的系數，再將每個特殊點取倒數后代入$B(x)$得到的點值表示的每個數除以$n$，得到的是$A(x)$的系數。
證明嘛，本人由於太弱，並不會證，引用一位大佬博客的證明：
至此，我們學會了DFT（啊哈哈哈哈哈哈哈哈哈猖狂地大笑）
不明情況的吃瓜群眾：可它還是$O(n^{2})$的啊
呃，不然我們要FFT有什么用呢？
FFT
終於，終於到了FFT。FFT其實就是在DFT的基礎上加上了一個分治。
我們先將多項式$A(x)=a_{0}+a_{1}x+\cdots+a_{n-1}x^{n-1}$拆開變成
$A(x)=(a_{0}+a_{2}x^{2}+\cdots+a_{n-2}x^{n-2})+(a_{1}+a_{3}x^{3}+\cdots+a_{n-1}x^{n-1})$
設兩個多項式
$A1(x)=a_{0}+a_{2}x+\cdots+a_{n-2}x^{\frac{n}{2}-1}$
$A2(x)=a_{1}+a_{3}x+\cdots+a_{n-1}x^{\frac{n}{2}-1}$
可得$A(x)=A1(x^{2})+xA2(x^{2})$
代入一個特殊點$w^{k}_{n} (k<\frac{n}{2})$:
$A(w^{k}_{n})$
$=A1(w^{2k}_{n})+w^{k}_{n}A2(w^{2k}_{n})$
$=A1(w^{k}_{\frac{n}{2}})+w^{k}_{n}A2(w^{k}_{\frac{n}{2}})$
再考慮$A(k+\frac{n}{2})$:
$A(w^{k+\frac{n}{2}}_{n})$
$=A1(w^{2k+n}_{n})+w^{k+\frac{n}{2}}_{n}A2(w^{2k+n}_{n})$
$=A1(w^{k+\frac{n}{2}}_{\frac{n}{2}})+w^{k+\frac{n}{2}}_{n}A2(w^{k+\frac{n}{2}}_{\frac{n}{2}})$
$=A1(-w^{k}_{\frac{n}{2}})-w^{k}_{n}A2(-w^{k}_{\frac{n}{2}})$
$=A1(w^{k}_{\frac{n}{2}})-w^{k}_{n}A2(w^{k}_{\frac{n}{2}})$
一番魔改之后，我們終於可以分治了，邊界為$n=1$，遞歸求解$A1(x)$和$A2(x)$。
代碼部分
遞歸
//並不是本人的代碼，因為遞歸實現常數過大，一般不會使用，給一份大佬的代碼（部分）
//本代碼中出現的complex類型是c++STL自帶的，不建議使用，最好自己手打
typedef complex<double> cp
cp omega(int n, int k){
return cp(cos(2 * PI * k / n), sin(2 * PI * k / n));
}
void fft(cp *a, int n, bool inv){
if(n == 1) return;
static cp buf[N];
int m = n / 2;
for(int i = 0; i < m; i++){ //將每一項按照奇偶分為兩組
buf[i] = a[2 * i];
buf[i + m] = a[2 * i + 1];
}
for(int i = 0; i < n; i++)
a[i] = buf[i];
fft(a, m, inv); //遞歸處理兩個子問題
fft(a + m, m, inv);
for(int i = 0; i < m; i++){ //枚舉x，計算A(x)
cp x = omega(n, i);
if(inv) x = conj(x);
//conj是一個自帶的求共軛復數的函數，精度較高。當復數模為1時，共軛復數等於倒數
buf[i] = a[i] + x * a[i + m]; //根據之前推出的結論計算
buf[i + m] = a[i] - x * a[i + m];
}
for(int i = 0; i < n; i++)
a[i] = buf[i];
}
遞推
遞推FFT需要經過仔細觀察才能發現（反正我沒有發現）。
遞歸FFT中，我們將每一項不斷分組遞歸求解子問題，而遞推我們需要預處理出遞歸最后每個初始位置分治后的位置。（有點繞口）
初始位置：0 1 2 3 4 5 6 7
第一輪后：0 2 4 6|1 3 5 7
第二輪后：0 4|2 6|1 5|3 7
第三輪后：0|4|2|6|1|5|3|7
“|”代表分組界限。
結論：一個位置a上的數，最后所在的位置是“a二進制翻轉得到的數”，例如6(011)最后到了3(110)，1(001)最后到了4(100)。
所以，遞推FFT要我們一個位置a上的數，最后所在的位置是“a二進制翻轉得到的數”，例如6(011)最后到了3(110)，1(001)最后到了4(100)。先把每個數放到最后的位置上，然后不斷向上還原，同時求出點值表示。
//洛谷P3803 會T
#include<iostream>
#include<cmath>
using namespace std;
const int N=1e7+5;
const double PI=acos(-1.0);
int lena,lenb,n=1,lim,r[N];
struct cp{//手打加速
double x,y;
cp(double _x=0,double _y=0){
x=_x;y=_y;
}
cp operator*(cp b){
return cp(x*b.x-y*b.y,x*b.y+y*b.x);
}
cp operator+(cp b){
return cp(x+b.x,y+b.y);
}
cp operator-(cp b){
return cp(x-b.x,y-b.y);
}
}a[N],b[N],buf[N];
inline int read(){//本題要卡常
int x=0,f=1;
char ch=getchar();
while(ch<'0'||ch>'9'){
if(ch=='-')f=-1;
ch=getchar();
}
while(ch>='0'&&ch<='9'){
x=x*10+ch-'0';
ch=getchar();
}
return x*f;
}
void FFT(cp *A,int tp){
for(int i=0;i<n;i++)if(i<r[i])swap(A[i],A[r[i]]);
for(int i=1;i<n;i<<=1){
cp W(cos(PI/i),tp*sin(PI/i));//單位根
for(int j=i<<1,k=0;k<n;k+=j){
cp w(1,0);//n等分
for(int l=0;l<i;l++,w=w*W){//buf充當一個緩存數組
buf[k+l]=A[k+l]+w*A[k+i+l];
buf[k+i+l]=A[k+l]-w*A[k+i+l];
}
}
for(int j=0;j<n;j++)A[j]=buf[j];
}
}
int main(){
scanf("%d%d",&lena,&lenb);
while(n<=lena+lenb)n<<=1,lim++;
for(int i=0;i<=lena;i++)a[i].x=read();
for(int i=0;i<=lenb;i++)b[i].x=read();
for(int i=0;i<n;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(lim-1));//預處理，r[i]表示初始位置為i的數最后分治到的位置
FFT(a,1);//計算點值表示
FFT(b,1);
for(int i=0;i<=n;i++)a[i]=a[i]*b[i];
FFT(a,-1);//計算系數
for(int i=0;i<=lena+lenb;i++)printf("%d ",(int)(a[i].x/n+0.5));
}
考慮優化buf數組。
我們發現buf數組存在的意義只是為了更新答案並避免更新時的錯誤而建立的一個緩存數組，所以我們可以將運算內容先記錄下來，舍棄buf數組。
這個優化還有個高大上的名字：蝴蝶操作（雖然我並不知道為啥要這么叫）
//洛谷P3803
#include<iostream> #include<cmath> using namespace std; const int N=1e7+5; const double PI=acos(-1.0); int lena,lenb,n=1,lim,r[N]; struct cp{ double x,y; cp(double _x=0,double _y=0){ x=_x;y=_y; } cp operator*(cp b){ return cp(x*b.x-y*b.y,x*b.y+y*b.x); } cp operator+(cp b){ return cp(x+b.x,y+b.y); } cp operator-(cp b){ return cp(x-b.x,y-b.y); } }a[N],b[N]; inline int read(){ int x=0,f=1; char ch=getchar(); while(ch<'0'||ch>'9'){ if(ch=='-')f=-1; ch=getchar(); } while(ch>='0'&&ch<='9'){ x=x*10+ch-'0'; ch=getchar(); } return x*f; } void FFT(cp *A,int tp){ for(int i=0;i<n;i++)if(i<r[i])swap(A[i],A[r[i]]); for(int i=1;i<n;i<<=1){ cp W(cos(PI/i),tp*sin(PI/i)); for(int j=i<<1,k=0;k<n;k+=j){ cp w(1,0); for(int l=0;l<i;l++,w=w*W){ cp x=A[k+l],y=w*A[k+i+l];//替代buf A[k+l]=x+y; A[k+i+l]=x-y; } } } } int main(){ lena=read();lenb=read(); while(n<=lena+lenb)n<<=1,lim++; for(int i=0;i<=lena;i++)a[i].x=read(); for(int i=0;i<=lenb;i++)b[i].x=read(); for(int i=0;i<n;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(lim-1)); FFT(a,1); FFT(b,1); for(int i=0;i<=n;i++)a[i]=a[i]*b[i]; FFT(a,-1); for(int i=0;i<=lena+lenb;i++)printf("%d ",(int)(a[i].x/n+0.5)); }
例題及題解
https://www.luogu.org/problem/P3803 洛谷P3803 https://www.cnblogs.com/Evan704/p/11406308.html