#include <stdio.h>
#include <math.h>
#include <time.h>
#include <float.h>
#define N 1048576
#define N2 N/2
main()
{
/* 
   Example of Apple Altivec coded binary radix FFT
   using intrinsics from Petersen and Arbenz "Intro. 
   to Parallel Computing," Section 3.6
 
   This is an expanded version of a generic work-space 
   FFT: steps are in-line. cfft2(n,x,y,w,sign) takes complex
   n-array "x" (Fortran real,aimag,real,aimag,... order) 
   and writes its DFT in "y". Both input "x" and the 
   original contents of "y" are destroyed. Initialization
   for array "w" (size n/2 complex of twiddle factors
   (exp(twopi*i*k/n), for k=0..n/2-1)) is computed once
   by cffti(n,w).

                      WPP, SAM. Math. ETHZ, 1 June, 2002 
*/

   int first,i,icase,it,ln2,n;
   int nits=1000000;
   static float seed = 331.0;
   float error,fnm1,sign,z0,z1,ggl();
   float *x,*y,*z,*w;
   double t1,mflops;
   void cffti(),cfft2();
/* allocate storage for x,y,z,w on 4-word bndr. */
   x = (float *) malloc(8*N);
   y = (float *) malloc(8*N);
   z = (float *) malloc(8*N);
   w = (float *) malloc(4*N);
   n     = 2;
   for(ln2=1;ln2<21;ln2++){
      first = 1;
      for(icase=0;icase<2;icase++){
         if(first){
            for(i=0;i<2*n;i+=2){
               z0 = ggl(&seed);     /* real part of array */
               z1 = ggl(&seed);     /* imaginary part of array */
               x[i] = z0;
               z[i] = z0;           /* copy of initial real data */
               x[i+1] = z1;
               z[i+1] = z1;         /* copy of initial imag. data */
            }
         } else {
            for(i=0;i<2*n;i+=2){
               z0 = 0;              /* real part of array */
               z1 = 0;              /* imaginary part of array */
               x[i] = z0;
               z[i] = z0;           /* copy of initial real data */
               x[i+1] = z1;
               z[i+1] = z1;         /* copy of initial imag. data */
            }
         }
/* initialize sine/cosine tables */
         cffti(n,w);
/* transform forward, back */
         if(first){
            sign = 1.0;
            cfft2(n,x,y,w,sign);
            sign = -1.0;
            cfft2(n,y,x,w,sign);
/* results should be same as initial multiplied by n */
            fnm1 = 1.0/((float) n);
            error = 0.0;
            for(i=0;i<2*n;i+=2){
               error += (z[i] - fnm1*x[i])*(z[i] - fnm1*x[i]) +
                   (z[i+1] - fnm1*x[i+1])*(z[i+1] - fnm1*x[i+1]);
            }
            error = sqrt(fnm1*error);
            printf(" for n=%d, fwd/bck error=%e\n",n,error);
            first = 0;
         } else {
            t1   = ((double)clock())/((double) CLOCKS_PER_SEC);
            for(it=0;it<nits;it++){
               sign = +1.0;
               cfft2(n,x,y,w,sign);
               sign = -1.0;
               cfft2(n,y,x,w,sign);
            }
            t1   = ((double)clock())/((double) CLOCKS_PER_SEC) - t1;
            t1   = 0.5*t1/((double) nits);
            mflops = 5.0*((double) n)*((double) ln2)/((1.e+6)*t1);
            printf(" for n=%d, t1=%e, mflops=%e\n",n,t1,mflops);
         }
      }
      if((ln2%4)==0) nits /= 10;
      n *= 2;
   }
}
void cfft2(unsigned int n,float x[][2],float y[][2],float w[][2], float sign)
{

/* 
   altivec version of cfft2 from Petersen and Arbenz book, "Intro.
   to Parallel Computing", Oxford Univ. Press, 2003, Section 3.6
                                        wpp 14. Dec. 2003
*/

   int jb,jc,jd,jw,k,k2,k4,lj,m,j,mj,mj2,pass,tgle;
   float rp,up,wr[4],wu[4];
   float *a,*b,*c,*d;
   const vector float vminus = (vector float)(-0.,0.,-0.,0.);
   const vector float vzero  = (vector float)(0.,0.,0.,0.);
   const vector unsigned char pv3201 =
   (vector unsigned char)(4,5,6,7,0,1,2,3,12,13,14,15,8,9,10,11);
   vector float V0,V1,V2,V3,V4,V5,V6,V7;
   vector float V8,V9,V10,V11,V12,V13,V14,V15;

   if(n<=1){
      y[0][0] = x[0][0];
      y[0][1] = x[0][1];
      return;
   }
   m    = (int) (log((float) n)/log(1.99));
   mj   = 1;
   mj2  = 2;
   lj   = n/2;
/* first pass thru data: x -> y */
   for(j=0;j<lj;j++){
      jb = n/2+j; jc  = j*mj2; jd = jc + 1;
      rp = w[j][0]; up = w[j][1];
      if(sign<0.0) up = -up;
      y[jd][0] = rp*(x[j][0] - x[jb][0]) - up*(x[j][1] - x[jb][1]);
      y[jd][1] = up*(x[j][0] - x[jb][0]) + rp*(x[j][1] - x[jb][1]);
      y[jc][0] = x[j][0] + x[jb][0];
      y[jc][1] = x[j][1] + x[jb][1];
   }
   if(n==2) return;
/* next pass is mj = 2 */
   mj  = 2;
   mj2 = 4;
   lj  = n/4;
   a = (float *)&y[0][0];
   b = (float *)&y[n/2][0];
   c = (float *)&x[0][0];
   d = (float *)&x[mj][0];
   if(n==4){
      c = (float *)&y[0][0];
      d = (float *)&y[mj][0];
   }
   for(j=0;j<lj;j++){
      jw = j*mj; jc = j*mj2; jd = 2*jc;
      rp = w[jw][0]; up = w[jw][1];
      if(sign<0.0) up = -up;
      wr[0] = rp; wr[1] = rp; wr[2] = rp; wr[3] = rp;
      wu[0] = up; wu[1] = up; wu[2] = up; wu[3] = up;
      V6 = vec_ld(0,wr);
      V7 = vec_ld(0,wu);
      V7 = vec_xor(V7,vminus);
      V0 = vec_ld(0,(vector float *) (a+jc));
      V1 = vec_ld(0,(vector float *) (b+jc));
      V2 = vec_add(V0,V1);                         /* a + b */
      vec_st(V2,0,(vector float *) (c+jd));     /* store c */
      V3 = vec_sub(V0,V1);                         /* a - b */
      V4 = vec_perm(V3,V3,pv3201);
      V0 = vec_madd(V6,V3,vzero);
      V1 = vec_madd(V7,V4,vzero);
      V2 = vec_add(V0,V1);                         /* w*(a - b) */
      vec_st(V2,0,(vector float*) (d+jd));         /* store d */
   }
   if(n==4) return;
   mj  *= 2;
   mj2  = 2*mj;
   lj   = n/mj2; 
   tgle = 0;  
   for(pass=2;pass<m-1;pass++){
      if(tgle){
         a = (float *)&y[0][0];
         b = (float *)&y[n/2][0];
         c = (float *)&x[0][0];
         d = (float *)&x[mj][0];
 	 tgle = 0;
      } else {
         a = (float *)&x[0][0];
         b = (float *)&x[n/2][0];
         c = (float *)&y[0][0];
         d = (float *)&y[mj][0];
	 tgle = 1;
      }
      for(j=0; j<lj; j++){
         jw = j*mj; jc = j*mj2; jd = 2*jc;
         rp = w[jw][0];
         up = w[jw][1];
         if(sign<0.0) up = -up;
         wr[0] = rp; wr[1] = rp; wr[2] = rp; wr[3] = rp;
         wu[0] = up; wu[1] = up; wu[2] = up; wu[3] = up;
         V6 = vec_ld(0,wr);
         V7 = vec_ld(0,wu);
         V7 = vec_xor(V7,vminus);
         for(k=0; k<mj; k+=4){
            k2   = 2*k; k4 = k2+4;
            V0 = vec_ld(0,(vector float *) (a+jc+k2));
            V1 = vec_ld(0,(vector float *) (b+jc+k2));
            V2 = vec_add(V0,V1);                        /* a + b */
            vec_st(V2,0,(vector float*) (c+jd+k2));   /* store c */
            V3 = vec_sub(V0,V1);                        /* a - b */
            V4 = vec_perm(V3,V3,pv3201);
            V0 = vec_madd(V6,V3,vzero);
            V1 = vec_madd(V7,V4,vzero);
            V2 = vec_add(V0,V1);                        /* w*(a - b) */
            vec_st(V2,0,(vector float *) (d+jd+k2));    /* store d */
	    V8 = vec_ld(0,(vector float *) (a+jc+k4));
	    V9 = vec_ld(0,(vector float *) (b+jc+k4));
	    V10 = vec_add(V8,V9);		        /* a + b */
	    vec_st(V10,0,(vector float *) (c+jd+k4));   /* store c */
	    V11 = vec_sub(V8,V9);		        /* a - b */
	    V12 = vec_perm(V11,V11,pv3201);
	    V8  = vec_madd(V6,V11,vzero);
	    V9  = vec_madd(V7,V12,vzero);
	    V10 = vec_add(V8,V9);		        /* w*(a - b) */
	    vec_st(V10,0,(vector float *) (d+jd+k4));   /* store d */
         }
      }
      mj  *= 2;
      mj2  = 2*mj;
      lj   = n/mj2;
   }
/* last pass thru data: in-place if previous in y */
   c = (float *)&y[0][0];
   d = (float *)&y[n/2][0];
   if(tgle) {
      a = (float *)&y[0][0];
      b = (float *)&y[n/2][0];
   } else {
      a = (float *)&x[0][0];
      b = (float *)&x[n/2][0];
   }
   for(k=0; k<(n/2); k+=4){
      k2 = 2*k; k4 = k2+4;
      V0 = vec_ld(0,(vector float *) (a+k2));
      V1 = vec_ld(0,(vector float *) (b+k2));
      V2 = vec_add(V0,V1);                      /* a + b */
      vec_st(V2,0,(vector float*) (c+k2));      /* store c */
      V3 = vec_sub(V0,V1);                      /* a - b */
      vec_st(V3,0,(vector float *) (d+k2));     /* store d */
      V4 = vec_ld(0,(vector float *) (a+k4));
      V5 = vec_ld(0,(vector float *) (b+k4));
      V6 = vec_add(V4,V5);                      /* a + b */
      vec_st(V6,0,(vector float *) (c+k4));     /* store c */
      V7 = vec_sub(V4,V5);                      /* a - b */
      vec_st(V7,0,(vector float *) (d+k4));     /* store d */
   }
}
void cffti(int n, float w[][2])
{

/* initialization routine for cfft2: computes
        cos(twopi*k),sin(twopi*k) for k=0..n/2-1
   - the "twiddle factors" for a binary radix FFT */

   int i,n2;
   float aw,arg,pi;
   pi = 3.141592653589793;
   n2 = n/2;
   aw = 2.0*pi/((float)n);
   for(i=0;i<n2;i++){
      arg   = aw*((float)i);
      w[i][0] = cos(arg);
      w[i][1] = sin(arg);
   }
}
#include <math.h>
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)));
}