快速傅里叶变换裸体,总结一下快速傅里叶变换的核心思想。
首先我们要搞明白的就是快速傅里叶变换的思想就是先将一个多项式转化为点-值表达式(DFT),然后再将点值表达式逆变换为多项式(IDFT),快速傅里叶变换就是用来优化这个过程的,就是让我们快速求出来一个点值表达式,然后再快速的变换回去。那么我们看如何做。首先我们要带入一些x才能解出y得到点值表达式,那么我们就要想一些合适的x能快速的计算,这里我们引入单位根。单位根就是w^n=1的n个不同的解,不懂得读者可以百度一下复数,这里不再赘述。它有很多非常有用的性质,由于不好打所以这里也不再赘述……这里往下我可能就说不明白了,可以参考picks的博客。我们首先要明确的是我们要求的是什么,我们要求的是一个n项序列,第i项代表多项式在第i个单位根下的解,那么我们要求的东西可以根据奇数项和偶数项划分出来算出来在n/2下的单位根的解,然后只要最后将答案综合一下就可以了。(好像说的不是很明白,不懂的还是找picks的博客吧)这样的话我们就能得到一个递归的版本来解决快速傅里叶变换,但是递归的常数巨大,我们这里给出一个非递归的版本,答题思想就是我们通过数学办法找到层层递归后到达最后一层时所有元素的位置,然后一层一层的向上反,经研究我们发现一个元素新的位置恰好就是将它的二进制倒过来。这样DFT就没什么问题了,然后就是IDFT的问题了,我们这里引入一个结论就是将求出的点值表达式重做一遍快速傅里叶变换,但是每次乘的复数变为原来的共轭复数,重新得到系数表达式,但是每一项重复计算了len次,所以每一项要除以len,这样就可以输出了。
#include<cstdio>
#include<cstdlib>
#include<iostream>
#include<iomanip>
#include<cstring>
#include<string>
#include<ctime>
#include<cmath>
#include<algorithm>
using namespace std;
struct cpx
{
double a,b;
cpx(double _=0.0,double __=0.0):a(_),b(__){}
cpx operator +(cpx x) {return cpx(a+x.a,b+x.b);}
cpx operator -(cpx x) {return cpx(a-x.a,b-x.b);}
cpx operator *(cpx x) {return cpx(a*x.a-b*x.b,a*x.b+b*x.a);}
}a[200000],b[200000],c[200000];
const double DFT=2.0;
const double IDFT=-2.0;
double trans_form;
const double pi=acos(-1);
int len;
int pos[200000];
void init()
{
for(int i=0;i<len;i++)
{
pos[i]=pos[i>>1]>>1;
if(i&1) pos[i]|=(len>>1);
}
}
void trans(cpx x[])
{
for(int i=0;i<len;i++) if(i<pos[i]) swap(x[i],x[pos[i]]);
for(int i=2;i<=len;i<<=1)
{
int step=i>>1;
cpx wm(cos(2*pi/(double)i),sin(trans_form*pi/(double)i));
for(int j=0;j<len;j+=i)
{
int limit=j+step;
cpx ww(1,0);
for(int k=j;k<limit;k++)
{
cpx a=x[k];
cpx b=x[k+step]*ww;
ww=ww*wm;
x[k]=a+b;
x[k+step]=a-b;
}
}
}
if(trans_form==IDFT) for(int i=0;i<len;i++) x[i].a/=(double)len;
}
char s[200000];
int ans[200000];
int main()
{
int n;
scanf("%d",&n);
len=1;
while(len<(n<<1)) len<<=1;
init();
scanf("%s",s);
for(int i=0;i<n;i++) a[i].a=s[i]-'0';
scanf("%s",s);
for(int i=0;i<n;i++) b[i].a=s[i]-'0';
trans_form=DFT;
trans(a);
trans(b);
for(int i=0;i<len;i++) c[i]=a[i]*b[i];
trans_form=IDFT;
trans(c);
int top=0;
for(int i=(n-1)<<1;~i;i--)
{
ans[top++]+=int(c[i].a+0.5);
ans[top]+=ans[top-1]/10;
ans[top-1]=ans[top-1]%10;
}
if(ans[top]) cout<<ans[top];
while(top)
{
printf("%d",ans[top-1]);
top--;
}
return 0;
}