Fafo.cin

// Discrete Approcimation of Fourier operator
 /* //Example of headers
#include<math.h>
#include<stdio.h>
#include <stdlib.h>
#include <complex>
using namespace std;
#define z_type complex<double>
*/
// The "miminal" Fast Fourier Transform by Dmitrii Kouznetsov
// do not forget to #define z_type complex<double>
// n should be 2^m ; o should be 1 or -1 ;

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;}
}

References

The routine is copypasted from http://tori.ils.uec.ac.jp/TORI/index.php/Fafo.cin