#include <math.h>
#define PI 3.14159265358979323846


double mult_real(double re_a, double im_a, double re_b, double im_b){
	return (re_a*re_b)-(im_a*im_b);
}

double mult_im(double re_a, double im_a, double re_b, double im_b){
	return (re_a*im_b)+(re_b*im_a);
}


void FFT(double *q, double *w, int n, int m){
	//USED FOR ORDERING
	int i, j, c, r;
	int b=1;
	int a=n/2;
		
	double new[256];
	double new_im[256];
	
	// ORdering algorithm
	for(i=0; i<(m-1); i++){
		int d=0;
		for (j=0; j<b; j++){
			for (c=0; c<a; c++){	
				int e=c+d;
				new[e]=q[(c*2)+d];
				new_im[e]=w[(c*2)+d];
				new[e+a]=q[2*c+1+d];
				new_im[e+a]=w[2*c+1+d];
			}
			d+=(n/b);
		}
		for (r=0; r<n;r++){
			q[r]=new[r];
			w[r]=new_im[r];
		}
		b*=2;
		a=n/(2*b);
	}
	//end ordering algorithm

	b=1;
	int k=0;
	for (j=0; j<m; j++){	
	//MATH
		for(i=0; i<n; i+=2){
			if (i%(n/b)==0 && i!=0)
				k++;
			//printf("%f\n", -PI*k/b);
			//printf ("cos: %f , sin: %f\n", cos(-PI*k/b), sin(-PI*k/b));
			double real=mult_real(q[i+1], w[i+1], cos(-PI*k/b), sin(-PI*k/b));	
			double imagine=mult_im(q[i+1], w[i+1], cos(-PI*k/b), sin(-PI*k/b));
			//printf ("q: %f , w: %f\n\n", real, imagine);	
			new[i]=q[i]+real;
			new_im[i]=w[i]+imagine;
			new[i+1]=q[i]-real;
			new_im[i+1]=w[i]-imagine;

		}
		for (i=0; i<n; i++){
			q[i]=new[i];
			w[i]=new_im[i];
		}
	//END MATH
	
	//REORDER
		for (i=0; i<n/2; i++){
			new[i]=q[2*i];
			new[i+(n/2)]=q[2*i+1];
			new_im[i]=w[2*i];
			new_im[i+(n/2)]=w[2*i+1];
		}
		for (i=0; i<n; i++){
			q[i]=new[i];
			w[i]=new_im[i];
		}
	//END REORDER	
		b*=2;
		k=0;
	}
	for (i=0; i<n; i++)
		printf("%f\t%f\n", q[i], w[i]);


}


int main(){
	//double a[2]={0, 1};
	//double a[4]={1, 2, 3, 4};
	//double a[8]={1, 2, 3, 4, 5, 6, 7, 8};
	double a[16]={0, 1, 2, 3 ,4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7};
	double j[16]={0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; 
	//double a[8]={4, 3, 5, 7, 2, 4, 1, 6};
	FFT(a, j, 16, 4);	
}