#include <stdio.h>
#include <math.h>
typedef double __complex__ c;
void pc(c a){printf("%18.11e + %18.11ei ",
    __real__ a, __imag__ a);}
c Sqrt(c x){
double A = __real__ x, B = __imag__ x;
if(A < 0) return 1.i* Sqrt(-x);
{double r = sqrt((A + sqrt(A*A + B*B))/2);
return r + 1i*(B/(2*r));}}
static double const er = 3.0e-14;

void fft(c a[], const int 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 = Sqrt(b);}}}

int main(){c w[1024];
// Initialize array w, thereafter a constant.
{c W[10];
  {int j=10; W[9] = -1; while(--j) W[j-1] = Sqrt(W[j]);}
  {int j=1024; while(j--) {int k=10; 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]; pc(r); printf("is r\n");
  {int j = 1024; while(j--) {
     c e = w[j]*r - w[(j+1)%1024];
     if(fabs(__real__ e) > er || fabs(__imag__ e) > er) {
         pc(w[j]); pc(w[(j+1)%1024]);
         printf("Roots: %d\n", j);}}}
 {int j = 10; while(j--) r = r*r; pc(r); printf("is r^1024\n");}
}
// Evaluate transform for each basis vector and
// demonstrate that the answer is correct.
{int k; for(k=0; k<1024; ++k) {
   c A[1024]; int j; for(j=0; j<1024; ++j) A[j]=0;
   A[k]=1;
   fft(A, 1024); {
     int n; for(n=0; n<1024; ++n) {
        c p = w[n*k % 1024];
        if(fabs(__real__ A[n] - __real__ p) > er
         || fabs(__imag__ A[n] - __imag__ p) > er)
          {pc(A[n]); pc(p); printf("%6d %6d\n", k, n);}
}}}}
return 0;}