【LuoguP3803】多项式乘法-FFT/NTT模板题(附带FFT/NTT简单介绍)
2018-02-14 17:29
549 查看
测试地址:多项式乘法
做法1:本题需要用到大名鼎鼎的FFT——快速傅里叶变换(Fast Fourier Transform)。
FFT是一个……计算向量卷积的东西(大概吧,才疏学浅,还请各位大佬见谅),而向量卷积,说白了就是多项式乘法。在学习FFT之前,先简单介绍DFT(离散傅里叶变换)和IDFT(逆离散傅里叶变换)的概念,简单来说,DFT是将多项式由系数表示法变为点值表示法的过程,而IDFT则相反,是将多项式由点值表示法变为系数表示法的过程。
什么是点值表示法?对于一个nn次多项式,我们既可以用n+1n+1个系数(分别为00~nn次项的系数)来表示(系数表示法),也可以用其图象上的n+1n+1个点来表示,这就是点值表示法。为什么我们要得到点值表示?因为如果我们知道了两个多项式的点值表示,那么我们只用将横坐标相同的点的纵坐标相乘,就可以得到这两个多项式卷积的点值表示,这样显然是O(n)O(n)的。可是我们注意到,如果我们用最简单的方法实现两种表示法之间的转换,时间复杂度将不低于O(n2)O(n2),而这题的nn达到106106,我们需要更加快速的方法。
这个时候我们需要向多项式中代入一些具有特殊性质的变量来方便我们计算,这时候我们选用nn次单位根(从这里开始,为了讨论方便,nn都代表进行变换的向量的维数)。
什么是nn次单位根?nn次单位根指的是方程xn=1xn=1在复数集中的解,这样的解有nn个:cos2kπn+isin2kπn(k=0,1,...,n−1)cos2kπn+isin2kπn(k=0,1,...,n−1)(也可表示为e2kπne2kπn),记kk取11时的nn次单位根为ωnωn,这里有一个结论:kk取tt时的nn次单位根就是ωnωn的tt次方,记作ωtnωnt。
我们选用了nn次单位根为要选取的变量,它们有什么有用的性质呢?明确我们现在的任务,对于一个多项式AA,我们要求的是A(ω0n),A(ω1n),...,A(ωn−1n)A(ωn0),A(ωn1),...,A(ωnn−1)的值。我们直接对A(ωtn)A(ωnt)进行讨论,我们将所有奇数次项提出来,然后提出一个ωtnωnt,我们得到了一个这样的形式:
A(ωtn)=A[0](ω2tn)+ωtnA[1](ω2tn)A(ωnt)=A[0](ωn2t)+ωntA[1](ωn2t)。
看起来分成了A[0]A[0]和A[1]A[1]两个子问题,但是按照目前的形式看好像并没有减少需要运算的次数,我们此时需要用到另一个结论:
ω2tn=ωtn/2ωn2t=ωn/2t
这个应该比较好证明,直接带回原来的形式中就会发现两个复数的实部和虚部都相等。因此我们将上面一个形式变换为:
A(ωtn)=A[0](ωtn/2)+ωtnA[1](ωtn/2)A(ωnt)=A[0](ωn/2t)+ωntA[1](ωn/2t)
因为里面的变量的取值从nn个变成了n/2n/2个,所以问题规模确确实实的减半了,而且子问题也拥有和原问题一样的形式,所以可以递归计算,时间复杂度为O(nlogn)O(nlogn),这就是FFT了。为了保证复杂度,我们必须一开始就将向量的维数补成2k2k这种形式,方法是在向量后面添加00。
然而上述方法还有一个问题:不够快……
我们知道递归会增大算法的常数,所以我们尝试用迭代的方法解决这个问题。我们尝试得到算法进行到最底层时,系数的编号发生了什么变化。我们发现,第i(0≤i<n)i(0≤i<n)个位置上的是原来rev(i)rev(i)次项的系数,其中rev(i)rev(i)为ii的kk位二进制反转之后得到的数。因此我们先预处理出rev(i)rev(i),然后把算法当做自底向上的区间合并,这就是迭代的写法了。
然而我们还有另一个优化可以减小算法的常数。注意到上面的式子中:
A(ωtn)=A[0](ωtn/2)+ωtnA[1](ωtn/2)A(ωnt)=A[0](ωn/2t)+ωntA[1](ωn/2t)
我们枚举tt来计算A(ωtn)A(ωnt),然而我们知道ωkn=ωk mod nnωnk=ωnk mod n(由ωnn=1ωnn=1证得),因此我们就重算了两遍A[0/1]A[0/1]这个东西,而我们又知道ωn/2n=−1ωnn/2=−1,因此:
A(ωn/2+tn)=A[0](ωtn/2)−ωtnA[1](ωtn/2)A(ωnn/2+t)=A[0](ωn/2t)−ωntA[1](ωn/2t)
于是我们在算A(ωtn)A(ωnt)的同时算出了A(ωn/2+tn)A(ωnn/2+t),成功减小了一半的算法常数。事实证明,添加了上述优化的FFT算法已经足够优秀了。
然而,我们注意到上述的算法解决了DFT的过程,那IDFT怎么办?通过一些推理,可以得知:只要将向量通过上述方法得到的点值表示当做新的向量,再做一次上面的算法,只不过这次我们选择代入变量ω−tnωn−t而非原来的ωtnωnt,然后把结果都除以nn就行了。详情请看代码。
以下是本人代码:
做法2(2018.2.20更新):本题还可以用快速数论变换(NTT,Number Theory Transform)来解决。
快速数论变换很像上面的快速傅里叶变换,只不过FFT是取复数域的单位根,而NTT是取模某大质数PP意义域的单位根。根据数论中的一些知识,可以得到:ωn=gP−1nωn=gP−1n,其中gg为PP的一个原根,它保证在nn为P−1P−1的因子时,ωknωnk各不相同。而为了进行离散傅里叶变换,我们又必须保证nn是22的幂,所以就必须保证PP为形如c×2k+1c×2k+1的素数。常用的P,gP,g有:
P=1004535809=479×221+1,g=3P=1004535809=479×221+1,g=3
P=998244353=119×223+1,g=3P=998244353=119×223+1,g=3
于是将FFT稍加修改后得到了NTT的代码。
以下是本人代码:
做法1:本题需要用到大名鼎鼎的FFT——快速傅里叶变换(Fast Fourier Transform)。
FFT是一个……计算向量卷积的东西(大概吧,才疏学浅,还请各位大佬见谅),而向量卷积,说白了就是多项式乘法。在学习FFT之前,先简单介绍DFT(离散傅里叶变换)和IDFT(逆离散傅里叶变换)的概念,简单来说,DFT是将多项式由系数表示法变为点值表示法的过程,而IDFT则相反,是将多项式由点值表示法变为系数表示法的过程。
什么是点值表示法?对于一个nn次多项式,我们既可以用n+1n+1个系数(分别为00~nn次项的系数)来表示(系数表示法),也可以用其图象上的n+1n+1个点来表示,这就是点值表示法。为什么我们要得到点值表示?因为如果我们知道了两个多项式的点值表示,那么我们只用将横坐标相同的点的纵坐标相乘,就可以得到这两个多项式卷积的点值表示,这样显然是O(n)O(n)的。可是我们注意到,如果我们用最简单的方法实现两种表示法之间的转换,时间复杂度将不低于O(n2)O(n2),而这题的nn达到106106,我们需要更加快速的方法。
这个时候我们需要向多项式中代入一些具有特殊性质的变量来方便我们计算,这时候我们选用nn次单位根(从这里开始,为了讨论方便,nn都代表进行变换的向量的维数)。
什么是nn次单位根?nn次单位根指的是方程xn=1xn=1在复数集中的解,这样的解有nn个:cos2kπn+isin2kπn(k=0,1,...,n−1)cos2kπn+isin2kπn(k=0,1,...,n−1)(也可表示为e2kπne2kπn),记kk取11时的nn次单位根为ωnωn,这里有一个结论:kk取tt时的nn次单位根就是ωnωn的tt次方,记作ωtnωnt。
我们选用了nn次单位根为要选取的变量,它们有什么有用的性质呢?明确我们现在的任务,对于一个多项式AA,我们要求的是A(ω0n),A(ω1n),...,A(ωn−1n)A(ωn0),A(ωn1),...,A(ωnn−1)的值。我们直接对A(ωtn)A(ωnt)进行讨论,我们将所有奇数次项提出来,然后提出一个ωtnωnt,我们得到了一个这样的形式:
A(ωtn)=A[0](ω2tn)+ωtnA[1](ω2tn)A(ωnt)=A[0](ωn2t)+ωntA[1](ωn2t)。
看起来分成了A[0]A[0]和A[1]A[1]两个子问题,但是按照目前的形式看好像并没有减少需要运算的次数,我们此时需要用到另一个结论:
ω2tn=ωtn/2ωn2t=ωn/2t
这个应该比较好证明,直接带回原来的形式中就会发现两个复数的实部和虚部都相等。因此我们将上面一个形式变换为:
A(ωtn)=A[0](ωtn/2)+ωtnA[1](ωtn/2)A(ωnt)=A[0](ωn/2t)+ωntA[1](ωn/2t)
因为里面的变量的取值从nn个变成了n/2n/2个,所以问题规模确确实实的减半了,而且子问题也拥有和原问题一样的形式,所以可以递归计算,时间复杂度为O(nlogn)O(nlogn),这就是FFT了。为了保证复杂度,我们必须一开始就将向量的维数补成2k2k这种形式,方法是在向量后面添加00。
然而上述方法还有一个问题:不够快……
我们知道递归会增大算法的常数,所以我们尝试用迭代的方法解决这个问题。我们尝试得到算法进行到最底层时,系数的编号发生了什么变化。我们发现,第i(0≤i<n)i(0≤i<n)个位置上的是原来rev(i)rev(i)次项的系数,其中rev(i)rev(i)为ii的kk位二进制反转之后得到的数。因此我们先预处理出rev(i)rev(i),然后把算法当做自底向上的区间合并,这就是迭代的写法了。
然而我们还有另一个优化可以减小算法的常数。注意到上面的式子中:
A(ωtn)=A[0](ωtn/2)+ωtnA[1](ωtn/2)A(ωnt)=A[0](ωn/2t)+ωntA[1](ωn/2t)
我们枚举tt来计算A(ωtn)A(ωnt),然而我们知道ωkn=ωk mod nnωnk=ωnk mod n(由ωnn=1ωnn=1证得),因此我们就重算了两遍A[0/1]A[0/1]这个东西,而我们又知道ωn/2n=−1ωnn/2=−1,因此:
A(ωn/2+tn)=A[0](ωtn/2)−ωtnA[1](ωtn/2)A(ωnn/2+t)=A[0](ωn/2t)−ωntA[1](ωn/2t)
于是我们在算A(ωtn)A(ωnt)的同时算出了A(ωn/2+tn)A(ωnn/2+t),成功减小了一半的算法常数。事实证明,添加了上述优化的FFT算法已经足够优秀了。
然而,我们注意到上述的算法解决了DFT的过程,那IDFT怎么办?通过一些推理,可以得知:只要将向量通过上述方法得到的点值表示当做新的向量,再做一次上面的算法,只不过这次我们选择代入变量ω−tnωn−t而非原来的ωtnωnt,然后把结果都除以nn就行了。详情请看代码。
以下是本人代码:
#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是取模某大质数PP意义域的单位根。根据数论中的一些知识,可以得到:ωn=gP−1nωn=gP−1n,其中gg为PP的一个原根,它保证在nn为P−1P−1的因子时,ωknωnk各不相同。而为了进行离散傅里叶变换,我们又必须保证nn是22的幂,所以就必须保证PP为形如c×2k+1c×2k+1的素数。常用的P,gP,g有:
P=1004535809=479×221+1,g=3P=1004535809=479×221+1,g=3
P=998244353=119×223+1,g=3P=998244353=119×223+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&g 10874 t;>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; }
相关文章推荐
- [uoj 34 多项式乘法] FFT&NTT 模板
- 洛谷.3803.[模板]多项式乘法(NTT)
- [uoj 34 多项式乘法] FFT&NTT 模板
- 洛谷P3803 【模板】多项式乘法(FFT和NTT)
- 洛谷.3803.[模板]多项式乘法(FFT)
- FFT/NTT总结+洛谷P3803 【模板】多项式乘法(FFT)(FFT/NTT)
- 洛谷P3803 【模板】多项式乘法(FFT) 【fft】
- 多项式乘法,FFT与NTT
- uoj34 多项式乘法【FFT or NTT】
- 洛谷 P3803 【模板】多项式乘法(FFT)
- 【UOJ#34】 多项式乘法(FFT && NTT)
- 洛谷 P3803 【模板】多项式乘法(FFT)
- 【NTT模板】[UOJ#34]多项式乘法
- 【Uoj34】多项式乘法(NTT,FFT)
- 【Uoj34】多项式乘法(NTT,FFT)
- 洛谷P3803 【模板】多项式乘法(FFT)
- UOJ#34 多项式乘法(FFT模板题)
- UVALive - 6886 Golf Bot 多项式乘法(FFT)
- 【UOJ 34】 #34. 多项式乘法 (FFT)
- [UOJ 34]多项式乘法(FFT)