很详细的FFT（快速傅里叶变换）概念与实现
FFT
首先要说明一个误区，很多人认为FFT只是用来处理多项式乘的，其实FFT是用来实现多项式的系数表示法和点值表示法的快速转换的，所以FFT的用处远不止多项式乘。
FFT的前置知识：点值表示法，复数运算，三角函数。
多项式的系数表示法和点值表示法
系数表示法
点值表示法
不妨将A视为关于x的函数，点值表示法就是在A的图像上取n个点，则该多项式可以被这n个点唯一确定。
点值表示法有什么好处呢？
我们知道系数表示法下两多项式相乘是\(O(n^2)\)，但在点值表示法下奇迹出现了：
显然这个是可以O(n)实现的。虽然但是，我们几乎不会在计算中用到点值表示法，但这也给了我们一个解决多项式乘的思路。系数转点值，相乘，点值转系数。
又很显然，我们可以随便取n个数往函数里带，可惜这样又使复杂度回到了\(O(n^2)\)。
于是FFT出现了，FFT使我们可以用\(O(n\log n)\)的复杂度将系数转换成一组特殊的点值，并再把点值转回系数。
复数的计算
简单的理解，复数就是实数加虚数，多少都知道点虚数吧，没错，知道点就够了。
单位根
记住以下性质（感兴趣可以自己推推，就是基础的三角函数）
好了，现在你已经掌握了所有FFT的前置知识了，自己来推推FFT吧。
正式开始FFT
将\(\omega_n^0,\dots,\omega_n^{n-1}\)这n个数带入得到点值表示，于是：
\[A(x)=a_0+a_1*x+a_2*{x^2}+a_3*{x^3}+a_4*{x^4}+a_5*{x^5}+ \dots+a_{n-2}*x^{n-2}+a_{n-1}*x^{n-1}\\ A(x)=(a_0+a_2*{x^2}+a_4*{x^4}+\dots+a_{n-2}*x^{n-2})+(a_1*x+a_3*{x^3}+a_5*{x^5}+ \dots+a_{n-1}*x^{n-1})\\ A_1(x)=a_0+a_2*{x}+a_4*{x^2}+\dots+a_{n-2}*x^{\frac{n}{2}-1}\\ A_2(x)=a_1+a_3*{x}+a_5*{x^2}+ \dots+a_{n-1}*x^{\frac{n}{2}-1}\\ A(x)=A_1(x^2)+xA_2(x^2)\\ \]我们将\(ωkn(k<n2)ωnk(k<n2)\)代入得
\[A(\omega_n^k)=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})\\ A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}(\omega_n^{2k+n})=A_1(\omega_n^{2k}*\omega_n^n)-\omega_n^kA_2(\omega_n^{2k}*\omega_n^n)=A_1(\omega_n^{2k})-\omega_n^kA_2(\omega_n^{2k}) \]摘自attak的blog
观察这两个式子
显然（怎么又是显然）只要求出\(A_1(\omega_n^{2k})\)和\(\omega_n^kA_2(\omega_n^{2k})\)就可以得出\(A(\omega_n^k)\)和\(A(\omega_n^{k+\frac{n}{2}})\)，而\(A_1(\omega_n^{2k})\)和\(\omega_n^kA_2(\omega_n^{2k})\)又可以进一步递归，正是这一点性质使FFT的复杂度达到了优秀的\(O(n\log n)\)
IFFT
前面提到了，系数转点值，相乘，点值转系数，所以我们还需要能在\(O(n\log n)\)复杂度下完成点值转系数，这就是IFFT，快速傅利叶逆变换。
可以将傅里叶变换和傅立叶逆变换的公式表示写出来（我也不会推逆变换）。
A(x) 表示多项式的系数表示法 B(x)表示多项式的点值表示法
IFFT就是在把FFT的\(\omega_n^{ik}改成\omega_n^{-ik}\),然后再乘个\(\frac 1 N\)
考虑\(\omega_n^k\)的几何意义（高一三角函数）可以得到
所以IFFT只需要在FFT上做一点改动。
千万别忘了最后乘个\(\frac 1 N\)
实现
个人认为这是所有算法最难也最重要的部分，然而很多blog都是将这部分一笔带过，所以我决定来详细的讲讲。
递归写法
竟然是递归，那必然有递归极限。每次递归多项式的项数剩下一半，只剩一项时\(A(x)=a_1\)与\(x\)无关，所以直接返回就好了。
在求出了\(A_1(\omega_n^{2k})\)和\(\omega_n^kA_2(\omega_n^{2k})\)后做两次多项式加减可以得出\(A(\omega_n^k)\)和\(A(\omega_n^{k+\frac{n}{2}})\)
由于多次的递归和数组新建和赋值使递归写法的常数大的出奇，所以我们需要更好的写法。
重难点来了！
递推做法
首先观察这张图。
我们的递归做法就是从上向下将原数列对半拆，再合并。
我们的合并顺序就是从下到上合并，详细的说先合并\(a_0和a_4,a_2和a_6,a1和a_5,a_3和a_7,然后a_0,a_4合并好的整体与a_2和a_6合并好的整体再合并\dots\)
观察最上和最下层的数列的二进制表示，发现就是将二进制翻转了。我们有这个神仙操作可以快速的翻转二进制。
rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1)); ///rev[i]保存i二进制翻转后的数
于是我们可以\(O(n)\)得到最下层的数列，再依次往上推，即不用递归也不用开大量数组，让代码变得飞快！
void fft(cp *a,int n,int inv)
{
int bit=0;
while((1<<bit)<n)bit++;
for(int i=0;i<=n-1;i++)
{
rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1)); ///翻转二进制 3(2)=011 to 6(2)=011
if(i<rev[i]) swap(a[i],a[rev[i]]);
}
for(int mid=1;mid<n;mid*=2) ///mid表示当前合并的行中每一块的长度的一半
{
cp ur(cos(pi/mid),inv*sin(pi/mid)); ///ur代表单位根
for(int i=0;i<n;i+=mid*2) ///i表示当前枚举到这一横行的哪一块的开头下标
{
cp tmp(1,0);
for(int j=0;j<mid;j++,tmp*=ur) ///与递归的写法一样，合并这一块和后面一块。
{
cp x=a[i+j],y=tmp*a[i+j+mid];
a[i+j]=x+y,a[i+j+mid]=x-y;
}
}
}
}
#include<bits/stdc++.h>
using namespace std;
typedef complex<double> cp;
const int N=(1<<21)+10;
const double pi=acos(-1);
cp a[N],b[N];
int n,m,bit,lim,rev[N];
void fft(cp *a,int type)
{
for(int i=0;i<lim;i++)
if(i<rev[i]) swap(a[i],a[rev[i]]);
for(int mid=1;mid<lim;mid*=2)
{
cp ur(cos(pi/mid),type*sin(pi/mid));
for(int i=0;i<lim;i+=mid*2)
{
cp tmp(1,0);
for(int j=0;j<mid;j++,tmp*=ur)
{
cp x=a[i+j],y=tmp*a[i+j+mid];
a[i+j]=x+y,a[i+j+mid]=x-y;
}
}
}
}
int main()
{
cin>>n>>m;
for(int i=0;i<=n;i++)
cin>>a[i];
for(int i=0;i<=m;i++)
cin>>b[i];
for(lim=1;lim<=n+m;lim<<=1)
bit++;
for(int i=0;i<lim;i++)
rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
fft(a,1),fft(b,1);
for(int i=0;i<=lim;i++)
a[i]*=b[i];
fft(a,-1);
for(int i=0;i<=n+m;i++)
cout<<(int)(0.5+a[i].real()/lim)<<' ';
return 0;
}
