#include <complex.h>
#define LN 13
#define N (1<<LN)
typedef double _Complex c;
typedef long double _Complex cx;
void pc(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; cx b=-1; while(j<n) {cx 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 = csqrtl(b);}}}

// tests
#include <stdio.h>
#include <math.h>
void pc(c a){printf("%20.17lf + %20.17lfi\n",
    creal(a),  cimag(a));}
int main(){c A[N];
typedef double R;
long double pi = 4.*atanl(1);
R ep = 6.7e-16;;
for(int s=0; s<N; ++s) {for(int j=0; j<N; ++j) A[j]=s==j;
fft(N, A);
if(0) {int j; for(j = 0; j<N; ++j) pc(A[j]);}
for(int j=0; j<N; ++j){R a = ((s*j)&(N-1))*(2*pi/N);
  if(fabs(creal(A[j]) - cos(a)) > ep
  || fabs(cimag(A[j]) - sin(a)) > ep) {
    printf("er: %20.17f + %20.17fi %4x %4x\n    ", cos(a), sin(a), s, j);
     pc(A[j]);}}}}

/* ep = 6.7e-16
er:  0.99883279685352799 + -0.04830159344948058i  37f 1cbf
     0.99883279685352788 + -0.04830159344947991i
er:  0.99883279685352799 + -0.04830159344948058i  b7f 14bf
     0.99883279685352788 + -0.04830159344947991i
er:  0.99883279685352799 + -0.04830159344948058i 137f  cbf
     0.99883279685352788 + -0.04830159344947991i
er:  0.99883279685352799 + -0.04830159344948058i 1b7f  4bf
     0.99883279685352788 + -0.04830159344947991i
*/