快速傅里叶变换 FFT
前置知识：离散时间傅里叶变换
和一些复变函数的小知识
离散时间傅里叶变换 & 离散傅里叶级数
离散时间傅里叶变换定义式：\(X(j\Omega)=\sum\limits^{\infty}_{n=-\infty}x_n e^{-j\Omega n}\)
性质：时域上的卷积等价于经过傅里叶变换后在频域上做乘法然后反变换回来
考虑到一个离散信号可以写成一个形式幂级数，那么两个信号的卷积自然对应多项式乘法
那么就可以通过傅里叶变换的方式算多项式乘法
离散时间傅里叶级数定义式：\(\tilde{X_k}=\sum\limits^{N}_{n=0}x_n e^{-j\frac{2\pi}{N} k}\)
离散傅里叶变换
离散傅里叶级数要求有周期并且长度有限
一般的离散信号没有周期，但是长度有限是满足的，但是值域是连续的
考虑抽样的过程，一个域抽样在另一个域表示周期延拓
所以我们可以很自然的把有限信号周期延拓一下，这样可以用离散傅里叶级数的计算公式了
这么做就相当于是在连续的频域上抽样，
所以我们接下来就可以利用抽样和恢复的技术把这个信号的频域离散下来存起来了
解释完原理之后我们直接照搬离散傅里叶级数的计算公式就可以
\(\tilde{X_k}=\sum\limits^{N}_{n=0}x_n e^{-j\frac{2\pi}{N} k}\)
这样卷积就可以变成乘法，做乘法的时间显然是 \(O(n)\) 的
但是这样朴素的过程，时间复杂度显然是 \(O(n^2)\) 考虑优化
快速傅里叶变换（递归）
注意到 \(e^{-j\frac{2\pi}{N} k}=-\omega ^{k}_{n}=\cos {\frac{2\pi}{N} k}-j\frac{2\pi}{N} k\)
而且还有一个 \(\omega^{k}_{2n}=-\omega^{k+n}_{2n}\)
那么我们其实悄无声息的凑出了一个 \((-1)^n\) 这个项
对于高中数列学的比较好的估计想到了凑奇偶项
假设所有偶数项构成 \(G(x)=a_0+a_2x^2+a_4x^4+\cdots\) 奇数项构成 \(H(x)=a_1x+a_3x^3+a_5x^5+\cdots\)
则有 \(F(x)=G(x)+xH(x)\)
对于 \(F(\omega^{k}_{2n})=G(\omega^{k}_{2n})+\omega^{k}_{2n}H(\omega^{k}_{2n})\)
对于 \(F(\omega^{k}_{2n})=G(\omega^{k+n}_{2n})+\omega^{k+n}_{2n}H(\omega^{k+n}_{2n})=G(\omega^{k}_{2n})-\omega^{k}_{2n}H(\omega^{k}_{2n})\)
并且还有 \(G(x),H(x)\) 也可以依照 \(F(x)\) 做处理
而且性质更好的是，\(G(x)\) 和 \(H(x)\) 显然对应 \(F(x^2)\)，也就是下角标除了个 \(2\)
也就是至多会递归 \(\log n\) 层
时间复杂度 \(O(n\log n)\)
代码直接模拟上述过程即可，
由于每一层递归处理的项数不变，所以每一层看似函数的数量增长很多，实际上是 \(O(n)\) 的
然后要注意的是，这个分治其实一直要求下角标是偶数，所以一定要把长度补成 \(2\) 的幂次
递归版FFT实现多项式乘法
#include <bits/stdc++.h>
using namespace std;
#define fu complex<double>
const fu I=(0,1);
const double pi=acos(-1),e=exp(1);
vector<fu> FFT(vector<fu>a,int val){
int N=a.size();
vector<fu> F(N);
if(a.size()==1)return a;
vector<fu> g,h;
for(int i=0;i<N;i++){
if(i&1)h.push_back(a[i]);
else g.push_back(a[i]);
}
fu now={1,0},omega={cos(2*pi/N),val*sin(2*pi/N)};
vector<fu> G=FFT(g,val),H=FFT(h,val);
for(int i=0;i<N/2;i++,now*=omega){
F[i]=G[i]+now*H[i];
F[i+N/2]=G[i]-now*H[i];
}
return F;
}
int read(){
int i=1,j=0;
char ch=getchar();
while(ch>'9'||ch<'0'){
if(ch=='-')i=-1;
ch=getchar();
}
while(ch>='0'&&ch<='9'){
j=j*10+ch-48;
ch=getchar();
}
return i*j;
}
int lg(int x){
if(x==1)return 0;
return lg(x/2)+1;
}
int main(){
int n=read(),m=read();
int len=(1<<lg(n+m+2))==(n+m+2)?(1<<lg(n+m+2)):(1<<(lg(n+m+2)+1));
vector<fu>a(len),b(len),c;
for(int i=0;i<=n;i++)a[i]=1.0*read();
for(int i=0;i<=m;i++)b[i]=1.0*read();
vector<fu> F=FFT(a,-1),G=FFT(b,-1),H(len);
for(int i=0;i<len;i++)H[i]=F[i]*G[i];
c=FFT(H,1);
for(int i=0;i<n+m+1;i++)cout<<(int)(c[i].real()/len+0.5)<<" ";
return 0;
}
优化
第一个优化是位逆序排列优化
在递归的过程，涉及到 \(G\) 和 \(H\) 的重新排列和向下递归
这个优化主要是针对重排到单点时的系数
具体而言：递归到底时，该位的系数在原序列的位置，是这个位的二进制倒写
证明考虑奇偶项分类的过程，实际上是将最后一位分开了
那么末尾取 \(0\) 向左递归 末位取 \(1\) 向右递归的过程，
实际上是根据当前最后一位的结果确定最高位的结果，这两位是完全相同的
把这一位删掉，递归地考虑这个过程即得。
所以可以实现预处理每一位最后的位置，然后向上合并的过程中更新
这个预处理的过程可以直接 \(O(n)\) 做完
然后是额外的 \(O(1)\) 空间更新
位逆序排列优化完了之后，相邻的两项就是一个 \(G\) 和 一个 \(H\)
同时这个 \(G\) 和 \(H\) 我们接下来也用不到
同时这两个位置恰好也是 \(F\) 的位置
所以实际上没有必要额外开空间就能做了
非递归版FFT实现多项式乘法
#include <bits/stdc++.h>
using namespace std;
#define fu complex<double>
const int o=2222222;
const double pi=acos(-1);
int r[o];
fu a[o],b[o];
void pre(int n){
int x=n/2,p=0;
r[p++]=0,r[p++]=x;
for(int i=2;i<=n;i*=2){
x/=2;
for(int j=0;j<i;j++)r[p++]=r[j]|x;
}
}
void FFT(fu *a,int n,int val){
for(int i=0;i<n;i++){
if(i<r[i])swap(a[i],a[r[i]]);
}
for(int i=2;i<=n;i*=2){
for(int j=0;j<n;j+=i){
fu now(1,0),omega(cos(2*pi/i),val*sin(2*pi/i));
for(int k=j;k<j+i/2;k++,now*=omega){
fu x=a[k],y=a[k+i/2];
a[k]=x+now*y;
a[k+i/2]=x-now*y;
}
}
}
}
int read(){
int i=1,j=0;
char ch=getchar();
while(ch>'9'||ch<'0'){
if(ch=='-')i=-1;
ch=getchar();
}
while(ch>='0'&&ch<='9'){
j=j*10+ch-48;
ch=getchar();
}
return i*j;
}
int lg(int x){
if(x==1)return 0;
return lg(x/2)+1;
}
int main(){
int n=read(),m=read();
int len=(1<<lg(n+m+2))==(n+m+2)?(1<<lg(n+m+2)):(1<<(lg(n+m+2)+1));
pre(len);
for(int i=0;i<=n;i++)a[i]={1.0*read(),0};
for(int i=0;i<=m;i++)b[i]={1.0*read(),0};
FFT(a,len,-1),FFT(b,len,-1);
for(int i=0;i<len;i++)a[i]=a[i]*b[i];
FFT(a,len,1);
for(int i=0;i<n+m+1;i++)printf("%d ",(int)(a[i].real()/len+0.5));
return 0;
}
浙公网安备 33010602011771号