//抄算导大法好

多项式A乘多项式B，求它们的乘积C，如果朴素地相乘，那么时间复杂度为O(n^2)，其中n为次数界。于是我们思考它的高效算法。

首先，从某种意义上，多项式的系数表达与点值表达是等价的，用点值形式表示的多项式都对应唯一系数形式的多项式（就像两点确定一条直线，三个不共线的点确定一条抛物线），选取n个点后通过霍纳法则，可以在O(n^2)的时间复杂度内求得。知道了A、B的点值后，因为对于任意点xk,C(xk)=A(xk)*B(xk)，就可以在O(n)的时间内求出C的点值（注意degree(C)=degree(A)+degree(B)，所以如果A、B的次数界为n，那么C的次数界为2n）。又可以证明得到插值多项式的唯一性，就可以运用拉格朗日公式在O(n^2)的时间复杂度内求出C的系数。这样看起来绕了一大圈什么收益都没有，于是我们对其进行优化。

我们可以巧妙地选取那n个点。我们断言：如果使用单位复数根，可以在O(nlgn)的时间内完成求值与插值运算。

n次单位负数根：满足ω^n=1的复数ω。n次单位负数根恰好有n个：对于k=0,1,...,n-1，这些根是e^(2πik/n) （e^(iu)=cos(u)+i*sin(u)），它们均匀地分布在以复平面的原点为圆心的单位半径的圆周上。n个n次单位负数根ωn^0,ωn^1,...,ωn^(n-1)在乘法意义下形成一个群，所以ωn^j*ωn^k=ωn^(j+k)=ωn^((j+k)mod n)。

然后我们将这n个单位复数根作为那n个点，于是有了DFT。但这是暴力计算，没有本质区别，所以我们可以利用单位复数根的特点，在O(nlgn)的时间内计算出DFT（注意，假设n是2的幂次方，如果不是就将它补到2的幂次方）。我们重新定义两个新的次数界为n/2的多项式A[0](x)和A[1](x):A[0](x)=a0+a2x+a4x^2+...+an-2x^(n/2-1)，A[1](x)=a1+a3x+a5x^2+...+an-1x^(n/2-1) ,于是有A(x)=A[0](x^2)+x*A[1](x^2)。那么我们就要求(ωn^0)^2,(ωn^1)^2,...,(ωn^(n-1))^2这n个点。根据折半引理：如果n>0为偶数，那么n个n次单位复数根的平方的集合就是n/2个n/2单位负数根的集合。可以得知那n个点的值并不是由n个不同值组成，而是仅由n/2个n/2次单位复数根所组成，每个根正好出现2次，那么就可以递归地求解，将原问题分解为2个形式相同，规模为一半的子问题。那么它的运行时间为T(n)=2T(n/2)+O(n)=O(nlgn)。

再考虑到在单位复数根处插值。既然点值的计算可以用范德蒙德矩阵来表示，那么我们只要给定一个逆矩阵，就可以推出IDFT，只需用ωn^(-1)替换ωn,将计算结果的每个元素除以n。证明？有个卷积定理：对任意两个长度为n的向量a和b，其中n是2的幂，a\otimesb=DFT2n^-1(DFT2n(a)·DFT2n(b))，其中向量a和b用0填充，使其长度达到2n，并用“·”表示2个2n个元素组成向量的点乘。

虽然这样的FFT时间复杂度是O(nlgn)，但是因为要递归等原因，并不够高效。所以，我们将递归的FFT改为迭代的FFT，只要将递归树直接从下往上推就行了。注意，叶子的顺序并不是a0,a1,...,an-1，它的顺序是一个位逆序置换。例如当n=8时，叶出现的次序为0,4,2,6,1,5,3,7，它们的二进制表示为000,100,010,110,001,101,011,111，逆序后就成了000,001,010,011,100,101,110,111，正好是按顺序的！

 

typedef complex<double> E;//复平面
//求逆序（朴素版，反正对时间复杂度影响小~）
int rev(int x)
{
    int a[30],i,tmp=0,ans=0;
    for(;x;x>>=1)
        a[++tmp]=x&1;
    for(i=1;i<=tmp;i++)
        ans+=a[i]<<(l-i);//l是lgn
    return ans;
}
void fft(E *a,int p)//p=1时是fft，p=-1时是ifft
{
    int i,j,k;
    for(i=0;i<n;i++)
        if(i<rev(i))swap(a[i],a[rev(i)]);
    for(i=1;i<n;i<<=1)//自底向上计算层次
    {
        E wn(cos(pi/i),sin(pi/i)*p);//每一层对应一个ωn
        for(j=0;j<n;j+=i<<1)//枚举每对DFT
        {
            E w(1,0);//初始化
            for(k=0;k<i;k++,w*=wn)//对两个n/2个元素的DFT进行组合
            {
                E x=a[j+k],y=a[j+k+i]*w;
                a[j+k]=x+y;a[j+k+i]=x-y;//蝴蝶操作
            }
        }
    }
    if(p==-1)//如果是ifft要对每一个元素除以n
        for(i=0;i<n;i++)
            a[i]/=n;
}