#include <c_FFT.h>
#include <math.h>
#include <string.h>

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

/* Private Helper: Verifies if an integer is a strict power of 2 */
C_STATIC_FORCE_INLINE
bool c_FFT_IsPowerOfTwo(c_size_t n) {
    return (n > 0) && ((n & (n - 1)) == 0);
}

/* Private Helper: In-place O(N) bit-reversal permutation shuffling */
static void c_FFT_BitReversePermutation(c_Complex_t* signal, c_size_t n) {
    c_size_t j = 0;
    for (c_size_t i = 0; i < n; ++i) {
        if (i < j) {
            c_Complex_t temp = signal[i];
            signal[i] = signal[j];
            signal[j] = temp;
        }
        c_size_t bit = n >> 1;
        while (j & bit) {
            j ^= bit;
            bit >>= 1;
        }
        j ^= bit;
    }
}

/* ------------------------------------------------------------------------------------------------------------------ */
/*  */


c_err_t c_FFT_Init(c_FFT_t* self, c_size_t N, c_Allocator_t* allocator) {
    if (!self || N == 0) return C_ERR_PARAM;
    if (!c_FFT_IsPowerOfTwo(N)) return C_ERR_PARAM; /* Algorithm constraint check */

    self->allocator = allocator ? *allocator : c_DefaultAllocator;
    self->N = N;

    return C_SUCCESS;
}

void c_FFT_Destroy(c_FFT_t* self) {
    if (!self) return;
    self->N = 0;
}

/* ================================================================================================================== */
/* Non-Recursive In-Place FFT Engine paired with your formal c_Complex_t APIs */

c_err_t c_FFT_Execute(const c_FFT_t* self, c_Complex_t* signal, c_bool_t inverse) {
    if (!self || !signal) return C_ERR_PARAM;
    if (self->N == 0) return C_SUCCESS;

    c_size_t n = self->N;

    /* 1. Execute the in-place bit-reversal shuffle up-front */
    c_FFT_BitReversePermutation(signal, n);

    /* 2. Bottom-up butterfly merging loops */
    /* len represents the current sub-problem size (2, 4, 8, ..., N) */
    for (c_size_t len = 2; len <= n; len <<= 1) {
        double angle = 2.0 * M_PI / (double)len * (inverse ? 1.0 : -1.0);

        /* Calculate principal twiddle factor step block using your standard constructor factory hook */
        c_Complex_t wlen = c_Complex_Make(cos(angle), sin(angle));

        for (c_size_t i = 0; i < n; i += len) {
            c_Complex_t w = c_Complex_Make(1.0, 0.0);
            c_size_t half_len = len >> 1;

            for (c_size_t j = 0; j < half_len; ++j) {
                c_size_t idx_u = i + j;
                c_size_t idx_v = i + j + half_len;

                c_Complex_t u = signal[idx_u];

                /* High-performance complex multiplication butterfly product segment: t = v * w */
                c_Complex_t t = c_Complex_Mul(signal[idx_v], w);

                /* Core butterfly update mappings leveraging clean addition and subtraction interfaces */
                signal[idx_u] = c_Complex_Add(u, t);
                signal[idx_v] = c_Complex_Sub(u, t);

                /* Advance intermediate rotation matrix step: w = w * wlen */
                w = c_Complex_Mul(w, wlen);
            }
        }
    }

    /* 3. Scale Inverse IFFT vectors uniformly by 1/N using scalar operations helper */
    if (inverse) {
        double scale = 1.0 / (double)n;
        for (c_size_t i = 0; i < n; ++i) {
            signal[i] = c_Complex_MulScalar(signal[i], scale);
        }
    }

    return C_SUCCESS;
}