#include <cmath>
#include <complex>
#include <testsuite_performance.h>

const double pi = 4 * std::atan(1.0);

void DFT(std::complex<double>* f, int N, bool inverse=false) 
{
  const std::complex<double> i(0.0, 1.0);
  
  std::complex<double>* g = new std::complex<double>[N];
  for (int n = 0; n < N; n++)
    g[n] = f[n];
  
  for (int k = 0; k < N; k++) {
    f[k] = 0.0;
    for (int n = 0; n < N; n++) {
      double theta = 2 * pi * n * k / double(N);
      if (inverse)
	theta = - theta;
      f[k] += g[n] * std::exp(- i * theta);
    }
    f[k] /= std::sqrt(double(N));
  }
  delete [] g;
}

int main() 
{
  /*
    std::cout << " Discrete Fourier Transform of f(t)=t(1-t) in [0,1]\n"
              << " --------------------------------------------------\n"
              << " Enter number of points N: ";
    int N;
    std::cin >> N;
  */

  int N = 5000;

  __gnu_test::time_counter time;
  __gnu_test::resource_counter resource;
  __gnu_test::start_counters(time, resource);

  std::complex<double> *f = new std::complex<double> [N];
  double dt = 1 / double(N - 1);
  for (int n = 0; n < N; n++) {
    double t = n * dt;
    f[n] = t * (1 - t);
  }

  for (int n = 0; n < N; n++) 
    {
      double t = n * dt;
      t = std::real(f[n]);
    }

  DFT(f, N);                     // compute DFT        

  for (int k = 0; k < N; k++) {
    double omega = 2 * k * pi / (N * dt);
    omega = f[k].real();
  }
  for (int k = 0; k < N; k++) {
    double omega = 2 * k * pi / (N * dt);
    omega = f[k].imag();
  }

  bool inverse = true;           // compute inverse DFT
  DFT(f, N, inverse);

  for (int n = 0; n < N; n++) {
    double t = n * dt;
    t = std::real(f[n]);
  }
  delete [] f;

  __gnu_test::stop_counters(time, resource);
  __gnu_test::report_performance(__FILE__, "", time, resource);
}