【LuoguP3803】多项式乘法-FFT/NTT模板题（附带FFT/NTT简单介绍）

测试地址：多项式乘法
做法1：本题需要用到大名鼎鼎的FFT——快速傅里叶变换（Fast Fourier Transform）。
FFT是一个……计算向量卷积的东西（大概吧，才疏学浅，还请各位大佬见谅），而向量卷积，说白了就是多项式乘法。在学习FFT之前，先简单介绍DFT（离散傅里叶变换）和IDFT（逆离散傅里叶变换）的概念，简单来说，DFT是将多项式由系数表示法变为点值表示法的过程，而IDFT则相反，是将多项式由点值表示法变为系数表示法的过程。
什么是点值表示法？对于一个n$n$次多项式，我们既可以用n+1$n+1$个系数（分别为0$0$~n$n$次项的系数）来表示（系数表示法），也可以用其图象上的n+1$n+1$个点来表示，这就是点值表示法。为什么我们要得到点值表示？因为如果我们知道了两个多项式的点值表示，那么我们只用将横坐标相同的点的纵坐标相乘，就可以得到这两个多项式卷积的点值表示，这样显然是O(n)$O(n)$的。可是我们注意到，如果我们用最简单的方法实现两种表示法之间的转换，时间复杂度将不低于O(n2)$O(n^2)$，而这题的n$n$达到106${10}^6$，我们需要更加快速的方法。
这个时候我们需要向多项式中代入一些具有特殊性质的变量来方便我们计算，这时候我们选用n$n$次单位根（从这里开始，为了讨论方便，n$n$都代表进行变换的向量的维数）。
什么是n$n$次单位根？n$n$次单位根指的是方程xn=1$x^n=1$在复数集中的解，这样的解有n$n$个：cos2kπn+isin2kπn(k=0,1,...,n−1)$\cos \frac{2k\pi }{n}+i\sin \frac{2k\pi }{n}(k=0,1,...,n-1)$（也可表示为e2kπn$e^{\frac{2k\pi }{n}}$），记k$k$取1$1$时的n$n$次单位根为ωn${\omega }_n$，这里有一个结论：k$k$取t$t$时的n$n$次单位根就是ωn${\omega }_n$的t$t$次方，记作ωtn${\omega }_n^t$。
我们选用了n$n$次单位根为要选取的变量，它们有什么有用的性质呢？明确我们现在的任务，对于一个多项式A$A$，我们要求的是A(ω0n),A(ω1n),...,A(ωn−1n)$A({\omega }_n^0),A({\omega }_n^1),...,A({\omega }_n^{n-1})$的值。我们直接对A(ωtn)$A({\omega }_n^t)$进行讨论，我们将所有奇数次项提出来，然后提出一个ωtn${\omega }_n^t$，我们得到了一个这样的形式：
A(ωtn)=A[0](ω2tn)+ωtnA[1](ω2tn)$A({\omega }_n^t)=A^{[0]}({\omega }_n^{2t})+{\omega }_n^t{A}^{[1]}({\omega }_n^{2t})$。
看起来分成了A[0]$A^{[0]}$和A[1]$A^{[1]}$两个子问题，但是按照目前的形式看好像并没有减少需要运算的次数，我们此时需要用到另一个结论：
ω2tn=ωtn/2${\omega }_n^{2t}={\omega }_{n/2}^t$
这个应该比较好证明，直接带回原来的形式中就会发现两个复数的实部和虚部都相等。因此我们将上面一个形式变换为：
A(ωtn)=A[0](ωtn/2)+ωtnA[1](ωtn/2)$A({\omega }_n^t)=A^{[0]}({\omega }_{n/2}^t)+{\omega }_n^t{A}^{[1]}({\omega }_{n/2}^t)$
因为里面的变量的取值从n$n$个变成了n/2$n/2$个，所以问题规模确确实实的减半了，而且子问题也拥有和原问题一样的形式，所以可以递归计算，时间复杂度为O(nlogn)$O(n\log n)$，这就是FFT了。为了保证复杂度，我们必须一开始就将向量的维数补成2k$2^k$这种形式，方法是在向量后面添加0$0$。
然而上述方法还有一个问题：不够快……
我们知道递归会增大算法的常数，所以我们尝试用迭代的方法解决这个问题。我们尝试得到算法进行到最底层时，系数的编号发生了什么变化。我们发现，第i(0≤i<n)$i(0\le i<n)$个位置上的是原来rev(i)$rev(i)$次项的系数，其中rev(i)$rev(i)$为i$i$的k$k$位二进制反转之后得到的数。因此我们先预处理出rev(i)$rev(i)$，然后把算法当做自底向上的区间合并，这就是迭代的写法了。
然而我们还有另一个优化可以减小算法的常数。注意到上面的式子中：
A(ωtn)=A[0](ωtn/2)+ωtnA[1](ωtn/2)$A({\omega }_n^t)=A^{[0]}({\omega }_{n/2}^t)+{\omega }_n^t{A}^{[1]}({\omega }_{n/2}^t)$
我们枚举t$t$来计算A(ωtn)$A({\omega }_n^t)$，然而我们知道ωkn=ωk mod nn${\omega }_n^k={\omega }_n^{kmodn}$（由ωnn=1${\omega }_n^n=1$证得），因此我们就重算了两遍A[0/1]$A^{[0/1]}$这个东西，而我们又知道ωn/2n=−1${\omega }_n^{n/2}=-1$，因此：
A(ωn/2+tn)=A[0](ωtn/2)−ωtnA[1](ωtn/2)$A({\omega }_n^{n/2+t})=A^{[0]}({\omega }_{n/2}^t)-{\omega }_n^t{A}^{[1]}({\omega }_{n/2}^t)$
于是我们在算A(ωtn)$A({\omega }_n^t)$的同时算出了A(ωn/2+tn)$A({\omega }_n^{n/2+t})$，成功减小了一半的算法常数。事实证明，添加了上述优化的FFT算法已经足够优秀了。
然而，我们注意到上述的算法解决了DFT的过程，那IDFT怎么办？通过一些推理，可以得知：只要将向量通过上述方法得到的点值表示当做新的向量，再做一次上面的算法，只不过这次我们选择代入变量ω−tn${\omega }_n^{-t}$而非原来的ωtn${\omega }_n^t$，然后把结果都除以n$n$就行了。详情请看代码。
以下是本人代码：

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <algorithm>
#include <cmath>
using namespace std;
int n,m,r[3000010];
const double pi=acos(-1.0);

