// fft_verify.c - correctness check for FFT(): compare vs naive DFT, and check
// the forward->inverse round-trip recovers the input.
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include "../src/fft.h"

// Reference DFT accumulated in double precision so the *reference* isn't the
// source of error when comparing a float FFT at large n.
static void naive_dft(const float complex* x, float complex* X, int n) {
    for (int k = 0; k < n; k++) {
        double complex s = 0;
        for (int t = 0; t < n; t++)
            s += x[t] * cexp(-2.0 * M_PI * I * k * t / n);
        X[k] = (float complex)s;
    }
}

int main(void) {
    int sizes[] = {2, 4, 8, 16, 64, 256, 1024, 2048};
    double worst_dft = 0, worst_rt = 0;
    for (unsigned si = 0; si < sizeof(sizes)/sizeof(sizes[0]); si++) {
        int n = sizes[si];
        float complex* x   = malloc(n*sizeof(float complex));
        float complex* X   = malloc(n*sizeof(float complex));
        float complex* Xn  = malloc(n*sizeof(float complex));
        float complex* xrt = malloc(n*sizeof(float complex));
        srand(1234 + n);
        for (int i = 0; i < n; i++)
            x[i] = (rand()/(float)RAND_MAX - 0.5f) + (rand()/(float)RAND_MAX - 0.5f)*I;

        FFT(x, X, n, false);
        naive_dft(x, Xn, n);
        FFT(X, xrt, n, true);

        double norm = 0;
        for (int i = 0; i < n; i++) norm += cabsf(Xn[i]);
        norm /= n; // mean magnitude, for a relative tolerance
        for (int i = 0; i < n; i++) {
            double e1 = cabsf(X[i] - Xn[i]) / norm; // relative to spectrum scale
            double e2 = cabsf(xrt[i] - x[i]);
            if (e1 > worst_dft) worst_dft = e1;
            if (e2 > worst_rt)  worst_rt = e2;
        }
        free(x); free(X); free(Xn); free(xrt);
    }
    printf("max relative |FFT - naive DFT| = %.3e\n", worst_dft);
    printf("max round-trip error           = %.3e\n", worst_rt);
    printf("%s\n", (worst_dft < 1e-4 && worst_rt < 1e-4) ? "PASS" : "FAIL");
    return (worst_dft < 1e-4 && worst_rt < 1e-4) ? 0 : 1;
}