FFT/NTT总结+洛谷P3803 【模板】多项式乘法（FFT）（FFT/NTT）

前言

众所周知，这两个东西都是用来算多项式乘法的。

对于这种常人思惟难以理解的东西，就少些理解，多背板子吧！

所以只总结一下思路和代码，什么概念和推式子就靠巨佬们吧

推荐自为风月马前卒巨佬的概念和定理都很是到位的总结

推荐ppl巨佬的简明易懂的总结

FFT

多项式乘法的蹊径——点值表示法

通常咱们把两个长度为\(n\)的多项式乘起来，就相似于作竖式乘法，一位一位地乘再加到对应位上，是\(O(n^2)\)的

如何优化？直接看是没有思路的，只好另辟蹊径了。

多项式除了咱们经常使用的系数表示法\(y=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}\)之外，还能够用点值表示法。

所谓点值表示法，就是至关于用函数图像上\(n\)个点的坐标\((x_i,y_i)\)表示一个\(n\)次多项式，显然这个表示是惟一的。

咱们能够把系数表示转化成点值表示。点值表示下的多项式怎么相乘呢？就是选相同的\(x_i\)，把对应的\(y_i\)相乘。

固然，两个长度为\(n\)的多项式相乘获得的是长度为\(2n-1\)的多项式，须要\(2n-1\)个点值才能惟一表示，所以一开始两个多项式也要选\(2n-1\)个点表示。

举一个例子

\((x+1)(x+2)\)→ → → → \(x^2+3x+2\)
↓\(\qquad\qquad\qquad\qquad\qquad\qquad\)↑
↓（点值）\(\qquad\qquad\qquad\) （系数）↑
↓\(\qquad\qquad\qquad\qquad\qquad\qquad\)↑
\((0,1)(1,2)(2,3)\) （相乘）\(\qquad\quad\) ↑
\((0,2)(1,3)(2,4)\)→ → → →\((0,2)(1,6)(2,12)\)
但是，随意选\(O(n)\)个点，计算\(y\)是\(O(n)\)的，总时间仍是\(O(n^2)\)的，把它还原成系数表达式更不轻松。

因此，若是选的点很普通，这只是一条蹊径，并非一条捷径。

神奇的单位根

介绍一个神奇的东西——\(n\)次单位根（记为\(\omega\)）。

它是\(n\)个复数的集合，每个的\(n\)次方都等于\(1\)。其中的一个是\(e^{{2\pi i}\over n}\)，记为\(\omega_n\)。

普及一下欧拉公式，\(e^{\theta i}=\cos\theta+(\sin\theta)i\)，\(\theta\)就是这个复数向量的旋转角。显然知足\((\omega_n)^n=e^{2\pi i}=1\)，那么它的\(k\)次方\((\omega_n^k)^n=e^{2k\pi i}\)也等于\(1\)。

因而能够看出，知足\(n\)次方等于\(1\)的复数的取值只会有\(n\)个，为\(\omega_n^k(k\in[0,n-1])\)，由于会有\(\omega_n^{n+k}=\omega_n^k\)。

这\(n\)个复数向量在单位圆上呈放射状。下面是算导的图片：

介绍消去引理\(ω^{dk}_{dn}=ω^k_n\)，证实很容易的。

DFT&IDFT

单位根有什么用呢？

看看咱们把\(\omega_n^k(k\in[0,n-1])\)分别带入多项式求点值会发生什么就知道了。这个过程叫DFT。

假设\(n\)为偶数，那么咱们能够把它的奇偶项分开，用两个多项式表示它

\(A^{[0]}(x)=a_0+a_2x+a_4x^2+...+a_{n−2}x^{{\frac n 2}−1}\)

\(A^{[1]}(x)=a_1+a_3x+a_5x^2+...+a_{n−1}x^{{\frac n 2}−1}\)

那么\(A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)\)。

注意，下面的变化用到了\(ω^{2k}_n=ω^k_{\frac n 2}\)，\(ω^n_n=1\)和\(ω^{\frac n 2}_n=-1\)。

