/*
  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
*/

// Complex must be a complex type

template<typename Complex, typename Real> void fft( Complex *X, unsigned long N, int inv)
{
 // precompute W to improve speed -- but requires extra memory
 //
  int pre = 1; Complex c = Complex( 0, -2 * state<Complex,Real>::pi / N); if( inv) c = conj(c);

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

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

    if( !Wpre_new)
    {
      std::cerr << "fft: WARNING: realloc " << N*sizeof(Complex) << " bytes failed." << std::endl;
      pre = 0;
    }
    else
    {
      #if 0 // debug
        std::cout << "fft: computing Wpre" << std::endl;
      #endif

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

      for( unsigned long wi = 0; wi < N; ++wi) Wpre[wi] = exp(c*(Real)wi);
    }
  }
  else if( Wpre_inv != inv) // replace with conjugate values
  {
    #if 0 // debug
      std::cout << "fft: conjugate Wpre" << std::endl;
    #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)
      {
        Complex *XJ  = X + J;
        Complex *XJH = X + J + HalfSize;
        Complex W = pre ? Wpre[wi] : exp(c*(Real)wi);
        Complex t = *XJ, xjh = *XJH, xj = t + xjh;
        *XJ = xj; xjh = W*(t - xjh); *XJH = xjh;
        wi += NumOfProblems;
      }
    }
    NumOfProblems *= 2;
    ProblemSize = HalfSize;
  }
}