#include <stdio.h>
#include <stdlib.h>
#include <strings.h>
#include <math.h>
#include <complex.h>
#include <cilk/cilk.h>
#include "matrixlib.h"
#include "tictoc.h"

void dft(int N,complex x[N],complex b[N]){
	bzero(b,sizeof(complex)*N);
	cilk_for(int k=0;k<N;k++){
		for(int l=0;l<N;l++){
			b[k]+=cexp(-I*2*M_PI*k*l/N)*x[l];
		}
	}
}
void sfft(int N,int s,complex x[N],complex b[N]){
	if(N==1) {
		b[0]=x[0];
		return;
	}
	if(N%2){
		printf("N not divisible by 2!\n");
		exit(1);
	}
	int K=N/2;
	sfft(K,2*s,&x[0],&b[0]);
	sfft(K,2*s,&x[s],&b[K]);
	for(int k=0;k<K;k++){
		complex even=b[k],odd=b[k+K];
		complex w=cexp(-I*2*M_PI*k/N);
		b[k]=even+w*odd;
		b[k+K]=even-w*odd;
	}
}
void fft(int N,int s,complex x[N],complex b[N]){
	if(s>8){
		sfft(N,s,x,b);
		return;
	}
	if(N==1) {
		b[0]=x[0];
		return;
	}
	if(N%2){
		printf("N not divisible by 2!\n");
		exit(1);
	}
	int K=N/2;
	cilk_spawn fft(K,2*s,&x[0],&b[0]);
	fft(K,2*s,&x[s],&b[K]);
	cilk_sync;
	cilk_for(int k=0;k<K;k++){
		complex even=b[k],odd=b[k+K];
		complex w=cexp(-I*2*M_PI*k/N);
		b[k]=even+w*odd;
		b[k+K]=even-w*odd;
	}
}

#define FFTSIZE 16384
int main(){
	complex X[FFTSIZE],B[FFTSIZE];
	for(int i=0;i<FFTSIZE;i++) X[i]=i;
	tic();
	fft(FFTSIZE,1,X,B);
	double t=toc();
	printf("FFT time %g seconds.\n",t);
	tic();
	dft(FFTSIZE,X,B);
	t=toc();
	printf("Slow FT time %g seconds.\n",t);
//	cvecprint(FFTSIZE,B);
	return 0;
}