typedef long clock_t;
     const double RESOLUTION_DEFAULT = 2.0;
     const int RANDOM_SEED = 101010;



     const int FFT_SIZE = 1024;
     const int SOR_SIZE =100;
     const int SPARSE_SIZE_M = 1000;
     const int SPARSE_SIZE_nz = 5000;
     const int LU_SIZE = 100;



     const int LG_FFT_SIZE = 1048576;
     const int LG_SOR_SIZE =1000;
     const int LG_SPARSE_SIZE_M = 100000;
     const int LG_SPARSE_SIZE_nz =1000000;
     const int LG_LU_SIZE = 1000;





     const int TINY_FFT_SIZE = 16;
     const int TINY_SOR_SIZE =10;
     const int TINY_SPARSE_SIZE_M = 10;
     const int TINY_SPARSE_SIZE_N = 10;
     const int TINY_SPARSE_SIZE_nz = 50;
     const int TINY_LU_SIZE = 10;




typedef struct
{
  int m[17];
  int seed;
  int i;
  int j;
  int haveRange;
  double left;
  double right;
  double width;
}
Random_struct, *Random;

Random new_Random_seed(int seed);
double Random_nextDouble(Random R);
void Random_delete(Random R);
double *RandomVector(int N, Random R);
double **RandomMatrix(int M, int N, Random R);






typedef struct{
    int running;
    double last_time;
    double total;

} *Stopwatch, Stopwatch_struct;



double seconds();

void Stopwtach_reset(Stopwatch Q);

Stopwatch new_Stopwatch(void);
void Stopwatch_delete(Stopwatch S);
void Stopwatch_start(Stopwatch Q);
void Stopwatch_resume(Stopwatch Q);
void Stopwatch_stop(Stopwatch Q);
double Stopwatch_read(Stopwatch Q);



void FFT_transform(int N, double *data);
void FFT_inverse(int N, double *data);
void FFT_bitreverse(int N, double *data);
double FFT_num_flops(int N);
double seconds()
{
        return ((double) clock()) / (double) ((clock_t)1000);
}

void Stopwtach_reset(Stopwatch Q)
{
    Q->running = 0;
    Q->last_time = 0.0;
    Q->total= 0.0;
}


Stopwatch new_Stopwatch(void)
{
    Stopwatch S = (Stopwatch) malloc(sizeof(Stopwatch_struct));
    if (S == 0)
        return 0;

    Stopwtach_reset(S);
    return S;
}

void Stopwatch_delete(Stopwatch S)
{
    if (S != 0)
        free(S);
}




void Stopwatch_start(Stopwatch Q)
{
    if (! (Q->running) )
    {
        Q->running = 1;
        Q->total = 0.0;
        Q->last_time = seconds();
    }
}







void Stopwatch_resume(Stopwatch Q)
{
    if (!(Q->running))
    {
        Q-> last_time = seconds();
        Q->running = 1;
    }
}

void Stopwatch_stop(Stopwatch Q)
{
    if (Q->running)
    {
        Q->total += seconds() - Q->last_time;
        Q->running = 0;
    }
}


double Stopwatch_read(Stopwatch Q)
{

    if (Q->running)
    {
        double t = seconds();
        Q->total += t - Q->last_time;
        Q->last_time = t;
    }
    return Q->total;
}



static const int m1 = (1 << (32 -2)) + ((1 << (32 -2) )-1);
static const int m2 = 1 << 32/2;
static double dm1;




static void initialize(Random R, int seed);

Random new_Random_seed(int seed)
{
    Random R = (Random) malloc(sizeof(Random_struct));

    initialize(R, seed);
    R->left = 0.0;
    R->right = 1.0;
    R->width = 1.0;
    R->haveRange = 0 ;

    return R;
}

Random new_Random(int seed, double left, double right)
{
    Random R = (Random) malloc(sizeof(Random_struct));

    initialize(R, seed);
    R->left = left;
    R->right = right;
    R->width = right - left;
    R->haveRange = 1;

    return R;
}

void Random_delete(Random R)
{
    free(R);
}





double Random_nextDouble(Random R)
{
    int k;

    int I = R->i;
    int J = R->j;
    int *m = R->m;

    k = m[I] - m[J];
    if (k < 0) k += m1;
    R->m[J] = k;

    if (I == 0)
        I = 16;
    else I--;
    R->i = I;

    if (J == 0)
        J = 16 ;
    else J--;
    R->j = J;

    if (R->haveRange)
        return R->left + dm1 * (double) k * R->width;
    else
        return dm1 * (double) k;

}
static void initialize(Random R, int seed)
{

    int jseed, k0, k1, j0, j1, iloop;

    dm1 = 1.0 / (double) m1;

    R->seed = seed;

    if (seed < 0 ) seed = -seed;
    jseed = (seed < m1 ? seed : m1);
    if (jseed % 2 == 0) --jseed;
    k0 = 9069 % m2;
    k1 = 9069 / m2;
    j0 = jseed % m2;
    j1 = jseed / m2;
    for (iloop = 0; iloop < 17; ++iloop)
    {
        jseed = j0 * k0;
        j1 = (jseed / m2 + j0 * k1 + j1 * k0) % (m2 / 2);
        j0 = jseed % m2;
        R->m[iloop] = j0 + m2 * j1;
    }
    R->i = 4;
    R->j = 16;

}

