多项式乘法与FFT——学习笔记

前言
这篇博客很长，因为考虑到读者很多都是对FFT的了解除了名字再无其它的大白，所以（并不）详细地解释了许多，只要你耐心地把它看完，你一定会对FFT有更深刻的理解。

---

准备工作
首先我们要从泰勒公式说起。
泰勒公式是指若一个函数在x=x0$x=x_0$处有n阶导数，我们可以用含这n阶导数和(x−x0)$(x-x_0)$的多项式来模拟函数在(x−Δx,x+Δx)（其中Δx无限接近于0）$(x-{\mathrm{\Delta }}_x,x+{\mathrm{\Delta }}_x)（其中{\mathrm{\Delta }}_x无限接近于0）$的取值。用公式表示就是：

（此图来源于百度百科)
然后，我们可以根据泰勒公式推导出欧拉公式，即

（证明请自行查阅百度百科）
现在我们来考虑方程xn=1$x^n=1$的解，由此引出n次单位根的概念，即e2πi/n$e^{2\pi i/n}$，单位根的0——n-1次方都是这个方程的解

---

多项式
多项式有两种常用的表示方法，系数表示法和点值表示法，系数表示法就是我们常用的y=kx+b$y=kx+b$，点值表示法是说我们可以用至少n个点来表示n-1次多项式（求出各个系数）。
而这两种表示法是可以互相转化的，从系数到点值我们称为离散傅里叶变换（DFT），从点值到系数我们称为逆离散傅里叶变换（IDFT）。

---

快速傅里叶变换（FFT）
结合上面两点，我们取单位根的0——n-1次方为取样点，即可在nlogn$nlogn$的时间内快速完成离散傅里叶变换。（在操作之前我们需要先将阶数补为2的t次方）

对于一个多项式：
A(x)=a0+a1∗x+a2∗x2+...+an−1∗xn$A(x)=a_0+a_1\ast x+a_2\ast x^2+...+a_{n-1}\ast x^n$
我们取e2πi/n$e^{2\pi i/n}$的0到n-1次方为取样点，用A[k]表示第k个取样点的结果
然后我们采取分支的思路，将A分为O和E
O(x)=a1+a3∗x+a5∗x2+...+an−1∗xn/2−1$O(x)=a_1+a_3\ast x+a_5\ast x^2+...+a_{n-1}\ast x^{n/2-1}$
E(x)=a0+a2∗x+a4∗x2+...+an−2∗xn/2−1$E(x)=a_0+a_2\ast x+a_4\ast x^2+...+a_{n-2}\ast x^{n/2-1}$
我们会发现A(x)=x∗O(x2)+E(x2)$A(x)=x\ast O(x^2)+E(x^2)$
即A[k]=x∗O[k]+E[k]$A[k]=x\ast O[k]+E[k]$（因为A[k]代表x$x$取e2π(k−1)/n$e^{2\pi (k-1)/n}$，O[k],E[k]代表x2$x^2$取e2π(k−1)/(n/2)=e4π(k−1)/n$e^{2\pi (k-1)/(n/2)}=e^{4\pi (k-1)/n}$，两者刚好相同）
这样我们就成功地把问题的规模缩小到了1/2，即在nlogn$nlogn$的时间内完成了DFT。

在DFT的过程中，我们相当于给初始系数矩阵A乘上了一个n*n的矩阵，如下图

那么我们做IDFT的话，只需要求这个矩阵的逆矩阵即可。
可以证明，对于该矩阵每一项取倒数再除以n就是该矩阵的逆矩阵。
这样时间复杂度就和DFT完全一样，是nlogn$nlogn$

至此我们完成了FFT

---

FFT的写法
从上面的叙述中，FFT显然可以通过递归来实现分治，但这样会花费大量空间用于创建和维护数组，而我们可以使用迭代法（循环）避免这一部分的空间使用。（NOI实际评测时理论上会开无限栈，即栈的大小算在总内存里面，所以递归写法应该也能过）

我们以0——7为例来研究这样分治的特点
0,1,2,3,4,5,6,7
(0,2,4,6)(1,3,5,7)
(0,4)(2,6)(1,5)(3,7)
0,4,2,6,1,5,3,7
写成二进制：000,100,010,110,001,101,011,111
如果每个数都从右往左读的话，我们就会（惊讶地）发现，这是按从小往大的顺序排序的，现在我们只要能将数按照这种方式排序，就能够用循环来写FFT了。
于是我们递推预处理出rev数组

void get_rev(int bit)
{
    for(int i=0;i<(1<<bit);i++)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
}

然后根据rev对a重新排序，就可以进行迭代FFT了，一个元素的DFT为其本身，两两合并得到两个元素的DFT，两两合并得到四个元素的DFT，进行下去直到得到整串的DFT。

---

例题
洛谷板子题
两多项式相乘，先各做一遍DFT，每个点值相乘，再做一遍IDFT即可。
代码：

#include <bits/stdc++.h>
using namespace std;
typedef complex<double> cpx;
const double pi=acos(-1);

cpx a[3000000],b[3000000];
int rev[3000000];
int n,m;

void get_rev(int bit)
{
    for(int i=0;i<(1<<bit);i++)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
}

void fft(cpx *a,int len,int dft)
{
    for(int i=0;i<len;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
    for(int step=1;step<len;step<<=1)
    {
        cpx wn=exp(cpx(0,pi*dft/step));
        for(int j=0;j<len;j+=(step<<1))
        {
            cpx wnk(1,0);
            for(int k=j;k<j+step;k++)
            {
                cpx x=a[k];
                cpx y=wnk*a[k+step];
                a[k]=x+y;
                a[k+step]=x-y;
                wnk*=wn;
            }
        }
    }
    if(dft==-1) for(int i=0;i<len;i++) a[i]/=len;
}

int main()
{
    scanf("%d%d",&n,&m);
    n++;m++;
    for(int i=0;i<n;i++)
    {
        double kk;
        scanf("%lf",&kk);
        a[i]=cpx(kk,0);
    }
    for(int i=0;i<m;i++)
    {
        double kk;
        scanf("%lf",&kk);
        b[i]=cpx(kk,0);
    }
    int bit=1,s=2;
    for(;(1<<bit)<n+m;bit++) s<<=1;
    get_rev(bit);
    fft(a,s,1);
    fft(b,s,1);
    for(int i=0;i<s;i++) a[i]*=b[i];
    fft(a,s,-1);
    for(int i=0;i<n+m-1;i++) printf("%d ",(int)(a[i].real()+0.5)); printf("\n");
    return 0;
}