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)
(qquadqquadqquadqquadqquadqquad)
↓(点值)(qquadqquadqquad) (系数)↑
(qquadqquadqquadqquadqquadqquad)
((0,1)(1,2)(2,3)) (相乘)(qquadquad)
((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^{{2pi i}over n}),记为(omega_n)

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

于是可以看出,满足(n)次方等于(1)的复数的取值只会有(n)个,为(omega_n^k(kin[0,n-1])),因为会有(omega_n^{n+k}=omega_n^k)

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

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

DFT&IDFT

单位根有什么用呢?

看看我们把(omega_n^k(kin[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))

(kgeqfrac 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)\qquadquad=A^{[0]}(ω^{2k}_n)+ω^{frac n 2}_nω^k_nA^{[1]}(ω^{2k}_n)\qquadquad =A^{[0]}(ω^k_{frac n 2})-ω^k_nA^{[1]}(ω^k_{frac n 2}))

也就是说,如果我们要求一个长度为(n)的多项式取(omega_n^k(kin[0,n-1]))(n)个点值,我们只需要求出两个长度为(n/2)的多项式取(omega_{frac n 2}^k(kin[0,frac n 2-1]))(frac n 2)个点值,再通过上述两个式子合并结果。这完全就是个递归过程,很容易写一个函数来计算。

可以由算法写出DFT的复杂度(T(n)=2T(frac n 2)+O(n)=O(nlog n))


当然,别忘了还原成系数表示,这个过程叫做IDFT。

蒟蒻觉得这里理解太麻烦了,因此不再证明IDFT的过程,想了解的参见其它的总结。

只需要记住,IDFT的(omega_n)(e^{-frac{2pi 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]的值。我们把这样的一次运算叫做“蝴蝶操作”。

这样的操作有什么特点呢?我们试着将a0a1合并成一个数组a,每一次蝴蝶操作后,取出了a[k]a[k+n]的值,又求出了a[k]a[k+n]的值。最后,整个数组都完成了求值,而并没有用到两个数组!

(n=8)为例,看看递归过程的结构

其实,我们完全可以递推求解!假设(a)数组已经变成了第四层的样子,先对a0和a4、a2和a6、a1和a5、a3和a7蝴蝶操作,变成第三层;再对a0和a2、a4和a6、a1和a3、a5和a7蝴蝶操作,变成第二层;最后对a0和a1、a2和a3、a4和a5、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;
}
原文地址:https://www.cnblogs.com/flashhu/p/9323521.html