【快速傅里叶变换】FFT
参考：从多项式乘法到快速傅里叶变换 by miskcoo
FFT 学习笔记 by Menci
（一）多项式的表示法
系数表示法：f(x)=a[n-1]*x^(n-1)+...+a[0]，称为n-1次多项式。
点值表示法：一个n-1次多项式在复数域中有n个根，即n个(x,y)可以唯一确定一个n-1次多项式。
对于一个多项式，从其系数表示法到其点值表示法的变换称为离散傅里叶变换（DFT），反之称为傅里叶逆变换（IDFT）。
朴素的离散傅里叶变换，枚举实现的复杂度为O(n^2)。
快速傅里叶变换是指以O(n log n)的复杂度实现IDF和IDFT的算法，常用Cooley-Tukey算法。
（二）复数
复数是形如a+bi的数，当b=0时为实数。
定义一个平面为复平面，那么平面内的每个点(a,b)唯一对应一个复数a+bi，i可以理解为y轴上的单位长度，正如1是x轴上的单位长度。
i的本质是在数轴上定义旋转变换，i是1逆时针旋转90°，那么i^2=-1。
复数相加，遵循平行四边形定则。
复数相乘，模长相乘，幅角相加。
（三）单位根
以圆点为起点，以复平面单位圆的n等分点为终点，作n个向量，设所得幅角为正且最小的向量对应的复数为ω(1,n)，即n次单位根。（括号左为上标，右为下标）。
图片来源：zball
其中B点是单位根ω(1,n)，逆时针依次为ω(2,n)，ω(3,n)...，ω(n,n)=ω(n,0)=1。
计算公式：ω(k,n)=cos ( 2kπ/n ) + i*sin ( 2kπ/n )
单位根的性质：
(1)消去：ω(2n,2k)=ω(n,k)
(2)折半：ω(n,k+n/2)=-ω(n,k)
将ω(n,0)~ω(n,n-1)这n个单位根作为代表n-1次多项式的n个点的横坐标，可以得到很好的性质。
（四）快速傅里叶变换（FFT解决DFA）
这部分因为不会操作数学公式，直接粘贴Menci博客QAQ。
将n-1次多项式A(x)的系数奇偶分成两个多项式A1(x)和A2(x)，则A(x)=A1(x^2)+x*A2(x^2)。
对于k<n/2，有A(ω(n,k))=A1(ω(n/2,k)) + ω(n,k)*A2(ω(n/2,k))
同时，有A(ω(n,k+n/2))=A1(ω(n/2,k)) - ω(n,k)*A2(ω(n/2,k))
对于一个k次多项式，通过奇偶分项得到两个k/2次多项式，分别计算后再调用其值解决k次多项式，即分治解决。
（五）傅里叶逆变换（IDFA）
对于n-1次多项式，其n-1维系数向量{a0,a1...an-1}通过DFA得到点值向量{b0,b1...bn-1}，反之操作称为IDFA。
将点值向量作为系数，以单位根的倒数进行FFT，得到的每个数除以n，就是IDFA的结果。
（六）迭代实现FFT
对于多项式A(x)，已知系数向量a[]，求横坐标为ω(n,0)~ω(n,n-1)的点值向量b[]。
将多项式奇偶分项后，对于k<n/2，有A(ω(n,k))=A1(ω(n/2,k)) + ω(n,k)*A2(ω(n/2,k))，同时有A(ω(n,k+n/2))=A1(ω(n/2,k)) - ω(n,k)*A2(ω(n/2,k))，分治边界是a[i]*ω(0,0)即a[i]。
边界元素：FFT递归边界的数组排布恰好是原数组每个位置二进制反转后的数字，例如：
原：00 01 10 11
终：00 10 01 11
蝴蝶操作：为了在合并时不引入新数组，进行一下操作。ω(l,k)=ω(n,n/l*k)，预处理以n为底的ω[]，IDFT时预处理倒数。
t=ω(n/l*k)*a[l/2+k]
a[l/2+k]=a[k]-t
a[k]+=t
（七）多项式乘法
多项式的点值表示法易于进行乘法，因为对于fc(x)=fa(x)*fb(x)，每个点x在多项式A,B中的点值相乘即可得到在多项式C中的点值。
将n-1次多项式A和m-1次多项式B均视为n+m-2次多项式（高位补0），进行DFT后相乘再通过IDFT即可得到多项式C。
#include<cstdio> #include<cstring> #include<cctype> #include<cmath> #include<queue> #include<stack> #include<set> #include<vector> #include<algorithm> #define ll long long #define lowbit(x) x&-x using namespace std; int read(){ char c;int s=0,t=1; while(!isdigit(c=getchar()))if(c=='-')t=-1; do{s=s*10+c-'0';}while(isdigit(c=getchar())); return s*t; } int min(int a,int b){return a<b?a:b;} int max(int a,int b){return a<b?b:a;} int ab(int x){return x>0?x:-x;} //int MO(int x){return x>=MOD?x-MOD:x;} //void insert(int u,int v){tot++;e[tot].v=v;e[tot].from=first[u];first[u]=tot;} /*------------------------------------------------------------*/ const int inf=0x3f3f3f3f; const int maxn=300010;//2^18!!! const double PI=acos(-1); int n,m; struct cp{ double x,y; cp(double a,double b){x=a;y=b;} cp(){x=y=0;}; cp operator + (cp a){return cp(x+a.x,y+a.y);} cp operator - (cp a){return cp(x-a.x,y-a.y);} cp operator * (cp a){return cp(x*a.x-y*a.y,x*a.y+y*a.x);} }a[maxn],b[maxn]; void fft(cp *a,int n,int f){ int k=0; for(int i=0;i<n;i++){ if(i>k)swap(a[i],a[k]); for(int j=n>>1;(k^=j)<j;j>>=1); } for(int l=2;l<=n;l<<=1){ int m=l/2; cp wn(cos(2*PI*f/l),sin(2*PI*f/l)); for(cp *p=a;p!=a+n;p+=l){ cp w(1,0); for(int i=0;i<l/2;i++){ cp t=w*p[i+m]; p[i+m]=p[i]-t; p[i]=p[i]+t; w=w*wn; } } } if(f==-1){for(int i=0;i<n;i++)a[i].x/=n;}// } int main(){ scanf("%d%d",&n,&m);n++;m++; for(int i=0;i<n;i++){int u;scanf("%d",&u);a[i]=cp(u,0);} for(int i=0;i<m;i++){int u;scanf("%d",&u);b[i]=cp(u,0);} int N=1;while(N<n+m)N*=2; fft(a,N,1);fft(b,N,1); for(int i=0;i<N;i++)a[i]=a[i]*b[i]; fft(a,N,-1); for(int i=0;i<n+m-1;i++)printf("%d ",(int)(a[i].x+0.1)); return 0; }