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

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