#include <complex.h>
typedef double _Complex c;

void fft(const int n, c a[n]){
  {int j=0, i; for(i=0; i<n-1; ++i) {
    if(i<j) {c t=a[j]; a[j]=a[i]; a[i]=t;}
    {int m=n/2; while(m<=j) {j -= m; m /= 2;} j += m;}}} 
  {int j=1, m; c b=-1; while(j<n) {c w = 1;
    for(m=0; m<j; ++m) {int i;
      for(i=m; i<n; i+=j+j)
        {c t = w*a[i+j]; a[i+j]=a[i]-t; a[i] += t;}
      w *= b;}
    j *= 2; b = csqrt(b);}}}

// tests
#include <stdio.h>
#include <math.h>
void pc(c a){printf("%19.16f + %19.16fi\n",
    creal(a),  cimag(a));}
int main(){c A[1024]; int j;
typedef double R;
R ep = 2.8e-14;
for(j=0; j<1024; ++j) A[j]=0;
A[11]=1;
fft(1024, A);
if(0) {int j; for(j = 0; j<1024; ++j) pc(A[j]);}
{ int j=1024; while(j--) {R a = 2*(11*j&(1024-1))*3.14159265358979324/1024;
  if(fabs(creal(A[j]) - cos(a)) > ep
  || fabs(cimag(A[j]) - sin(a)) > ep) {
    printf("er: %19.16f + %19.16fi %4d\n    ", cos(a), sin(a), j);
     pc(A[j]);}}}
return 0;}
/* ep=2.8e-14
er:  0.9977230666441916 + -0.0674439195636643i 1023
     0.9977230666441635 + -0.0674439195636614i
er: -0.9977230666441916 +  0.0674439195636642i  511
    -0.9977230666441635 +  0.0674439195636614i
*/