#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>

#define DATA_L 1024
#define FP_acc 1000

#define PI 3.141592653589793238462643383279502884197169399375105820974944



int cos_delta_FP_arr[DATA_L] = {0,};
int sin_delta_FP_arr[DATA_L] = {0,};

void FFT_fp(uint32_t* inp, uint32_t inp_L, uint32_t* buf){

	    // buf имеет длину inp_L * 2 (Re, Im, Re, Im, ...)
	    // inp содержит inp_L значений uint32_t

	    uint32_t i, j, bit;
//	    uint32_t N = inp_L / 2;  // длина комплексного массива (inp содержит только Re-часть)
	    uint32_t N = inp_L;  // длина комплексного массива (inp содержит только Re-часть)

	    // --- копирование входных данных в буфер ---
	    for (i = 0; i < inp_L; i++) {
	        buf[i * 2]     = (int)inp[i];  // Re
	        buf[i * 2 + 1] = 0;            // Im = 0
	    }

	    // --- bit-reversal перестановка ---
	    j = 0;
	    for (i = 1; i < N; i++) {
	        bit = N >> 1;
	        while (j & bit) {
	            j ^= bit;
	            bit >>= 1;
	        }
	        j |= bit;
	        if (i < j) {
	            int tmp_re = buf[i * 2];
	            int tmp_im = buf[i * 2 + 1];
	            buf[i * 2]     = buf[j * 2];
	            buf[i * 2 + 1] = buf[j * 2 + 1];
	            buf[j * 2]     = tmp_re;
	            buf[j * 2 + 1] = tmp_im;
	        }
	    }

	    // --- уровни бабочек ---
	    uint32_t m = 2;
	    while (m <= N) {
	        uint32_t half = m >> 1;

	        int c_delta_FP = cos_delta_FP_arr[m - 2];
	        int s_delta_FP = sin_delta_FP_arr[m - 2];

	        for (uint32_t start = 0; start < N; start += m) {
	            int wr = FP_acc;
	            int wi = 0;

	            for (uint32_t k = 0; k < half; k++) {
	                int u_re = buf[(start + k) * 2];
	                int u_im = buf[(start + k) * 2 + 1];
	                int v_re = buf[(start + k + half) * 2];
	                int v_im = buf[(start + k + half) * 2 + 1];

	                // t = w * v  (в фиксированной точке)
	                int64_t t_re = ((int64_t)wr * v_re - (int64_t)wi * v_im) / FP_acc;
	                int64_t t_im = ((int64_t)wr * v_im + (int64_t)wi * v_re) / FP_acc;

	                // верх/низ
	                buf[(start + k) * 2]              = u_re + (int)t_re;
	                buf[(start + k) * 2 + 1]          = u_im + (int)t_im;
	                buf[(start + k + half) * 2]       = u_re - (int)t_re;
	                buf[(start + k + half) * 2 + 1]   = u_im - (int)t_im;

	                // w *= w_m (обновление twiddle)
	                int64_t wr_new = ((int64_t)wr * c_delta_FP - (int64_t)wi * s_delta_FP) / FP_acc;
	                int64_t wi_new = ((int64_t)wr * s_delta_FP + (int64_t)wi * c_delta_FP) / FP_acc;
	                wr = (int)wr_new;
	                wi = (int)wi_new;
	            }
	        }
	        m <<= 1;
	    }
	    for (uint32_t i = 0; i < inp_L; ++i){
	    	//buf[i*2] /= inp_L;
	    	//buf[i*2+1] /= inp_L;
	    	printf("re,im: %d, %d\n", buf[i*2], buf[i*2 +1]);
	    }

}

void imagery2abs(uint32_t* inp, uint32_t inp_L, uint32_t *out){
	printf("+++++++++++++++++++++\n");
	printf("calculating abs^2\n");
	for (uint32_t i = 0; i < inp_L; ++i){
		int val = (inp[i*2]*inp[i*2] + inp[i*2 +1]*inp[i*2 +1])/FP_acc;
		printf("%d\n", val);
		out[i] = val;
	}
}



void sin_cos_filler(int* sin_arr, int* cos_arr, uint32_t L){
	//[int(-sin(2.0 * pi / m) * FP_acc) for m in range(2,inp_L//2 + 2)]
	for (uint32_t i = 0; i < L; ++i){
		int m = i + 2;
		double angle = 2.0* PI/m;
		sin_arr[i] = lround(-sin(angle) * FP_acc);
		cos_arr[i] = lround(cos(angle) * FP_acc);
		printf("i, angle, sin, cos: %d %g %d %d\n", i, angle, sin_arr[i], cos_arr[i]);
	}
}

void data_generator(uint32_t* X, uint32_t* Y, uint32_t N){
	for (int i = 0; i < N; ++i){
		X[i] = (int) i;
		Y[i] = lround(sin(((double)i)*2.0* PI/(N/10))*FP_acc);
//		Y[I] += 10*FP_acc;
	}
}

void main(){
	char* logfilename = "tmp";
	uint32_t data_X[DATA_L] = {0,};
	uint32_t data_Y[DATA_L] = {0,};
	uint32_t data_res_FFT_imag[DATA_L*2] = {0,};
	uint32_t data_res_FFT_abs[DATA_L] = {0,};

	sin_cos_filler(sin_delta_FP_arr, cos_delta_FP_arr, DATA_L);

	data_generator(data_X, data_Y, DATA_L);


	FFT_fp(data_Y, DATA_L, data_res_FFT_imag);


	imagery2abs(data_res_FFT_imag, DATA_L, data_res_FFT_abs);

	FILE *logfile_ptr;
    logfile_ptr = fopen(logfilename, "w");
    fprintf(logfile_ptr, "X, Y, FFT\n");
	for (uint32_t i = 0; i < DATA_L; ++i){
		fprintf(logfile_ptr, "%d, %d, %d \n", data_X[i], data_Y[i],data_res_FFT_abs[i]);
	}
	fclose(logfile_ptr);
}