Fast Fourier Transform

/*
 * 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