首先带入单位根

\(A(ω^k_n)=A^{[0]}(ω^{2k}_n)+ω^k_nA^{[1]}(ω^{2k}_n)\\\qquad=A^{[0]}(ω^k_{\frac n 2})+ω^k_nA^{[1]}(ω^k_{\frac n 2})(k<\frac n 2)\)

那\(k\geq\frac n 2\)时又会发生什么呢？把它变成\(\frac n 2+k\)

\(A(ω^{\frac n 2+k}_n)=A^{[0]}(ω^{n+2k}_n)+ω^{\frac n 2+k}_nA^{[1]}(ω^{n+2k}_n)\\\qquad\quad=A^{[0]}(ω^{2k}_n)+ω^{\frac n 2}_nω^k_nA^{[1]}(ω^{2k}_n)\\\qquad\quad =A^{[0]}(ω^k_{\frac n 2})-ω^k_nA^{[1]}(ω^k_{\frac n 2})\)

也就是说，若是咱们要求一个长度为\(n\)的多项式取\(\omega_n^k(k\in[0,n-1])\)的\(n\)个点值，咱们只须要求出两个长度为\(n/2\)的多项式取\(\omega_{\frac n 2}^k(k\in[0,\frac n 2-1])\)的\(\frac n 2\)个点值，再经过上述两个式子合并结果。这彻底就是个递归过程，很容易写一个函数来计算。

能够由算法写出DFT的复杂度\(T(n)=2T(\frac n 2)+O(n)=O(n\log n)\)

---

固然，别忘了还原成系数表示，这个过程叫作IDFT。

蒟蒻以为这里理解太麻烦了，所以再也不证实IDFT的过程，想了解的参见其它的总结。

只须要记住，IDFT的\(\omega_n\)是\(e^{-\frac{2\pi i}n}\)，最后的结果除以\(n\)，其它过程同DFT，能够写在一个函数里。具体看下面的代码：

void FFT(R complex<double>*a,R int n,R int op){//op=1为DFT，-1为IDFT
    if(!n)return;//为了方便，n的意义与上面不同，这里的n是a0、a1的长度
    complex<double>a0[n],a1[n];
    for(R int k=0;k<n;++k)
        a0[k]=a[k<<1],a1[k]=a[k<<1|1];//奇偶项分离
    FFT(a0,n>>1,op);FFT(a1,n>>1,op);//递归处理
    R complex<double>wn(cos(PI/n),sin(PI/n)*op),w(1,0);//单位根
    for(R int k=0;k<n;++k,w*=wn)//k从到n/2
        a[k]=a0[k]+w*a1[k],a[k+n]=a0[k]-w*a1[k];
}

递归版过程

引入例题：洛谷P3803 【模板】多项式乘法（FFT）

因为系数很小，咱们没必要担忧精度的问题了（固然float是使不得的）

递归版代码：

#include<cstdio>
#include<cmath>
#include<complex>
#include<iostream>
#define R register
#define G c=getchar()
using namespace std;
const int N=4.2e6;
const double PI=acos(-1);//自定义π
complex<double>f[N],g[N];
inline int in(){
    R char G;
    while(c<'-')G;
    return c&15;
}
void FFT(R complex<double>*a,R int n,R int op){//同上
    if(!n)return;
    complex<double>a0[n],a1[n];
    for(R int i=0;i<n;++i)
        a0[i]=a[i<<1],a1[i]=a[i<<1|1];
    FFT(a0,n>>1,op);FFT(a1,n>>1,op);
    R complex<double>W(cos(PI/n),sin(PI/n)*op),w(1,0);
    for(R int i=0;i<n;++i,w*=W)
        a[i]=a0[i]+w*a1[i],a[i+n]=a0[i]-w*a1[i];
}
int main(){
    R int n,m,i;
    scanf("%d%d",&n,&m);
    for(i=0;i<=n;++i)f[i]=in();
    for(i=0;i<=m;++i)g[i]=in();
    for(m+=n,n=1;n<=m;n<<=1);//把长度补到2的幂，没必要担忧高次项的系数，由于默认为0
    FFT(f,n>>1,1);FFT(g,n>>1,1);//DFT
    for(i=0;i<n;++i)f[i]*=g[i];//点值直接乘
    FFT(f,n>>1,-1);//IDFT
    for(i=0;i<=m;++i)printf("%.0f ",fabs(f[i].real())/n);//注意结果除以n，当心“-0”
    puts("");
    return 0;
}

