#include <stdio.h>
#include <math.h>
#include <time.h>
/* timing test program for normal RNG generators */
#define N 2000000
static unsigned long jz,jsr=123456789;
#define SHR3 (jz=jsr, jsr^=(jsr<<13), jsr^=(jsr>>17), jsr^=(jsr<<5),jz+jsr)
#define UNI (.5 + (signed) SHR3 * .2328306e-9)
#define IUNI SHR3
static long hz;
static int seed=331;
static unsigned long iz, kn[128], ke[256];
static float wn[128],fn[128], we[256],fe[256];
#define RNOR (hz=SHR3, iz=hz&127, (abs(hz)<kn[iz])? hz*wn[iz] : nfix())
#define REXP (jz=SHR3, iz=jz&255, (    jz <ke[iz])? jz*we[iz] : efix())
/* nfix() generates variates from the residue when rejection in RNOR occurs. */
/* The starting of the right tail */
static float x, y;
float nfix(void) 
{   
/* The starting of the right tail */
   const float r = 3.442620f; 
   static float x, y;
   float ggl();
   for(;;){
      x=hz*wn[iz];
      if(iz==0){   /* iz==0, handle the base strip */
         do{   
            x=-log(ggl(&seed))*0.2904764;   /* .2904764 is 1/r */  
            y=-log(ggl(&seed)); 
         } while(y+y<x*x);
         return (hz>0)? r+x : -r-x;      
      }
      /* iz>0, handle the wedges of other strips */ 
      if( fn[iz]+ggl(&seed)*(fn[iz-1]-fn[iz]) < exp(-.5*x*x) ) return x;
      /* start all over */ 
      hz=SHR3;
      iz=hz&127;
      if(abs(hz)<kn[iz]) return (hz*wn[iz]);   
   }
}
/* efix() generates variates from the residue when rejection in REXP occurs. */
float efix(void)
{
   float x,ggl();
   for(;;){
      /* iz==0 */
      if(iz==0) return (7.69711-log(ggl(&seed)));
      x=jz*we[iz];
      if( fe[iz]+ggl(&seed)*(fe[iz-1]-fe[iz]) < exp(-x) ) return (x);
      /* Start all over again */
      jz=SHR3;
      iz=(jz&255);
      if(jz<ke[iz]) return (jz*we[iz]);   
   }
}
/*--------This procedure sets the seed and creates the tables------*/
void zigset(unsigned long jsrseed) 
{
   const double m1 = 2147483648.0, m2 = 4294967296.;
   double dn=3.442619855899,tn=dn,vn=9.91256303526217e-3, q;
   double de=7.697117470131487, te=de, ve=3.949659822581572e-3;
   int i;
   jsr^=jsrseed;
/* Set up tables for RNOR */
   q     = vn/exp(-.5*dn*dn);
   kn[0]   = (dn/q)*m1;
   kn[1]   = 0;
   wn[0]   = q/m1;
   wn[127] = dn/m1;
   fn[0]   = 1.;
   fn[127] = exp(-.5*dn*dn);      
   for(i=126;i>=1;i--){
       dn      = sqrt(-2.*log(vn/dn+exp(-.5*dn*dn)));
       kn[i+1] = (dn/tn)*m1;
       tn      = dn;
       fn[i]   = exp(-.5*dn*dn);
       wn[i]   = dn/m1;
   }
/* Set up tables for REXP */
   q       = ve/exp(-de);
   ke[0]   = (de/q) * m2;
   ke[1]   = 0;
   we[0]   = q/m2;
   we[255] = de/m2;
   fe[0]   = 1.;
   fe[255] = exp(-de);      
   for(i=254;i>=1;i--){
       de      = -log(ve/de+exp(-de));
       ke[i+1] = (de/te)*m2;
       te      = de;
       fe[i]   = exp(-de);
       we[i]   = de/m2;
   }
}
float xv[N],yv[N],zv[N];
main()
{
/* 
   compares performance of Box-Muller, Ziggurat
   vs. the Polar methods for generating normals 
*/
   float r1,r2,s,x1,x2,y1,y2,z1,z2,z3,z4,twopi;
   float t1,t2,u1,u2,ts,te;
   float ggl();
   void zigset();
   static float seed=331.0;
   int i,iter,n,nv;
   twopi = 6.2831853071795862;
   n     = 10000;
   for(nv=0;nv<4;nv++){
/* Box-Muller method */
   ts = ((float)clock())/((float) CLOCKS_PER_SEC);
   for(iter=0;iter<10;iter++){
#pragma ivdep
   for(i=0;i<n;i+=2){
      u1 = ggl(&seed);
      u2 = ggl(&seed);
      r1 = twopi*u1;
      t1 = sin(r1);
      t2 = cos(r1);
      r2 = sqrt(-2.0*log(u2));
      xv[i] = t1*r2;
      xv[i+1] = t2*r2;
   }
   }
   te = ((float)clock())/((float) CLOCKS_PER_SEC);
   te = te - ts;
   te *= 0.1;
   printf(" Box-Muller time for n=%d is %e\n",n,te);
   x1 = 0.0;
   x2 = 0.0;
   for(i=0;i<n;i++){
      x1 += xv[i];
      x2 += xv[i]*xv[i];
   }
   x1 /= (float)n;
   x2 /= (float)n;
   x2 -= x1*x1;
   printf("    - mean and variance of B-M method: %e, %e\n",x1,x2);
/* polar method */
   ts = ((float)clock())/((float) CLOCKS_PER_SEC);
   for(iter=0;iter<10;iter++){
   for(i=0;i<n;i+=2){
new:
      u1 = ggl(&seed);
      u2 = ggl(&seed);
      u1 = 1.0 - 2.0*u1;
      u2 = 1.0 - 2.0*u2;
      s  = u1*u1 + u2*u2;
      if(s < 1.0){
          yv[i  ] = u1*sqrt(-2.0*log(s)/s);
          yv[i+1] = u2*sqrt(-2.0*log(s)/s);
      } else {
          goto new;
      }
   }
   }
   te = ((float)clock())/((float) CLOCKS_PER_SEC);
   te = te - ts;
   te *= 0.1;
   y1 = 0.0;
   y2 = 0.0;
   for(i=0;i<n;i++){
      y1 += yv[i];
      y2 += yv[i]*yv[i];
   }
   y1 /= (float)n;
   y2 /= (float)n;
   y2 -= y1*y1;
   printf(" Polar method time for n=%d is %e\n",n,te);
   printf("    - mean and variance of polar method: %e, %e\n",y1,y2);
/* ziggurat method */
   zigset(331);
   ts = ((float)clock())/((float) CLOCKS_PER_SEC);
   for(iter=0;iter<10;iter++){
   for(i=0;i<n;i+=2){
       zv[i  ] = RNOR;
       zv[i+1] = RNOR;
   }
   }
   te = ((float)clock())/((float) CLOCKS_PER_SEC);
   te = te - ts;
   te *= 0.1;
   z1 = 0.0;
   z2 = 0.0;
   z3 = 0.0;
   z4 = 0.0;
   for(i=0;i<n;i++){
      z1 += zv[i];
      z2 += zv[i]*zv[i];
      z3 += zv[i]*zv[i]*zv[i];
      z4 += zv[i]*zv[i]*zv[i]*zv[i];
   }
   z1 /= (float)n;
   z2 /= (float)n;
   z3 /= (float)n;
   z4 /= (float)n;
   z4 += -4.0*z3*z1 + 6.0*z2*z1*z1 - 3.0*z1*z1*z1*z1;
   z3 += -3.0*z2*z1 + 2.0*z1*z1*z1;
   z2 -= z1*z1;
   printf(" ziggurat method time for n=%d is %e\n",n,te);
   printf("    - mean and variance of ziggurat method: %e, %e\n",z1,z2);
   printf("    - skewness and kurtosis of ziggurat method: %e, %e\n",z3,z4);
   n *= 5;
   }
}
#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 */

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