#include <math.h>
#include <stdio.h>
#include <complex.h>
typedef float _Complex c;
typedef double _Complex cx;
#define C const
static void pc(c C a){printf("%9.2f + %9.2fi\n", __real__ a, __imag__ a);}

static void fft(c a[], C 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; 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 = csqrt(b);}}}

#if 0
//Inverse transform
int main(){int C z = 5, sz = 1<<z;
c a[sz]; {int j=sz; while(j--) a[j]=0;}
a[3] = 1.+2i;
fft(a, sz);
{int k=sz; while(k--) a[k] = conj(a[k])/sz;}
fft(a, sz);
{int i; for(i=0; i<sz; ++i) {printf("%2d ", i); pc(a[i]);}}

return 0;}
#endif

// Convolution demo
int main(){int C z = 5; int C sz = 1<<z;
c a[sz], b[sz]; {int j=sz; while(j--) {a[j]=0; b[j]=0;}}
a[1] = 1; a[2] = 3; b[1] = 1; b[3] = 2;
fft(a, sz); fft(b, sz);
{int j=sz; while(j--) a[j] = conj(a[j]*b[j])/sz;}
fft(a, sz);
{int i; for(i=0; i<sz; ++i) {printf("%2d ", i); pc(a[i]);}}
return 0;}