前置芝士：
欧拉说：$e^{i\theta}=cos(\theta)+isin(\theta)$
定义单位根$\omega_{n}{k}=e{n}}$
记$\theta=2\pi/n$ ,则$\omega_{n}{k}=e=cos(k\theta)+isin(k\theta)$
$1.\omega_{n}{0}=e=1,\omega_{n}{n}=e=cos(2\pi)+isin(2\pi)=1$
$2.(\omega_{n}{k})2=(e^{i 2\pi \frac{k}{n}})2=e{n/2}}=\omega_{n/2}^{k}$
$3.\omega_{n}{k}=cos(k\theta)+isin(k\theta)=-cos(k\theta+\pi)-isin(k\theta+\pi)=-e=-e^{i(\frac{2\pi k}{n}+\pi)}= -e{i2\pi(\frac{k+n/2}{n})}=-\omega_{n}$
$4.\omega_{n}{a}*\omega_{n}=e^{i2\pi \frac{a+b}{n}}=\omega_{n}^{a+b}$
基本思路
对于一个最高次为$n$次的多项式，我们可以通过$n+1$组$(x_i,F(x_i))$来唯一确定这个多项式的系数，这种表示方式叫做点值表示法。
$FFT$：把一个多项式从系数表示转化为点值表示法
$IFFT$：把一个多项式从点值表示法转为系数表示法
为求解$G(x)=H(x)F(x)$ 的系数，我们的思路是先把$H(x)$和$F(x)$通过$FFT$转为点值表示，分别为$(x_i,H(x_i))和(x_i,F(x_i))$，则$(x_i,H(x_i)F(x_i))$为$G(x)$的点值表示，再通过$IFFT$求出$G(x)$ 的系数
FFT
设$F(x)=a_0+a_1x+a_2x2...+a_{n-1}x$
我们可以利用单位根的一些性质进行快速计算，不妨将点值表示定为$(\omega_{n}{0}，F(\omega_{n})),(\omega_{n}{1}，F(\omega_{n})),...(\omega_{n}{n-1}，F(\omega_{n}))$
问题转化为快速求解这$n$个点值
将$F(x)$的奇偶项分开（我们暂把$n$当作偶数）
$F(x)=(a_0+a_2x2+...+a_{n-2}x)+(a_1x+a_3x3+...+a_{n-1}x)$
$F(x)=(a_0+a_2x2+...+a_{n-2}x)+x(a_1+a_3x2+...+a_{n-1}x)$
记为$F(x)=A_0(x2)+xA_1(x2)$
代入 $x=\omega_{n}^{k}$
① $k<=n/2$
$F(\omega_{n}{k})=A_0(\omega_{n/2})+\omega_{n}{k}A_1(\omega_{n/2})$
②$k>n/2$
记$k'=k-n/2$
$F(\omega_{n}{k})=A_0((-\omega_{n})2)+\omega_{n}A_1((-\omega_{n}{k'})2)$
$F(\omega_{n}{k})=A_0(\omega_{n/2})-\omega_{n}{k'}A_1(\omega_{n/2})$
那么对于$k<=n/2$
$F(\omega_{n}{k})=A_0(\omega_{n/2})+\omega_{n}{k}A_1(\omega_{n/2})$
$F(\omega_{n}{k+n/2})=A_0(\omega_{n/2})-\omega_{n}{k}A_1(\omega_{n/2})$
发现求出$k<=n/2$部分的所有$A_0(x)$和$A_1(x)$，则可同时得到$k>n/2$部分的$F(x)$
那么可以分治：
我们最开始要求解$n$个长度为$n$的多项式$F(x)$
而我们只要求$\frac{n}{2}$ 个长度为$\frac{n}{2}$ 的多项式$A_0(x)$和另$\frac{n}{2}$个长度为$\frac{n}{2}$的$A_1(x)$
......以此类推
最终只要求解共$n$ 个长度为 1 的多项式即可，长度为1意味着最高次是0，直接返回系数即可
则每一层的计算复杂度为O(n)，共log(n)层，复杂度为O(nlogn)
为了避免一些不必要的讨论，且为了后续的一个优化，我们始终希望计算的每一个多项式的长度$n$是一个偶数，因此我们常把多项式的长度填充到$2^k$次为止，这也是前面直接把$n$当作偶数的原因。
IFFT
假设$G(x)=b_0+b_1x+b_2x2...+b_{n-1}x$
而我们现在已经求出了$G(x)$的点值表示$(x_0,y_0),(x_1,y_1),...,(x_{n-1},y_{n-1})$其中$x_i=\omega_{n}^{i}$
类似于中国剩余定理，我们可以构造出一个多项式
令$\phi(x)=y_0+y_1x+y_2x2...+y_{n-1}x$
则$\phi(\omega_{n}{-i})=\displaystyle\sum_{j=0} y_j *(\omega_{n}{-i})j$
打开$y_i$有：$\phi(\omega_{n}{-i})=\displaystyle\sum_{j=0}\displaystyle\sum_{k=0}{n-1}b_k*(\omega_{n})k*(\omega_{n})j=\displaystyle\sum_{j=0}\displaystyle\sum_{k=0}{n-1}b_k*\omega_{n}=\displaystyle\sum_{k=0}{n-1}\displaystyle\sum_{j=0}b_k\omega_{n}{j(k-i)}=\displaystyle\sum_{k=0}(b_k\displaystyle\sum_{j=0}{n-1}\omega_{n})$
这个形式跟中国剩余定理的答案构造异曲同工，那我们进行分类
① $k=i$
$b_i\displaystyle\sum_{j=0}{n-1}\omega_{n}=b_in$
②$k\neq i$
$\displaystyle\sum_{j=0}{n-1}\omega_{n}=\displaystyle\sum_{j=0}{n-1}(\omega_{n})j=(\omega_{n})0+(\omega_{n})1...+(\omega_{n}){n-1}=\frac{(\omega_{n})n-1}{\omega_{n}-1}=\frac{(\omega_{n}{n})-1} {\omega_{n}^{k-i}-1} =\frac{1-1}{\omega_{n}^{k-i}-1}=0$
$\phi(\omega_{n}^{-i})=n*b_i$
$b_i=\frac{\phi(\omega_{n}^{-i})}{n}$
所以我们只需要对$\phi(x)$代入$\omega_{n}{-0},\omega_{n},...,\omega_{n}^{-(n-1)}$做一遍FFT，把他的每个点值除以$n$就得到了$G(x)$的系数
优化
直接进行分治需要预开大量空间，而且常数大，我们尝试不要递归分治，而从最底下往上进行多项式的计算
所以我们需要把多项式长度扩充到$2^k$整，并且因此最底层的序数十分有规律,我们模拟一次最底层系数的分配过程
$(a_0,a_1,a_2,a_3, a_4,a_5,a_6,a_7)$
$(a_0,a_2,a_4,a_6),(a_1,a_3,a_5,a_7)$
$(a_0,a_4),(a_2,a_6),(a_1,a_5),(a_3,a_7)$
$(a_0),(a_4),(a_2),(a_6),(a_1),(a_5),(a_3),(a_7)$
如果把最高层当作第一层，则第$i$层相当于按照第$i$位的二进制对序数进行分类，为$0$放左，为$1$放右
所以整个过程相当于从低位到高位进行了一次高位到低位的排序
其结果就是最底层的二进制是反过来的
正常序$(000),(001),(010),(011),(100),(101),(110),(111)$
最底层$(000),(100),(010),(110),(001),(101),(011),(111)$
记$rev(i)$为$i$的逆序
显然有递推关系 $i=rev(rev(i/2))*2+i&1$
则有$rev(i)=rev(rev(rev(i/2))*2+i&1)$
注意到$rev(rev(i/2))*2$是一个偶数，则最低位是空出来的，与$i&1$不产生进位，所以我们可以直接拆开括号，有
$rev(i)=rev(rev(rev(i/2))*2)+rev(i&1)$
由于rev本身代表了逆序操作，打开$rev$时要把位数倒过来，如下
$rev(i)=rev(rev(rev(i/2)))/2+(i&1)*2^{k-1}$
于是有$rev(i)=rev(i/2)/2+(i&1)*2^{k-1}$
#include<bits/stdc++.h>
using namespace std;const int N = 10000000 + 5;
const double Pi = acos(-1.0); inline int read(){char ch=getchar();int x=0;while(ch<'0'||ch>'9') ch=getchar();for(;ch>='0'&&ch<='9';ch=getchar())x=(x<<3)+(x<<1)+(ch^48);return x;
}complex <double> a[N],b[N];
int n,m,lim,l,rev[N];void FFT(complex<double> *A,int opt){for(int i=0;i<lim;++i) if(i<rev[i]) swap(A[i],A[rev[i]]);for(int m=2;m<=lim;m<<=1){int mid=m>>1;complex<double> Wn(cos(2.0*Pi/m),opt*sin(2.0*Pi/m));for(int i=0;i<lim;i+=m){complex <double> w(1,0);for(int k=0;k<mid;++k,w=w*Wn){complex <double> x=A[i+k],y=w*A[i+k+mid];A[i+k]=x+y;A[i+k+mid]=x-y;}}}
}int main(){n=read(),m=read();for(int i=0;i<=n;++i) a[i].real(read());for(int i=0;i<=m;++i) b[i].real(read());lim=1,l=0;while(lim<n+m+1) lim<<=1,++l;for(int i=0;i<lim;++i)rev[i]=(rev[i>>1]>>1)+((i&1)<<(l-1));FFT(a,1);FFT(b,1);for(int i=0;i<lim;++i) a[i]=a[i]*b[i];FFT(a,-1);for(int i=0;i<=n+m;++i)printf("%d ",int(a[i].real()/lim+0.5));return 0;
}