好记又好写的递归版结果怎样呢？

Fast Fast TLE!只有77分。

最主要的缘由在于，空间调用太多了。

蝴蝶操做和Rader排序

针对效率过低的问题，咱们继续观察FFT实现过程当中的更多特色。

观察这一句代码a[k]=a0[k]+w*a1[k],a[k+n]=a0[k]-w*a1[k];，在操做过程当中，取出了a0[k]和a1[k]的值，经过计算又求出了a[k]和a[k+n]的值。咱们把这样的一次运算叫作“蝴蝶操做”。

这样的操做有什么特色呢？咱们试着将a0和a1合并成一个数组a，每一次蝴蝶操做后，取出了a[k]和a[k+n]的值，又求出了a[k]和a[k+n]的值。最后，整个数组都完成了求值，而并无用到两个数组！

以\(n=8\)为例，看看递归过程的结构

其实，咱们彻底能够递推求解！假设\(a\)数组已经变成了第四层的样子，先对a0和a四、a2和a六、a1和a五、a3和a7蝴蝶操做，变成第三层；再对a0和a二、a4和a六、a1和a三、a5和a7蝴蝶操做，变成第二层；最后对a0和a一、a2和a三、a4和a五、a6和a7蝴蝶操做，变成第一层，FFT就完成了，而空间复杂度仅为\(O(n)\)。这个过程能够用循环来控制。

剩下的问题就是把初始的数组变成最后一层的样子了。先别急着写一个递归函数暴力把位置换过去。来观察一下最后序列的编号的二进制表示000,100,010,110,001,101,011,111，是否是与原来000,001,010,011,100,101,110,111相比，每一个位置上的二进制位都反过来了？这样的变化叫作Rader排序。

咱们能够\(O(n)\)将Rader排序的映射关系求出。设\(i\)Rader排序后的数为\(r_i\)，咱们能够经过\(r_{\frac i 2}\)递推求出，具体方法能够看看代码&动动脑筋qwq

#include<cmath>
#include<cstdio>
#define R register
#define I inline
using namespace std;
const int N=4.2e6;
const double PI=acos(-1);
int n,r[N];
struct C{//手写complex，比STL快一点点
    double r,i;
    I C(){r=i=0;}
    I C(R double x,R double y){r=x;i=y;}
    I C operator+(R C&x){return C(r+x.r,i+x.i);}
    I C operator-(R C&x){return C(r-x.r,i-x.i);}
    I C operator*(R C&x){return C(r*x.r-i*x.i,r*x.i+i*x.r);}
    I void operator+=(R C&x){r+=x.r;i+=x.i;}
    I void operator*=(R C&x){R double t=r;r=r*x.r-i*x.i;i=t*x.i+i*x.r;}
}f[N],g[N];
I int in(){
    R char c=getchar();
    while(c<'-')c=getchar();
    return c&15;
}
I void FFT(R C*a,R int op){
    R C W,w,t,*a0,*a1;
    R int i,j,k;
    for(i=0;i<n;++i)//根据映射关系交换元素
        if(i<r[i])t=a[i],a[i]=a[r[i]],a[r[i]]=t;
    for(i=1;i<n;i<<=1)//控制层数
        for(W=C(cos(PI/i),sin(PI/i)*op),j=0;j<n;j+=i<<1)//控制一层中的子问题个数
            for(w=C(1,0),a1=i+(a0=a+j),k=0;k<i;++k,++a0,++a1,w*=W)
                t=*a1*w,*a1=*a0-t,*a0+=t;//蝴蝶操做
}
int main(){
    R int m,i,l=0;
    scanf("%d%d",&n,&m);
    for(i=0;i<=n;++i)f[i].r=in();
    for(i=0;i<=m;++i)g[i].r=in();
    for(m+=n,n=1;n<=m;n<<=1,++l);
    for(i=0;i<n;++i)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));//递推求r
    FFT(f,1);FFT(g,1);
    for(i=0;i<n;++i)f[i]*=g[i];
    FFT(f,-1);
    for(i=0;i<=m;++i)printf("%.0f ",fabs(f[i].r)/n);
    puts("");
    return 0;
}

