一、引入
首先，定義多項式的形式為 \(f(x)=\sum_{i=0}^n a_ix^i\)，其中 \(a_i\) 為系數，\(n\) 為次數，這種表示方法稱為“系數表示法”，一個多項式是由其系數確定的。
可以證明，\(n+1\) 個點可以唯一確定一個 \(n\) 次多項式。對於 \(f(x)\)，代入 \(n+1\) 個不同的 \(x\)，得到 \(n+1\) 個不同的 \(y\)。一個 \(n\) 次的多項式就可以等價地換成 \(n+1\) 個等式，相當於平面上的 \(n+1\) 組坐標 \((x_i,y_i)\)，這種表示方法稱為“點值表示法”。
多項式乘法 (卷積)：設 \(C(x)=A(x)\cdot B(x)\)，\(A,B,C\) 的系數構成的數列分別為 \(a,b,c\)，則 \(c_k=\sum_{i=0}^ka_ib_{k-i}\)。理解：因為 \(x^i\times x^{k-i}=x^k\)，\(a_i\) 和 \(b_{k-i}\) 相乘后，它們后面的未知數就變成了 \(x^k\)，對 \(c_k\) 產生貢獻。
暴力求解：兩個 \(n\) 次多項式相乘，時間復雜度 \(\mathcal{O}(n^2)\)。
快速傅里葉變換（Fast Fourier Transform，簡稱 FFT）可以 \(\mathcal{O}(n\log n)\) 求解。
二、基本步驟
利用點值表示法，分三步快速求出多項式乘積：
- 由系數表示法轉換成點值表示法。（DFT）
- 利用點值表示法，求兩個多項式的乘積。
- 再將點值表示法轉化成系數表示法。（IDFT）
對於步驟二，給出兩個 \(n\) 次多項式 \(A\) 和 \(B\) 的點值表達式，我們可以 \(\mathcal{O}(n)\) 求出其乘積 \(C\) 的點值表達式。顯然 \(C\) 是 \(2n\) 次的，取 \(2n+1\) 個不同的 \(x_i\)。
- 代入 \(A\)：\(\{(x_0,y_0),(x_1,y_1),\cdots,(x_{2n},y_{2n})\}\)。
- 代入 \(B\)：\(\{(x_0,y'_0),(x_1,y'_1),\cdots,(x_{2n},y'_{2n})\}\)。
- 則 \(C\) 的點值表達式為：\(\{(x_0,y_0y'_0),(x_1,y_1y'_1),\cdots,(x_{2n},y_{2n}y'_{2n})\}\)。
下面重點分析步驟一。有一些前置概念。
三、復數
我們把形如 \(z=a+bi\)（\(a,b\) 為實數）的數稱為 復數，其中，\(a,b\) 分別叫做復數 \(z\) 的實部與虛部，\(i\) 為虛數單位，\(i^2=-1\)。
我們把復數表示在 復平面 上（\(x\) 軸叫實軸，\(y\) 軸叫虛軸），就像把實數表示在數軸上。如圖，復數 \(z=a+bi\) 與復平面內的點 \(Z(a,b)\) 一一對應。
連接 \(OZ\)。復數 \(z=a+bi\) 與平面向量 \(\vec{OZ}\) 一一對應。
復數的四則運算：
- 加減法：\((a+bi)\pm (c+di)=(a\pm c)+(b\pm d)i\)。
- 乘法：\((a+bi)(c+di)=ac+adi+bci+bdi^2=(ac-bd)+(bc+ad)i\)。
- 除法：\(\frac{a+bi}{c+di}=\frac{(a+bi)\times (c-di)}{(c+di)\times(c-di)}=\frac{(ac+bd)+(bc-ad)i}{c^2+d^2}=\frac{ac+bd}{c^2+d^2}+\frac{bc-ad}{c^2+d^2}i\)。
complex<double>a; //STL。a.real() 返回復數 a 的實部，a.imag() 返回復數 a 的虛部 struct cp{ //手寫，比 STL 快一點 double a,b; cp operator+(cp &x){return (cp){a+x.a,b+x.b};} cp operator-(cp &x){return (cp){a-x.a,b-x.b};} cp operator*(cp &x){return (cp){a*x.a-b*x.b,a*x.b+b*x.a};} cp operator/(cp &x){double v=x.a*x.a+x.b*x.b; return (cp){(a*x.a+b*x.b)/v,(b*x.a-a*x.b)/v};} //除法一般不用 qwq };
我們稱 \(\overline{z}=a-bi\) 為復數 \(z=a+bi\) 的 共軛復數。
對於復數 \(z=a+bi\)，模長 \(|z|\) 為點 \(Z(a,b)\) 到原點的距離，幅角 為對應的向量 \(\vec{OZ}\) 與橫軸正半軸的夾角。
復數的乘法可以表達為：模長相乘，輻角相加。
四、單位根
1. 定義
\(n\) 次單位根是滿足 \(x^n=1\) 的復數 \(x\)。
首先，單位根的模長必然為 \(1\)。因為若 \(|x|>1\)，則 \(|x^n|=|x|^n>1\)；若 \(|x|<1\)，則 \(|x^n|=|x|^n<1\)。
所以單位根表示的點一定在 單位圓（圓心為原點，半徑為 \(1\)）上。
其次，單位根的輻角 \(\theta\) 一定滿足 \(\frac{n\theta}{2\pi}\in\mathbb{Z}\)。也就是一個向量從 \((1,0)\) 開始，每次旋轉 \(\theta\) 的角度，旋轉 \(n\) 次后還落在 \((1,0)\) 上，那么它一定旋轉了整數圈。
然后發現 \(n\) 次單位根正好是模長為 \(1\)，輻角為 \(\frac{2k\pi}{n}\) 的向量對應的復數。
記模長為 \(1\)，輻角為 \(\frac{2k\pi}{n}\) 的向量對應的 \(n\) 次單位根為 \(\omega_n^k\)，稱為第 \(k\) 個 \(n\) 次單位根。
還能發現，\(\omega_n^k=\omega_n^{k\bmod n}\)，所以一般情況下，我們認為 \(n\) 次單位根有 \(n\) 個，即 \(\omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1}\)。
2. 性質
單位根的性質：（這些性質在后文會被用到）
-
性質 1：\(n\) 次單位根對應的向量將單位圓 \(n\) 等分。
兩個相鄰的 \(n\) 次單位根對應的向量的夾角相等。單位根的輻角是周角的 \(\frac{1}{n}\)。 -
性質 2：\(\omega_n^k=\omega_n^{k\bmod n}\)。
在弧度制下，任意弧度 \(\theta\) 與 \(\theta+2k\pi\,(k\in\mathbb{Z})\) 表示相同的角。 -
性質 3：\({(\omega_n^k)}^p=\omega_n^{kp}\)。
第 \(k\) 個 \(n\) 次單位根對應向量的輻角變為原來的 \(p\) 倍，相當於 \(\omega_n^{kp}\) 對應的向量。 -
性質 4：\(\omega_{dn}^{dk}=\omega_{n}^k\)。
考慮兩者對應向量的輻角，\(\frac{2dk\pi}{dn}=\frac{2k\pi}{n}\)。也可以這樣理解：\(dn\) 次單位根對應的向量將單位圓 \(dn\) 等分，取第 \(dk\) 個。\(n\) 次單位根對應的向量將單位圓 \(n\) 等分，取第 \(k\) 個。兩者等價。 -
性質 5：\(\omega_n^{k+n/2}=-\omega_n^k\)。其中 \(n\) 為偶數。
相當於一個復數對應的向量進行一次中心對稱，\(a+bi\) 變為 \(-a-bi\)。
根據性質 3 有，\({(\omega_n^k)}^2=\omega_{n}^{2k}\)。根據性質 4 有，\(\omega_n^{2k}=\omega_{n/2}^k\)，其中 \(n\) 為偶數。
3. 求法
根據性質 3，有 \(\omega_n^k=(\omega_n^1)^k\)。也就是說，只要求出 \(\omega_n^1\)，就能得到 \(\omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1}\)。
\(\omega_n^1\) 所對應的向量模長為 \(1\)，輻角為 \(\frac{2\pi}{n}\)，得到 \(\omega_n^1\) 所對應的點為 \((\cos(\frac{2\pi}{n}),\sin(\frac{2\pi}{n}))\)。
求 \(\pi\)：double pi=acos(-1)。
補充：\(\omega_n^k=e^{\frac{2\pi ik}{n}}=\cos(\frac{2\pi k}{n})+i\cdot \sin(\frac{2\pi k}{n})\)。其中 \(i\) 為虛數單位。
五、DFT
將系數表示法轉換成點值表示法。
1. 基本思路
對於 \(n-1\) 次多項式（也就是有 \(n\) 項） \(f(x)=\sum_{i=0}^{n-1} a_ix^i\)，將奇偶次數分離。
\(f(x)=(a_0+a_2x^2+a_4x^4+\cdots+a_{n-2}x^{n-2})+(a_1x+a_3x^3+a_5x^5+\cdots+a_{n-1}x^{n-1})\)
定義兩個新的多項式 \(f_1(x)\) 和 \(f_2(x)\)：
- \(f_1(x)=a_0+a_2x+a_4x^2+\cdots+a_{n-2}x^{{n/2-1}}\)。
- \(f_2(x)=a_1+a_3x+a_5x^2+\cdots+a_{n-1}x^{n/2-1}\)。
於是有 \(f(x)=f_1(x^2)+xf_2(x^2)\)。
將 \(\omega_n^k\,(k<\frac{n}{2})\) 代入得：\(f(\omega_n^k)=f_1(\omega_n^{2k})+\omega_n^kf_2(\omega_n^{2k})=f_1(\omega_{n/2}^k)+\omega_n^kf_2(\omega_{n/2}^k)\)。
同理，將 \(\omega_n^{k+n/2}\,(k<\frac{n}{2})\) 代入得：\(f(\omega_n^{k+n/2})=f_1(\omega_n^{2k+n})+\omega_n^{k+n/2}f_2(\omega_n^{2k+n})\)
\(=f_1(\omega_n^{2k})+\omega_n^{k+n/2}f_2(\omega_n^{2k})=f_1(\omega_{n/2}^k)+\omega_n^{k+n/2}f_2(\omega_{n/2}^k)=f_1(\omega_{n/2}^k)-\omega_n^kf_2(\omega_{n/2}^k)\)。
發現兩者的右邊只有正負號的區別。
第一個式子的 \(k\) 取遍 \([0,\frac{n}{2}-1]\) 時，\(k+\frac{n}{2}\) 取遍 \([\frac{n}{2},n-1]\)。
如果我們知道 \(f_1(x),f_2(x)\) 分別在 \(x=\omega_{n/2}^0,\omega_{n/2}^1,\cdots,\omega_{n/2}^{n/2-1}\) 的點值表示，就可以 \(\mathcal{O}(n)\) 求出 \(f(x)\) 在 \(x=\omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1}\) 的點值表示。
然后，發現 \(f_1(x),f_2(x)\) 和 \(f(x)\) 的性質完全相同，這樣就把問題分成了兩個子問題，對於這兩個子問題再進行遞歸求解。這樣就可以在 \(\mathcal{O}(n\log n)\) 的時間復雜度內求出點值表達式。
2. 代碼實現
DFT 是利用單位根的特殊性質進行分治。
考慮到它能處理的多項式長度只能為 \(2^k\)（\(k\) 為整數），否則在分治時左右的項數就會不同。我們可以在最高次補一些系數為 \(0\) 的項，把原來的 \(n\) 補到 \(2^k\,(2^k\geq n)\)。這樣不影響計算結果。
別忘了：\(f(\omega_n^k)=f_1(\omega_{n/2}^k)+\omega_n^kf_2(\omega_{n/2}^k),f(\omega_n^{k+n/2})=f_1(\omega_{n/2}^k)-\omega_n^kf_2(\omega_{n/2}^k)\)。
遞歸實現 DFT：
#include<bits/stdc++.h> #define int long long using namespace std; const int N=1e6+5; int n,len; double x,pi=acos(-1); complex<double>a[N]; void DFT(complex<double>*a,int n){ if(n==1) return ; //邊界 int m=n/2; complex<double>a1[m],a2[m]; for(int i=0;i<m;i++) a1[i]=a[i<<1],a2[i]=a[i<<1|1]; //按奇偶分類 DFT(a1,m),DFT(a2,m); //處理子問題 complex<double>x(cos(2*pi/n),sin(2*pi/n)),w(1,0); //x: w_n^1。w: 當前的 n 次單位根，初始值為 w_n^0，即 1。 for(int i=0;i<m;i++) a[i]=a1[i]+w*a2[i],a[i+m]=a1[i]-w*a2[i],w*=x; //w*=x: 得到下一個單位根。 } signed main(){ scanf("%lld",&n); for(int i=0;i<n;i++) scanf("%lf",&x),a[i]=x; //a[i] 的實部賦值為 x for(len=1;len<n;len<<=1); //把項數補到 2 的整冪，高次項的系數默認為 0 DFT(a,len); for(int i=0;i<len;i++) printf("(%.4lf,%.4lf)\n",a[i].real(),a[i].imag()); return 0; }
六、IDFT
DFT 的逆運算。將點值表示法轉化成系數表示法。
結論：把 DFT 中的 \(\omega_n^1\) 換成它的共軛復數，即 \((\cos(\frac{2\pi}{n}),-\sin(\frac{2\pi}{n}))\)，得到的系數再除以 \(n\) 即可。證明略。
可以把 IDFT 和 DFT 放在一起寫。
void FFT(complex<double>*a,int n,int opt){ //opt=1 為 DFT，opt=-1 為 IDFT if(n==1) return ; int m=n/2; complex<double>a1[m],a2[m]; for(int i=0;i<m;i++) a1[i]=a[i<<1],a2[i]=a[i<<1|1]; FFT(a1,m,opt),FFT(a2,m,opt); complex<double>x(cos(2*pi/n),sin(2*pi/n)*opt),w(1,0); for(int i=0;i<m;i++) a[i]=a1[i]+w*a2[i],a[i+m]=a1[i]-w*a2[i],w*=x; }
七、迭代實現
目前的代碼如下。可以發現它的效率並不是很高。
//Luogu P3803 #include<bits/stdc++.h> using namespace std; const int N=3e6+5; int n,m,len; double x,pi=acos(-1); complex<double>a[N],b[N]; void FFT(complex<double>*a,int n,int opt){ //opt=1/-1: DFT/IDFT if(n==1) return ; int m=n/2; complex<double>a1[m],a2[m]; for(int i=0;i<m;i++) a1[i]=a[i<<1],a2[i]=a[i<<1|1]; FFT(a1,m,opt),FFT(a2,m,opt); complex<double>x(cos(2*pi/n),sin(2*pi/n)*opt),w(1,0); for(int i=0;i<m;i++) a[i]=a1[i]+w*a2[i],a[i+m]=a1[i]-w*a2[i],w*=x; //蝴蝶操作（只是一個名字 qwq）。這里 w*a2[i] 算了兩次，先記錄下來再算可以減小常數。 } signed main(){ scanf("%d%d",&n,&m); for(int i=0;i<=n;i++) scanf("%lf",&x),a[i]=x; for(int i=0;i<=m;i++) scanf("%lf",&x),b[i]=x; n=n+m+1; //一個 n 次多項式和一個 m 次多項式的乘積是一個 n+m 次多項式 (有 n+m+1 項) for(len=1;len<n;len<<=1); FFT(a,len,1),FFT(b,len,1); for(int i=0;i<len;i++) a[i]*=b[i]; //點值直接乘 FFT(a,len,-1); for(int i=0;i<n;i++) printf("%d%c",(int)(a[i].real()/len+0.5),i==n-1?'\n':' '); //注意這里是除以 len return 0; }
如圖，考慮遞歸的結構：
求出第 \(4\) 層的狀態，就能往上合並求出前 \(3\) 層的狀態。
觀察發現，最后的數組下標的序列是原序列的 二進制翻轉。比如 \(6=(110)_2\)，反過來就是 \((011)_2=3\)。而第 \(4\) 層 \(a_3\) 的位置就是原來 \(a_6\) 的位置。
//r[i] 表示 i 二進制翻轉后的結果。求 0~len-1 在二進制位數為 log2(len)-1 意義下的二進制翻轉。 for(int i=0;i<len;i++) //len 是 2 的冪次 r[i]=(r[i>>1]>>1)|((i&1)?len>>1:0);
理解：考慮 \(i\) 與 \(\frac{i}{2}\) 在二進制下的關系。\(i\) 可以看作是 \(\frac{i}{2}\) 在二進制下的每一位左移一位得到。翻轉后，\(i\) 是 \(\frac{i}{2}\) 在二進制下的每一位右移一位得到，然后判一下最后一位即可。
迭代實現：
//Luogu P3803 #include<bits/stdc++.h> using namespace std; const int N=3e6+5; int n,m,len,r[N]; double x,pi=acos(-1); complex<double>a[N],b[N]; void FFT(complex<double>*a,int n,int opt){ //opt=1/-1: DFT/IDFT for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]); //求出最后一層的序列 for(int k=2;k<=n;k<<=1){ //枚舉區間長度 int m=k>>1; //待合並的長度 complex<double>x(cos(2*pi/k),sin(2*pi/k)*opt),w(1,0),v; for(int i=0;i<n;i+=k,w=1) //枚舉起始點 for(int j=i;j<i+m;j++) //遍歷區間 v=w*a[j+m],a[j+m]=a[j]-v,a[j]=a[j]+v,w*=x; //蝴蝶操作。注意先 a[j+m]=a[j]-v 再 a[j]=a[j]+v。 } } signed main(){ scanf("%d%d",&n,&m); for(int i=0;i<=n;i++) scanf("%lf",&x),a[i]=x; for(int i=0;i<=m;i++) scanf("%lf",&x),b[i]=x; n=n+m+1; for(len=1;len<n;len<<=1); for(int i=0;i<len;i++) //二進制翻轉 r[i]=(r[i>>1]>>1)|((i&1)?len>>1:0); FFT(a,len,1),FFT(b,len,1); for(int i=0;i<len;i++) a[i]*=b[i]; FFT(a,len,-1); for(int i=0;i<n;i++) printf("%d%c",(int)(a[i].real()/len+0.5),i==n-1?'\n':' '); return 0; }
八、模板
以 P3803 【模板】多項式乘法（FFT） 為例。
遞歸實現：
#include<bits/stdc++.h> using namespace std; const int N=3e6+5; int n,m,len; double x,pi=acos(-1); complex<double>a[N],b[N]; void FFT(complex<double>*a,int n,int opt){ //opt=1/-1: DFT/IDFT if(n==1) return ; int m=n/2; complex<double>a1[m],a2[m]; for(int i=0;i<m;i++) a1[i]=a[i<<1],a2[i]=a[i<<1|1]; FFT(a1,m,opt),FFT(a2,m,opt); complex<double>x(cos(2*pi/n),sin(2*pi/n)*opt),w(1,0); for(int i=0;i<m;i++) a[i]=a1[i]+w*a2[i],a[i+m]=a1[i]-w*a2[i],w*=x; } signed main(){ scanf("%d%d",&n,&m); for(int i=0;i<=n;i++) scanf("%lf",&x),a[i]=x; for(int i=0;i<=m;i++) scanf("%lf",&x),b[i]=x; n=n+m+1; for(len=1;len<n;len<<=1); FFT(a,len,1),FFT(b,len,1); for(int i=0;i<len;i++) a[i]*=b[i]; FFT(a,len,-1); for(int i=0;i<n;i++) printf("%d%c",(int)(a[i].real()/len+0.5),i==n-1?'\n':' '); return 0; }
迭代實現：（比遞歸快）
#include<bits/stdc++.h> using namespace std; const int N=3e6+5; int n,m,len,r[N]; double x,pi=acos(-1); complex<double>a[N],b[N]; void FFT(complex<double>*a,int n,int opt){ //opt=1/-1: DFT/IDFT for(int i=0;i<n;i++) if(i<r[i]) swap(a[i],a[r[i]]); for(int k=2;k<=n;k<<=1){ int m=k>>1; complex<double>x(cos(2*pi/k),sin(2*pi/k)*opt),w(1,0),v; for(int i=0;i<n;i+=k,w=1) for(int j=i;j<i+m;j++) v=w*a[j+m],a[j+m]=a[j]-v,a[j]=a[j]+v,w*=x; } } signed main(){ scanf("%d%d",&n,&m); for(int i=0;i<=n;i++) scanf("%lf",&x),a[i]=x; for(int i=0;i<=m;i++) scanf("%lf",&x),b[i]=x; n=n+m+1; for(len=1;len<n;len<<=1); for(int i=0;i<len;i++) r[i]=(r[i>>1]>>1)|((i&1)?len>>1:0); FFT(a,len,1),FFT(b,len,1); for(int i=0;i<len;i++) a[i]*=b[i]; FFT(a,len,-1); for(int i=0;i<n;i++) printf("%d%c",(int)(a[i].real()/len+0.5),i==n-1?'\n':' '); return 0; }
記憶：
- \(f(x)=f_1(x^2)+xf_2(x^2)\)。
- \(f(\omega_n^k)=f_1(\omega_{n/2}^k)+\omega_n^kf_2(\omega_{n/2}^k)\)。
- \(f(\omega_n^{k+n/2})=f_1(\omega_{n/2}^k)-\omega_n^kf_2(\omega_{n/2}^k)\)。