FFT 相关

FFT 相关

FFT 数组记得开两倍!

参考博客

快速傅里叶变换(一)FFT

快速傅里叶变换(二)NTT

快速傅里叶变换(三)

超详细傅里叶变换

FFT自己卷自己

FFT

简介

用于求卷积(\(a,b\) 已知):

\[\sum_{i=0}^n a_ib_{n-i} \]

或者多项式乘法(\(A,B\) 均是多项式且已知):

\[C=A\cdot 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\) 看成多项式的系数,卷积就变成求:

\[c_n = \sum_{i+j=n} a_i b_j \]

求卷积或者多项式乘法的时间复杂度是 \(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\) 按奇偶分为两部分:

\[A_0(x)=a_0+a_2 x+a_4 x^2+\dots a_{n-2} x^{\frac{n-2}{2}}\\ A_1(x)=a_1+a_3x^2+a_5x^4+\dots a_{n-1}x^{\frac{n-2}{2}}\]

因此有:

\[A(x)=A_0(x^2)+xA_1(x^2) \]

代入 \(x=\omega_n^k\),有:

\[\begin{aligned} A(\omega_n^k)& =A_0(\omega_n^{2k})+\omega_n^k A_1(\omega_n^{2k})\\ & = A_0(\omega_{\frac{n}{2}}^k)+\omega_n^k A_1(\omega_{\frac{n}{2}}^k) \end{aligned} \]

我们先求出 \(k< \frac{n}{2}\)\(A_0,A_1\) 的点值。这个范围是缩小了一半的。然后把它们相加就得到了 \(A\)\(k< \frac{n}{2}\) 的点值。

然后我们求剩下一半的 \(A_1,A_2\) 的点值。仍然设 \(k< \frac{n}{2}\) 发现:

\[\begin{aligned} A(\omega_n^{k+\frac{n}{2}}) & =A_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}} A_2(\omega_n^{2k+n})\\ &=A_1(\omega_n^{2k})-\omega_n^{k} A_2(\omega_n^{2k})\\ &=A_1(\omega_{\frac{n}{2}}^k)-\omega_n^{k} A_2(\omega_{\frac{n}{2}}^k) \end{aligned} \]

因此你发现,这俩十分地相似,因此你求出 \(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)\)。有

\[dft(A)_i = \sum_{j=0}^{n-1} \omega_{n}^{ij} a_j= \sum_{j=0}^{n-1} p_{i,j} a_j \]

\(p_{i,j}\) 写成矩阵就是:

\[\begin{bmatrix} 1 & 1 & \cdots & 1\\ 1 & \omega_{n}^1 & \cdots & \omega_{n}^{n-1}\\ \vdots & \vdots & \ddots & \vdots\\ 1 & \omega_{n}^{n-1} & \cdots & \omega_{n}^{(n-1)(n-1)} \end{bmatrix} \]

只要把参数改成逆矩阵,其他不变,就可以由 \(dft(A)\) 求出 \(A\) 了。

根据单位根的性质,矩阵的逆恰好是指数取反:

\[\begin{bmatrix} 1 & 1 & \cdots & 1\\ 1 & \omega_{n}^{-1} & \cdots & \omega_{n}^{-(n-1)}\\ \vdots & \vdots & \ddots & \vdots\\ 1 & \omega_{n}^{-(n-1)} & \cdots & \omega_{n}^{-(n-1)(n-1)} \end{bmatrix} \]

改进

code

然后你会发现过不了板子……

因为这个递归的过程常数很大,假设 FFT 的常数本身就大,然后就超时了。

考虑从下往上递推分治的过程。

\[\left(a_{0}, a_{1}, a_{2}, a_{3}, a_{4}, a_{5}, a_{6}, a_{7}\right)\\ \left(a_{0}, a_{2}, a_{4}, a_{6}\right)\left(a_{1}, a_{3}, a_{5}, a_{7}\right)\\ \left(a_{0}, a_{4}\right)\left(a_{2}, a_{6}\right)\left(a_{1}, a_{5}\right)\left(a_{3}, a_{7}\right)\\ \left(a_{0}\right)\left(a_{4}\right)\left(a_{2}\right)\left(a_{6}\right)\left(a_{1}\right)\left(a_{5}\right)\left(a_{3}\right)\left(a_{7}\right) \]

然后你惊奇地发现最下层的顺序是 \(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

Luogu 模板

形如 \(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)\)

posted @   wing_heart  阅读(11)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· DeepSeek “源神”启动!「GitHub 热点速览」
· 我与微信审核的“相爱相杀”看个人小程序副业
· 上周热点回顾(2.17-2.23)
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· 如何使用 Uni-app 实现视频聊天(源码,支持安卓、iOS)
点击右上角即可分享
微信分享提示