知識點：FFT詳解
目錄
- 前言
- 前置知識
- 知識點講解
- 概要
- 多項式相乘的樸素算法
- 系數表示法與點值表示法
- 復數的引入
- 單位復根
- 有關定理的證明
- DFT
- DFT的優化
- IDFT
- AC代碼（luogu3803）
前言
FFT其實在很早的時候就已經接觸到了，但是那個時候學起來有點仙，感覺這東西離實際解題的距離有點遠，不如那些其他的數據結構那麽直接。但是半年多下來的做題，發現FFT其實應用的十分廣泛，並且很多數學題推出公式之後就可以套用FFT進行計算。所以對於FFT的理解也不能僅僅只是停留於背板子的階段了，而應該更加深入的去理解它。
前置知識
復數的運算及性質，多項式
知識點講解
概要
首先講什麽是FFT，FFT的全稱為快速傅裏葉變換（Fast Fourier Transformation），是基於離散傅裏葉變換的快速求法。最基礎的運用就是解決多項式相乘的問題，可以將樸素算法的\(O(n^2)\)
多項式相乘的樸素算法
我們假設我們現在有兩個多項式：
那麽我們令這兩個的多項式相乘為\(h(x)\)，即：
那麽我們可以得到：
所以我們就可以得到一個復雜度為\(O(n^2)\)的解法了，即枚舉兩個多項式的每一位的系數，那麽第一個多項式第\(i\)
系數表示法與點值表示法
首先，我們需要知道，多項式的一個表示方法為系數表示法，即一個n次的多項式可以表示為\(\sum_{i=0}^na_ix^i\)。這裏的每一個\(a_i\)表示每一項的系數。然後我們發現，對於一個多項式，我們只會關心這個多項式的系數，而並不需要真正的記錄下它的指數，因為在系數表示法中，\(a_i\)的指數就是\(i\)。所以一個多項式可以表示成這個樣子：
然後我們就可以引出點值表示法了。點值表示法顧名思義就是用幾個點來表示一個多項式。兩點確定一條直線，三點可以確定一條拋物線，所以一個\(n\)
然後多項式的乘積就可以這樣表示：
這樣子我們就把一個多項式轉化成了一些離散的點，這樣的過程就叫做離散傅裏葉變換（DFT）。
然後將一些離散的點重新轉化成一個多項式，這個過程就叫做離散傅裏葉反變換（IDFT）。
這樣子我們對於FFT就會有一個大致的思路了：先將原來的兩個多項式進行DFT，轉化成點值之後再進行相乘，最後做IDFT重新變成答案的多項式。這樣在相乘過程中的復雜度就可以被優化成\(O(n)\)了。但是由於我們在相乘時用的是點值表達式，而解題的時候一般不會給你點值表達式，所以我們還需要將系數表示轉化為點值表示，於是下面開始進入DFT與IDFT。
復數的引入
我們在一般的計算當中都不會對\(\sqrt{-1}\)進行定義，然而在復數中，\(\sqrt{-1}\)等於一個神奇的數:\(i\)，這個數在復數的定義下相當於\(1\)的作用。下面列舉一些有關\(i\)的計算：
復數分為實部和虛部，所以一個復數可以表示為\(a+bi\)，這裏的\(a\)表示的是這個復數的實部，\(bi\)表示的是這個復數的虛部。同時，一個復數也是可以用坐標系來表示的，只不過這個坐標系\(y\)軸的單位不再是1，而是\(i\)，這樣我們就可以引入極角與復數的乘法了。
一個復數一共有三種表示方法（在坐標系中）：\(a+bi\)，\((a,b)\)，\((r,\theta)\)（這裏的r表示到圓心的距離，\(\theta\)表示極角）。於是復數的乘法也可以定義了：
通過這個，我們對於另一種表示法也可以定義出它的乘法了：
對於一個復數\(a+bi\)，\(r=\sqrt{a^2+b^2}，\theta=arctan(\frac{b}{a})\)
復數\(a+bi\)與\(c+di\)相乘，\(r_1=\sqrt{a^2+b^2}\)，\(r_2=\sqrt{c^2+d^2}\)，\(\theta _1=arctan(\frac{b}{a})\)，\(\theta _2=arctan(\frac{d}{c})\)
那麽乘積為\((ac-bd)+(ad+bc)i\)
\(R=\sqrt{(ac-bd)^2+(ad+bc)^2}=\sqrt{a^2c^2+b^2d^2-2abcd+b^2c^2+a^2d^2+2abcd}=\sqrt{(a^2+b^2)(b^2+d^2)}=r1\times r2\)
\(\Theta=arctan(\frac{bc+ad}{ac-bd})=arctan(\frac{\frac{bc+ad}{ac}}{\frac{ac-bd}{ac}})=arctan(\frac{\frac{b}{a}+\frac{d}{c}}{1-\frac{bd}{ac}})=arctan(\frac{\frac{b}{a}+\frac{d}{c}}{1-\frac{b}{a}\cdot \frac{d}{c}})=arctan(\frac{b}{a})+arctan(\frac{d}{c})=\theta_1+\theta_2\)
所以\((r1,\theta1)\times(r2,\theta_2)=(r1\times r2,\theta_1+\theta_2)\)，總結下來就是兩個復數相乘，長度相乘，極角相加。由於c++自帶的\(complex\)函數有點菜菜的，所以我們可以自己手寫\(complex\)函數。然後我們可以想象一下，如果兩個復數到原點的距離都為1的話，那麽這兩個復數相乘就相當於繞著以原點為圓心的，半徑為1的圓進行旋轉。
struct comp {
double r,i;
comp() {}
comp(double r,double i):r(r),i(i) {}
};
comp operator + (comp a,comp b) {
return comp(a.r+b.r,a.i+b.i);
}
comp operator - (comp a,comp b) {
return comp(a.r-b.r,a.i-b.i);
}
comp operator * (comp a,comp b) {
return comp(a.r*b.r-a.i*b.i,a.i*b.r+a.r*b.i);
}
單位復根
我們現在已經知道了答案多項式的點值表達式，我們現在所需要做的任務就是把這個點值表達式轉變成系數表達式。最明顯的一種方法就是用高斯消元法，暴力去解除這個方程的解，可是這樣的復雜度又會變成\(O(n^2)\)了，因為我們在計算\(x_0,x_0^2,x_0^3\cdots x_0^n\)的時候會重復計算很多次，而這個計算在實數域中似乎是不可避免的，所以我們就可以回到復數域中。在復數中，我們需要的數是\(\omega^k=1\)的數，而由於復數的乘法的定義，我們可以發現，對於所有的\((1,\theta)\)的復數，他的\(i\)次方總是在同一個圓上，並且經過一定次數的乘方之後一定是可以等於1的，所以我們就可以有效的利用這種數，將這個數帶入，就可以很妙妙的求值了。
我們定義這種\(\omega^k=1\)的數為\(k\)次單位復根，計作\(\omega_k^n\)，而這個\(n\)實際上是一個序號，表示的是將所有的\(k\)次單位根按極角序進行排序之後從零開始編號的單位根。
而由於\(\omega_k^0=1\)，所以我們只要知道\(\omega_k^1\)，就可以計算出所有的\(k\)次單位根了。
有關定理的證明
有關單位復根有很多比較有用的性質，這裏對這些性質進行介紹並證明。
基本性質
既然單位復根相乘滿足“長度相乘，極角相加”，而單位復根的長度又是1，所以我們可以發現所有的\(k\)次單位根是均勻的排布在一個半徑為1，以原點為圓心的圓上的。根據歐拉公式\(e^{\pi i}=-1\)，而\(\omega_n^n=1\)，所以\(\omega_n^n=(-1)^2=(e^{\pi i})^2=e^{2\pi i}\)，所以：
這樣我們就可以快速的算出我們需要的\(k\)次單位根了。
消去引理
先給出引理：
這個其實很好理解的，就居\(\omega_8^2\)來說吧，根據引理，我們可以得出\(\omega_4^1\)，這就相當於我們原本把一個圓劃分成了8份，然後取的是其中的2份，而這又與將圓劃分成4份，取其中的1份是等價的，由此可以發現這個消去引理是正確的。
由此我們可以推出\(\omega_{2n}^n=\omega_2^1=-1\)。
折半引理
引理：
這個定理看起來是比較的玄，首先前半部分是毋庸置疑的，因為\(\omega_n^k=-\omega_n^{k+\frac{1}{2}n}\)，所以這兩個的平方也是相同的。然後考慮後半部分，我們還是以\(n=8\)為例，然後畫出圓被分成八份的圖像：
然後我們可以想象一下，如果我們將任意一個單位根進行平方，那麽相當於他的極角擴大了一倍，而擴大後的單位根一定落在4次單位根上，並且單位根的編號與原來單位根的編號是相同的（這裏只考慮\(k\leq \frac{1}{2}n\)）。
求和引理
引理：
這個定理看起來很高妙，仔細看會發現這是一個等比數列求和，於是我們套用等比數列求和公式，那麽對於第一個公式：
\(S=\frac{1-(\omega_n^k)^n}{1-\omega_n^k}=\frac{1-\omega_n^{nk}}{1-\omega_n^k}=\frac{1-1^k}{1-\omega_n^k}=0\)
對於第二個公式，就更加簡單了：
\(S=\sum_{j=0}^{n-1}(\omega_m^{dn})^j=\sum_{j=0}^{n-1}1^j=n\)
這樣這兩個公式就推完了。
DFT
首先我們先回到原本的系數表示法：
將我們的單位復根帶入之後，就會得到：
然後我們得到的\(\lbrace y_0,y_1,\cdots \rbrace\)就是\(\lbrace a_0,a_1,\cdots \rbrace\)的點值了。
現在知道了求DFT的方法了，我們就要開始想能夠讓DFT的復雜度在\(O(nlogn)\)的做法了。
這裏用到的是分治的思想，分治計算所有出\(k\)次單位根所對應的點值，而由於我們一共需要\(n+1\)個點值，所以對於一個\(n\)次的多項式，我們需要計算的就是它的\(n\)次單位根對應的點值。分治的思想主要體現在將多項式分為奇數項與偶數項，那麽一個多項式可以表示為：
令\(G(x)=(a_0+a_2x+a_4x^2+a_6x^3)\)，\(H(x)=(a_1+a_3x+a_5x^2+a_7x^3)\)
那麽\(f(x)=G(x^2)\times +x\times H(x^2)\)
當\(x=\omega^k\)時，\(DFT(f(x))_k=DFT(G(x^2)\times +\omega^k \times H(x^2))\)
而由於折半引理，我麽可以得出：
\(DFT(f(x+\frac{n}{2}))_k=DFT(G(x^2)\times(-\omega^k)\times H(x^2))\)
所以我們發現，對於一個多項式，我們只要求出了前一半部分的值，我們就可以推出後一半的值，這樣我們利用了分治在DFT上做到了復雜度為\(O(nlogn)\)了。但是這樣也是有一個問題的，就是這樣子只能處理出多項式系數的個數為2的整數次冪的多項式，所以我們在做DFT的時候，還需要將不滿2的整數次冪的多項式在前面補零，這樣就能做到不影響結果並且正確計算出答案。
DFT的優化
在這個分治上，我們采用的是遞歸的形式，然而由於每一次遞歸都要下傳數組，導致這樣的常數非常大，空間也有了一些無意義上的消耗，所以我們希望的是能夠采用叠代的形式。然後我們開始思考，這個分治的實際結果是什麽。我們每次將數組中的數按奇偶分組，實際上就是對於二進制下的這一位進行1與0的分組。並且我們可以觀察一下最後分組的結果，以長度為8為例：
或許這樣還看不出來什麽，但是我們將每個元素的下表用二進制表示之後，就可以發現一些有趣的事了：
我們可以發現，每一位的在分治之後的序號正好與分治之前的序號在二進制表示上是翻轉過來的。於是我們就可以考慮預處理出每一位在分治之後的序號，這樣就不用進行遞歸的過程，而是用叠代來代替了。對於預處理這個序號，我們可以類似DP一樣的進行：
void GetRev() {
for(int i=0;i<lim;i++) {
rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
}
}這樣就得出了我們在翻轉之後，每一位的序號了，接下來就是簡單的叠代了。
IDFT
處理到這裏，我們的FFT就只剩下最後的一步了，就是用IDFT把點值表達式轉化成系數表達式。首先我們先考慮我們是如何從系數表達式變成點值表達式的。我們是將單位復根帶入原來的系數表達式來計算點值的，我們可以將單位復根帶入之後的n個多項式的系數用一個矩陣來表示：
\[ \left[ \begin{matrix} 1 & 1 & 1 & \cdots & 1 \ 1 & \omega_n^1 & \omega_n^2 & \cdots & \omega_n^{n-1} \ 1 & \omega_n^2 & \omega_n^4 & \cdots & \omega_n^{2(n-1)} \ \cdots & \cdots & \cdots & \cdots & \cdots \ 1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \cdots & \omega_n^{(n-1)^2}\\end{matrix} \right] \]
現在我們想要讓這所有的多項式的系數還原成用\(a\)表示的多項式，我們就必須在這個矩陣之後乘上一個這個矩陣的逆矩陣。首先我們要先了解逆矩陣的定義。一個矩陣的逆矩陣\(V\)定義：
其中\(I\)為單位矩陣，單位矩陣的樣子是這樣的：
\[
\left[
\begin{matrix}
1 & 0 & 0 & \cdots & 0 \ 0 & 1 & 0 & \cdots & 0 \ \cdots & \cdots & \cdots & \cdots & \cdots \ 0 & 0 & 0 & \cdots & 1 \\end{matrix}
\right]
\]
也就是說這個單位矩陣的主對角線都是1，其余全是0。
在這裏我先給出這個矩陣的逆矩陣長什麽樣：
\[
\left[
\begin{matrix}
\frac{1}{n} & \frac{1}{n} & \frac{1}{n} & \cdots & \frac{1}{n} \ \frac{1}{n} & \frac{\omega_n^{-1}}{n} & \frac{\omega_n^{-2}}{n} & \cdots & \frac{\omega_n^{-(n-1)}}{n} \ \frac{1}{n} & \frac{\omega_n^{-2}}{n} & \frac{\omega_n^{-4}}{n} & \cdots & \frac{\omega_n^{-2(n-1)}}{n} \ \cdots & \cdots & \cdots & \cdots & \cdots \ \frac{1}{n} & \frac{\omega_n^{-(n-1)}}{n} & \frac{\omega_n^{-2(n-1)}}{n} & \cdots & \frac{\omega_n^{-(n-1)^2}}{n}\\end{matrix}
\right]
\]
然後我們開始證明這個矩陣是原矩陣的逆矩陣：
首先我們可以知道原矩陣位於\((j,k)\)位置的值為\(\omega_n^{jk}\)，而位於逆矩陣中的\((j,k)\)處的值為\(\omega_n^{-jk}\)，於是我們可以計算出\([V\cdot V^{-1}]\)的第\((j_1,j_2)\)位置上的值為：
然後根據求和引理，我們發現只有當\(j_1-j_2=0\)時，\([V\cdot V^{-1}]_{j_1j_2}=1\)，其他時候\([V\cdot V^{-1}]_{j_1j_2}=0\)，這樣就可以證明上圖的矩陣就是我們要乘的逆矩陣了。
我們將矩陣帶入公式，就可以得到\(a_j=\frac{1}{n}\sum_{k=0}^{n-1}y_k\omega_n^{-jk}\)（這裏的\(y_k\)是DFT中的\(y_k\)）。然後我們就會驚奇的發現IDFT與DFT在公式上實際上是差不多的：
所以我麽完全可以用同一個函數傳不同的參數進去進行計算：
void FFT(comp *a,int IDFT) {//IDFT傳進來時為-1，DFT傳進來時為1
for(int i=0;i<lim;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
for(int mid=1;mid<lim;mid<<=1) {
comp w=comp(cos(PI/mid),IDFT*sin(PI/mid));//IDFT是單位復根應該取負，而DFT時單位復根應該取正的。
for(int l=mid<<1,j=0;j<lim;j+=l) {
comp wn=comp(1.0,0.0);
for(int k=0;k<mid;k++) {
comp x=a[k+j];
comp y=a[k+j+mid]*wn;
a[k+j]=x+y;
a[k+j+mid]=x-y;
wn=wn*w;
}
}
}//這裏沒有除以n是因為在最後除了
}到這裏FFT就已經講的差不多了，如果仍然沒有聽明白的話，那麽可以移步這裏。
AC代碼（luogu3803）
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
bool Finish_read;
template<class T>inline void read(T &x){Finish_read=0;x=0;int f=1;char ch=getchar();while(!isdigit(ch)){if(ch==‘-‘)f=-1;if(ch==EOF)return;ch=getchar();}while(isdigit(ch))x=x*10+ch-‘0‘,ch=getchar();x*=f;Finish_read=1;}
template<class T>inline void print(T x){if(x/10!=0)print(x/10);putchar(x%10+‘0‘);}
template<class T>inline void writeln(T x){if(x<0)putchar(‘-‘);x=abs(x);print(x);putchar(‘\n‘);}
template<class T>inline void write(T x){if(x<0)putchar(‘-‘);x=abs(x);print(x);}
/*================Header Template==============*/
const int maxn=5e6+500;
const double PI=acos(-1);
int n,m;
int rev[maxn];
int lim=1,len;
/*==================Define Area================*/
struct comp {
double r,i;
comp() {}
comp(double r,double i):r(r),i(i) {}
}a[maxn],b[maxn];
comp operator + (comp a,comp b) {
return comp(a.r+b.r,a.i+b.i);
}
comp operator - (comp a,comp b) {
return comp(a.r-b.r,a.i-b.i);
}
comp operator * (comp a,comp b) {
return comp(a.r*b.r-a.i*b.i,a.i*b.r+a.r*b.i);
}
void GetRev() {
for(int i=0;i<lim;i++) {
rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
}
}
void FFT(comp *a,int IDFT) {
for(int i=0;i<lim;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
for(int mid=1;mid<lim;mid<<=1) {
comp w=comp(cos(PI/mid),IDFT*sin(PI/mid));
for(int l=mid<<1,j=0;j<lim;j+=l) {
comp wn=comp(1.0,0.0);
for(int k=0;k<mid;k++) {
comp x=a[k+j];
comp y=a[k+j+mid]*wn;
a[k+j]=x+y;
a[k+j+mid]=x-y;
wn=wn*w;
}
}
}
}
int main() {
read(n);read(m);
while(lim<=n+m) lim<<=1,len++;
GetRev();
for(int i=0;i<=n;i++) scanf("%lf",&a[i].r);
for(int i=0;i<=m;i++) scanf("%lf",&b[i].r);
FFT(a,1);
FFT(b,1);
for(int i=0;i<lim;i++) {
a[i]=a[i]*b[i];
}
FFT(a,-1);
for(int i=0;i<=m+n;i++) {
printf("%d ",(int)(a[i].r/lim+0.5));
}
return 0;
}
知識點：FFT詳解