#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <ctime>
#include <iostream>
#include <iomanip>
using namespace std;
const double PI = 3.14159265359;


void DirectFourierTransformation(double* const massiveF, double** const matrixA, const int SIZE)
{
    for(int k = 0; k < SIZE; ++k)
    {
        for(int j = 0; j < SIZE; ++j)
        {
            matrixA[0][k] = matrixA[0][k] + cos(-2.0 * PI * k * j / SIZE) * massiveF[j];
            matrixA[1][k] = matrixA[1][k] + sin(-2.0 * PI * k * j / SIZE) * massiveF[j];
        }
        matrixA[0][k] = matrixA[0][k] / SIZE;
        matrixA[1][k] = matrixA[1][k] / SIZE;
    }

}

void InverseFourierTransform(double ** const matrixA, double ** const matrixF, const int SIZE)
{
    for(int k = 0; k < SIZE; ++k)
    {
        for(int j = 0; j < SIZE; ++j)
        {
            matrixF[0][k] = matrixF[0][k] + cos(2.0 * PI * k * j / SIZE) * matrixA[0][j] - sin(2.0 * PI * k * j / SIZE) * matrixA[1][j];
            matrixF[1][k] = matrixF[1][k] +  cos(2.0 * PI * k * j / SIZE) * matrixA[1][j] + sin(2.0 * PI * k * j / SIZE) * matrixA[0][j];
        }
    }
}

void Semi_FastFourierTransformation(double* const massivef, double** const matrixA, const int SIZE)
{
    int p1, p2;
    p1 = (int) sqrt(SIZE);
    while (SIZE % p1)
        --p1;
    p2 = SIZE / p1;

    for(int k = 0; k < SIZE; ++k)
    {
        int k1 = k % p1, k2 = k / p1;
        for(int j2 = 0; j2 < p2; ++j2)
        {
            double A[2]{};
            for(int j1 = 0; j1 < p1; ++j1)
            {
                A[0] += cos(-2.0 * PI * j1 * k1 / p1) * massivef[j2 + p2 * j1];
                A[1] += sin(-2.0 * PI * j1 * k1 / p1) * massivef[j2 + p2 * j1];
            }
            matrixA[0][k] += (A[0] * cos(-2.0 * PI * j2 / (p1 * p2) * (k1 + p1 * k2)) - A[1] * sin(-2.0 * PI * j2 / (p1 * p2) * (k1 + p1 * k2))) / p1;
            matrixA[1][k] += (A[0] * sin(-2.0 * PI * j2 / (p1 * p2) * (k1 + p1 * k2)) + A[1] * cos(-2.0 * PI * j2 / (p1 * p2) * (k1 + p1 * k2))) / p1;
        }
        matrixA[0][k] /= p2;
        matrixA[1][k] /= p2;
    }
}

void FastFourierTransform(double* massivef, double** const matrixA, const int SIZE)
{
    int r = (int)log2(SIZE);
    int* k = new int[r];
    for(int km = 0; km < SIZE; ++km)
    {

    }
}


int main()
{
    srand(time(NULL));
    int SIZE;
    cout << "INPUT SIZE: ";
    cin >> SIZE;
    double* massivef = new double[SIZE];
    double** matrixA = new double*[2];
    double** matrixF = new double*[2];
    for(int i = 0; i < SIZE; i++)
    {
        massivef[i] = rand()%100;
        if(i < 2)
        {
            matrixA[i] = new double[SIZE];
            matrixF[i] = new double[SIZE];
        }
    }
    for(int i = 0; i < 2; i++)
    {
        for(int j = 0; j < SIZE; j++)
        {
            matrixF[i][j] = 0;
            matrixA[i][j] = 0;
        }
    }

    cout << endl << " Massive: ";
    for (int i = 0; i < SIZE; ++i)
        cout << massivef[i] << setw(5);
    cout << endl << endl;

    DirectFourierTransformation(massivef, matrixA, SIZE);

    cout << "DirectFourierTransformation: " << endl;
    for (int i = 0; i < SIZE; ++i)
        cout << "( " << matrixA[0][i] << setw(5) << " , i" << matrixA[1][i] << " )" << endl;

    InverseFourierTransform(matrixA, matrixF, SIZE);

    cout << endl << "InverseFourierTransform: " << endl;
    for (int i = 0; i < SIZE; ++i)
        cout << "( " << matrixF[0][i]  << setw(5) << " , i" << matrixF[1][i] << " )" << endl;
    cout << endl;

    for(int i = 0; i < 2; i++)
        for(int j = 0; j < SIZE; j++)
            matrixA[i][j] = 0;

    Semi_FastFourierTransformation(massivef, matrixA, SIZE);

    cout << endl << "Semi_FastFourierTransformation: " << endl;
    for (int i = 0; i < SIZE; ++i)
        cout << "( " << matrixA[0][i] << setw(5) << " , i" << matrixA[1][i] << " )" << endl;
    cout << endl;

    for(int i = 0; i < 2; i++)
        for(int j = 0; j < SIZE; j++)
            matrixA[i][j] = 0;

    for(int i = 0; i < 2; i++)
    {
        delete[] matrixA[i];
        delete[] matrixF[i];
    }
    delete [] matrixA;
    delete [] matrixF;
    delete [] massivef;
    system("pause");
}