/*
  in-place FFT and inverse FFT, 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 len(X) = 2**n

  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 int n, int inv)
{
  unsigned long N = 1UL << n;

 // 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 = std::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*) std::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] = std::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] = std::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] : std::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;
  }

 // bit-reverse and scale
 // see https://stackoverflow.com/questions/932079/in-place-bit-reversed-shuffle-on-an-array
 //
  Real s = 1/std::sqrt(N);
  for( unsigned long I = 0, J = 0; I < N; ++I) {
    J = bitrev( I, n); Complex t, *XI = X + I, *XJ  = X + J;
    if( I < J) { t = *XI; *XI = *XJ; *XJ = t; *XI *= s; *XJ *= s; }
    else if( I == J) *XI *= s;
  }
}