Fast Fourier Transform

来源：互联网发布：人工智能外文参考文献编辑：程序博客网时间：2024/06/05 23:01

/* * Fast Fourier Transform * By Pan Yuchong */#include <cstdio>#include <cmath>#include <complex>#define maxN 1000000#define pi 3.14159265358979323846#define re real()using namespace std;int n,m,l,N,L,rev[maxN+10],dig[maxN+10];double x;complex<double> a[maxN+10],b[maxN+10],c[maxN+10],A[maxN+10],B[maxN+10],C[maxN+10];void FFT(complex<double> A[],int f){for (int i=2;i<=N;i<<=1){complex<double> Wn(cos(2*pi/i),f*sin(2*pi/i));for (int k=0;k<N;k+=i){complex<double> W(1);for (int j=0;j<i/2;++j){complex<double> x=W*A[k+j+i/2],y=A[k+j];A[k+j]=y+x,A[k+j+i/2]=y-x;W*=Wn;}}}if (f==-1) for (int i=0;i<N;++i) A[i].re/=N;}int main(){scanf("%d%d",&n,&m);for (int i=0;i<n;++i) scanf("%lf",&x),a[i].re=x;for (int i=0;i<m;++i) scanf("%lf",&x),b[i].re=x;for (N=1;N<n||N<m;N<<=1,L++);N<<=1,L++;for (int i=0;i<N;++i){l=0;for (int t=i;t;t>>=1) dig[l++]=t&1;for (int j=0;j<L;++j) rev[i]=rev[i]*2+dig[j];A[i]=a[rev[i]],B[i]=b[rev[i]];}FFT(A,1),FFT(B,1);for (int i=0;i<N;++i) c[i]=A[i]*B[i];for (int i=0;i<N;++i) C[i]=c[rev[i]];FFT(C,-1);for (int i=0;i<N;++i) printf("%lf\n",C[i].re);return 0;}

By Charlie Pan

Apr 30,2014

0 0