#include "PGFFT.h"
#include <sys/time.h>
#include <sys/resource.h>
#include <iostream>
#include <algorithm>
#include <limits>
#include <cstdlib>
#include <cstring>

#include <fftw3.h>

using namespace std;


// RandomBnd(n) returns a random number in [0..n).
// Assumes n > 0.
// FIXME: uses brain-dead rand() function
long
RandomBnd(long n)
{
   const int BPL = std::numeric_limits<unsigned long>::digits;
   const int ROTAMT = 7;
   unsigned long x = 0;
   for (long i = 0;  i < 12; i++) {
      unsigned long x1 = std::rand();
      // rotate x ROTAMT bits
      x = (x << ROTAMT) | (x >> (BPL-ROTAMT));
      x = x ^ x1;
   }

   return long(x % ((unsigned long) n));
}

void
SetSeed()
{
   srand(time(0));
}

double GetTime()
{
   struct rusage used;

   getrusage(RUSAGE_SELF, &used);
   return (used.ru_utime.tv_sec + used.ru_stime.tv_sec +
      (used.ru_utime.tv_usec + used.ru_stime.tv_usec) / 1e6);
}

#define TIME_IT(t, action) \
do { \
   double _t0, _t1; \
   long _iter = 1; \
   long _cnt = 0; \
   do { \
      _t0 = GetTime(); \
      for (long _i = 0; _i < _iter; _i++) { action; _cnt++; } \
      _t1 = GetTime(); \
   } while ( _t1 - _t0 < 3 && (_iter *= 2)); \
   t = (_t1 - _t0)/_iter; \
} while(0)


typedef complex<double> cmplx_t;





int main()
{
   SetSeed();

   vector<long> nvec = { 1*1024, 4*1024, 16*1024, 64*1024, 
                         1409,   6007,   24001, 96001 };

   long N = nvec.size();
   vector<double> avec(N), bvec(N);


   for (long i = 0; i < N; i++) {
      long n = nvec[i];

      std::cout << "n=" << n << "\n";

      vector<cmplx_t> v(n);
      for (long i = 0; i < n; i++)
	 v[i] = RandomBnd(10) - 5;

      vector<cmplx_t> w(n);

      PGFFT pgfft(n);

      double t, t_min;
      TIME_IT(t, pgfft.apply(&v[0], &w[0])); 
      t_min = t;
      TIME_IT(t, pgfft.apply(&v[0], &w[0])); 
      t_min = min(t_min, t);
      TIME_IT(t, pgfft.apply(&v[0], &w[0])); 
      t_min = min(t_min, t);

      std::cout << t_min << "\n";
      avec[i] = t_min;

      fftw_complex *in, *out;
      fftw_plan p;
      in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n);
      out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n);
      p = fftw_plan_dft_1d(n, in, out, FFTW_FORWARD, FFTW_ESTIMATE);

      std::memcpy(in, &v[0], sizeof(fftw_complex) * n);

      TIME_IT(t, fftw_execute(p));
      t_min = t;
      TIME_IT(t, fftw_execute(p));
      t_min = min(t_min, t);
      TIME_IT(t, fftw_execute(p));
      t_min = min(t_min, t);

      std::cout << t_min << "\n";
      bvec[i] = t_min;

      fftw_destroy_plan(p);
      fftw_free(in);
      fftw_free(out);
   }

   printf("n          ");
   for (long i = 0; i < N; i++) printf("%9d", nvec[i]);
   printf("\n");

   printf("PGFFT (ms) ");
   for (long i = 0; i < N; i++) printf("%9.1e", avec[i]*1000);
   printf("\n");

   printf("FFTW  (ms) ");
   for (long i = 0; i < N; i++) printf("%9.1e", bvec[i]*1000);
   printf("\n");





}
