FFT
简介
用于求卷积(\(a,b\) 已知):
\[\sum_{i=0}^n a_ib_{n-i} \]或者多项式乘法(\(A(x),B(x)\) 已知):
\[C(x)=A(x)B(x) \]\(A(x)=\sum_{i=0}^{n} a_i x^i\\ B(x)=\sum_{i=0}^{m} b_i x^i\)
可见 \(C(x)\) 是 \(n+m\) 次多项式。
如果我们把卷积的 \(a_i,b_i\) 看成多项式的系数,卷积就变成求:
\[[x^n] C(x)=A(x)B(x) \]其中 \([x^n]\) 表示 \(x^n\) 的系数。
求卷积或者多项式乘法的时间复杂度是 \(O(n^2)\) 的。使用 FFT 可以做到 \(O(n\log n)\)。
大体思路
设 \(C(x)\) 的项数为 \(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)。
单位根
写作 \(\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_1(x)=a_0+a_2 x+a_4 x^2+\dots a_{n-2} x^{\frac{n-2}{2}}\\ A_2(x)=a_1+a_3x^2+a_5x^4+\dots a_{n-1}x^{\frac{n-2}{2}}\]因此有:
\[A(x)=A_1(x^2)+xA_2(x^2) \]代入 \(x=\omega_n^k\),有:
\[A(\omega_n^k)=A_1(\omega_n^{2k})+\omega_n^k A_2(\omega_n^{2k})\\ A(\omega_n^k)=A_1(\omega_{\frac{n}{2}}^k)+\omega_n^k A_2(\omega_{\frac{n}{2}}^k)\]我们先求出 \(k< \frac{n}{2}\) 的 \(A_1,A_2\) 的点值。这个范围是缩小了一半的。然后把它们相加就得到了 \(A\) 在 \(k< \frac{n}{2}\) 的点值。
然后我们求剩下一半的 \(A_1,A_2\) 的点值。仍然设 \(k< \frac{n}{2}\) 发现:
\[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(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k})-\omega_n^{k} A_2(\omega_n^{2k})\\ A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_{\frac{n}{2}}^k)-\omega_n^{k} A_2(\omega_{\frac{n}{2}}^k)\]因此你发现,这俩十分地相似,因此你求出 \(k< \frac{n}{2}\) 的 \(A_1,A_2\) 的点值之后,可以直接 \(O(n)\) 求出 \(A\) 的 \(n\) 个点值了。然后就这样分治下去,点值时间复杂度为 \(O(n\log n)\)。
插值
求插值的过程是类似的,改几个参数就行了。我还没搞懂是怎么推出来的
然后我们就这么求两次点值,再求一次插值就做完了。
改进
然后你会发现过不了板子……
因为这个递归的过程常数很大,假设 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
#include<bits/stdc++.h>
#define sf scanf
#define pf printf
using namespace std;
typedef long long ll;
const int N=4e6+3;
int n,m;
struct fushu{
double x,y;
fushu (double xx=0,double yy=0){x=xx,y=yy;}
}a[N],b[N];
fushu operator + (fushu a,fushu b){ return (fushu){a.x+b.x,a.y+b.y}; };
fushu operator - (fushu a,fushu b){ return (fushu){a.x-b.x,a.y-b.y}; };
fushu operator * (fushu a,fushu b){ return (fushu){a.x*b.x-a.y*b.y,a.x*b.y+b.x*a.y}; }
int len;
const double pi=acos(-1.0);
int re[N];
void FFT(fushu *c,int type){
for(int i=0;i<(1<<len);i++){
if(i<re[i]) swap(c[i],c[re[i]]);
}
for(int mid=1;mid<(1<<len);mid<<=1){
fushu Wn( cos(pi/mid) , type*sin(pi/mid) );
for(int r=mid<<1,j=0;j<(1<<len);j+=r){
fushu w(1,0);
for(int i=0;i<mid;i++,w=w*Wn){
fushu x=c[j+i],y=w*c[j+mid+i];
c[j+i]=x+y;
c[j+mid+i]=x-y;
}
}
}
}
int main(){
sf("%d%d",&n,&m);
for(int i=0;i<=n;i++) sf("%lf",&a[i].x);
for(int i=0;i<=m;i++) sf("%lf",&b[i].x);
while((1<<len)<=n+m) len++;
for(int i=0;i<(1<<len);i++){
re[i]=(re[i>>1]>>1)|((i&1)<<(len-1));
}
FFT(a,1);
FFT(b,1);
for(int i=0;i<(1<<len);i++) a[i]=a[i]*b[i];
FFT(a,-1);
for(int i=0;i<=n+m;i++){
pf("%d ",(int)(a[i].x/(1<<len)+0.5));
}
}
算法缺陷
由于复数域是用浮点数计算的,所以会存在掉精度问题。如果答案是对一个特别的指数取模,如著名的 \(998244353\),可以使用原根代替单位根计算,在剩余系里计算而不是在复数域计算。详见 NTT。
标签:right,frac,FFT,点值,omega,fushu,left From: https://www.cnblogs.com/liyixin0514/p/18412826