#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#define M_PI 3.14159265358979323846                                             
                                                                                
struct complex {                                                                  
  double re,im;                                                                 
};                                                                              
                                                                                
void fft(int dir, int e, struct compl *x)                                       
{                                                                               
  unsigned  k, l;                                                                  
  unsigned  i, j, ll, n, nn;                                                  
  double c, p, tr, ti, wr, wi, ur, ui;                                       
  struct complex w, *a, *b;                                                   
                                                                                
  dir&=1;                                                                       
  n=1<<e;                                                                       
  nn=n/2;                                                                       
  for (k = ll = 0; k < (n - 1); k++) {                                                     
    if (ll > k) w = x[k], x[k] = x[ll], x[ll] = w;                                    
    l = nn;                                                                     
    while (l <= ll) ll -= l, l >>= 1;                                               
    ll += l;                                                                    
  }                                                                          
  p = (dir == 0) ? (-M_PI) : (M_PI);                                                
  for (i = l = 1; i <= e; i++) {                                                      
    ll = l;  
    l <<= 1;                                                             
    ur = 1.0;
    ui = 0.0;                                                   
    wr = cos(p / ll);
    wi = sin(p / ll);                                             
    for (j = 0; j < ll; j++) {                                                     
      for (k=j, b = x + k, a = x + k + ll; k < n; k += l, b += l, a += l) {                     
        tr = ur * a->re - ui * a->im;                                           
        ti = ur * a->im + ui * a->re;                                           
        a->re = b->re - tr;                                                 
        a->im = b->im - ti;                                                 
        b->re += tr;                                                        
        b->im += ti;                                                        
      } /* for k */                                              
      tr = ur * wr - ui * wi;                                                        
      ui = ui * wr + ur * wi;                                                        
      ur = tr;                                                                 
    } /* for j */                                                         
  } /* for i */                                                             
  if (dir != 0) { 
    c = (1.0 / n);
    for(i = 0; i < n; i++) { 
      x[i].re *= c;
      x[i].im *= c;
    } 
  }              
}