double *RandomVector(int N, Random R)
{
    int i;
    double *x = (double *) malloc(sizeof(double)*N);

    for (i=0; i<N; i++)
        x[i] = Random_nextDouble(R);

    return x;
}


double **RandomMatrix(int M, int N, Random R)
{
    int i;
    int j;



    double **A = (double **) malloc(sizeof(double*)*M);

    if (A == 0) return 0;

    for (i=0; i<M; i++)
    {
        A[i] = (double *) malloc(sizeof(double)*N);
        if (A[i] == 0)
        {
            free(A);
            return 0;
        }
        for (j=0; j<N; j++)
            A[i][j] = Random_nextDouble(R);
    }
    return A;
}




static int int_log2(int n);

double FFT_num_flops(int N)
{

     double Nd = (double) N;
     double logN = (double) int_log2(N);

     return (5.0*Nd-2)*logN + 2*(Nd+1);
}

static int int_log2 (int n)
{
    int k = 1;
    int log = 0;
    for(k = 1; k < n; k *= 2, log++);
    if (n != (1 << log))
    {
      printf("FFT: Data length is not a power of 2!: %d ",n);
      exit(1);
    }
    return log;
}

static void FFT_transform_internal (int N, double *data, int direction) {
    int n = N/2;
    int bit = 0;
    int logn;
    int dual = 1;

    if (n == 1) return;
    logn = int_log2(n);


    if (N == 0) return;


    FFT_bitreverse(N, data) ;



    for (bit = 0; bit < logn; bit++, dual *= 2) {
      double w_real = 1.0;
      double w_imag = 0.0;
      int a;
      int b;

      double theta = 2.0 * direction * 3.1415926535897932 / (2.0 * (double) dual);
      double s = sin(theta);
      double t = sin(theta / 2.0);
      double s2 = 2.0 * t * t;

      for (a=0, b = 0; b < n; b += 2 * dual) {
        int i = 2*b ;
        int j = 2*(b + dual);

        double wd_real = data[j] ;
        double wd_imag = data[j+1] ;

        data[j] = data[i] - wd_real;
        data[j+1] = data[i+1] - wd_imag;
        data[i] += wd_real;
        data[i+1]+= wd_imag;
      }


      for (a = 1; a < dual; a++) {

        {
          double tmp_real = w_real - s * w_imag - s2 * w_real;
          double tmp_imag = w_imag + s * w_real - s2 * w_imag;
          w_real = tmp_real;
          w_imag = tmp_imag;
        }
        for (b = 0; b < n; b += 2 * dual) {
          int i = 2*(b + a);
          int j = 2*(b + a + dual);

          double z1_real = data[j];
          double z1_imag = data[j+1];

          double wd_real = w_real * z1_real - w_imag * z1_imag;
          double wd_imag = w_real * z1_imag + w_imag * z1_real;

          data[j] = data[i] - wd_real;
          data[j+1] = data[i+1] - wd_imag;
          data[i] += wd_real;
          data[i+1]+= wd_imag;
        }
      }
    }
  }


void FFT_bitreverse(int N, double *data) {

        int n=N/2;
        int nm1 = n-1;
        int i=0;
        int j=0;
        for (i = 0; i < nm1; i++) {


                int ii = i << 1;


                int jj = j << 1;


                int k = n >> 1;

                if (i < j) {
                        double tmp_real = data[ii];
                        double tmp_imag = data[ii+1];
                        data[ii] = data[jj];
                        data[ii+1] = data[jj+1];
                        data[jj] = tmp_real;
                        data[jj+1] = tmp_imag; }

                while ( (k < j) || (k == j) )
                {

                        j -= k;


                        k >>= 1 ;
                }
                j += k ;
        }
}


void FFT_transform(int N, double *data)
{
    FFT_transform_internal(N, data, -1);
}


void FFT_inverse(int N, double *data)
{
    int n = N/2;
    double norm = 0.0;
    int i=0;
    FFT_transform_internal(N, data, +1);




    norm=1/((double) n);
    for(i=0; i<N; i++)
      data[i] *= norm;

}


double kernel_measureFFT(
                                                 int N,
                                                 Random R,
                                                 long cycles,
                                                 double * time,
                                                 double * mflops
                                                 )
{


        int twoN = 2*N;
        double *x = RandomVector(twoN, R);
        Stopwatch Q = new_Stopwatch();
        int i=0;
        double result = 0.0;

        Stopwatch_start(Q);

        for (i=0; i<cycles; i++)
        {
                FFT_transform(twoN, x);
                FFT_inverse(twoN, x);
        }

        Stopwatch_stop(Q);

        *time = Stopwatch_read(Q);
        *mflops = FFT_num_flops(N)*cycles/ Stopwatch_read(Q) * 1.0e-6;
        Stopwatch_delete(Q);
        free(x);
        return result;
}

int main(int argc, char** argv)
{
        double time = 0.0;
        double mflops = 0.0;
        int FFT_size = LG_FFT_SIZE;
        long cycles = 10;
        Random R = new_Random_seed(RANDOM_SEED);

        kernel_measureFFT( FFT_size, R, cycles, &time, &mflops );

        printf("cciftFFT (size %d) (cycles %d)\n", FFT_size, cycles);
        printf("\t time:  %8.2f seconds\n", time);
        printf("\t perf:  %8.2f mflops\n", mflops);

}
