「学习笔记」FFT 快速傅里叶变换
「学习笔记」FFT 快速傅里叶变换
几个星期之后,继 扩展欧拉定理 之后, \(lj\) 大佬又给我们来了一发数论...
虽然听得心态爆炸, 但是还好的是没有 \(ymx\) 大佬的飞机开得好...
至少我还没有坐飞机...
啥是 FFT 呀?它可以干什么?
首先,你需要知道 矩阵乘法 的相关知识。
通过 矩阵乘法 的知识,我们知道,对于一个 \(f(x)\) 与 \(g(x)\) ,如果我们用朴素的矩阵乘法计算 \(f(x)\cdot g(x)\) ,时间复杂度是 \(O(N^2M)\) ,但是对于这是函数,复杂度中的 \(M=1\) ,那么时间复杂度为 \(O(N^2)\) 。
似乎这是比较优秀的,但是看看这道 板题 :
嗯? \(n,m\le 10^6\) ,似乎过于太小了?(不要管这病句...)
而 FFT 就是用来解决多项式乘法的问题的,可以把它的时间复杂度优化到 \(O(nlog^n_2)\) 。
这里引用一句话:
实际上,FFT 并不是直接计算多项式乘法,而是把原来的多项式 \(f(x),g(x)\) 在 \(O(nlog^n_2)\) 的复杂度内转换为它的 点值表示(后面会讲),而点值表示的多项式相乘的时间复杂度是 \(O(n)\) 的。最后再用 \(O(nlog_2^n)\) 的时间复杂度把所得多项式的点值表示转化为一般形式。
必备芝士
在学习 FFT 之前,我们需要知道 点值表示 和 复数 。
点值表示
什么是点值表示
对于一个函数 \(f(x)\) ,我们可以用一个多项式表示它:
这种表达方法我们叫做 一般形式 。
如果熟悉,那么可以确定一个东西,如果我们有在这个函数上的 \(n+1\) 个点的具体的坐标 \((x_i,y_i)\) ,那么我们就可以唯一确定这个函数的解析式,那么,我们也可以用这 \(n+1\) 个点的具体坐标来表示这个多项式。
当然,如果知道 拉格朗日插值法 的大佬,或许可以更好地理解这句话。
也就是说一个多项式与一个点值表示是一一对应的。
那么 FFT 完成的操作就是:
- 把已知的一个多项式转化成对应的点值表示;
- 把已知的点值表示转换成对应的多项式。
复杂度都是 \(O(nlog_2^n)\) 。
点值表示的乘法
那么,假设我们用 \(\{x_0,x_1,x_2,...,x_n\}\) 来表示 \(f(x)\) 为 \(\{(x_0,f(x_0)),(x_1,f(x_1)),...,(x_n,f(x_n))\}\) ,表示 \(g(x)\) 为 \(\{(x_0,g(x_0)),(x_1,g(x_1)),...,(x_n,g(x_n))\}\) 。
那么我们就知道 \(f(x)\cdot g(x)\) 的点值表达式为 \(\{(x_0,f(x_0)\cdot g(x_0)),...,(x_n,f(x_n)\cdot g(x_n))\}\) ,其实就是对应项相乘,那么这个的时间复杂度为 \(O(n)\) ,似乎可以接受。
复数
复数的定义
\(Hint.\) 是 复数(complex) 而不是 负数(negative) 。
首先我们要知道 虚数 ,即有一个很特别的数 \(i\) ,满足 \(i^2=-1\) ,即 \(i=\sqrt{-1}\) 。
不用理解为什么,因为他本来就是特殊定义的,只需要记住就好。
然后一个复数 \(a\) 可以表示为实部 \(x\) 和虚部 \(y\times i\) 之和—— \(a=x+y\times i\) 。对于一个实数,我们可以把它放在“一维数轴”上,即数轴;那么对于一个复数,我们需要把它放在“二维数轴”,也就是直角坐标系上。
其实,我们可以把 \(a\) 表示成 向量(vector) 的形式,即将 \(a=x+y\times i\) 表示为向量 \((x,y)\) ,但是对于后面虚数的乘法,其实是不满足向量乘法的,但是这可以作为一个理解方式。
复数的乘法运算
对于两个虚数 \(a,b\) ,其中 \(a=x+y\times i,b=x'+y'\times i\) ,计算 \(ab\) 。
我们有两种方法:
- 可以直接展开相乘: \(ab=(x+yi)(x'+y'i)=(xx'-yy')+(xy'+x'y)i\) 。
- 可以从几何角度理解:两个复数相乘等于“幅角相加,模长相乘”(但并不满足向量的运算)
更推荐就用方法一来理解。
单位复数
\(Hint.\) 复数的模长 \(|a|\) (这个不是绝对值),表示把它表示成向量后的长度,即 \(|a|=\sqrt{x^2+y^2}\) 。
我们定义 \(n\) 次单位复根 \(\omega_n^i\) 为满足 \(\text{z}^n=1\) 的复数 \(\text{z}\) 。
本人的理解方法:单位复根 \(\omega_n^i\) 即为将一个单位圆分成 \(n\) 分,而 \(\omega_n^i\) 就是 \(i\) 个这样的角拼起来。特别地, \(\omega_n^0=1\) 。
其实,之后的理解都可以通过 三角函数 有更好的理解
可以发现 \(n\) 次单位复根满足以下性质( \(\omega_n^i\) 表示第 \(i\) 个单位复根):
- \(|\omega_n^i|=1\) (都在单位圆上);
- \(\omega_n^i\) 的幅角为 \(\frac{2\pi i}{n}\) ;\(Tab.\) 这里的 \(i\) 是下标不是虚数,角度为弧度制( \(360^{\circ}\) 弧度角为 \(2\pi\) )
- 总共有 \(n\) 个单位复根,记为 \(\omega_n^0,\omega_n^1,\omega_n^2,...,\omega_n^{n-1}\)
另外,还有一些比较明显的性质:
- \(\omega_n^i\cdot\omega_n^j=\omega_n^{i+j}\) ;
- \(\omega_n^{i+n}=\omega_n^i\) ;
- \(\omega_n^i=\omega_{kn}^{ki}\) ;
- \(\omega_n^i=-\omega_n^{i+n/2}\) (这里的“负”表示大小相等、方向相反)
以上这些结论都可以通过在单位圆上画出单位根来证明。
单位复根有什么用呢?因为 \(n\) 次单位复根恰好有 \(n\) 个,就可以把这 \(n\) 个单位复根分别代入 \(n-1\) 次多项式 \(f(x)\) —— 得到点值表达 \(\{(\omega_n^0,f(\omega_n^0)),(\omega_n^1,f(\omega_n^1)),...,(\omega_n^{n-1},f(\omega_n^{n-1}))\}\) 。
傅立叶正变换
有了这些辅助知识 没错,虽然你可能已经晕了,但是他们真的只是辅助知识 ,我们终于可以进行正题了。
所谓变换,那么一定有正也有逆,现在我我们先来掌握它的正变换。
FFT 的正变换实现,是基于对多项式进行奇偶项分开递归再合并的分治进行的。
对于 n-1 次多项式,我们选择插入 n 次单位根求出其点值表达式。
记多项式 \(A(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}\) 。
再记 \(A_o(x)=a_1+a_3x+a_5x^2+...\) 。
再记 \(A_e(x)=a_0+a_2x+a_4x^2+...\) 。
有 \(A(x)=x*A_o(x^2)+A_e(x^2)\) 。
令 \(n = 2*p\) 。则有:
\(A(w_n^k)=w_n^k*A_o[(w_{n/2}^{k/2})^2]+A_e[(w_{n/2}^{k/2})^2]=w_n^k*A_o(w_p^k)+A_e(w_p^k)\) ;
\(A(w_n^{k+p})=w_n^{k+p}*A_o(w_p^{k+p})+A_e(w_p^{k+p})=-w_n^k*A_o(w_p^k)+A_e(w_p^k)\)
在已知 \(A_o(w_p^k)\) 与 \(A_e(w_p^k)\) 的前提下,可以 \(O(1)\) 算出 \(A(w_n^k)\) 与 \(A(w_n^{k+p})\) 。
因此,假如我们递归求解 \(A_o(x),A_e(x)\) 两个多项式 \(p\) 次单位根的插值,就可以在 \(O(n)\) 的时间内算出 \(A(x)\) n 次单位根的插值。
时间复杂度是经典的 \(T(n)=2*T(n/2)+O(n)=O(n\log n)\)。
傅里叶逆变换
刚刚研究完正的,现在我们来研究逆变换,其实也比较好理解。
观察我们刚刚的插值过程,实际上就是进行了如下的矩阵乘法。
我们记上面的系数矩阵为 \(V\) 。
对于下面定义的 \(D\) :
考虑 \(D*V\) 的结果:
当 \(i = j\) 时, \((D*V)_{ij}=n\) ;
当 \(i ≠ j\) 时, \((D*V)_{ij}=1+w_n^{j-i}+(w_n^{j-i})^2+...=\frac{1-(w_n^{j-i})^n}{1-w_n^{j-i}}=0\) ;
根据定义, \(n\) 次单位根的 \(n\) 次方都等于 \(1\)
所以:\(\frac1n*D=V^{-1}\)
因此将这个结果代入最上面那个公式里面,有:
“这样,逆变换 就相当于把 正变换 过程中的 \(w_n^k\) 换成 \(w_n^{-k}\),之后结果除以 \(n\) 就可以了。”——摘自某博客。
……
还是有点难理解。比如为什么我们不直接把 \(w_n^k\) 换成 \(\frac1n*w_n^{-k}\) 算了。
实际上,因为 \(w_n^{-k}=w_n^{n-k}\) ,也就是说它 TM 还是一个 \(n\) 次单位根。所以我们插值还是正常的该怎么插怎么插。如果换成 \(\frac1n*w_n^{-k}\) 它就不是一个单位根,以上性质就不满足了。
FFT 的代码实现
我们有两个版本——递归、迭代,相信大家都也想到了吧?
毋庸置疑的,递归版本确实很好写,将 \(FFT\) 的思路全部放到上面即可,但是方便也有他的坏处——递归实现常数较大,这样对于一些题就会 \(T\) 飞,但是我还是写了,这里上代码了
#include<cstdio>
#include<cmath>
#define rep(i,__l,__r) for(register int i=__l,i##_end_=__r;i<=i##_end_;++i)
#define fep(i,__l,__r) for(register int i=__l,i##_end_=__r;i>=i##_end_;--i)
#define writc(a,b) fwrit(a),putchar(b)
#define mp(a,b) make_pair(a,b)
#define ft first
#define sd second
#define LL long long
#define ull unsigned long long
#define pii pair<int,int>
#define Endl putchar('\n')
// #define FILEOI
// #define int long long
#ifdef FILEOI
#define MAXBUFFERSIZE 500000
inline char fgetc(){
static char buf[MAXBUFFERSIZE+5],*p1=buf,*p2=buf;
return p1==p2&&(p2=(p1=buf)+fread(buf,1,MAXBUFFERSIZE,stdin),p1==p2)?EOF:*p1++;
}
#undef MAXBUFFERSIZE
#define cg (c=fgetc())
#else
#define cg (c=getchar())
#endif
template<class T>inline void qread(T& x){
char c;bool f=0;
while(cg<'0'||'9'<c)f|=(c=='-');
for(x=(c^48);'0'<=cg&&c<='9';x=(x<<1)+(x<<3)+(c^48));
if(f)x=-x;
}
inline int qread(){
int x=0;char c;bool f=0;
while(cg<'0'||'9'<c)f|=(c=='-');
for(x=(c^48);'0'<=cg&&c<='9';x=(x<<1)+(x<<3)+(c^48));
return f?-x:x;
}
template<class T,class... Args>inline void qread(T& x,Args&... args){qread(x),qread(args...);}
template<class T>inline T Max(const T x,const T y){return x>y?x:y;}
template<class T>inline T Min(const T x,const T y){return x<y?x:y;}
template<class T>inline T fab(const T x){return x>0?x:-x;}
inline int gcd(const int a,const int b){return b?gcd(b,a%b):a;}
inline void getInv(int inv[],const int lim,const int MOD){
inv[0]=inv[1]=1;for(int i=2;i<=lim;++i)inv[i]=1ll*inv[MOD%i]*(MOD-MOD/i)%MOD;
}
template<class T>void fwrit(const T x){
if(x<0)return (void)(putchar('-'),fwrit(-x));
if(x>9)fwrit(x/10);putchar(x%10^48);
}
inline LL mulMod(const LL a,const LL b,const LL mod){//long long multiplie_mod
return ((a*b-(LL)((long double)a/mod*b+1e-8)*mod)%mod+mod)%mod;
}
const int MAXN=3e6;
const double Pi=acos(-1.0);
class task{
private:
struct cplx{
double vr,vi;//实部和虚部
cplx(const double R=0,const double I=0):vr(R),vi(I){}//构造函数
//------------------overload----------------//
cplx operator + (const cplx a)const{return cplx(vr+a.vr,vi+a.vi);}//重载加法
cplx operator - (const cplx a)const{return cplx(vr-a.vr,vi-a.vi);}
cplx operator * (const cplx a)const{return cplx(vr*a.vr-vi*a.vi,vr*a.vi+a.vr*vi);}
cplx operator / (const double var)const{return cplx(vr/var,vi/var);}
};
int n,m;
cplx a[MAXN+5],b[MAXN+5];
void fft(cplx* f,const int len,const short opt=1){
//opt==-1 : FFT 的逆变换
if(!len)return;
cplx f0[len+5],f1[len+5];
for(int i=0;i<len;++i)
f0[i]=f[i<<1],f1[i]=f[i<<1|1];
fft(f0,len>>1,opt);
fft(f1,len>>1,opt);
cplx w=cplx(cos(Pi/len),opt*sin(Pi/len)),buf=cplx(1,0);
for(int i=0;i<len;++i,buf=buf*w){
f[i]=f0[i]+buf*f1[i];
f[i+len]=f0[i]-buf*f1[i];
}
}
public:
inline void launch(){
qread(n,m);
rep(i,0,n)scanf("%lf",&a[i].vr);
rep(i,0,m)scanf("%lf",&b[i].vr);
for(m+=n,n=1;n<=m;n<<=1);
fft(a,n>>1);
fft(b,n>>1);
rep(i,0,n-1)a[i]=a[i]*b[i];
fft(a,n>>1,-1);
rep(i,0,m)writc((int)((a[i].vr)/n+0.5),' ');
Endl;
}
}This;
signed main(){
#ifdef FILEOI
freopen("file.in","r",stdin);
freopen("file.out","w",stdout);
#endif
This.launch();
return 0;
}
似乎递归版本比较好写,现在我们来看一下迭代(递推)版本应该怎么做:
原序列: \(0,1,2,3,4,5,6,7\)
终序列: \(0,4,2,6,1,5,3,7\)
转换为二进制再来看看。
原序列: \(000,001,010,011,100,101,110,111\)
终序列: \(000,100,010,110,001,101,011,111\)
可以发现终序列是原序列每个元素的二进制翻转。
于是我们可以先把要变换的系数排在相邻位置,从下往上迭代。
这个二进制翻转过程可以自己脑补方法,只要保证时间复杂度 \(O(n\log n)\) ,代码简洁就可以了。
在这里给出一个参考的方法:
我们对于每个 \(i\) ,假设已知 \(i-1\) 的翻转为 \(j\) 。考虑不进行翻转的二进制加法怎么进行:从最低位开始,找到第一个为 \(0\) 的二进制位,将它之前的 \(1\) 变为 \(0\) ,将它自己变为 \(1\) 。因此我们可以从 \(j\) 的最高位开始,倒过来进行这个过程。
这是迭代版本:
#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;
#define rep(i,__l,__r) for(register int i=__l,i##_end_=__r;i<=i##_end_;++i)
#define fep(i,__l,__r) for(register int i=__l,i##_end_=__r;i>=i##_end_;--i)
#define writc(a,b) fwrit(a),putchar(b)
#define mp(a,b) make_pair(a,b)
#define ft first
#define sd second
#define LL long long
#define ull unsigned long long
#define pii pair<int,int>
#define Endl putchar('\n')
// #define FILEOI
// #define int long long
#ifdef FILEOI
#define MAXBUFFERSIZE 500000
inline char fgetc(){
static char buf[MAXBUFFERSIZE+5],*p1=buf,*p2=buf;
return p1==p2&&(p2=(p1=buf)+fread(buf,1,MAXBUFFERSIZE,stdin),p1==p2)?EOF:*p1++;
}
#undef MAXBUFFERSIZE
#define cg (c=fgetc())
#else
#define cg (c=getchar())
#endif
template<class T>inline void qread(T& x){
char c;bool f=0;
while(cg<'0'||'9'<c)f|=(c=='-');
for(x=(c^48);'0'<=cg&&c<='9';x=(x<<1)+(x<<3)+(c^48));
if(f)x=-x;
}
inline int qread(){
int x=0;char c;bool f=0;
while(cg<'0'||'9'<c)f|=(c=='-');
for(x=(c^48);'0'<=cg&&c<='9';x=(x<<1)+(x<<3)+(c^48));
return f?-x:x;
}
template<class T,class... Args>inline void qread(T& x,Args&... args){qread(x),qread(args...);}
template<class T>inline T Max(const T x,const T y){return x>y?x:y;}
template<class T>inline T Min(const T x,const T y){return x<y?x:y;}
template<class T>inline T fab(const T x){return x>0?x:-x;}
inline int gcd(const int a,const int b){return b?gcd(b,a%b):a;}
inline void getInv(int inv[],const int lim,const int MOD){
inv[0]=inv[1]=1;for(int i=2;i<=lim;++i)inv[i]=1ll*inv[MOD%i]*(MOD-MOD/i)%MOD;
}
template<class T>void fwrit(const T x){
if(x<0)return (void)(putchar('-'),fwrit(-x));
if(x>9)fwrit(x/10);putchar(x%10^48);
}
inline LL mulMod(const LL a,const LL b,const LL mod){//long long multiplie_mod
return ((a*b-(LL)((long double)a/mod*b+1e-8)*mod)%mod+mod)%mod;
}
const int MAXN=3e6;
const double Pi=acos(-1.0);
class task{
private:
struct cplx{
double vr,vi;//实部和虚部
cplx(const double R=0,const double I=0):vr(R),vi(I){}//构造函数
//------------------overload----------------//
cplx operator + (const cplx a)const{return cplx(vr+a.vr,vi+a.vi);}//重载加法
cplx operator - (const cplx a)const{return cplx(vr-a.vr,vi-a.vi);}
cplx operator * (const cplx a)const{return cplx(vr*a.vr-vi*a.vi,vr*a.vi+a.vr*vi);}
cplx operator / (const double var)const{return cplx(vr/var,vi/var);}
};
int n,m;
cplx a[MAXN+5],b[MAXN+5];
int revi[MAXN+5];
/*
f(w^x)
=f0(w^{2x})+w^x*f1(w^{2x})
=a0+a2+a4+...,a1+a3+a5+...
=a0+a4+a8+...+a2+a6+a10...+a1+a5+a9+...+a3+a7+a11...
=f00(w^{4x})+w^{2x}*f01(w^{4x})+w^x*f10(w^{4x})+w^3x*f11(w^{4x})
=a0+a4+a8+...,w^{2x}*(a2+a6+a10...),w^x*(a1+a5+a9+...),w^{3x}*(a3+a7+a11...)
f_s -> 下标将 s 反过来之后, a_i 的 i 的二进制反过来与 s 相同
.
.
.
s
000 001 010 011
|反过来
v
000 100 010 110
a0 a{k/2} a{3*k/4} a{k}
a0 a4 a2 a6
*/
/*
void fft(cplx* f,const int len,const short opt=1){
if(!len)return;
cplx f0[len+5],f1[len+5];
for(int i=0;i<len;++i)
f0[i]=f[i<<1],f1[i]=f[i<<1|1];
fft(f0,len>>1,opt);
fft(f1,len>>1,opt);
cplx w=cplx(cos(Pi/len),opt*sin(Pi/len)),buf=cplx(1,0);
for(int i=0;i<len;++i,buf=buf*w){
f[i]=f0[i]+buf*f1[i];
f[i+len]=f0[i]-buf*f1[i];
}
}
*/
inline void fft(cplx* f,const short opt=1){
for(int i=0;i<n;++i)if(i<revi[i])
swap(f[i],f[revi[i]]);
for(int p=2;p<=n;p<<=1){
//枚举层数
int len=p/2;//上一层的一半的长度
cplx tmp(cos(Pi/len),opt*sin(Pi/len));
//单位复根
for(int k=0;k<n;k+=p){
cplx buf(1,0);//记录 omega 的次方
for(int l=k;l<k+len;++l,buf=buf*tmp){
//每次 buf 累成单位 omega
cplx tt=buf*f[len+l];
f[len+l]=f[l]-tt;
f[l]=f[l]+tt;
/*
此处与递归版本的 FFT 中这一段是一样的:
f[i]=f0[i]+buf*f1[i];
f[i+len]=f0[i]-buf*f1[i];
*/
}
}
}
if(opt==-1)for(int i=0;i<n;++i)f[i]=f[i]/n;
//如果是 逆变换 , 那么需要全部 /n
}
public:
inline void launch(){
qread(n,m);
rep(i,0,n)scanf("%lf",&a[i].vr);
rep(i,0,m)scanf("%lf",&b[i].vr);
for(m+=n,n=1;n<=m;n<<=1);
rep(i,0,n-1)revi[i]=(revi[i>>1]>>1)|((i&1)?n>>1:0);//处理反转
fft(a);fft(b);
rep(i,0,n-1)a[i]=a[i]*b[i];
fft(a,-1);
rep(i,0,m)writc((int)(a[i].vr+0.5),' ');
Endl;
}
}This;
signed main(){
#ifdef FILEOI
freopen("file.in","r",stdin);
freopen("file.out","w",stdout);
#endif
This.launch();
return 0;
}
如果我可以回到过去,我一定会去当杀手。
为什么?因为我要去干翻欧某和傅某某...