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

void fft(const int n, c a[n]){
  {int j=0; for(int 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; cx b=-1; while(j<n) {cx w = 1;
    for(int m=0; m<j; ++m) {
      for(int 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>
#define p 17
void pc(c a){printf("%19.16f + %19.16fi\n", creal(a),  cimag(a));}
void mx(int m){c A[1<<p];
typedef double R;
R pi = 4.*atan(1), ep = 2.55351296e-15;
for(int j=0; j<(1<<p); ++j) A[j]=j==m;
fft((1<<p), A);
if(0) {for(int j = 0; j<(1<<p); ++j) pc(A[j]);}
printf("m = %d\n", m);
for(int j=0; j<(1<<p); ++j){R a = 2*(m*j&((1<<p)-1))*pi/(1<<p);
  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]);}}}
int main(){mx(0); mx(1); mx(11); mx(244); mx(833); mx(1023); 
  mx(1024); mx(1025); mx(14572); return 0;}
/* ep = 2.55351296e-15
m = 0
m = 1
m = 11
m = 244
m = 833
m = 1023
er: -0.1843592794914510 + -0.9828589197160460i 65308
    -0.1843592794914510 + -0.9828589197160487i
m = 1024
m = 1025
m = 14572
*/