#include <stdio.h>
#include <math.h>
#include <time.h>
#include <float.h>
#define NZ 1048576
#define N2 NZ/2
#define FLT_REL 0.11920928955078e-06
static float xx[NZ][2],zz[NZ][2];
/* test driver for complex FFT routine ccfft2. WPP 28 June 1996 */
main()
{
   int n,n2exp,i,nit,kk,nexec;
/*
  repetition count:
                  16,   32,   64,  128,  256, 512, 1024, 2048,
                   |     |     |     |     |    |     |     | */
   int NX[17]={1000 ,1000 ,1000 ,1000 ,1000 ,100 , 100 , 100 ,
           100,  100,  10 ,  10 ,    1,     1,     1,     1,      1};
/*           |     |     |     |     |      |      |      |       |
           4096,8192,16384,32768,65536,131072,262144,524288,1048576 */

   float eps,z0,z1,fnm1,errCC,ggl(float*);
   float *x;
   double ln2,mflops,t0,t1,t2,t3,t4,walltime(double*);
   static float seed=331.0;
   int sign;
   void ccfft2(), rcfft2(), crfft2(), refft2(), rofft2();

/* print headers */

   printf("  ----------------------------------- \n");
   printf("     loop inversion variant of CFFT2   \n");
   printf("  ------------------------------------ \n");

   for(n2exp=4;n2exp<21;n2exp++){
      n     = pow(2,n2exp);
      eps   = ((float)n)*FLT_REL;
      printf("\n N=%d tests, rms error criterion=%8.2e: \n",n,eps);
      printf(" ------------------------------------------\n",n);
      t3    = 100000.;
      t4    = 100000.;
      nexec = 1;
      for(nit=0;nit<10;nit++){
         if(nit==0){
            for(i=0;i<n;i++){
               z0       = ggl(&seed);
               z1       = ggl(&seed);
               xx[i][0] = z0;
               xx[i][1] = z1;
               zz[i][0] = z0;
               zz[i][1] = z1;
            }
         } else {
            if(nit==1){
               for(i=0;i<n;i++){
                  z0       = 0.;
                  z1       = 0.;
                  xx[i][0] = z0;
                  xx[i][1] = z1;
                  zz[i][0] = z0;
                  zz[i][1] = z1;
               }
            }
         }
/* forward transform */

         x    = (float *) xx;
         sign = 1;
         t0   = 0.0;
         t1 = walltime(&t0);
         for(kk=0;kk<nexec;kk++){
             ccfft2(n,x,sign);
         }
         t2 = walltime(&t1);
         if(nit != 0) t3 = (t3<t2)?t3:t2;

/* backward transform */ 

         sign = -1;
         t0   = 0.0;
         t1   = walltime(&t0);
         for(kk=0;kk<nexec;kk++){
             ccfft2(n,x,sign);
         }
         t2   = walltime(&t1);
         if(nit != 0) t4 = (t4<t2)?t4:t2;

/* first pass through, check errors */ 

         if(nit == 0){
            fnm1  = 1./((float) n);
            for(i=0;i<n;i++){
                zz[i][0] = fnm1*xx[i][0] - zz[i][0];
                zz[i][1] = fnm1*xx[i][1] - zz[i][1];
            }
            errCC = 0.0;
            for(i=0;i<n;i++){
                errCC = errCC + zz[i][0]*zz[i][0]+ zz[i][1]*zz[i][1];
            }
            errCC = sqrt(fnm1*errCC);
            if(errCC>eps){
               printf("\n error in CCFFT F/B: errCC=%e\n",errCC);
            } else {
               printf("\n     CCFFT F/B test OK\n");
            }
         }
         nexec = NX[n2exp-4];
      }
      t1     = 0.5*(t3+t4)/((float) nexec);
      t1     = (t1>0.0)?t1:0.0;
      ln2    = (double) n2exp;
      if(t1>0.0){
         mflops = 5.0*((double) n)*ln2/((1.e+6)*t1);
         printf("        Time=%8.2e, Mflops=%6.1f for ccfft2 on %d points\n",
         t1,mflops,n);
         printf("\n");
      }
   }
}
float ggl(float *ds)
{

/* generate u(0,1) distributed random numbers. 
   Seed ds must be saved between calls. ggl is 
   essentially the same as the IMSL routine RNUM. 

   W. Petersen and M. Troyer, 24 Oct. 2002, ETHZ: 
   a modification of a fortran version from 
   I. Vattulainen, Tampere Univ. of Technology, 
   Finland, 1992 */

   double t,d2=0.2147483647e10;
   t   = (float) *ds;
   t   = fmod(0.16807e5*t,d2);
   *ds = (float) t;
   return((float) ((t-1.0e0)/(d2-1.0e0)));
}

