#include <stdio.h>
#include <string.h>
#include <math.h>
#include "zhufft.h"


void add( COMPLEX *a, COMPLEX *b, COMPLEX *ret)
{
	ret->re = a->re + b->re;
	ret->im = a->im + b->im;
}

void sub(COMPLEX *a,COMPLEX *b,COMPLEX *ret)
{
	ret->re = a->re - b->re;
	ret->im = a->im - b->im;
}

void mul(COMPLEX *m,COMPLEX *n,COMPLEX *ret)
{
	double a,b,c;

	a = ( m->re - m->im ) * n->im;
	b = m->re * ( n->re - n->im );
	c = m->im * ( n->re + n->im );
	ret->re = a + b;
	ret->im = a + c;
}

int reverse(int t,int k)
{
	int i,j,m,s;

	m = t;
	s = 0;
	for ( i = 0 ; i < k ; i++ )
	{
		j = m / 2;
		s = s * 2 + ( m - j * 2 );
		m = j;
	}

	return s;
}

void root( COMPLEX omega[],int n)
{
	int j;

	omega->re = 1;
	omega->im = 0;
	omega[1].re = (double)cos(2 * PI / n);
	omega[1].im = -(double)sin(2 * PI / n);

	for ( j = 2 ; j < n/2 ; j++ )
		mul(&omega[1],&omega[j-1],&omega[j]);

	for ( j = n/2 ; j < n ; j++ )
	{
		omega[j].re = -omega[j-n/2].re;
		omega[j].im = -omega[j-n/2].im;
	}
}

void fft(COMPLEX in[],COMPLEX out[], COMPLEX omega[],int n)
{
	int s,k,m,l,nv,t,j;
	COMPLEX podd,ret;

	k = (int)(log(n) / log(2) + 0.5);
	nv = n;
	m = 1;
	for ( l = k-1 ; l >= 0 ; l-- )
	{
		for ( t = 0 ; t < m * nv ; t+=nv )
			for ( j = 0 ; j < nv/2 ; j++ )
			{
				s = (t+j) / (int)pow(2,l);
				s = reverse(s,k);

				ret = omega[s];
				mul(&ret,&in[t+j+nv/2],&podd);
				sub(&in[t+j],&podd,&in[t+j+nv/2]);
				add(&in[t+j],&podd,&in[t+j]);
			}
		m *= 2;
		nv /= 2;
	}

	for ( t = 0 ; t < n ; t++ )
	{
		s = reverse(t,k);
		out[t] = in[s];
	}
}