struct Complex
{
    double x,y;
}a[3000010],b[3000010];
Complex operator + (Complex a,Complex b) {Complex s;s.x=a.x+b.x,s.y=a.y+b.y;return s;}
Complex operator - (Complex a,Complex b) {Complex s;s.x=a.x-b.x,s.y=a.y-b.y;return s;}
Complex operator * (Complex a,Complex b) {Complex s;s.x=a.x*b.x-a.y*b.y,s.y=a.x*b.y+a.y*b.x;return s;}

void FFT(Complex *a,int type)
{
    for(int i=0;i<n;i++)
        if (i<r[i]) swap(a[i],a[r[i]]);
    for(int mid=1;mid<n;mid<<=1) //枚举当前区间的半长mid
    {
        Complex W={cos(pi/mid),(double)type*sin(pi/mid)}; //W=omega_n
        for(int l=0,g=mid<<1;l<n;l+=g) //枚举区间的左端点l，g为区间的长度
        {
            Complex w={1.0,0.0}; //w从omega_n^0开始枚举
            for(int k=0;k<mid;k++,w=w*W)
            {
                Complex x=a[l+k],y=w*a[l+mid+k];
                a[l+k]=x+y;
                a[l+mid+k]=x-y;
            }
        }
    }
    if (type==-1)
    {
        for(int i=0;i<n;i++)
            a[i].x/=n;
    }
}

int read()
{
    int s=0;
    char c;
    c=getchar();
    while(c<'0'||c>'9') c=getchar();
    while(c>='0'&&c<='9') s=s*10+c-'0',c=getchar();
    return s;
}

