0XFF---FFT是啥?
FFT是一種DFT的高效算法，稱為快速傅立葉變換（fast Fourier transform），它根據離散傅氏變換的奇、偶、虛、實等 特性，對離散傅立葉變換的算法進行改進獲得的。 ---百度百科
對於兩個多項式 \(F(x)\) 和 \(G(x)\) ，要求你將他們乘起來。
那還不簡單？直接暴力相乘啊：
設 \(F(x)\) 的系數數列為 \(C\)。
\(F(x) \times G(x) = C_nx^nG(x) + C_{n-1}x^{n-1}G(x) + C_{n-2}x^{n-2}G(x) \cdots C_2x^2G(x) + C_1x^1G(x) + C_0G(x)\)
這樣下來需要做做 \(n\) 次單項式乘多項式，每次的時間復雜度 \(O(n)\) ，則總復雜度高達 \(O(n^2)\)
基本上 \(n\) 上了\(4000\) 就會被卡吧......那怎么提速呢？
這就需要我們偉大而又神奇的神器：\(FFT\) (快速博立葉變換)
復雜度就只有 \(O(nlogn)\) 了。
0X1F---FFT的前置知識.
1.復數是什么?
我們把形如 \(z=a+bi\)（ \(a,b\) 均為實數）的數稱為復數，其中 \(a\) 稱為實部， \(b\) 稱為虛部， \(i\) 稱為虛數單位。當虛部等於零時，這個復數可以視為實數；當z的虛部不等於零時，實部等於零時，常稱z為純虛數。復數域是實數域的代數閉包，即任何復系數多項式在復數域中總有根。 復數是由意大利米蘭學者卡當在十六世紀首次引入，經過達朗貝爾、棣莫弗、歐拉、高斯等人 的工作，此概念逐漸為數學家所接受。 ---百度百科
想必大家都知道實數是啥(不知道重讀幼兒園吧......)，實數位於數軸上，就像下圖這樣：
我們稍微觀察一下，\(1\) 是怎么變到 \(-1\) 的呢？
在數軸上轉了 \(180°\)。
如果，是 \(90°\) 的話，會發生什么呢？
這個時候，會轉到 \(0\) 上面的位置，但是那里，好像沒有數啊!
不對，其實是有的，只不過這個數不在實數軸上，而是在虛數軸上！
虛數軸的單位是 \(i\) ，我們可以這么表示：
嗯，對。這顯然是一個平面坐標系。現在我們的數僅限於數軸上，如果是這個平面坐標系上的一個點怎么表達呢？
對於下面的紅色點：
這個點的坐標很容易的可以得到：\((2,i)\) ，也可以表示成 \(2+i\) .
你沒猜錯！這個就叫復數！
一個很重要的結論：復數相乘時，模長相乘，幅角相加！
2.點值表示法是什么?
我們用一個二維平面坐標系，在上面畫 \(N+1\) 個點，最終可以解出一個 \(n\) 元的函數。證明略。
同樣，我們可以用 \(N-1\) 個點來表達一個多項式。
因為點值相乘的復雜度只有 \(O(n)\) 顯然優秀許多。
3.單位根是什么?
*n次單位根(n為正整數)是n次冪為1的復數!
*n次單位根(n為正整數)是n次冪為1的復數!
*n次單位根(n為正整數)是n次冪為1的復數!
我們先在復平面上畫個點，就像這樣：
Ta叫做單位圓。
圓邊上的任意一點的模長都是 \(1\).
只有單位圓上的點表示的復數才有可能成為\(n\)次單位根!
單位根的基本符號：\(ω\)
一個單位圓，我們將它切成 \(n\) 份，從 \((1,0)\) 開始旋轉，每次旋轉 \(\frac{1}{n} \times 360\) 度，每次旋轉后的點都記為 \(ω_{n}^{k}\)，特別的，\(ω_{n}^{0}\) 和 \(ω_{n}^{n}\) 都是 \((1,0)\) 點。
還有，當 \(k>=n\) 或者 \(k<0\) 時，\(ω_{n}^{k}\) 也是合法的。
單位根的性質：
\(1.\) 對於任意的 \(n\) , \(ω_{n}^{0}\) 都為 \((1,0)\) 點。
\(2.\) $ω_{n}^{a} \times ω_{n}^{b} = ω_{n}^{a+b} $
\(3.\) $ω_{an}^{ak} = ω_{n}^{k} $
\(4.\) $(ω_{n}^{x})^y = (ω_{n}^{y})^x $
\(5.\) $ω_{n}^{k+n/2} = -ω_{n}^{k} $ if(n%2==0)
0X2F---FFT的求解過程.
- 分治思想很重要！
我們將多項式 \(F(x)\) 按位置分成兩塊。
那么變成了(保證n是2的正整數次冪)：
\(F(x) = (C_0+C_2x^2+C_4x^4+ \cdots +C_{n-2}x^{n-2}) + (C_1x+C_3x^3+C_5x^5+ \cdots +C_{n-1}x^{n-1})\)
設兩個多項式 \(F1(x),F2(x)\)。
\(F1(x) = C_0+C_2x+C_4x^2+ \cdots +C_{n-2}x^{n/2-1}\)
\(F2(x) = C_1x+C_3x+C_5x^2+ \cdots +C_{n-1}x^{n/2-1}\)
則我們可以得出：
\(F(x) = F1(x^2) + F2(x^2) \times x\)
設 \(k<n/2\) ， 將 \(ω_{n}^{k}\) 帶入多項式 \(F(x)\).
\(F(ω_{n}^{k}) = F1((ω_{n}^{k})2) + F2((ω_{n}^{k})^2) \times ω_{n}^{k}\)
簡化得： \(F(ω_{n}^{k}) = F1(ω_{n/2}^{k}) + F2(ω_{n/2}^{k}) \times ω_{n}^{k}\)
再假設 \(k<n/2\) ，將 \(ω_{n}^{k+n/2}\) 帶入多項式 \(F(x)\).
\(F(ω_{n}^{k+n/2}) = F1((ω_{n}^{k+n/2})2) + F2((ω_{n}^{k+n/2})^2) \times ω_{n}^{k}\)
\(F(ω_{n}^{k+n/2}) = F1(ω_{n}^{2k+n}) + F2(ω_{n}^{2k+n}) \times ω_{n}^{k+n/2}\)
\(F(ω_{n}^{k+n/2}) = F1(ω_{n}^{2k}) + F2(ω_{n}^{2k}) \times ω_{n}^{k+n/2}\)
\(F(ω_{n}^{k+n/2}) = F1(ω_{n/2}^{k}) + F2(ω_{n/2}^{k}) \times ω_{n}^{k+n/2}\)
\(F(ω_{n}^{k+n/2}) = F1(ω_{n/2}^{k}) - F2(ω_{n/2}^{k}) \times ω_{n}^{k}\)
比較一下兩個式子：
- \(F(ω_{n}^{k}) = F1(ω_{n/2}^{k}) + F2(ω_{n/2}^{k}) \times ω_{n}^{k}\)
- \(F(ω_{n}^{k+n/2}) = F1(ω_{n/2}^{k}) - F2(ω_{n/2}^{k}) \times ω_{n}^{k}\)
等式右邊只有一個負號的差別!
這兩個式子很關鍵！
0X3F---FFT的代碼實現.
對於復數的使用
雖然 \(C++ STL\) 里面有復數 \((complex)\) 但是太慢不建議大家使用。
你可以自己手打 \(complex\)
- 手打的 \(complex\) :
struct complex{complex(double a=0,double b=0){x=a,y=b;}double x,y;};
complex operator + (complex a,complex b){return complex(a.x+b.x,a.y+b.y);}
complex operator - (complex a,complex b){return complex(a.x-b.x,a.y-b.y);}
complex operator * (complex a,complex b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
- FFT:
complex a[N],b[N];
inline void FFT(complex *f,int len,short inv){
if(!len)return;complex f1[len+1],f2[len+2];
for(int k=0;k<len;++k)f1[k]=f[k<<1],f2[k]=f[k<<1|1];//按位置分
FFT(f1,len>>1,inv);FFT(f2,len>>1,inv);//遞歸處理子問題
complex tmp=complex(cos(PI/len),inv*sin(PI/len)),buf=complex(1,0);
/*tmp:做一次平方后坐標的變換*/ /*buf:初始位置*/
for(RI k=0;k<len;++k){
complex t=buf*f2[k];
f[k]=f1[k]+t,f[k+len]=f1[k]-t;buf=buf*tmp;//按照公式還原
}return;
}
//注意，inv的作用是判斷是 "系數轉點值" 還是 "點值轉系數"
\(Code\) 中提到的公式是這兩項：
- \(F(ω_{n}^{k}) = F1(ω_{n/2}^{k}) + F2(ω_{n/2}^{k}) \times ω_{n}^{k}\)
- \(F(ω_{n}^{k+n/2}) = F1(ω_{n/2}^{k}) - F2(ω_{n/2}^{k}) \times ω_{n}^{k}\)
對於文中的"坐標的變換"：
我們依舊來看單位圓：
實際上，這個坐標的變換，直接用園中的三角形，運用三角函數就可以得出解了。
過程略.
最后我們得到的結果是：\(ω_{n}^{1} = (cos(\frac{2π}{n}),sin(\frac{2π}{n}))\)
求出 \(ω_{n}^{1}\) 后將它乘 \(n\) 次，可以得到：\({ω_{n}^{0},ω_{n}^{1},ω_{n}^{2},ω_{n}^{3}，ω_{n}^{4}，ω_{n}^{5} \cdots ω_{n}^{n-1}}\)
貼出最終的代碼：
#include<cstdio>
#include<cmath>
#include<string>
#define ll long long
#define RI register int
#define inf 0x3f3f3f3f
#define PI 3.1415926535898
using namespace std;
const int N=6e4+2;
template <typename _Tp> inline _Tp max(const _Tp&x,const _Tp&y){return x>y?x:y;}
template <typename _Tp> inline _Tp min(const _Tp&x,const _Tp&y){return x<y?x:y;}
template <typename _Tp> inline void IN(_Tp&x){
char ch;bool flag=0;x=0;
while(ch=getchar(),!isdigit(ch))if(ch=='-')flag=1;
while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
if(flag)x=-x;
}
struct complex{complex(double a=0,double b=0){x=a,y=b;}double x,y;};
complex operator + (complex a,complex b){return complex(a.x+b.x,a.y+b.y);}
complex operator - (complex a,complex b){return complex(a.x-b.x,a.y-b.y);}
complex operator * (complex a,complex b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
complex a[N],b[N];
inline void FFT(complex *f,int len,short inv){
if(!len)return;complex f1[len+1],f2[len+2];
for(int k=0;k<len;++k)f1[k]=f[k<<1],f2[k]=f[k<<1|1];
FFT(f1,len>>1,inv);FFT(f2,len>>1,inv);
complex tmp=complex(cos(PI/len),inv*sin(PI/len)),buf=complex(1,0);
for(RI k=0;k<len;++k){
complex t=buf*f2[k];
f[k]=f1[k]+t,f[k+len]=f1[k]-t;buf=buf*tmp;
}return;
}
int n,m;
int main(){
scanf("%d%d",&n,&m);
for(RI i=0;i<=n;++i)scanf("%lf",&a[i].x);
for(RI i=0;i<=m;++i)scanf("%lf",&b[i].x);
for(m+=n,n=1;n<=m;n<<=1);
FFT(a,n>>1,1);FFT(b,n>>1,1);
for(int i=0;i<n;++i)a[i]=a[i]*b[i];
FFT(a,n>>1,-1);
for(int i=0;i<=m;++i)printf("%.0f ",fabs(a[i].x)/n);
putchar('\n');
return 0;
}
聽說可以優化，那啥的我還不會，就到這吧.
過了一會兒......
"原來FFT小優化這么簡單啊!"
0X4F---FFT的一些小優化.
不用遞歸：
遞歸版(數組下標，先偶后奇，從0開始)：
0 1 2 3 4 5 6 7 --第1層
0 2 4 6 |1 3 5 7 --第2層
0 4 |2 6 |1 5 |3 7 --第3層
0 |4 |2 |6 |1 |5 |3| 7 --第4層
發現了什么嗎？
最后的序列是原序列的二進制反轉!
比如： \(6 = (110)_2\) 反過來變成了 \((011)_2 = 3\) ！
如何得到二進制翻轉后的數列？遞推即可!
for(RI i=0;i<n;++i)filp[i]=(filp[i>>1]>>1)|((i&1)?n>>1:0);
//filp[i] 即為 i 的二進制位翻轉
Code:
#include<cstdio>
#include<cmath>
#include<string>
#define ll long long
#define RI register int
#define inf 0x3f3f3f3f
#define PI 3.1415926535898
using namespace std;
const int N=3e6+2;
int n,m,filp[N];
template <typename _Tp> inline _Tp max(const _Tp&x,const _Tp&y){return x>y?x:y;}
template <typename _Tp> inline _Tp min(const _Tp&x,const _Tp&y){return x<y?x:y;}
template <typename _Tp> inline void IN(_Tp&x){
char ch;bool flag=0;x=0;
while(ch=getchar(),!isdigit(ch))if(ch=='-')flag=1;
while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
if(flag)x=-x;
}
struct complex{complex(double a=0,double b=0){x=a,y=b;}double x,y;};
complex operator + (complex a,complex b){return complex(a.x+b.x,a.y+b.y);}
complex operator - (complex a,complex b){return complex(a.x-b.x,a.y-b.y);}
complex operator * (complex a,complex b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
complex a[N],b[N];
inline void FFT(complex *f,short inv){
for(RI i=0;i<n;++i)if(i<filp[i]){complex tmp=f[i];f[i]=f[filp[i]];f[filp[i]]=tmp;}
/*換位置*/
for(RI p=2;p<=n;p<<=1){//每局區間長度
RI len=p/2;//合並子區間的長度(所以是p/2)
complex tmp=complex(cos(PI/len),inv*sin(PI/len));
for(RI k=0;k<n;k+=p){//每局左端點
complex buf=complex(1,0);
for(RI l=k;l<k+len;++l){//遍歷區間
complex t=buf*f[len+l];
f[len+l]=f[l]-t,f[l]=f[l]+t,buf=buf*tmp;//賦值有微小的變化，注意順序！
}
}
}return;
}
/*主程序不變*/
int main(){
scanf("%d%d",&n,&m);
for(RI i=0;i<=n;++i)scanf("%lf",&a[i].x);
for(RI i=0;i<=m;++i)scanf("%lf",&b[i].x);
for(m+=n,n=1;n<=m;n<<=1);
for(RI i=0;i<n;++i)filp[i]=(filp[i>>1]>>1)|((i&1)?n>>1:0);
FFT(a,1);FFT(b,1);
for(RI i=0;i<n;++i)a[i]=a[i]*b[i];
FFT(a,-1);
for(RI i=0;i<=m;++i)printf("%.0f ",fabs(a[i].x)/n);
putchar('\n');
return 0;
}
luogu上的題，遞歸的總是T最后一個點，改成非遞歸版的就A了？
emmmmmmmmmmmmmm
所有優化全開：
很作死，建議不要輕易嘗試[滑稽]
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#pragma GCC optimize("Ofast")
#pragma GCC optimize("inline")
#pragma GCC optimize("-fgcse")
#pragma GCC optimize("-fgcse-lm")
#pragma GCC optimize("-fipa-sra")
#pragma GCC optimize("-ftree-pre")
#pragma GCC optimize("-ftree-vrp")
#pragma GCC optimize("-fpeephole2")
#pragma GCC optimize("-ffast-math")
#pragma GCC optimize("-fsched-spec")
#pragma GCC optimize("unroll-loops")
#pragma GCC optimize("-falign-jumps")
#pragma GCC optimize("-falign-loops")
#pragma GCC optimize("-falign-labels")
#pragma GCC optimize("-fdevirtualize")
#pragma GCC optimize("-fcaller-saves")
#pragma GCC optimize("-fcrossjumping")
#pragma GCC optimize("-fthread-jumps")
#pragma GCC optimize("-funroll-loops")
#pragma GCC optimize("-fwhole-program")
#pragma GCC optimize("-freorder-blocks")
#pragma GCC optimize("-fschedule-insns")
#pragma GCC optimize("inline-functions")
#pragma GCC optimize("-ftree-tail-merge")
#pragma GCC optimize("-fschedule-insns2")
#pragma GCC optimize("-fstrict-aliasing")
#pragma GCC optimize("-fstrict-overflow")
#pragma GCC optimize("-falign-functions")
#pragma GCC optimize("-fcse-skip-blocks")
#pragma GCC optimize("-fcse-follow-jumps")
#pragma GCC optimize("-fsched-interblock")
#pragma GCC optimize("-fpartial-inlining")
#pragma GCC optimize("no-stack-protector")
#pragma GCC optimize("-freorder-functions")
#pragma GCC optimize("-findirect-inlining")
#pragma GCC optimize("-fhoist-adjacent-loads")
#pragma GCC optimize("-frerun-cse-after-loop")
#pragma GCC optimize("inline-small-functions")
#pragma GCC optimize("-finline-small-functions")
#pragma GCC optimize("-ftree-switch-conversion")
#pragma GCC optimize("-foptimize-sibling-calls")
#pragma GCC optimize("-fexpensive-optimizations")
#pragma GCC optimize("-funsafe-loop-optimizations")
#pragma GCC optimize("inline-functions-called-once")
#pragma GCC optimize("-fdelete-null-pointer-checks")
[滑稽][滑稽][滑稽][滑稽][滑稽][滑稽][滑稽][滑稽][滑稽][滑稽][滑稽]
最后，因為本人實在太弱了，太蒟了，所以實在寫不出啥了。
\(by Qiuly\)