快速傅里叶变换（FFT）
写在前面的话：
感觉最近每次考试都是同样的一种情形：考试的时候绞尽脑汁地去思考正解无果，在Solution发下来以后，却发现这是一种我没有学习过的算法......并且感觉这些东西貌似就我一个人不会了，赶紧来学习了一下，首先，便是FFT
（一）预备知识
要学习FFT，首先得学习一些关于复数的东西
这一篇写得不错，理解起来很容易，具体看here
（二）开始学习
附上我在学习FFT的时候用的3篇博客，很多内容都是一样的，但是中间有些知识可以互补
首先学习博客1，这篇博客很多东西的证明写得很清楚，有利于初学者
然后再看一看博客2，这篇博客中有两副图有助于理解，并且也详细解释了使用矩阵来计算IDFT的过程
全部理解了以后，自然是要写代码的呀，看一看博客3，这篇代码写得还是蛮详细的。
（三）重要知识
为了方便我以后复习，先将3篇博客的重要部分浓缩至此，略去所有证明以及基础概念性的东西，为什么？因为我懒得打字啊~~~
（1）多项式乘法
定义$ A(x)=\sum_{i=0}^{n-1} a_{i}x^{i} $ 与 \(B(x)=\sum_{i=0}^{n-1} b_{i}x^{i}\) 相乘的结果为\(C(x)\)，那么\(C(x)\)的值为
$$C(x)=\sum_{k=0}^{2n-2}(\sum_{k=i+j} a_{i}b_{j})x^{k}$$
这样做的时间复杂度无疑是\(O(n^{2})\)，显然不行
但是如果我们知道了两个多项式在2n-1个点处所取得的点值表示，即将那2n-1个点作为自变量x带入A与B之后所得的结果，那么我们就可以在\(O(n)\)时间之内计算出\(C(x)\)
那么我们该如何做到这一点呢？这就是为什么要用FFT的原因了。
（2）单位根
在复平面上，以原点为圆心，1为半径作圆，所得的圆叫做单位圆。以原点为起点，单位圆的n等分点为终点，作n个向量。设所得的幅角为正且最小的向量对应的复数为\(\omega_{n}\)，称为n次单位根。
在代数中，若\(z^{n}=1\)，我们把z称为n次单位根
性质一 ：\(\omega_{2n}^{2k}=\omega_{n}^{k}\)
性质二 ：\(\omega_{n}^{k+ \frac{n}{2}}=-\omega_{n}^{k}\)
（3）快速傅里叶变换
考虑多项式\(A(x)\)的表示。将n次单位根的0到n-1次幂带入多项式的系数表示，所得点值向量\((A(\omega_{n}^{0}),...,A(\omega_{n}^{n-1}))\)称为其系数向量的离散傅里叶变换。
Cooley-Tukey 算法是FFT的最常用算法，具体见3个博客
（4）傅里叶逆变换
将点值表示的多项式转化为系数表示，同样可以使用快速傅里叶变换，这个过程叫做傅里叶逆变换。
设 \((y_{0},...,y{n-1})\) 为 \((a_{0},...,a{n-1})\) 的傅里叶变换。考虑另一个向量 \((c_{0},...,c_{n-1})\) 满足
即多项式 \(B(x)=y_{0}+y_{1}x+y_{2}x^{2}+...+y_{n-1}x^{n-1}\) 在 \(\omega_{n}^{0},\omega_{n}^{-1},...,\omega_{n}^{-(n-1)}\)处的点值表示。
。。。。。。（此处省略一大堆证明）
然后，我们可以得到：使用单位根的倒数代替单位根，做一次类似快速傅里叶变换的过程，再将结果每个数除以 n，即为傅里叶逆变换的结果。
（5）流程图
总得来说，FFT的程序流程图长这样
图还是蛮好的，很形象
（6）迭代实现
图像更清晰
但是细心的话可以发现，其实这个图右下角有些问题，画反了，反正不是我画的，所以呢，这锅我就不背了，哈哈。
（7）蝴蝶操作
懒得写了，看博客1，写得很详细
（8）代码
哈哈，终于可以上代码了！
原题，洛谷P3803
#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
#include<cmath>
typedef long long ll;
typedef double dd;
#define For(i,j,k) for (int i=j;i<=k;++i)
#define Forr(i,j,k) for (int i=j;i>=k;--i)
#define Set(a,p) memset(a,p,sizeof(a))
using namespace std;
template<typename T>bool chkmax(T& a,T b) {return a<b?a=b,1:0;}
template<typename T>bool chkmin(T& a,T b) {return a>b?a=b,1:0;}
const int maxn=5e6+1e2;
const dd Pi=acos(-1.0);
struct Complex {
dd x,y;
}a[maxn],b[maxn];
int n,m,N,cnt;
int p[maxn];
inline int read() {
int x=0,p=1;
char c=getchar();
while (!isdigit(c)) {if (c=='-') p=-1; c=getchar();}
while (isdigit(c)) {x=(x<<1)+(x<<3)+(c-'0'); c=getchar();}
return x*p;
}
Complex operator + (const Complex &aa,const Complex &bb) {
return (Complex){aa.x+bb.x,aa.y+bb.y};
}
Complex operator - (const Complex &aa,const Complex &bb) {
return (Complex){aa.x-bb.x,aa.y-bb.y};
}
Complex operator * (const Complex &aa,const Complex &bb) {
return (Complex){aa.x*bb.x-aa.y*bb.y,aa.x*bb.y+aa.y*bb.x};
}
inline void FFT(Complex *s,int type) {
For (i,0,N-1)
if (i<p[i]) swap(s[i],s[p[i]]);//得到新序列
for (int mid=1;mid<N;mid<<=1) { //枚举合并区间的中点
Complex Wn=(Complex){cos(Pi/mid),type*sin(Pi/mid)}; //单位根
for (int len=mid<<1,j=0;j<N;j+=len) {
//len是区间的长度，j是区间的左端点
Complex w=(Complex){1,0}; //幂
for (int k=0;k<mid;++k,w=w*Wn) {
//枚举左区间
Complex x=s[j+k],y=w*s[j+mid+k];
s[j+k]=x+y; s[j+mid+k]=x-y;
//蝴蝶操作
}
}
}
}
int main() {
n=read(); m=read();
For (i,0,n) a[i].x=read();
For (i,0,m) b[i].x=read();
for (N=1;N<=n+m;N<<=1,++cnt) ;
For (i,0,N-1)
p[i]=(p[i>>1]>>1) | ((i&1)<<(cnt-1));
/* 这里就是翻转，对于原序列的第i个数，将i在二进制下翻转之后就是原来元素在新序列
中的位置。这里预处理，因为要翻转i，其实就是先将i/2翻转之后，再看i的奇偶，判
断最高位。*/
FFT(a,1); FFT(b,1);
For (i,0,N) a[i]=a[i]*b[i];
FFT(a,-1);
For (i,0,n+m) printf("%d ",(int)(a[i].x/N+0.5));
return 0;
}
（8）注意
\(\Pi\) 是doule型的，不要习惯性的定义为int