int main()
{
    n=read(),m=read();
    n++,m++;
    for(int i=0;i<n;i++)
    {
        a[i].x=(double)read();
        a[i].y=0.0;
    }
    for(int i=0;i<m;i++)
    {
        b[i].x=(double)read();
        b[i].y=0.0;
    }

    int bit=0,x=1,y=n;
    while(x<n+m-1) bit++,x<<=1;
    for(int i=n;i<x;i++) a[i].x=a[i].y=0.0;
    for(int i=m;i<x;i++) b[i].x=b[i].y=0.0;
    for(int i=0;i<x;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(bit-1));
    n=x;

    FFT(a,1),FFT(b,1);
    for(int i=0;i<n;i++) a[i]=a[i]*b[i];
    FFT(a,-1);
    for(int i=0;i<y+m-1;i++)
        printf("%d ",(int)(a[i].x+0.5));

    return 0;
}

做法2（2018.2.20更新）：本题还可以用快速数论变换（NTT，Number Theory Transform）来解决。
快速数论变换很像上面的快速傅里叶变换，只不过FFT是取复数域的单位根，而NTT是取模某大质数P$P$意义域的单位根。根据数论中的一些知识，可以得到：ωn=gP−1n${\omega }_n=g^{\frac{P-1}{n}}$，其中g$g$为P$P$的一个原根，它保证在n$n$为P−1$P-1$的因子时，ωkn${\omega }_n^k$各不相同。而为了进行离散傅里叶变换，我们又必须保证n$n$是2$2$的幂，所以就必须保证P$P$为形如c×2k+1$c\times 2^k+1$的素数。常用的P,g$P,g$有：
P=1004535809=479×221+1,g=3$P=1004535809=479\times 2^{21}+1,g=3$
P=998244353=119×223+1,g=3$P=998244353=119\times 2^{23}+1,g=3$
于是将FFT稍加修改后得到了NTT的代码。
以下是本人代码：

#include <bits/stdc++.h>
#define ll long long
#define mod 998244353
#define g 3
using namespace std;
int n,m,r[3000010];
ll a[3000010]={0},b[3000010]={0};

ll power(ll a,ll b)
{
    ll s=1,ss=a;
    while(b)
    {
        if (b&1) s=(s*ss)%mod;
        ss=(ss*ss)%mod;b>>=1;
    }
    return s;
}

void NTT(ll *a,int type)
{
    for(int i=0;i<n;i++)
        if (i<r[i]) swap(a[i],a[r[i]]);
    for(int mid=1;mid<n;mid<<=1)
    {
        ll W=power(g,(mod-1)/(mid<<1));
        if (type==-1) W=power(W,mod-2);
        for(int l=0,r=mid<<1;l<n;l+=r)
        {
            ll w=1;
            for(int k=0;k<mid;k++,w=(w*W)%mod)
            {
                ll x=a[l+k],y=(w*a[l+mid+k])%mod;
                a[l+k]=(x+y)%mod;
                a[l+mid+k]=((x-y)%mod+mod)%mod;
            }
        }
    }
    if (type==-1)
    {
        ll inv=power(n,mod-2);
        for(int i=0;i<n;i++)
            a[i]=(a[i]*inv)%mod;
    }
}

int read()
{
    int s=0;
    char c;
    c=getchar();
    while(c<'0'||c>'9') c=getchar();
    while(c>='0'&&c<='9') s=s*10+c-'0',c=getchar();
    return s;
}

int main()
{
    n=read(),m=read();
    n++,m++;
    for(int i=0;i<n;i++) a[i]=(ll)read();
    for(int i=0;i<m;i++) b[i]=(ll)read();

    int bit=0,x=1,y=n;
    while(x<n+m-1) bit++,x<<=1;
    for(int i=0;i<x;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(bit-1));
    n=x;

    NTT(a,1),NTT(b,1);
    for(int i=0;i<n;i++) a[i]=(a[i]*b[i])%mod;
    NTT(a,-1);
    for(int i=0;i<y+m-1;i++)
        printf("%lld ",a[i]);

    return 0;
}