/* ======================================= */
/*    split here for library functions     */
/* ======================================= */

void cffti(int n,float *w)
/* 
   Initialization routine for binary radix FFTs.
 
   OUTPUT
         w = complex n/2 dimensional array,
           = exp(2*Pi*i*j/n), for j=0..n/2-1
 
         i.e. n/2 powers of the n-th root of unity
 
   W Petersen, SAM, Mathmatik, ETHZ, 25 Aug., 1999
*/
{
   int i,ii;
   float pi,arg,t;

   pi = 3.141592653589793;
   t  = 2.0*pi/((float) n);
   ii = 0;
   for(i=0;i<n;i+=2){
      arg     = t*((float) ii++);
      w[i  ] = cos(arg);
      w[i+1] = sin(arg);
   }
}
void ccfft2(int n,float *x, int sign)
/*
 
  complex, in-place, in-order binary radix FFT
 
   INPUT
 
        n = power of 2 integer > 1
 
        x = complex n dimensional array
 
        sign = integer sign of transform (+- below) = +1, or -1.
 
   OUTPUT (where w = n-th root of unity)
 
        x_j <- sum(k=0..n-1) w^{+- jk} x_k, for j=0..n-1
 
   W Petersen, SAM, Mathematik, ETHZ, 25 Aug., 1999
*/
 
{
    void cffti(), cfft_2();
    static float *w,*pw1,*pw2;
    static int n1=-33,n2=-33;
    int nt;
    if(n==n1){
       w = pw1;
    } else {
       if(n==n2){
          w   = pw2;
          pw2 = pw1;
          pw1 = w;
          nt  = n1;
          n1  = n2;
          n2  = nt;
       } else {
          if(n2>0) free(pw2);
          w   = (void *) malloc(n*sizeof(float));
          cffti(n,w);
          pw2 = pw1;
          pw1 = w;
          n2  = n1;
          n1  = n;
       }
    }
    cfft_2(n,x,w,1,sign);
}
#define max(a,b) (a>b)?a:b
void cfft_2(n,x,w,iw,sign)
int n, iw,sign;
float x[][2], w[];
/*
 
  in-place, in-order binary radix FFT from C. Temperton,
  SIAM J. on Sci. and Stat. Computing, vol 12, no. 4,
  July 1991. This is has a loop-inversion optimization: 
  first log2(n)/2 passes are Cooley-Tukey, last log2(n)/2 
  passes compute 'bugs' in pairs. Two variants of the
  self-ordering phase are possible - step1 and step2. A
  schedule for choosing step1 or step2 is given by the
  empirical choice given by the formula for various machines
  in the comments. Pick your favorite. Parameter iw=stride in 
  array w, used for real/complex routines (iw=2 in that case).
 
  W Petersen, SAM, Mathematik, ETHZ, 25 Aug., 1999

*/
 

