模板题。
大概理解了fft的过程,但是代码还不是很熟练,只能理解比较简单的模板,大佬们的代码都看不懂qaq。
具体过程这篇博客讲的还是挺好的。
http://blog.csdn.net/under_sky_dxj/article/details/52778350
学FFT之前首先得学会DFT,DFT的作用是将时域转换到频域,其实就是将数组x[m]转换到数组X[m]的一个公式。而IDFT可以将X[m]转换回x[m]。对于数组a,b,求出其对应的A,B数组后,将其对应项相乘就对应A和B的卷积,但是这时候还需要IDFT将卷积转换回来,就得到a*b的卷积了。但是这个过程是O(n^2)的,因为DFT的时候x[i]的转换需要n次的计算。
FFT就是在这基础加速求卷积的一个算法。由某些(并不知道)性质得
A[m]=A1[m]+W[m]*A2[m],(W[m]是螺旋因子)
A[m+N/2]=A1[m]-W[m]*A2[m].
A1,A2是将a按下标奇偶分组后DFT求出的,A1对应奇数。
于是我们只要求出A1,A2就可以在线性时间内求出A。
设给长度为n的序列左DFT的复杂度为T(n),则有:
T(n)=2T(n/2)+O(n);
根据主定理T(n)=O(nlogn)。
(这个过程貌似叫蝶形运算)
这样就做DFT的复杂度就从O(n^2)降到O(nlogn)了。
接下来只要相乘然后再做一次IDFT就可以求出卷积。
相乘的过程貌似有精度优化,不是很懂原理,所以目前看不懂,以后再补。
写出来就是给自己看看,大佬要是万一看见了欢迎吐槽。
代码:
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
#define N 1000005
#define pi acos(-1.0) // PI值
using namespace std;
struct complex
{
double r,i;
complex(double real=0.0,double image=0.0){
r=real; i=image;
}
// 以下为三种虚数运算的定义
complex operator + (const complex o){
return complex(r+o.r,i+o.i);
}
complex operator - (const complex o){
return complex(r-o.r,i-o.i);
}
complex operator * (const complex o){
return complex(r*o.r-i*o.i,r*o.i+i*o.r);
}
}x1[N],x2[N];
double z[N];
char a[N/2],b[N/2];
int sum[N]; // 结果存在sum里
void brc(complex *y,int l) // 二进制平摊反转置换 O(logn)
{
register int i,j,k;
for(i=1,j=l/2;i<l-1;i++)
{
if(i<j) swap(y[i],y[j]); // 交换互为下标反转的元素
// i<j保证只交换一次
k=l/2;
while(j>=k) // 由最高位检索,遇1变0,遇0变1,跳出
{
j-=k;
k/=2;
}
if(j<k) j+=k;
}
}
void fft(complex *y,int l,double on) // FFT O(nlogn)
// 其中on==1时为DFT,on==-1为IDFT
{
register int h,i,j,k;
complex u,t;
brc(y,l); // 调用反转置换
for(h=2;h<=l;h<<=1) // 控制层数
{
// 初始化单位复根
complex wn(cos(on*2*pi/h),sin(on*2*pi/h));
for(j=0;j<l;j+=h) // 控制起始下标
{
complex w(1,0); // 初始化螺旋因子
for(k=j;k<j+h/2;k++) // 配对
{
u=y[k];
t=w*y[k+h/2];
y[k]=u+t;
y[k+h/2]=u-t;
w=w*wn; // 更新螺旋因子
} // 据说上面的操作叫蝴蝶操作…
}
}
if(on==-1) for(i=0;i<l;i++) y[i].r/=l; // IDFT
}
int main(void)
{
int l1,l2,l;
register int i;
int n, m;
while(scanf("%d%d",&n,&m)!=EOF)
{
l1=n+1;
l2=m+1;
l=1;
while(l<l1*2 || l<l2*2) l<<=1; // 将次数界变成2^n
// 配合二分与反转置换
for(i=0;i<=n;i++) // 倒置存入
{
scanf("%lf", &x1[i].r);
x1[i].i=0.0;
}
for(;i<l;i++) x1[i].r=x1[i].i=0.0;
// 将多余次数界初始化为0
for(i=0;i<=m;i++)
{
scanf("%lf", &x2[i].r);
x2[i].i=0.0;
}
for(;i<l;i++) x2[i].r=x2[i].i=0.0;
fft(x1,l,1); // DFT(a)
fft(x2,l,1); // DFT(b)
for(i=0;i<l;i++) x1[i]=x1[i]*x2[i]; // 点乘结果存入a
fft(x1,l,-1); // IDFT(a*b)
printf("%.0f", x1[0].r+0.1);
for(i=1; i<l1+l2-1; i++)
{
printf(" %.0f", x1[i].r+0.1);
}
printf("\n");
/*
for(i=0;i<l;i++) sum[i]=x1[i].r+0.5; // 四舍五入
for(i=0;i<l;i++) // 进位
{
sum[i+1]+=sum[i]/10;
sum[i]%=10;
}
l=l1+l2-1;
while(sum[l]<=0 && l>0) l--; // 检索最高位
for(i=l;i>=0;i--) putchar(sum[i]+'0'); // 倒序输出
putchar('\n');
*/
}
return 0;
}