#include <mkl.h>
#include <math.h>
#include <time.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#ifdef windows
#include <windows.h>
#endif
#include <pmmintrin.h>
#include <immintrin.h>
#include "task.h"
#include "timer.h"
#include "mkl_dfti.h"
#include <mkl_service.h>
#include <mkl_dfti.h>

int task_FFT()
{
#define FORWARD 0

	int idxLoop;
	int idxData;
	int idxTimeLoop;
	int idxFFTSizeLoop;
	int FFT_sizeBuf[FFT_SIZE_NUM]={128, 256, 512, 1024, 2048};

	//FFT structure variables
	int FFTSize;
    int i, j;
	DFTI_DESCRIPTOR_HANDLE FFT_desc;
	MKL_LONG status;

	//pointer variables
	short *FFT_in_fixed;
	float *FFT_in_fixed_float;
	float *FFT_out_fixed_float;
	float *cs_float;
	float *FFT_in_singlePrecision;
	float *FFT_out_singlePrecision;

    FFTSize = FFT_sizeBuf[4];

    //allocate buffers
	cs_float = (float *)mkl_malloc(FFTSize * 2 * sizeof(float), 64);
    FFT_in_fixed = (short *)mkl_malloc(FFTSize * 2 * sizeof(short), 4);
	FFT_in_fixed_float = (float *)mkl_malloc(FFTSize * 2 * sizeof(float), 64);
	FFT_out_fixed_float = (float *)mkl_malloc(FFTSize * 2 * sizeof(float), 64);
    FFT_in_singlePrecision = (float *)mkl_malloc(FFTSize * 2 * sizeof(float), 64);
    FFT_out_singlePrecision = (float *)mkl_malloc(FFTSize * 2 * sizeof(float), 64);

    //generate the random input, uncomment below will make the random input change from running to running
    //secondAsSeed = time(0);
    //srand(secondAsSeed);

    for (idxData = 0; idxData < FFTSize; idxData++)
    {
        float randVal1 = sin(((double)rand()) / RAND_MAX * 2 * PI);
		float randVal2 = sin(((double)rand()) / RAND_MAX * 2 * PI);

		FFT_in_fixed_float[idxData * 2] = randVal1;
		FFT_in_fixed_float[idxData * 2 + 1] = randVal2;

		FFT_in_singlePrecision[idxData * 2] = randVal1;
		FFT_in_singlePrecision[idxData * 2 + 1] = randVal2;
    }

    for (i = 0; i < FFT_SIZE_NUM; i++)
    {
        FFTSize = FFT_sizeBuf[i];
        printf("Doing %d-point FFT!\n", FFTSize);
        fprintf(fp_timingLog, "Doing %d-point FFT:\n", FFTSize);
        fprintf(fp_PrecisionLog, "Doing %d-point FFT:\n", FFTSize);
		//mkl_set_num_threads(2);
        //configure and run single precision floating-point FFT
        //configure fft descriptor: single precision, 2048-point FFT
        //DFTI_SINGLE is single precision, DFTI_DOUBLE is double precision
        status = DftiCreateDescriptor(&FFT_desc, DFTI_SINGLE, DFTI_COMPLEX, 1, FFTSize);
        //DFTI_INPLACE is FFT output overwrites input, DFTI_NOT_INPLACE is FFT output does not overwrite input
        status = DftiSetValue(FFT_desc, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
		status = DftiSetValue(FFT_desc, DFTI_NUMBER_OF_USER_THREADS, 2);
#if FORWARD
        //status = DftiSetValue(FFT_desc, DFTI_ORDERING, DFTI_BACKWARD_SCRAMBLED);
#else
        status = DftiSetValue(FFT_desc, DFTI_ORDERING, DFTI_BACKWARD_SCRAMBLED);
#endif
        //frease FFT descriptor
        status = DftiCommitDescriptor(FFT_desc);


#if FORWARD
		//run fft with forward method
		status = DftiComputeForward(FFT_desc, FFT_in_singlePrecision, FFT_out_singlePrecision);
#else
		//run fft with backward method
		status = DftiComputeBackward(FFT_desc, FFT_in_singlePrecision, FFT_out_singlePrecision);
#endif


        status = DftiFreeDescriptor(&FFT_desc);

		for (idxData = 0; idxData < FFTSize; idxData++)
		{
			FFT_out_singlePrecision[idxData * 2] = FFT_out_singlePrecision[idxData * 2] / sqrt((double)FFTSize);
			FFT_out_singlePrecision[idxData * 2 + 1] = FFT_out_singlePrecision[idxData * 2 + 1] / sqrt((double)FFTSize);
		}
        {
#if FORWARD
            FILE *fp = (FILE *)fopen("FFT_MKL_for_data_Log.txt", "w");
#else
            FILE *fp = (FILE *)fopen("FFT_MKL_back_data_Log.txt", "w");

#endif
            for (idxData = 0; idxData < FFTSize; idxData++)
            {
                fprintf(fp, " %.7f, %.7f \n", FFT_out_singlePrecision[idxData * 2], FFT_out_singlePrecision[idxData * 2 + 1]);
            }
            fclose(fp);
        }

	//free allocated buffers
	mkl_free(cs_float);
	mkl_free(FFT_in_fixed);
	mkl_free(FFT_in_fixed_float);
	mkl_free(FFT_out_fixed_float);
	mkl_free(FFT_in_singlePrecision);
	mkl_free(FFT_out_singlePrecision);

	return 0;
}