{
   int n2,m,j,mj,p2,p3,p4,BK;
   void step0(), step1(), step2();

   m  = (int) (log((float) n)/log(1.99));
/* 
optimization point for Y-MP:     
   BK = max(m/3-2,0);
optimization point for J-90            
   BK = max((m+1)/6,0);
optimization point for IBM:    
   BK = max((m-3)/2,0);
optimization point for DEC Alpha:                             
   BK = max(m/2-2,0);
optimization point for HP9000                     
   BK = max(m/2-3,0);
optimization point for SGI Indigo
   BK = max((m-1)/4,0);
optimization point for SGI TFP (R8000)
   BK = max((m-1)/2,0);
optimization point for i860/XP                                  
   BK = max((20*(m-5)-(m-5)*(m-5))/25),0);
optimization point for Intel P-6
   BK = max((7+5*(m-4))/10+1,0);
optimization point for Sparc-10
   BK = max(0,2*m/3-4); 
optimization point for Mac Power-PC
   BK = 4*m/9
plain vanilla: */ 
   BK = 0; 

   mj = 1;
   n2 = n/2;

   for(j=0;j<m;j++){
      if(j < (m+1)/2){
         p2 = n2/mj;
         step0(n,mj,&x[0][0],&x[p2][0],w,iw,sign);
      }
      else{
         p2 = n2/mj;
         p3 = mj;
         p4 = p2+mj;
         if(j < (m/2+BK)){
            step1(n,mj,&x[0][0],&x[p2][0],&x[p3][0],&x[p4][0],w,iw,sign);
         } 
         else{
            step2(n,mj,&x[0][0],&x[p2][0],&x[p3][0],&x[p4][0],w,iw,sign);
         }
      }
      mj = 2*mj;
   }
}
void step0(n,mj,a,b,w,iw,sign)
int iw,n,mj;
float a[][2],b[][2],w[][2];
int sign;
{
   float wkr,wku,wambr,wambu;
   int i,k,ks,kw,lj,ii,ij;

   lj = n/(2*mj);
   ij = n/mj;
   ks = iw*mj;
#pragma omp parallel shared(w,a,b,ij,ks,lj,mj,sign) \
 private(i,k,ii,kw,wkr,wku,wambr,wambu)
#pragma omp for nowait
   for(i=0;i<mj;i++){
      ii = i*ij;
      if(sign > 0){
#pragma ivdep
         for(k=0;k<lj;k++){
             kw    = k*ks;
             wkr   = w[kw][0];
             wku   = w[kw][1];
             wambr = wkr*(a[ii+k][0]-b[ii+k][0]) - wku*(a[ii+k][1]-b[ii+k][1]);
             wambu = wku*(a[ii+k][0]-b[ii+k][0]) + wkr*(a[ii+k][1]-b[ii+k][1]);
             a[ii+k][0] = a[ii+k][0]+b[ii+k][0];
             a[ii+k][1] = a[ii+k][1]+b[ii+k][1];
             b[ii+k][0] = wambr;
             b[ii+k][1] = wambu;
         } 
      } else {
#pragma ivdep
         for(k=0;k<lj;k++){
             kw    = k*ks;
             wkr   = w[kw][0];
             wku   = -w[kw][1];
             wambr = wkr*(a[ii+k][0]-b[ii+k][0]) - wku*(a[ii+k][1]-b[ii+k][1]);
             wambu = wku*(a[ii+k][0]-b[ii+k][0]) + wkr*(a[ii+k][1]-b[ii+k][1]);
             a[ii+k][0] = a[ii+k][0]+b[ii+k][0];
             a[ii+k][1] = a[ii+k][1]+b[ii+k][1];
             b[ii+k][0] = wambr;
             b[ii+k][1] = wambu;
         } 
      }
   }
}
void step1(n,mj,a,b,c,d,w,iw,sign)
int iw,n,mj;
float a[][2],b[][2],c[][2],d[][2],w[][2];
int sign;
{
   float wkr,wku,wambr,wambu,wcmdr,wcmdu;
   int mj2,i,j,k,ks,kw,lj;

   mj2 = 2*mj;
   lj  = n/mj2;
   ks  = iw*mj;

#pragma omp parallel shared(w,a,b,c,d,ks,n,lj,mj,mj2,sign) \
 private(i,j,k,kw,wkr,wku,wambr,wambu,wcmdr,wcmdu)
#pragma omp for nowait
   for(j=0;j<mj;j+=n/mj){
      if(sign > 0){
         for(i=j;i<n;i+=mj2){
#pragma ivdep
            for(k=0;k<lj;k++){
                kw    = k*ks;
                wkr   = w[kw][0];
                wku   = w[kw][1];
                wambr = wkr*(a[i+k][0]-b[i+k][0]) - wku*(a[i+k][1]-b[i+k][1]);
                wambu = wku*(a[i+k][0]-b[i+k][0]) + wkr*(a[i+k][1]-b[i+k][1]);
                a[i+k][0] = a[i+k][0]+b[i+k][0];
                a[i+k][1] = a[i+k][1]+b[i+k][1];
                b[i+k][0] = c[i+k][0]+d[i+k][0];
                b[i+k][1] = c[i+k][1]+d[i+k][1];
                wcmdr = wkr*(c[i+k][0]-d[i+k][0]) - wku*(c[i+k][1]-d[i+k][1]);
                wcmdu = wku*(c[i+k][0]-d[i+k][0]) + wkr*(c[i+k][1]-d[i+k][1]);
                c[i+k][0] = wambr;
                c[i+k][1] = wambu;
                d[i+k][0] = wcmdr;
                d[i+k][1] = wcmdu;
            } 
         } 
      } else {
         for(i=j;i<n;i+=mj2){
#pragma ivdep
            for(k=0;k<lj;k++){
                kw    = k*ks;
                wkr   = w[kw][0];
                wku   = -w[kw][1];
                wambr = wkr*(a[i+k][0]-b[i+k][0]) - wku*(a[i+k][1]-b[i+k][1]);
                wambu = wku*(a[i+k][0]-b[i+k][0]) + wkr*(a[i+k][1]-b[i+k][1]);
                a[i+k][0] = a[i+k][0]+b[i+k][0];
                a[i+k][1] = a[i+k][1]+b[i+k][1];
                b[i+k][0] = c[i+k][0]+d[i+k][0];
                b[i+k][1] = c[i+k][1]+d[i+k][1];
                wcmdr = wkr*(c[i+k][0]-d[i+k][0]) - wku*(c[i+k][1]-d[i+k][1]);
                wcmdu = wku*(c[i+k][0]-d[i+k][0]) + wkr*(c[i+k][1]-d[i+k][1]);
                c[i+k][0] = wambr;
                c[i+k][1] = wambu;
                d[i+k][0] = wcmdr;
                d[i+k][1] = wcmdu;
            }
         }
      }
   }
}
void step2(n,mj,a,b,c,d,w,iw,sign)
int iw,n,mj;
float a[][2],b[][2],c[][2],d[][2],w[][2];
int sign;
{
   float wkr,wku,wambr,wambu,wcmdr,wcmdu;
   int mj2,i,j,k,ks,kw,lj,ii;

   mj2 = 2*mj;
   lj  = n/mj2;
   ks  = iw*mj;

   if(mj==(n/2)){
#pragma ivdep
      for(j=0;j<mj;j+=n/mj){
          wambr   = a[j][0]-b[j][0];
          wambu   = a[j][1]-b[j][1];
          wcmdr   = c[j][0]-d[j][0];
          wcmdu   = c[j][1]-d[j][1];
          a[j][0] = a[j][0]+b[j][0];
          a[j][1] = a[j][1]+b[j][1];
          b[j][0] = c[j][0]+d[j][0];
          b[j][1] = c[j][1]+d[j][1];
          c[j][0] = wambr;
          c[j][1] = wambu;
          d[j][0] = wcmdr;
          d[j][1] = wcmdu;
      }
      return; 
   } 
#pragma omp parallel shared(w,a,b,c,d,ks,lj,mj,mj2,n,sign) \
 private(i,j,k,ii,kw,wkr,wku,wambr,wambu,wcmdr,wcmdu)
#pragma omp for nowait
   for(k=0;k<lj;k++){
      kw    = k*ks;
      wkr   = w[kw][0];
      wku = (sign > 0)?w[kw][1]:(-w[kw][1]);
      for(i=0;i<lj;i++){
         ii = i*mj2;
#pragma ivdep
         for(j=k;j<mj;j+=n/mj){
             wambr = wkr*(a[ii+j][0]-b[ii+j][0]) - wku*(a[ii+j][1]-b[ii+j][1]);
             wambu = wku*(a[ii+j][0]-b[ii+j][0]) + wkr*(a[ii+j][1]-b[ii+j][1]);
             a[ii+j][0] = a[ii+j][0]+b[ii+j][0];
             a[ii+j][1] = a[ii+j][1]+b[ii+j][1];
             b[ii+j][0] = c[ii+j][0]+d[ii+j][0];
             b[ii+j][1] = c[ii+j][1]+d[ii+j][1];
             wcmdr = wkr*(c[ii+j][0]-d[ii+j][0]) - wku*(c[ii+j][1]-d[ii+j][1]);
             wcmdu = wku*(c[ii+j][0]-d[ii+j][0]) + wkr*(c[ii+j][1]-d[ii+j][1]);
             c[ii+j][0] = wambr;
             c[ii+j][1] = wambu;
             d[ii+j][0] = wcmdr;
             d[ii+j][1] = wcmdu;
         } 
      } 
   } 
}
void rcfft2(int n, float *x, int sign)
/*
 
   real -> complex FFT: x = real input array (n dimensional)
 
   INPUT
         n = power of 2 dimension of transform
         x = n/2 dimensional real array
 
   OUTPUT (where w = n-th root of unity)
 
         x_j <- sum(j=1,n) w^{jk} x_k, j=0..n/2
             = "half complex" output array, x_{n-j} = x_j
               hence, only n/2+1 unique elements.
 
         Since input is real, output x_{0} and x_{n/2} are
         purely real. To be done in-place, then, the data 
         structure of output has a packed first complex 
         element output: 

            x[0][0] = real(x_{0}),  and 
            x[0][1] = real(x_{n/2}).

         Otherwise, x[i][0] = real(x_{i}), and x[i][1] =
         imag(x_{i}), for i=1..n/2-1.

   W Petersen, IPS, ETH Zurich, 7 March, 1995
*/
{
   void cffti(),cfft_2();
   void post_p(int,float*,float*,int);
   static float *w,*pw1,*pw2;
   static int n1=-33,n2=-33;
   int nt;
   if(n==n1){
      w = pw1;
   } else {
      if(n==n2){
         w   = pw2;
         pw2 = pw1;
         pw1 = w;
         nt  = n1;
         n1  = n2;
         n2  = nt;
      } else {
         if(n2>0) free(pw2);
         w   = (void *) malloc(n*sizeof(float));
         cffti(n,w);
         pw2 = pw1;
         pw1 = w;
         n2  = n1;
         n1  = n;
      }
   }
   cfft_2(n/2,x,w,2,sign);
   post_p(n,x,w,sign);
}
void crfft2(int n, float *x, int sign)
/* 
 
   "Half complex" -> real FFT: x = Hermitian complex input
   assumed, such that for all j=1,n, x_{n-j} = conjg(x_j), 
   hence only n/2+1 unique elements. Furthermore, first and
   last are purely real. In-place computation has a packed 
   data structure as below.
 
   INPUT
         n = power of 2 dimension of transform

         x = n/2 dimensional COMPLEX array, where
             the first elements 

                 x[0][0] = real(x_{0}) and, 
                 x[0][1] = real(x_{n/2}).
 
             Otherwise, for i=1..n/2-1, 

                  x[i][0] = real(x_{i}) and
                  x[i][1] = imag(x_{i}), 

             treated as a complex array. 
 
   OUTPUT (where w = n-th root of unity) x[j], j=0..n-1 get
 
         x_j <- sum(k=1,n) w^{jk} x_k, j=0..n-1, are REAL.
 
   W Petersen, SAM, Mathematik, ETHZ, 25 Aug, 1999

*/
{
   void cfft_2(),pre_proc();
   static float *w,*pw1,*pw2;
   static int n1=-33,n2=-33;
   int nt;
   if(n==n1){
      w = pw1;
   } else {
      if(n==n2){
         w   = pw2;
         pw2 = pw1;
         pw1 = w;
         nt  = n1;
         n1  = n2;
         n2  = nt;
      } else {
         if(n2>0) free(pw2);
         w   = (void *) malloc(n*sizeof(float));
         cffti(n,w);
         pw2 = pw1;
         pw1 = w;
         n2  = n1;
         n1  = n;
      }
   }
   pre_proc(n,x,w,sign);
   cfft_2(n/2,x,w,2,sign);
}
void post_p(int n, float *z, float *w, int sign)
{
   float pr,pi,qr,qi,zr,zi;
   float ar,ai,br,bi,wr,wi;
   int j;

   zr     = *z + *(z+1);
   zi     = *z - *(z+1);
   *z     = zr; 
   *(z+1) = zi;
   if(sign<0){
      for(j=2;j<n/2-1;j+=2){
         pr    = 0.5*(*(z+j) + *(z+n-j));
         pi    = 0.5*(*(z+j+1) - *(z+n+1-j));
         wr    = *(w+j);
         wi    = - (*(w+j+1));
         ar    = *(z+j);
         ai    = *(z+j+1);
         br    = *(z+n-j);
         bi    = *(z+n+1-j);
         qr    = 0.5*(wr*(ar-br) - wi*(ai+bi));
         qi    = 0.5*(wi*(ar-br) + wr*(ai+bi));
         *(z+j)     = pr+qi;
         *(z+j+1)   = pi-qr;
         *(z+n-j)   = pr-qi;
         *(z+n+1-j) = -pi-qr;
      }
      *(z+n/2+1) = -(*(z+n/2+1));
   } else {
      for(j=2;j<n/2-1;j+=2){
         pr    = 0.5*(*(z+j) + *(z+n-j));
         pi    = 0.5*(*(z+j+1) - *(z+n+1-j));
         wr    = *(w+j);
         wi    = *(w+j+1);
         ar    = *(z+j);
         ai    = *(z+j+1);
         br    = *(z+n-j);
         bi    = *(z+n+1-j);
         qr    = 0.5*(wr*(ar-br) - wi*(ai+bi));
         qi    = 0.5*(wi*(ar-br) + wr*(ai+bi));
         *(z+j)     = pr+qi;
         *(z+j+1)   = pi-qr;
         *(z+n-j)   = pr-qi;
         *(z+n+1-j) = -pi-qr;
      }
   }
} 
void pre_proc(int n, float *z, float *w, int sign)
{
   float pr,pi,qr,qi,zr,zi;
   float ar,ai,br,bi,wr,wi;
   int j;

   zr     = *z + *(z+1);
   zi     = *z - *(z+1);
   *z     = zr; 
   *(z+1) = zi;
   if(sign<0){
      for(j=2;j<n/2-1;j+=2){
         pr    = *(z+j) + *(z+n-j);
         pi    = *(z+j+1) - *(z+n+1-j);
         wr    = *(w+j);
         wi    = - (*(w+j+1));
         ar    = *(z+j);
         ai    = *(z+j+1);
         br    = *(z+n-j);
         bi    = *(z+n+1-j);
         qr    = wr*(ar-br) - wi*(ai+bi);
         qi    = wi*(ar-br) + wr*(ai+bi);
         *(z+j)     = pr-qi;
         *(z+j+1)   = pi+qr;
         *(z+n-j)   = pr+qi;
         *(z+n+1-j) = -pi+qr;
      }
      zr         = *(z+n/2);
      zi         = *(z+n/2+1);
      *(z+n/2)   = 2.0*zr;
      *(z+n/2+1) = 2.0*zi;
   } else {
      for(j=2;j<n/2-1;j+=2){
         pr    = *(z+j) + *(z+n-j);
         pi    = *(z+j+1) - *(z+n+1-j);
         wr    = *(w+j);
         wi    = *(w+j+1);
         ar    = *(z+j);
         ai    = *(z+j+1);
         br    = *(z+n-j);
         bi    = *(z+n+1-j);
         qr    = wr*(ar-br) - wi*(ai+bi);
         qi    = wi*(ar-br) + wr*(ai+bi);
         *(z+j)     = pr-qi;
         *(z+j+1)   = pi+qr;
         *(z+n-j)   = pr+qi;
         *(z+n+1-j) = -pi+qr;
      }
      zr         = *(z+n/2);
      zi         = *(z+n/2+1);
      *(z+n/2)   = 2.0*zr;
      *(z+n/2+1) = -2.0*zi;
   }
} 
void refft2(int n,float x[])
/* 
   real/even FFT: x = real input array (n/2+1 dimensional)
   such that for all j=0,n-1, x_{n-j} = x_j, hence only n2/+1
   unique elements. This is a cosine transform. See:
   PN Swarztrauber, Math. Comp., 47, 174, pp 323-346, 1986.
   It uses Dollimore's method.
 
   INPUT
         n = power of 2 dimension of transform
         x = n/2+1 dimensional real array
 
   OUTPUT (where w = n-th root of unity)
 
         x_j <- sum(k=0,n-1) w^{jk} x_k, j=0..n/2
 
   W Petersen, SAM, Mathematik, ETHZ, 25 Aug, 1999
*/
{
    void cffti(int,float*),rcfft2(int,float*,int);
    static float *w,*pw1,*pw2;
    float x1,x2,xm,xp;
    static int n1=-33,n2=-33;
    int i,ii,nt;
    if(n==n1){
       w = pw1;
    } else {
       if(n==n2){
          w   = pw2;
          pw2 = pw1;
          pw1 = w;
          nt  = n1;
          n1  = n2;
          n2  = nt;
       } else {
          if(n2>0) free(pw2);
          w   = (void *) malloc(n*sizeof(float));
          cffti(n,w);
          pw2 = pw1;
          pw1 = w;
          n2  = n1;
          n1  = n;
       }
    }
    x1      = x[0] + x[n/2];
    x2      = x[0] - x[n/2];
    x[n/4]  = 2.0*x[n/4];
    x[0]    = x1;
    for(i=1;i<n/4;i++){
       ii       = 2*i;
       xp       = x[i] + x[n/2-i]; 
       xm       = x[i] - x[n/2-i];
       x2      += 2.0*(*(w+ii))*xm;
       x[i]     = xp + (*(w+ii+1))*xm; 
       x[n/2-i] = xp - (*(w+ii+1))*xm;
    }

    rcfft2(n/2,x,1);

    x1     = x[1];
    x[1]   = x2;
    x[n/2] = x1;
    for(i=3; i<n/2; i += 2){
       x[i]  = x[i-2] - 2.0*x[i];
    }
}
void rofft2(int n,float x[])
/* 
   real/odd FFT: x = real input array (n/2 dimensional)
   such that for all j=0,n-1, x_{n-j} = - x_j, hence only n2-1
   unique elements. This is a sine transform. See:
   PN Swarztrauber, Math. Comp., 47, 174, pp 323-346, 1986
 
   INPUT
         n = power of 2 dimension of transform
         x = n/2+1 dimensional real array
             Important: x(1) = 0
 
   OUTPUT (where w = n-th root of unity), 

         x[j], j=0..n/2-1 get
 
             x_j <- sum(k=0,n-1) w^{jk} x_k, j=0..n/2-1
                
         and x[0] = x_{0} = 0.
 
   W Petersen, SAM, Mathematik, ETHZ, 25 Aug, 1999
*/
{
    void cffti(), rcfft2();
    static float *w,*pw1,*pw2;
    static int n1=-33,n2=-33;
    float x1,xm,xp;
    int i,ii,nt;
    if(n==n1){
       w = pw1;
    } else {
       if(n==n2){
          w   = pw2;
          pw2 = pw1;
          pw1 = w;
          nt  = n1;
          n1  = n2;
          n2  = nt;
       } else {
          if(n2>0) free(pw2);
          w   = (void *) malloc(n*sizeof(float));
          cffti(n,w);
          pw2 = pw1;
          pw1 = w;
          n2  = n1;
          n1  = n;
       }
    }
    for(i=1;i<n/4;i++){
       ii       = 2*i;
       xp       = x[i] + x[n/2-i];
       xm       = x[i] - x[n/2-i];
       x[i]     =  xm + (*(w+ii+1))*xp;
       x[n/2-i] = -xm + (*(w+ii+1))*xp;
    }
    x[n/4] *= 2.0;

    rcfft2(n/2,x,1);

    for(i=0; i<n/2; i += 2){
       x1     = x[i];
       x[i]   = x[i+1];
       x[i+1] = x1;
    }
    x[0] = 0;
    for(i=3; i<n/2; i += 2){
       x[i]   = 2.0*x[i] + x[i-2];
    }
}
# include <sys/time.h>
double walltime(double *t0)
{
   double mic, time;
   double mega=0.000001;
   struct timeval tp;
   struct timezone tzp;
   static long base_sec = 0;
   static long base_usec = 0;

   (void) gettimeofday(&tp, &tzp);
   if (base_sec == 0) {
     base_sec  = tp.tv_sec;
     base_usec = tp.tv_usec;
   }
   time = (double)(tp.tv_sec - base_sec);
   mic = (double)(tp.tv_usec - base_usec);
   time = (time + mic * mega) - *t0;
   return(time);
}