文章目录
- 前言
- 系数表示法和点值表示法
- 单位根
- 离散傅立叶变换(DFT)
- 位逆序置换(蝴蝶变换)
- 离散傅立叶逆变换(IDFT)
- 代码
- NTT
- 代码
所谓快速傅立叶变换,就是傅立叶发明的一种快速的变换
(逃)
前言
FFT,用于在nlognnlognnlogn的复杂度内求两个多项式相乘。
算是多项式的入门吧。
十一被rabbit的数学打蒙了所以开始碰这些东西。
不知道以后会不会直接鸽掉。
(update on 2022.1.4:三个月后终于回来填坑了)。
很妙的算法。
代码写成迭代,精炼之后很优美。
系数表示法和点值表示法
系数表示法就是我们常用的函数的表示形式。
对于一个n次的函数,点值表示法则是给出了n+1个互异的点,通过这n+1个互异的点就可以唯一确定的描述出这个函数(可以想想高斯消元)。
这东西有啥用?
考虑如果两个多项式都表示成了点值形式,且选取的点的横坐标相同,我把它们的纵坐标对应乘起来,就能得到它们乘积的多项式的点值表示
乘起来的复杂度降到了O(n)O(n)O(n)
但是这个点值表示法求起来暴力似乎还是n方啊…
所以快速傅立叶变换其实解决的就是系数表示法和点值表示法快速互相转换的问题
单位根
规定:
若复数ω\omegaω满足ωn=1\omega^n=1ωn=1,就称ω\omegaω是n次单位根
对于一个n项的多项式,我们让它的点值表示取的横坐标为ωnk(0≤k<n)\omega_n^k(0\leq k <n)ωnk(0≤k<n)
然后就会有很多很妙的性质:
- ωnmkm=ωnk\omega_{nm}^{km}=\omega_n^kωnmkm=ωnk
- ωnk+n=ωnk\omega_{n}^{k+n}=\omega_n^kωnk+n=ωnk
- ωnk+n2=−wnk\omega_n^{k+\frac{n}{2}}=-w_n^kωnk+2n=−wnk
都可以结合几何意义较为显然的证明。
离散傅立叶变换(DFT)
那么回到问题。
首先,我们要把多项式填零直到 n=2kn=2^kn=2k 的长度。
接下来,我们要求:
f(x)=a0+a1x+a2x2+a3x3+...+an−1xn−1f(x)=a_0+a_1x+a_2x^2+a_3x^3+...+a_{n-1}x^{n-1}f(x)=a0+a1x+a2x2+a3x3+...+an−1xn−1
让我们把这个函数分奇偶拆成两个函数:
A(x)=a0+a2x+a4x2+a6x3+...+an−1xn2A(x)=a_0+a_2x+a_4x^2+a_6x^3+...+a_{n-1}x^{\frac{n}{2}}A(x)=a0+a2x+a4x2+a6x3+...+an−1x2n
B(x)=a1+a3x+a5x2+a7x3+...+an−2xn2B(x)=a_1+a_3x+a_5x^2+a_7x^3+...+a_{n-2}x^{\frac{n}{2}}B(x)=a1+a3x+a5x2+a7x3+...+an−2x2n
那么就有:
f(x)=A(x2)+xB(x2)f(x)=A(x^2)+xB(x^2)f(x)=A(x2)+xB(x2)
我们把 x=ωnkx=\omega_n^kx=ωnk 带入。
分情况讨论:
(设:0≤k<n20\le k<\dfrac{n}{2}0≤k<2n)
当指数 <n2< \dfrac{n}{2}<2n 时:
f(ωnk)=A((ωnk)2)+ωnkB((ωnk)2)f(\omega_n^k)=A((\omega_n^k)^2)+\omega_n^kB((\omega_n^k)^2)f(ωnk)=A((ωnk)2)+ωnkB((ωnk)2)
=A(ωn2k)+ωnkB(ωn2k)=A(\omega_n^{2k})+\omega_n^kB(\omega_n^{2k})=A(ωn2k)+ωnkB(ωn2k)
=A(ωn2k)+ωnkB(ωn2k)=A(\omega_{\frac{n}{2}}^{k})+\omega_n^kB(\omega_{\frac{n}{2}}^{k})=A(ω2nk)+ωnkB(ω2nk)
类似的,当指数 ≥n2\ge\dfrac{n}{2}≥2n 时:
f(ωnk+n2)=A(ωn2k+n)+ωnk+n2B(ωn2k+n)f(\omega_n^{k+\frac{n}{2}})=A(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}B(\omega_n^{2k+n})f(ωnk+2n)=A(ωn2k+n)+ωnk+2nB(ωn2k+n)
=A(ωn2k)−ωnkB(ωn2k)=A(\omega_n^{2k})-\omega_n^{k}B(\omega_n^{2k})=A(ωn2k)−ωnkB(ωn2k)
=A(ωn2k)−ωnkB(ωn2k)=A(\omega_{\frac{n}{2}}^{k})-\omega_n^{k}B(\omega_{\frac{n}{2}}^{k})=A(ω2nk)−ωnkB(ω2nk)
这样,我们就可以在 O(len)O(len)O(len) 的复杂度内合并两个长度为 lenlenlen 的多项式,就可以利用分治把复杂度降到 O(nlogn)O(n\log n)O(nlogn)。
位逆序置换(蝴蝶变换)
这个也很妙
有点线性求逆元那味了
通过观察发现,系数 i 调换位置之后会去往它的二进制表示反过来的位置(记为rir_iri)
而上面那个东西可以通过:
r[i]=r[i>>1]>>1+[i&1]∗(n>>1)r[i]=r[i>>1]>>1+[i\&1]*(n>>1) r[i]=r[i>>1]>>1+[i&1]∗(n>>1)
来线性的求
这个还是挺好理解的画一画就行了
离散傅立叶逆变换(IDFT)
我们变完之后当然还要变回来。
我会 n3n^3n3 高斯消元!
设之前DFT得到的点值表示 bk=f(ωnk)=∑i=0n−1ai(ωnk)ib_k=f(\omega_n^k)=\sum_{i=0}^{n-1}a_i(\omega_{n}^k)^ibk=f(ωnk)=∑i=0n−1ai(ωnk)i
构造一个函数:
g(x)=b0+b1x+b2x2+b3x3+...+bn−1xn−1g(x)=b_0+b_1x+b_2x^2+b_3x^3+...+b_{n-1}x_{n-1}g(x)=b0+b1x+b2x2+b3x3+...+bn−1xn−1
然后我们用 ωn0,ωn−1,ωn−2,...,ωn−n\omega_n^{0},\omega_n^{-1},\omega_n^{-2},...,\omega_n^{-n}ωn0,ωn−1,ωn−2,...,ωn−n 作为横坐标带入。(不难发现刚才 DFT 需要的性质依然成立,所以可以用同样的方法求解)
那么我们在第 kkk 位得到的新的点值就是:
g(ωn−k)=∑i=0n−1bi(ωn−k)ig(\omega_n^{-k})=\sum_{i=0}^{n-1}b_i(\omega_n^{-k})^ig(ωn−k)=i=0∑n−1bi(ωn−k)i
=∑i=0n−1∑j=0n−1aj(ωni)j(ωn−k)i=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(\omega_{n}^i)^j(\omega_n^{-k})^i=i=0∑n−1j=0∑n−1aj(ωni)j(ωn−k)i
=∑j=0j−1aj∑i=0n−1(ωnj−k)i=\sum_{j=0}^{j-1}a_j\sum_{i=0}^{n-1}(\omega_{n}^{j-k})^i=j=0∑j−1aji=0∑n−1(ωnj−k)i
考虑 ∑i=0n−1(ωnj−k)i\sum_{i=0}^{n-1}(\omega_{n}^{j-k})^i∑i=0n−1(ωnj−k)i,这一项,当 j=kj=kj=k 时,显然等于0;当 j≠kj\ne kj=k 时,根据等比公式,有:
∑i=0n−1(ωnj−k)i=(wnj−k)n−1wnj−k−1\sum_{i=0}^{n-1}(\omega_{n}^{j-k})^i=\frac{(w_n^{j-k})^n-1}{w_n^{j-k}-1}i=0∑n−1(ωnj−k)i=wnj−k−1(wnj−k)n−1
=(wnn)j−k−1wnj−k−1=1j−k−1wnj−k−1=0=\frac{(w_n^{n})^{j-k}-1}{w_n^{j-k}-1}=\frac{1^{j-k}-1}{w_n^{j-k}-1}=0=wnj−k−1(wnn)j−k−1=wnj−k−11j−k−1=0
所以原来的式子只有当 j=kj=kj=k 的时候有值,等于 nakna_knak。
所以我们直接 NTT 完除 nnn 即可。
代码
#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define il inline
const int N=5e6+100;
const int M=150;
const int mod=998244353;
const double pi=acos(-1.0);
inline ll read(){ll x=0,f=1;char c=getchar();while(!isdigit(c)){if(c=='-') f=-1;c=getchar();}while(isdigit(c)){x=x*10+c-'0';c=getchar();}return x*f;
}
int n,m,len,k;
struct node{double x,y;node(double a=0,double b=0){x=a;y=b;}
}A[N],B[N];
il node operator * (node a,node b){return (node){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};
}
il node operator + (node a,node b){return (node){a.x+b.x,a.y+b.y};
}
il node operator - (node a,node b){return (node){a.x-b.x,a.y-b.y};
}
int r[N];
il void fft(node *x,int flag){for(int i=0;i<=len;i++){if(i<r[i]) swap(x[i],x[r[i]]);}for(int l=1;l<len;l<<=1){node o(cos(pi/l),flag*sin(pi/l));for(int st=0;st<len;st+=l<<1){node t(1,0);for(int j=0;j<l;j++,t=t*o){node u=x[st+j],v=t*x[st+j+l];x[st+j]=u+v;x[st+j+l]=u-v;}}}return;
}
int main(){n=read();m=read();for(int i=0;i<=n;i++) A[i].x=read();for(int i=0;i<=m;i++) B[i].x=read();len=n+m;while((1<<k)<=len) k++;len=1<<k;for(int i=1;i<=len;i++) r[i]=(r[i>>1]>>1)|((i&1)<<(k-1));fft(A,1);fft(B,1);for(int i=0;i<=len;i++) A[i]=A[i]*B[i];fft(A,-1);for(int i=0;i<=n+m;i++) printf("%d ",(int)(A[i].x/(len)+0.5));return 0;
}
/**/
NTT
NTT说简单也简单说难也难。
简单是因为板子几乎和FFT一模一样。
难是因为那个群的概念完全看不懂qwq。
然后我就选择调过了证明部分。
背下来背下来!——宋队
总的来说,写一写对背板子有帮助的我可以理解的。
对于一个质数P,若其存在原根(基本全是3)且可以表示为:
P=k∗2n+1P=k*2^n+1P=k∗2n+1
它就是一个ntt模数。
ggg为什么叫原根?
因为它有一些关键性质:
1.gk∗n=1(modP)1. g^{k*n}=1(mod P)1.gk∗n=1(modP)
2.gi(modP)g^i(mod P)gi(modP) 在 0≤i≤p0\leq i\leq p0≤i≤p 的范围内取值两两不同。
设:gkg^kgk 为单位根,记为gng_ngn。
那么它就有很多熟悉的性质:
- gnn=1(modP)g_n^n=1(mod P)gnn=1(modP)
- g2k2k=gkkg_{2k}^{2k}=g_k^kg2k2k=gkk
- gnk+n/2=−gnkg_n^{k+n/2}=-g_n^kgnk+n/2=−gnk
没错,推FFT时 ω\omegaω 需要有的性质它全有!
所以就可以直接上FFT的板子啦!
然后烦人的复数、浮点运算和精度问题全都拜拜了。
针不戳。
宋队诚不欺我。
逆变换还是倒数,从FFT的求共轭改成求逆元。
就像喝水一样。
代码
il void ntt(ll *x,int flag){for(int i=0;i<lim;i++){if(i<r[i]) swap(x[i],x[r[i]]);}for(register int l=1;l<lim;l<<=1){register ll g=ksm(3,(mod-1)/(l<<1));if(flag==-1) g=ksm(g,mod-2);for(register int st=0;st<lim;st+=l<<1){register ll now=1;for(register int i=0;i<l;i++,now=now*g%mod){ll u=x[st+i],v=x[st+i+l]*now%mod;x[st+i]=(u+v)%mod;x[st+i+l]=(u-v+mod)%mod;}}}if(flag==-1){register ll ni=ksm(lim,mod-2);for(register int i=0;i<lim;i++) x[i]=x[i]*ni%mod;}
}