NTT

原根

FFT的全过程都是基于小数运算，一旦数据较大就会丢失精度。咱们如何在整数范围内寻找和单位根\(\omega_n^n=1\)性质同样的数呢？

显然，若是不在模意义下这样的式子是不成立的。

引入原根，设为\(g\)，若是\(g^n=1(\mod p)\)中\(n\)的最小正整数解为欧拉函数\(\phi(p)\)，则称\(g\)为\(p\)的一个原根。

为了理解NTT，了解定义就够了，固然更多内容能够参考YL的总结

用\(\omega_n\)替换\(g\)，根据定义显然也知足消去引理\(ω^{dk}_{dn}=ω^k_n\)。那么放到多项式乘法里没问题。

一个问题在于，递推过程当中如何根据求出不一样子问题的\(ω\)？显然由于\(g^{p-1}=\omega_n^n=1\)，有\(\omega_n=g^{\frac{p-1}n}\)

固然，多项式的长度都被补成了\(2\)的幂，这就要求\(\phi(p)\)的质因子中含有大量的\(2\)；另外，IDFT中的单位根要对原根求逆，最后除以\(n\)也须要求逆，为了方便还要求\(p\)是质数。这样，能被用做NTT的模数的数很是少，常见的就是\(998244353(998244352=2^{23}×7×17)\)，\(3\)是它的一个原根。

上面那题的NTT代码

#include<cstdio>
#define I inline
#define RG register
#define RL RG LL
#define R RG int
typedef long long LL;
const int N=4.2e6,YL=998244353;//一块儿来%YL
int n,f[N],g[N],r[N];
I int in(){
    RG char c=getchar();
    while(c<'-')c=getchar();
    return c&15;
}
I int qpow(RL b,R k){//快速幂
    RL a=1;
    for(;k;k>>=1,b=b*b%YL)
        if(k&1)a=a*b%YL;
    return a;
}
I void NTT(R*a,R op){
    R i,j,k,d,wn,w,t,*a0,*a1;
    for(i=0;i<n;++i)
        if(i<r[i])t=a[i],a[i]=a[r[i]],a[r[i]]=t;
    for(i=1;i<n;i<<=1){
        wn=qpow(3,(YL-1)/(d=i<<1));//原根变换
        if(op&2)wn=qpow(wn,YL-2);//注意要求逆
        for(j=0;j<n;j+=d)
            for(w=1,a1=(a0=a+j)+i,k=0;k<i;++k,++a0,++a1,w=(LL)w*wn%YL)
                t=(LL)*a1*w%YL,*a1=(*a0-t+YL)%YL,*a0=(*a0+t)%YL;
    }
}
int main(){
    R m,i,l=0;
    scanf("%d%d",&n,&m);
    for(i=0;i<=n;++i)f[i]=in();
    for(i=0;i<=m;++i)g[i]=in();
    for(m+=n,n=1;n<=m;n<<=1,++l);
    for(i=0;i<n;++i)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
    NTT(f,1);NTT(g,1);
    for(i=0;i<n;++i)f[i]=(LL)f[i]*g[i]%YL;
    NTT(f,-1);
    R inv=qpow(n,YL-2);//一样注意要求逆
    for(i=0;i<=m;++i)printf("%lld ",(LL)f[i]*inv%YL);
    puts("");
    return 0;
}