#include <math.h>

int cfft( float x[], float y[], int n, int m, int itype )
{
   const float TWO_PI = 6.283185307179586476925287;
   register int i, j, k, l, n1, n2, ie, ia;
   float c, s, xt, yt, p, a;

// Quick return for n=1

   if ( n == 1 )
   {
      return 0;
   }

// Conjugate if backward transform

   if ( itype == - 1 )
   {
      for ( i = 0; i < n; i++ )
      {
         y[i] = - y[i];
      }
   }

// Main loop

   p = TWO_PI / float( n );
   n2 = n;
   for ( k = 1; k <= m; k++ )
   {
      n1 = n2;
      n2 /= 2;
      ie = n / n1;
      ia = 1;
      for ( j = 0; j < n2; j++ )
      {
         a = float( ia - 1 ) * p;
         c = cos( a );
         s = sin( a );
         ia += ie;
         for ( i = j; i < n; i += n1 )
         {
	    l = i + n2;
	    xt = x[i] - x[l];
	    x[i] += x[l];
	    yt = y[i] - y[l];
	    y[i] += y[l];
	    x[l] = c * xt + s * yt;
	    y[l] = c * yt - s * xt;
	 }
      }
   }

// Bit reversal permutation

   j = 0;
   for ( i = 0; i < n-1; i++ )
   {
      if ( i < j )
      {
         xt = x[j];
         x[j] = x[i];
         x[i] = xt;
         yt = y[j];
         y[j] = y[i];
         y[i] = yt;
      }
      k = n / 2;
      while ( k < j+1 )
      {
         j -= k;
         k /= 2;
      }
      j += k;
   }

// Conjugate and normalize if backward transform

   if ( itype == - 1 )
   {
      p = float( n );
      {
         for ( i = 0; i < n; i++ )
         {
            y[i] /= - p;
            x[i] /= p;
         }
      }
   }

   return 0;
}