#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "fourier.h"

// RNG implemented localy to avoid library incongruences
#ifdef RAND_MAX
#undef RAND_MAX
#endif
#define RAND_MAX 32767
static unsigned long long int next = 1;

int rand( void ) {
    next = next * 1103515245 + 12345;
    return (unsigned int)(next / 65536) % RAND_MAX+1;
}

void srand( unsigned int seed ) {
    next = seed;
}
// End of RNG implementation

/* Return 0 when V1 and V2 do not match within the allowed FP_ABSTOLERANCE.  */
static inline int
check_FP(float V1, float V2) {
  double AbsTolerance = FP_ABSTOLERANCE;
  double Diff = fabs(V1 - V2);
  if (Diff > AbsTolerance) {
    fprintf(stderr, "A = %lf and B = %lf differ more than"
                  " FP_ABSTOLERANCE = %lf\n", V1, V2, AbsTolerance);
    return 0;
  }

  return 1;
}

int main(int argc, char *argv[]) {
	unsigned MAXSIZE;
	unsigned MAXWAVES;
	unsigned i,j;
	float *RealIn;
	float *ImagIn;
	float *RealOut;
	float *ImagOut;
	float *RealOut_StrictFP;
	float *ImagOut_StrictFP;
	float *coeff;
	float *amp;
	int invfft=0;

	if (argc<3)
	{
		printf("Usage: fft <waves> <length> -i\n");
		printf("-i performs an inverse fft\n");
		printf("make <waves> random sinusoids");
		printf("<length> is the number of samples\n");
		exit(-1);
	}
	else if (argc==4)
		invfft = !strncmp(argv[3],"-i",2);
	MAXSIZE=atoi(argv[2]);
	MAXWAVES=atoi(argv[1]);
  srand(1);
		
 RealIn=(float*)malloc(sizeof(float)*MAXSIZE);
 ImagIn=(float*)malloc(sizeof(float)*MAXSIZE);
 RealOut=(float*)malloc(sizeof(float)*MAXSIZE);
 ImagOut=(float*)malloc(sizeof(float)*MAXSIZE);
 RealOut_StrictFP=(float*)malloc(sizeof(float)*MAXSIZE);
 ImagOut_StrictFP=(float*)malloc(sizeof(float)*MAXSIZE);
 coeff=(float*)malloc(sizeof(float)*MAXWAVES);
 amp=(float*)malloc(sizeof(float)*MAXWAVES);

 /* Makes MAXWAVES waves of random amplitude and period */
	for(i=0;i<MAXWAVES;i++) 
	{
		coeff[i] = i%1000;
		amp[i] = i%1000;
	}
 for(i=0;i<MAXSIZE;i++) 
 {
   /*   RealIn[i]=rand();*/
	 RealIn[i]=0;
	 for(j=0;j<MAXWAVES;j++) 
	 {
		 /* randomly select sin or cos */
		 if (rand()%2)
		 {
		 		RealIn[i]+=coeff[j]*cos(amp[j]*i);
			}
		 else
		 {
		 	RealIn[i]+=coeff[j]*sin(amp[j]*i);
		 }
  	 ImagIn[i]=0;
	 }
 }

 /* regular*/
 fft_float (MAXSIZE,invfft,RealIn,ImagIn,RealOut,ImagOut);
 fft_float_StrictFP (MAXSIZE,invfft,RealIn,ImagIn,RealOut_StrictFP,ImagOut_StrictFP);
 
 printf("RealOut:\n");
 for (i=0;i<MAXSIZE;i++) {
   if (!check_FP(RealOut[i], RealOut_StrictFP[i]))
     return 1;
   printf("%f \t", RealOut_StrictFP[i]);
 }
 printf("\n");

printf("ImagOut:\n");
 for (i=0;i<MAXSIZE;i++) {
   if (!check_FP(ImagOut[i], ImagOut_StrictFP[i]))
     return 1;
   printf("%f \t", ImagOut_StrictFP[i]);
 }
 printf("\n");

 free(RealIn);
 free(ImagIn);
 free(RealOut);
 free(ImagOut);
 free(RealOut_StrictFP);
 free(ImagOut_StrictFP);
 free(coeff);
 free(amp);
 exit(0);


}