Fafo.cin

// Discrete Approcimation of Fourier operator /* //Example of headers using namespace std; */ // The "miminal" Fast Fourier Transform by Dmitrii Kouznetsov // do not forget to #define z_type complex // n should be 2^m ; o should be 1 or -1 ;
 * 1) include
 * 2) include
 * 3) include 
 * 4) include
 * 1) define z_type complex

void fft(z_type *a, int N, int o) // Arry is FIRST! {z_type u,w,t; int i,j,k,l,e=1,L,p,q,m; q=N/2; p=2; for(m=1;p<N;m++) p*=2; p=N-1; z_type y=z_type(0.,M_PI*o); j=0; for(i=0;i<p;i++){ if(i<j) { t=a[j]; a[j]=a[i]; a[i]=t;} k=q; while(k<=j){ j-=k; k/=2;} j+=k; } for(l=0;l<m;l++) { L=e; e*=2; u=1.; w=exp(y/((double)L)); for(j=0;j<L;j++) { for(i=j;i<N;i+=e){k=i+L; t=a[k]*u; a[k]=a[i]-t; a[i]+=t;} u*=w; } } }

void fafo(z_type *b,int N,int o) { int j; z_type c; double q=sqrt(1./N); for(j=0;j<N/2;j++) {c=b[j];b[j]=b[j+N/2];b[j+N/2]=c;} fft(b,N,o); for(j=0;j<N/2;j++) {c=b[j];b[j]=b[j+N/2]*q;b[j+N/2]=c*q;} }