#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>
#define LN 10
#define N (1<<LN)
#define m 11
static double const ep = 4.45e-16;
#include <math.h>
void pc(c a){printf("%14.7e + %14.7ei\n",
   creal(a), cimag(a));}

int main(){c w[N];
// Initialize array w, thereafter a constant.
{c W[LN];
  {W[LN-1] = -1; for(int j=LN-1; j; --j) W[j-1] = csqrtl(W[j]);}
  {int j=N; while(j--) {int k=LN; c a = 1;
     while(k--) if(j & (1 << k)) a *= W[k]; w[j] = a;}}}
// Demonstrate that w[1] is the correct root of unity
// and then that the other values of w are correct.
{c r = w[1]; printf("r = "); pc(r); 
  {int j = N; while(j--) {
     c e = w[j]*r - w[(j+1)&(N-1)];
     if(fabs(__real__ e) > ep || fabs(__imag__ e) > ep) {
         pc(w[j]); pc(w[(j+1)&(N-1)]);
         printf("Roots: %d\n", j);}}}
 {int j = LN; while(j--) r *= r; printf("r^%d-1 = ", N); pc(r-1);}
}
// Evaluate transform for each basis vector and
// demonstrate that the answer is correct.
{int k=N; while(k--) {c A[N]; int j=N; while(j--) A[j]=0;
   A[k]=1;
   fft(N, A); {int n=N; while(n--) {
        c p = w[n*k &(N-1)];
        if(fabs(creal(A[n]) - creal(p)) > ep
         || fabs(cimag(A[n]) - cimag(p)) > ep)
          {pc(A[n]); pc(p); printf("er: %6d %6d\n", k, n);}}}}}
return 0;}
/*
r =  9.9998118e-01 +  6.1358846e-03i
r^1024-1 = -2.5313085e-14 + -6.6613381e-16i
*/