/* fft_time.c  time a few fft's and inverses */

#include <time.h>
#include <stdio.h>
#include <math.h>
#include "fftc.h"

int main(void)
{
  float A[8193];
  float B[8192];
  float A1[8192];
  float block = -99.99; /* sentinal   */
  double t1, t2;
  float Q;
  int i, j;
  int N;
  int rep;

  printf("fft_time.c \n");
  for(i=0; i<8193; i++) A[i] = block;

  N = 16;
  Q = 0.00001;
  rep = 20000;
  for(i=0; i<N; i++)
  {
     A[2*i]   =  (float)i/(float)N;
     A[2*i+1] = -A[2*i];
  }
  A[2] = -0.25;
  A[3] =  0.75;
  A[6] =  0.875;
  A[7] =  0.375;
  for(i=0; i<2*N; i++) B[i] = A[i];

  fft16(A);
  if(A[2*N] != block) printf("BAD outside block (%g, %g)\n", A[2*N], A[2*N+1]);
  for(j=0; j<2*N; j++)A1[j]=A[j];
  ifft16(A);
  if(A[2*N] != block) printf("BADI outside block (%g, %g)\n", A[2*N], A[2*N+1]);
  for(i=0; i<N; i++)
  {
    if( fabs(A[2*i]-B[2*i])    > Q &&
        fabs(A[2*i+1]-B[2*i+1]) > Q)
         printf("BAD %d  A=(%g, %g)  B=(%g, %g)\n",
              i, A[2*i], A[2*i+1], B[2*i], B[2*i+1]);
  } 
  t1 = (double)clock()/(double)CLOCKS_PER_SEC;
  for(i=0; i<rep; i++)
  {
    for(j=0; j<2*N; j++)A[j]=B[j];
    fft16(A);
  }
  t2 = (double)clock()/(double)CLOCKS_PER_SEC;
  t2 = (t2-t1)/(double)rep;
  printf("N=%d FFT time=%g \n", N, t2);
  t1 = (double)clock()/(double)CLOCKS_PER_SEC;
  for(i=0; i<rep; i++)
  {
    for(j=0; j<2*N; j++)A[j]=A1[j];
    ifft16(A);
  }
  t2 = (double)clock()/(double)CLOCKS_PER_SEC;
  t2 = (t2-t1)/(double)rep;
  printf("     inv time=%g \n", t2);
  t1 = (double)clock()/(double)CLOCKS_PER_SEC;
  for(i=0; i<rep; i++)
  {
    for(j=0; j<2*N; j++)A[j]=B[j];
    fft16(A);
    ifft16(A);
  }
  t2 = (double)clock()/(double)CLOCKS_PER_SEC;
  t2 = (t2-t1)/(double)rep;
  printf("    both time=%g \n\n", t2);
  if(A[2*N] != block) printf("BADB outside block (%g, %g)\n", A[2*N], A[2*N+1]);

  N = 32;
  Q = 0.00001;
  rep = 10000;
  for(i=0; i<N; i++)
  {
     A[2*i]   =  (float)i/(float)N;
     A[2*i+1] = -A[2*i];
  }
  A[2] = -0.25;
  A[3] =  0.75;
  A[6] =  0.875;
  A[7] =  0.375;
  for(i=0; i<2*N; i++) B[i] = A[i];

  fft32(A);
  if(A[2*N] != block) printf("BAD outside block (%g, %g)\n", A[2*N], A[2*N+1]);
  for(j=0; j<2*N; j++)A1[j]=A[j];
  ifft32(A);
  if(A[2*N] != block) printf("BADI outside block (%g, %g)\n", A[2*N], A[2*N+1]);
  for(i=0; i<N; i++)
  {
    if( fabs(A[2*i]-B[2*i])    > Q &&
        fabs(A[2*i+1]-B[2*i+1]) > Q)
         printf("BAD %d  A=(%g, %g)  B=(%g, %g)\n",
              i, A[2*i], A[2*i+1], B[2*i], B[2*i+1]);
  } 
  t1 = (double)clock()/(double)CLOCKS_PER_SEC;
  for(i=0; i<rep; i++)
  {
    for(j=0; j<2*N; j++)A[j]=B[j];
    fft32(A);
  }
  t2 = (double)clock()/(double)CLOCKS_PER_SEC;
  t2 = (t2-t1)/(double)rep;
  printf("N=%d FFT time=%g \n", N, t2);
  t1 = (double)clock()/(double)CLOCKS_PER_SEC;
  for(i=0; i<rep; i++)
  {
    for(j=0; j<2*N; j++)A[j]=A1[j];
    ifft32(A);
  }
  t2 = (double)clock()/(double)CLOCKS_PER_SEC;
  t2 = (t2-t1)/(double)rep;
  printf("     inv time=%g \n", t2);
  t1 = (double)clock()/(double)CLOCKS_PER_SEC;
  for(i=0; i<rep; i++)
  {
    for(j=0; j<2*N; j++)A[j]=B[j];
    fft32(A);
    ifft32(A);
  }
  t2 = (double)clock()/(double)CLOCKS_PER_SEC;
  t2 = (t2-t1)/(double)rep;
  printf("    both time=%g \n\n", t2);
  if(A[2*N] != block) printf("BADB outside block (%g, %g)\n", A[2*N], A[2*N+1]);

  N = 64;
  Q = 0.0001;
  rep = 5000;
  for(i=0; i<N; i++)
  {
     A[2*i]   =  (float)i/(float)N;
     A[2*i+1] = -A[2*i];
  }
  A[2] = -0.25;
  A[3] =  0.75;
  A[6] =  0.875;
  A[7] =  0.375;
  for(i=0; i<2*N; i++) B[i] = A[i];

  fft64(A);
  if(A[2*N] != block) printf("BAD outside block (%g, %g)\n", A[2*N], A[2*N+1]);
  for(j=0; j<2*N; j++)A1[j]=A[j];
  ifft64(A);
  if(A[2*N] != block) printf("BADI outside block (%g, %g)\n", A[2*N], A[2*N+1]);
  for(i=0; i<N; i++)
  {
    if( fabs(A[2*i]-B[2*i])    > Q &&
        fabs(A[2*i+1]-B[2*i+1]) > Q)
         printf("BAD %d  A=(%g, %g)  B=(%g, %g)\n",
              i, A[2*i], A[2*i+1], B[2*i], B[2*i+1]);
  } 
  t1 = (double)clock()/(double)CLOCKS_PER_SEC;
  for(i=0; i<rep; i++)
  {
    for(j=0; j<2*N; j++)A[j]=B[j];
    fft64(A);
  }
  t2 = (double)clock()/(double)CLOCKS_PER_SEC;
  t2 = (t2-t1)/(double)rep;
  printf("N=%d FFT time=%g \n", N, t2);
  t1 = (double)clock()/(double)CLOCKS_PER_SEC;
  for(i=0; i<rep; i++)
  {
    for(j=0; j<2*N; j++)A[j]=A1[j];
    ifft64(A);
  }
  t2 = (double)clock()/(double)CLOCKS_PER_SEC;
  t2 = (t2-t1)/(double)rep;
  printf("     inv time=%g \n", t2);
  t1 = (double)clock()/(double)CLOCKS_PER_SEC;
  for(i=0; i<rep; i++)
  {
    for(j=0; j<2*N; j++)A[j]=B[j];
    fft64(A);
    ifft64(A);
  }
  t2 = (double)clock()/(double)CLOCKS_PER_SEC;
  t2 = (t2-t1)/(double)rep;
  printf("    both time=%g \n\n", t2);
  if(A[2*N] != block) printf("BADB outside block (%g, %g)\n", A[2*N], A[2*N+1]);

  N = 128;
  Q = 0.00001;
  rep = 2000;
  for(i=0; i<N; i++)
  {
     A[2*i]   =  (float)i/(float)N;
     A[2*i+1] = -A[2*i];
  }
  A[2] = -0.25;
  A[3] =  0.75;
  A[6] =  0.875;
  A[7] =  0.375;
  for(i=0; i<2*N; i++) B[i] = A[i];

  fft128(A);
  if(A[2*N] != block) printf("BAD outside block (%g, %g)\n", A[2*N], A[2*N+1]);
  for(j=0; j<2*N; j++)A1[j]=A[j];
  ifft128(A);
  if(A[2*N] != block) printf("BADI outside block (%g, %g)\n", A[2*N], A[2*N+1]);
  for(i=0; i<N; i++)
  {
    if( fabs(A[2*i]-B[2*i])    > Q &&
        fabs(A[2*i+1]-B[2*i+1]) > Q)
         printf("BAD %d  A=(%g, %g)  B=(%g, %g)\n",
              i, A[2*i], A[2*i+1], B[2*i], B[2*i+1]);
  } 
  t1 = (double)clock()/(double)CLOCKS_PER_SEC;
  for(i=0; i<rep; i++)
  {
    for(j=0; j<2*N; j++)A[j]=B[j];
    fft128(A);
  }
  t2 = (double)clock()/(double)CLOCKS_PER_SEC;
  t2 = (t2-t1)/(double)rep;
  printf("N=%d FFT time=%g \n", N, t2);
  t1 = (double)clock()/(double)CLOCKS_PER_SEC;
  for(i=0; i<rep; i++)
  {
    for(j=0; j<2*N; j++)A[j]=A1[j];
    ifft128(A);
  }
  t2 = (double)clock()/(double)CLOCKS_PER_SEC;
  t2 = (t2-t1)/(double)rep;
  printf("     inv time=%g \n", t2);
  t1 = (double)clock()/(double)CLOCKS_PER_SEC;
  for(i=0; i<rep; i++)
  {
    for(j=0; j<2*N; j++)A[j]=B[j];
    fft128(A);
    ifft128(A);
  }
  t2 = (double)clock()/(double)CLOCKS_PER_SEC;
  t2 = (t2-t1)/(double)rep;
  printf("    both time=%g \n\n", t2);
  if(A[2*N] != block) printf("BADB outside block (%g, %g)\n", A[2*N], A[2*N+1]);

  N = 256;
  Q = 0.0001;
  rep = 1000;
  for(i=0; i<N; i++)
  {
     A[2*i]   =  (float)i/(float)N;
     A[2*i+1] = -A[2*i];
  }
  A[2] = -0.25;
  A[3] =  0.75;
  A[6] =  0.875;
  A[7] =  0.375;
  for(i=0; i<2*N; i++) B[i] = A[i];

  fft256(A);
  if(A[2*N] != block) printf("BAD outside block (%g, %g)\n", A[2*N], A[2*N+1]);
  for(j=0; j<2*N; j++)A1[j]=A[j];
  ifft256(A);
  if(A[2*N] != block) printf("BADI outside block (%g, %g)\n", A[2*N], A[2*N+1]);
  for(i=0; i<N; i++)
  {
    if( fabs(A[2*i]-B[2*i])    > Q &&
        fabs(A[2*i+1]-B[2*i+1]) > Q)
         printf("BAD %d  A=(%g, %g)  B=(%g, %g)\n",
              i, A[2*i], A[2*i+1], B[2*i], B[2*i+1]);
  } 
  t1 = (double)clock()/(double)CLOCKS_PER_SEC;
  for(i=0; i<rep; i++)
  {
    for(j=0; j<2*N; j++)A[j]=B[j];
    fft256(A);
  }
  t2 = (double)clock()/(double)CLOCKS_PER_SEC;
  t2 = (t2-t1)/(double)rep;
  printf("N=%d FFT time=%g \n", N, t2);
  t1 = (double)clock()/(double)CLOCKS_PER_SEC;
  for(i=0; i<rep; i++)
  {
    for(j=0; j<2*N; j++)A[j]=A1[j];
    ifft256(A);
  }
  t2 = (double)clock()/(double)CLOCKS_PER_SEC;
  t2 = (t2-t1)/(double)rep;
  printf("     inv time=%g \n", t2);
  t1 = (double)clock()/(double)CLOCKS_PER_SEC;
  for(i=0; i<rep; i++)
  {
    for(j=0; j<2*N; j++)A[j]=B[j];
    fft256(A);
    ifft256(A);
  }
  t2 = (double)clock()/(double)CLOCKS_PER_SEC;
  t2 = (t2-t1)/(double)rep;
  printf("    both time=%g \n\n", t2);
  if(A[2*N] != block) printf("BADB outside block (%g, %g)\n", A[2*N], A[2*N+1]);

  N = 512;
  Q = 0.0001;
  rep = 500;
  for(i=0; i<N; i++)
  {
     A[2*i]   =  (float)i/(float)N;
     A[2*i+1] = -A[2*i];
  }
  A[2] = -0.25;
  A[3] =  0.75;
  A[6] =  0.875;
  A[7] =  0.375;
  for(i=0; i<2*N; i++) B[i] = A[i];

  fft512(A);
  if(A[2*N] != block) printf("BAD outside block (%g, %g)\n", A[2*N], A[2*N+1]);
  for(j=0; j<2*N; j++)A1[j]=A[j];
  ifft512(A);
  if(A[2*N] != block) printf("BADI outside block (%g, %g)\n", A[2*N], A[2*N+1]);
  for(i=0; i<N; i++)
  {
    if( fabs(A[2*i]-B[2*i])    > Q &&
        fabs(A[2*i+1]-B[2*i+1]) > Q)
         printf("BAD %d  A=(%g, %g)  B=(%g, %g)\n",
              i, A[2*i], A[2*i+1], B[2*i], B[2*i+1]);
  } 
  t1 = (double)clock()/(double)CLOCKS_PER_SEC;
  for(i=0; i<rep; i++)
  {
    for(j=0; j<2*N; j++)A[j]=B[j];
    fft512(A);
  }
  t2 = (double)clock()/(double)CLOCKS_PER_SEC;
  t2 = (t2-t1)/(double)rep;
  printf("N=%d FFT time=%g \n", N, t2);
  t1 = (double)clock()/(double)CLOCKS_PER_SEC;
  for(i=0; i<rep; i++)
  {
    for(j=0; j<2*N; j++)A[j]=A1[j];
    ifft512(A);
  }
  t2 = (double)clock()/(double)CLOCKS_PER_SEC;
  t2 = (t2-t1)/(double)rep;
  printf("     inv time=%g \n", t2);
  t1 = (double)clock()/(double)CLOCKS_PER_SEC;
  for(i=0; i<rep; i++)
  {
    for(j=0; j<2*N; j++)A[j]=B[j];
    fft512(A);
    ifft512(A);
  }
  t2 = (double)clock()/(double)CLOCKS_PER_SEC;
  t2 = (t2-t1)/(double)rep;
  printf("    both time=%g \n\n", t2);
  if(A[2*N] != block) printf("BADB outside block (%g, %g)\n", A[2*N], A[2*N+1]);

  N = 1024;
  Q = 0.0001;
  rep = 200;
  for(i=0; i<N; i++)
  {
     A[2*i]   =  (float)i/(float)N;
     A[2*i+1] = -A[2*i];
  }
  A[2] = -0.25;
  A[3] =  0.75;
  A[6] =  0.875;
  A[7] =  0.375;
  for(i=0; i<2*N; i++) B[i] = A[i];

  fft1024(A);
  if(A[2*N] != block) printf("BAD outside block (%g, %g)\n", A[2*N], A[2*N+1]);
  for(j=0; j<2*N; j++)A1[j]=A[j];
  ifft1024(A);
  if(A[2*N] != block) printf("BADI outside block (%g, %g)\n", A[2*N], A[2*N+1]);
  for(i=0; i<N; i++)
  {
    if( fabs(A[2*i]-B[2*i])    > Q &&
        fabs(A[2*i+1]-B[2*i+1]) > Q)
         printf("BAD %d  A=(%g, %g)  B=(%g, %g)\n",
              i, A[2*i], A[2*i+1], B[2*i], B[2*i+1]);
  } 
  t1 = (double)clock()/(double)CLOCKS_PER_SEC;
  for(i=0; i<rep; i++)
  {
    for(j=0; j<2*N; j++)A[j]=B[j];
    fft1024(A);
  }
  t2 = (double)clock()/(double)CLOCKS_PER_SEC;
  t2 = (t2-t1)/(double)rep;
  printf("N=%d FFT time=%g \n", N, t2);
  t1 = (double)clock()/(double)CLOCKS_PER_SEC;
  for(i=0; i<rep; i++)
  {
    for(j=0; j<2*N; j++)A[j]=A1[j];
    ifft1024(A);
  }
  t2 = (double)clock()/(double)CLOCKS_PER_SEC;
  t2 = (t2-t1)/(double)rep;
  printf("     inv time=%g \n", t2);
  t1 = (double)clock()/(double)CLOCKS_PER_SEC;
  for(i=0; i<rep; i++)
  {
    for(j=0; j<2*N; j++)A[j]=B[j];
    fft1024(A);
    ifft1024(A);
  }
  t2 = (double)clock()/(double)CLOCKS_PER_SEC;
  t2 = (t2-t1)/(double)rep;
  printf("    both time=%g \n\n", t2);
  if(A[2*N] != block) printf("BADB outside block (%g, %g)\n", A[2*N], A[2*N+1]);

  N = 2048;
  Q = 0.0001;
  rep = 100;
  for(i=0; i<N; i++)
  {
     A[2*i]   =  (float)i/(float)N;
     A[2*i+1] = -A[2*i];
  }
  A[2] = -0.25;
  A[3] =  0.75;
  A[6] =  0.875;
  A[7] =  0.375;
  for(i=0; i<2*N; i++) B[i] = A[i];

  fft2048(A);
  if(A[2*N] != block) printf("BAD outside block (%g, %g)\n", A[2*N], A[2*N+1]);
  for(j=0; j<2*N; j++)A1[j]=A[j];
  ifft2048(A);
  if(A[2*N] != block) printf("BADI outside block (%g, %g)\n", A[2*N], A[2*N+1]);
  for(i=0; i<N; i++)
  {
    if( fabs(A[2*i]-B[2*i])    > Q &&
        fabs(A[2*i+1]-B[2*i+1]) > Q)
         printf("BAD %d  A=(%g, %g)  B=(%g, %g)\n",
              i, A[2*i], A[2*i+1], B[2*i], B[2*i+1]);
  } 
  t1 = (double)clock()/(double)CLOCKS_PER_SEC;
  for(i=0; i<rep; i++)
  {
    for(j=0; j<2*N; j++)A[j]=B[j];
    fft2048(A);
  }
  t2 = (double)clock()/(double)CLOCKS_PER_SEC;
  t2 = (t2-t1)/(double)rep;
  printf("N=%d FFT time=%g \n", N, t2);
  t1 = (double)clock()/(double)CLOCKS_PER_SEC;
  for(i=0; i<rep; i++)
  {
    for(j=0; j<2*N; j++)A[j]=A1[j];
    ifft2048(A);
  }
  t2 = (double)clock()/(double)CLOCKS_PER_SEC;
  t2 = (t2-t1)/(double)rep;
  printf("     inv time=%g \n", t2);
  t1 = (double)clock()/(double)CLOCKS_PER_SEC;
  for(i=0; i<rep; i++)
  {
    for(j=0; j<2*N; j++)A[j]=B[j];
    fft2048(A);
    ifft2048(A);
  }
  t2 = (double)clock()/(double)CLOCKS_PER_SEC;
  t2 = (t2-t1)/(double)rep;
  printf("    both time=%g \n\n", t2);
  if(A[2*N] != block) printf("BADB outside block (%g, %g)\n", A[2*N], A[2*N+1]);

  N = 4096;
  Q = 0.0001;
  rep = 50;
  for(i=0; i<N; i++)
  {
     A[2*i]   =  (float)i/(float)N;
     A[2*i+1] = -A[2*i];
  }
  A[2] = -0.25;
  A[3] =  0.75;
  A[6] =  0.875;
  A[7] =  0.375;
  for(i=0; i<2*N; i++) B[i] = A[i];

  fft4096(A);
  if(A[2*N] != block) printf("BAD outside block (%g, %g)\n", A[2*N], A[2*N+1]);
  for(j=0; j<2*N; j++)A1[j]=A[j];
  ifft4096(A);
  if(A[2*N] != block) printf("BADI outside block (%g, %g)\n", A[2*N], A[2*N+1]);
  for(i=0; i<N; i++)
  {
    if( fabs(A[2*i]-B[2*i])    > Q &&
        fabs(A[2*i+1]-B[2*i+1]) > Q)
         printf("BAD %d  A=(%g, %g)  B=(%g, %g)\n",
              i, A[2*i], A[2*i+1], B[2*i], B[2*i+1]);
  } 
  t1 = (double)clock()/(double)CLOCKS_PER_SEC;
  for(i=0; i<rep; i++)
  {
    for(j=0; j<2*N; j++)A[j]=B[j];
    fft4096(A);
  }
  t2 = (double)clock()/(double)CLOCKS_PER_SEC;
  t2 = (t2-t1)/(double)rep;
  printf("N=%d FFT time=%g \n", N, t2);
  t1 = (double)clock()/(double)CLOCKS_PER_SEC;
  for(i=0; i<rep; i++)
  {
    for(j=0; j<2*N; j++)A[j]=A1[j];
    ifft4096(A);
  }
  t2 = (double)clock()/(double)CLOCKS_PER_SEC;
  t2 = (t2-t1)/(double)rep;
  printf("     inv time=%g \n", t2);
  t1 = (double)clock()/(double)CLOCKS_PER_SEC;
  for(i=0; i<rep; i++)
  {
    for(j=0; j<2*N; j++)A[j]=B[j];
    fft4096(A);
    ifft4096(A);
  }
  t2 = (double)clock()/(double)CLOCKS_PER_SEC;
  t2 = (t2-t1)/(double)rep;
  printf("    both time=%g \n\n", t2);
  if(A[2*N] != block) printf("BADB outside block (%g, %g)\n", A[2*N], A[2*N+1]);

} /* end fft_time.c */