/*
  in-place FFT and inverse FFT, with bit-reversed result locations

  R. Perry, Jan. 2018

  Algorithm 4.2 from "Inside the FFT Black Box: Serial and Parallel
  Fast Fourier Transform Algorithms", Eleanor Chu and Alan George,
  CRC Press, 1999

  requires N = len(X) = power of 2

  does not include the 1/sqrt(N) scale factor

  Adapted from the "FFT in C for Python" code originally written for use in
  the Bleichenbacher attack: http://fog.misty.com/perry/Bleichenbacher/notes.html
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>

void fft( double complex *X, unsigned long N, int inv)
{
 // precompute W to improve speed -- but requires extra memory
 //
  int pre = 1;  double pi = 4*atan(1.0);  double complex c = -2*I*pi/N;  if( inv) c = conj(c);

  static unsigned long Wpre_N = 0;  static int Wpre_inv = 0;  static double complex *Wpre = NULL;

  if( !Wpre || Wpre_N != N) // create or resize
  {
    double complex *Wpre_new = realloc( Wpre, N*sizeof(double complex));

    if( !Wpre_new)
    {
      printf( "fft: WARNING: realloc %lu bytes failed.\n", N*sizeof(double complex));
      pre = 0;
    }
    else
    {
      #if 0 // debug
        printf( "fft: computing Wpre\n");
      #endif

      Wpre = Wpre_new;  Wpre_N = N;  Wpre_inv = inv;

      for( unsigned long wi = 0; wi < N; ++wi) Wpre[wi] = cexp(c*wi);
    }
  }
  else if( Wpre_inv != inv) // replace with conjugate values
  {
    #if 0 // debug
      printf( "fft: conjugate Wpre\n");
    #endif

    for( unsigned long wi = 0; wi < N; ++wi) Wpre[wi] = conj(Wpre[wi]);

    Wpre_inv = inv;
  }

 // perform FFT
 //
  unsigned long ProblemSize = N, NumOfProblems = 1;

  while( ProblemSize > 1)
  {
    unsigned long HalfSize = ProblemSize/2;

    for( unsigned long K = 0; K < NumOfProblems; ++K)
    {
      unsigned long JFirst = K*ProblemSize, JLast = JFirst + HalfSize, wi = 0;

      for( unsigned long J = JFirst; J < JLast; ++J)
      {
        double complex *XJ  = X + J;
        double complex *XJH = X + J + HalfSize;
        double complex W = pre ? Wpre[wi] : cexp(c*wi);
        double complex t = *XJ, xjh = *XJH, xj = t + xjh;
        *XJ = xj; xjh = W*(t - xjh); *XJH = xjh;
        wi += NumOfProblems;
      }
    }
    NumOfProblems *= 2;
    ProblemSize = HalfSize;
  }
}