#include <math.h>
#include <dspc.h>

/* numeric C test code for 21K compiler */

void init_w(int n, complex float *w);
void fft_nc(int n, complex float *x, complex float *w);

/* COMPLEX STRUCTURE */

typedef struct {
    float real, imag;
} COMPLEX;

void fft_c(int n, COMPLEX *x, COMPLEX *w);

#define PI 3.14159265358979323846
int n = 16;
complex float x[1024],w[1024],y[1024];
void main()
{
  int i;
  COMPLEX *xc,*wc,*yc;
  xc = (COMPLEX *) x;
  wc = (COMPLEX *) w;
  yc = (COMPLEX *) y;
  x[0] = 1.0;
  for( ; n < 2048 ; n *= 2) {
    init_w(n,w);
    fft_nc(n,x,w);
    fft_c(n,xc,wc);
    {
        iter I = n;
        y[I]=x[I]*w[I];
    }

    for(i = 0 ; i < n ; i++) {
        yc[i].real = xc[i].real*wc[i].real - xc[i].imag*wc[i].imag;
        yc[i].imag = xc[i].real*wc[i].imag + xc[i].imag*wc[i].real;
    }
    {
        iter I = n;
        y[I]=x[I]+w[I];
    }

    for(i = 0 ; i < n ; i++) {
        yc[i].real = xc[i].real+wc[i].real;
        yc[i].imag = xc[i].imag+wc[i].imag;
    }
    {
        iter I = n;
        y[0]=sum(x[I]*w[I]);
    }

    for(i = 0 ; i < n ; i++) {
        yc[0].real += xc[i].real*wc[i].real - xc[i].imag*wc[i].imag;
        yc[0].imag += xc[i].real*wc[i].imag + xc[i].imag*wc[i].real;
    }
  }
}
void init_w(int n, complex float *w)
{
    iter I = n;
    float a = 2.0*PI/n;
    w[I] = cosf(I*a) + 1i*sinf(I*a);
}
void fft_nc(int n, complex float *x, complex float *w)
{
    int size,sect,deg = 1;
    for(size=n/2 ; size > 0 ; size/=2) {
        for(sect=0 ; sect < n ; sect += 2*size) {
            complex float *x1=x+sect;
            complex float *x2=x1+size;
            { iter I=size;
              for(I) {
                complex float temp;
                temp = x1[I] + x2[I];
                x2[I] = (x1[I] - x2[I]) * w[deg*I];
                x1[I] = temp;
              }
            }
        }
        deg *= 2;
    }
}

void fft_c(int n,COMPLEX *x,COMPLEX *w)
{
    COMPLEX u,temp,tm;
    COMPLEX *xi,*xip,*wptr;

    int i,j,le,windex;

/* start fft */

    windex = 1;
    for(le=n/2 ; le > 0 ; le/=2) {
        wptr = w;
        for (j = 0 ; j < le ; j++) {
            u = *wptr;
            for (i = j ; i < n ; i = i + 2*le) {
                xi = x + i;
                xip = xi + le;
                temp.real = xi->real + xip->real;
                temp.imag = xi->imag + xip->imag;
                tm.real = xi->real - xip->real;
                tm.imag = xi->imag - xip->imag;             
                xip->real = tm.real*u.real - tm.imag*u.imag;
                xip->imag = tm.real*u.imag + tm.imag*u.real;
                *xi = temp;
            }
            wptr = wptr + windex;
        }
        windex = 2*windex;
    }            
}