FFT 相关
FFT 相关
FFT 数组记得开两倍!
参考博客
FFT
简介
用于求卷积(\(a,b\) 已知):
或者多项式乘法(\(A,B\) 均是多项式且已知):
\(A=\sum_{i=0}^{n} a_i x^i\\ B=\sum_{i=0}^{m} b_i x^i\)
可见 \(C\) 是 \(n+m\) 次多项式。
如果我们把卷积的 \(a_i,b_i\) 看成多项式的系数,卷积就变成求:
求卷积或者多项式乘法的时间复杂度是 \(O(n^2)\) 的。使用 FFT 可以做到 \(O(n\log n)\)。
大体思路
设 \(C\) 的项数为 \(n\)(不是次数),若 \(A,B\) 不足 \(n\) 项就补系数 \(0\)。
显然这个过程很难优化,我们从另一个角度去想。
对于一个多项式,求其在 \(x\) 处的值的时间复杂度是 \(O(n)\) 的,我们把这个操作叫做点值(DFT)。
\(n\) 个点可以唯一确定一个 \(n\) 项多项式(即 \(n-1\) 次多项式),证明不显然但是略过我不会,大概感受一下吧。拉格朗日插值法给出了一种 \(O(n^2)\) 的求出这个多项式(即求出它的系数)的方法原理我不会,这个由点值求系数的过程叫做插值(IDFT)。
因此求多项式乘法,可以变成求任意 \(n\) 个点在两个多项式 \(A,B\) 的点值,然后由 \(C(x)=A(x)B(x)\),\(O(n)\) 求出多项式 \(C\) 在这 \(n\) 个点的点值,然后做一次插值求出 \(C\) 的系数。当然这个也是 \(O(n^2)\) 的,但是聪明的傅里叶给出了一种基于单位根的特殊性质的分治方法求 DFT 和 IDFT,成为快速傅里叶变换(FFT)。
就是先求出 \(A,B\) 的点值(DFT),进行点值相乘得到 \(C\) 的点值,然后插值(IDFT)得到 \(C\) 的系数。
单位根
写作 \(\omega_n^k\),读作 \(n\) 次单位根的 \(k\) 次方。
\(\omega_n^k\) 是在复数域上的向量,形如在复平面上分成 \(n\) 等分,其中 \(\omega_n^0=\omega_n^n=1\)。按逆时针分别为 \(\omega_n^0,\omega_n^1\dots \omega_n^{n-1}\)。
后面为了方便,我们会假设 \(n=2^k\)。
有几个重要的性质。
- \(\omega_n^n=1\) 正确性显然
- \(\omega_{an}^{ak}=\omega_n^k\) 在平面上想象一下,显然正确
- \(\omega_n^{k-\frac{n}{2}}=-\omega_n^k\) 相当于把向量转 \(180^。\)。
点值
基于这些性质,我们求 \(A,B\) 在 \(\omega_n^0,\omega_n^1\dots \omega_n^{n-1}\) 处的点值。
以求 \(A(x)\) 为例,我们要求所有 \(A(\omega_n^k),0\le k<n\)。
我们把 \(A\) 按奇偶分为两部分:
设
因此有:
代入 \(x=\omega_n^k\),有:
我们先求出 \(k< \frac{n}{2}\) 的 \(A_0,A_1\) 的点值。这个范围是缩小了一半的。然后把它们相加就得到了 \(A\) 在 \(k< \frac{n}{2}\) 的点值。
然后我们求剩下一半的 \(A_1,A_2\) 的点值。仍然设 \(k< \frac{n}{2}\) 发现:
因此你发现,这俩十分地相似,因此你求出 \(k< \frac{n}{2}\) 的 \(A_1,A_2\) 的点值之后,可以直接 \(O(n)\) 求出 \(A\) 的 \(n\) 个点值了。然后就这样分治下去,点值时间复杂度为 \(O(n\log n)\)。
因为要一直对 \(n\) 除以 \(2\) 所以令 \(n=2^k\) 意义就在此。(位数不足高位补系数 \(0\))
插值
求插值的过程是类似的,改几个参数就行了。
把求点值的过程写成矩阵。
由 \(A\) 得到 \(dft(A)_i = A(\omega_{n}^i)\)。有
把 \(p_{i,j}\) 写成矩阵就是:
只要把参数改成逆矩阵,其他不变,就可以由 \(dft(A)\) 求出 \(A\) 了。
根据单位根的性质,矩阵的逆恰好是指数取反:
改进
然后你会发现过不了板子……
因为这个递归的过程常数很大,假设 FFT 的常数本身就大,然后就超时了。
考虑从下往上递推分治的过程。
然后你惊奇地发现最下层的顺序是 \(000,100,010,110,001,101,011,111\),刚好是 \(0\sim 7\) 的二进制的 reverse。
然后就有了如下板子:
Code
struct fushu { double x,y; fushu (double _x=0,double _y=0):x(_x),y(_y){} }a[N],b[N]; fushu operator + (fushu a,fushu b) { return {a.x+b.x,a.y+b.y}; } fushu operator - (fushu a,fushu b) { return {a.x-b.x,a.y-b.y}; } fushu operator * (fushu a,fushu b) { return {a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x}; } int len; const double pi=acos(-1.0); int re[N]; void FFT(fushu *c,int type) { rep(i,0,(1<<len)-1) { if(i<re[i]) swap(c[i],c[re[i]]); } for(int k=1;k<(1<<len);k<<=1) { fushu wn(cos(pi/k),type*sin(pi/k)); for(int r=k<<1,j=0;j<(1<<len);j+=r) { fushu w(1,0); for(int i=0;i<k;i++,w=w*wn) { fushu x=c[j+i],y=w*c[j+k+i]; c[j+i]=x+y; c[j+k+i]=x-y; } } } } int n,m; int main(){ sf("%d%d",&n,&m); rep(i,0,n) sf("%lf",&a[i].x); rep(i,0,m) sf("%lf",&b[i].x); while((1<<len)<=n+m) len++; rep(i,0,(1<<len)-1) { re[i]=(re[i>>1]>>1)|((i&1)<<(len-1)); } FFT(a,1),FFT(b,1); rep(i,0,(1<<len)-1) a[i]=a[i]*b[i]; FFT(a,-1); rep(i,0,n+m) pf("%d ",(int)(a[i].x/(1<<len)+0.5)); }
算法缺陷
由于复数域是用浮点数计算的,所以会存在掉精度问题。如果答案是对一个特别的指数取模,如著名的 \(998244353\),可以使用原根代替单位根计算,在剩余系里计算而不是在复数域计算。详见 NTT。
NTT
阶和原根
欧拉定理:若 \(\gcd(a,n)=1\),则 \(a^{\varphi(n)} \equiv 1(\bmod n)\)。
阶:设 \(m>1\),且 \(\gcd(a,m)=1\),根据欧拉定理一定存在正整数 \(d<m\) 使得 \(a^d \equiv 1(\bmod m)\)。把满足该式子的最小的正整数定义为 \(a\) 对模 \(m\) 的阶(指数),写作 \(\text{ord}_m(a)\)。
阶的一个结论:对于正整数 \(d\),\(a^d \equiv 1 (\bmod m)\) 的一个充要条件是 \(d \mid \text{ord}_m(a)\)。
原根:设 \(m>1\),\(\gcd(a,m)=1\),若 \(\text{ord}_m(a)=\varphi(m)\),则 \(a\) 为 \(m\) 的原根,记做 \(g_m\)。
常见质数 \(998244353\) 的原根是 \(3\)。
定理 1:\(g_m,g_m^2,\dots g_m^{\varphi(m)-1}\) 两两模 \(m\) 不同余且均与 \(m\) 互质。(其实是两条定理)
定理 2:若 \(m\) 是质数,\(g_m,g_m^2,\dots,g_{m}^{m-2}\) 模 \(m\) 的余数恰好形成 \(1 \sim m-1\) 的排列,而 \(g_m^{m-1}\) 模 \(m\) 的余数则是 \(0\)。
由此定理可以发现原根有很多类似单位根的性质。
定理 3:\(m\) 的原根有 \(\varphi(\varphi(m))\) 个。
简介
NTT 的原理是用原根代替单位根。质数 \(p\) 的原根在模 \(p\) 剩余系意义下具有我们利用的单位根的性质的相同性质。因此如果要求的多项式的系数只需要模 \(p\) 意义下的,可以用原根代替单位根。因为原根是整型,所以可以避免精度问题。
Code
const int N=4e6+7,mod=998244353,G=3,invG=332748118; ll a[N],b[N]; int len; int re[N]; ll ksm(ll a,ll b=mod-2) { ll s=1; while(b) { if(b&1) s=s*a%mod; a=a*a%mod; b>>=1; } return s; } void NTT(ll *c,int type) { rep(i,0,(1<<len)-1) { if(i<re[i]) swap(c[i],c[re[i]]); } for(int k=1;k<(1<<len);k<<=1) { ll wn=ksm(type==1?G:invG,(mod-1)/(k<<1)); for(int r=k<<1,j=0;j<(1<<len);j+=r) { ll w=1; for(int i=0;i<k;i++,w=w*wn%mod) { ll x=c[j+i],y=w*c[j+k+i]%mod; c[j+i]=(x+y)%mod; c[j+k+i]=(x-y+mod)%mod; } } } } int n,m; int main(){ #ifdef LOCAL freopen("in.txt","r",stdin); freopen("my.out","w",stdout); #endif sf("%d%d",&n,&m); rep(i,0,n) sf("%lld",&a[i]),a[i]%=mod; rep(i,0,m) sf("%lld",&b[i]),b[i]%=mod; while((1<<len)<=n+m) len++; ll inv=ksm(1<<len); rep(i,0,(1<<len)-1) { re[i]=(re[i>>1]>>1)|((i&1)<<(len-1)); } NTT(a,1),NTT(b,1); rep(i,0,(1<<len)-1) a[i]=a[i]*b[i]%mod; NTT(a,-1); rep(i,0,n+m) pf("%lld ",a[i]*inv%mod); }
分治 FFT
形如 \(f_i=\sum_{j=1}^i f_{i-j}g_j\) 的卷积(\(f_0\) 需要初值)。
一般的卷积我们是知道卷起来的两个多项式的系数的,但是这个卷积显然我们只知道 \(g\) 的系数,却不知道 \(f\) 的系数,那么如何卷呢?可以分治 FFT 解决。
类似于 CDQ 分治,要求 $f_{1\sim n} $,先求 \(f_{1 \sim mid}\),然后计算左边对右边的贡献,即算 \(f_{1\sim n}=\sum_{j=1}^{mid} f_{j}g_{i-j}\),做一次 \(O(n\log n)\) 的 FTT,因为这个卷积不是标准卷积(指 \(j\) 不是枚举到 \(i\)),因此我们可以给 \(j>mid\) 的 \(f_j\) 当做 \(0\) 来做卷积,相当于求 \(\sum_{j=1}^i f_j g_{i-j}\)。然后我们再求 \(f_{mid+1\sim r}\),因为左边已经对右边贡献过了,所以我们递归右边的时候,计算贡献不需要再带上左边。
也就是说,求 \(f_{l \sim r}\),先算左边,然后算左边对右边的贡献,即 \(\sum_{j=1}^{r-l+1} f_{j+l-1}g_{(r-l+1)-j}\),这是一个标准卷积, FFT 时间复杂度是 \(O(len \log len)\)。这样分治下去,直到长度为 \(1\) 时,不用卷积了,直接返回值。分治一共 \(\log\) 层,每一层没个区间都要做 \(O(len\log len)\) 的 FTT,总时间复杂度为 \(O(n\log^2 n)\)。
本文来自博客园,作者:wing_heart,转载请注明原文链接:https://www.cnblogs.com/wingheart/p/18412826
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· DeepSeek “源神”启动!「GitHub 热点速览」
· 我与微信审核的“相爱相杀”看个人小程序副业
· 上周热点回顾(2.17-2.23)
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· 如何使用 Uni-app 实现视频聊天(源码,支持安卓、iOS)