快速傅里葉變換(FFT)蝴蝶迭代

#include<cstdio>
#include<cmath>
#define Maxn 1000500
const double Pi=acos(-1);
using namespace std;
int n,m,r[Maxn<<2];
struct complex
{complex(double xx=0,double yy=0){x=xx,y=yy;}
double x,y;};
complex operator + (complex a,complex b){
  return complex(a.x+b.x,a.y+b.y);
}complex operator - (complex a,complex b){
  return complex(a.x-b.x,a.y-b.y);
}complex operator * (complex a,complex b){
  return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);
}complex b[Maxn<<2],c[Maxn<<2];
void fft(complex *f,short op)
{
  for (int i=0;i<n;i++)
   if (i<r[i]){
     complex tmp=f[i];
     f[i]=f[r[i]];
     f[r[i]]=tmp;
   }
  for(int p=2;p<=n;p<<=1){
    int len=p/2;
    complex tmp(cos(Pi/len),op*sin(Pi/len));
    for(register int k=0;k<n;k+=p){
      register complex buf(1,0);
      for(register int l=k;l<k+len;l++){
        register complex tt=buf*f[len+l] ;
        f[len+l]=f[l]-tt;
        f[l]=f[l]+tt;
        buf=buf*tmp;
      }
    }
  }
}
int main()
{
  scanf("%d%d",&n,&m);
  for (int i=0;i<=n;i++)scanf("%lf",&b[i].x);
  for (int i=0;i<=m;i++)scanf("%lf",&c[i].x);
  for(m+=n,n=1;n<=m;n<<=1);
  for(int i=0;i<n;i++)
   r[i]=(r[i>>1]>>1)|((i&1)?n>>1:0);
  fft(b,1);
  fft(c,1);//DFT
  for(int i=0;i<n;++i)b[i]=b[i]*c[i];
  fft(b,-1);//IDFT
  for(int i=0;i<=m;++i)printf("%.0f ",fabs(b[i].x)/n);
  return 0;
}