#include <stdio.h>
#include <math.h>
#define N 1048576
int main()
{
/* 
   Test program for two-term lagged Fibonacci RNG

        x_k = x_{k-p} + z_{k-q} mod 1.0

   levels=0..9 refer to which lags p,q in ptab, qtab
   are used. The longer the lags, the more memory 
   required for the buffer (size p, where p > q). 
   For information printing, these tables are repeated 
   in this main routine.

   Extract the library for zufall beginning at
        "cut here for library"
   where the code above that cutting point has only
   testing routines. test "nphere" generates uniformly
   distributed points either in/on an n-sphere. Only
   n=2 and n=3 spheres are tested, but the routine is
   designed for large n. save/restore sequences 
   are also tested. This is an extended and much
   modified version of zufall distributed on NETLIB.

                                     wpp 9/1/2004                  
*/
   float xstep,xs,x[N];
   float *tant,*vek,*vptr;
   float f2m1,eps;
   double fnm1,xx2,xx4,x1,x2,x3,x4,x5,x6;
   double d1,d2,d3,d4,d5,d6;
   int i,ii,lev,m,ibin[21];
   int tw2m;
   int iprime[]={331,337,347,349,353,359,367,373,379};
/* these are copies of library globals for printing */
   int pp[]={607,1279,2281,3217,4423,9689,19937,23209,44497};
   int qp[]={273, 418,1252, 576,2098,5502, 9842,13470,21034};
/*  */
   char *bewahr[]={"save331","save337","save347",
                   "save349","save353","save359",
                   "save367","save373","save379"};
   void zufall_init(),zufall_spar(),zufall_rest();
   float zufall(),zuf_g(),*nsphere();

   xstep = 21.0;
   for(lev=0;lev<9;lev++){
      printf("\n ===========================================\n");
      printf(" Begin lev=%d tests: lags=%d, %d, seed=%d\n",
               lev,pp[lev],qp[lev],iprime[lev]);
      printf(" ===========================================\n");
      zufall_init(lev,iprime[lev]);
      for(i=0;i<N/2;i++){
         x[i] = zufall();
      }
      zufall_spar(bewahr[lev]);
      for(i=N/2;i<N;i++){
         x[i] = zufall();
      }
      zufall_rest(bewahr[lev]);
      xs = zufall();
      if(fabsf(xs-x[N/2])>1.e-6){
         printf("\n save/restore error: xs=%e, x[%d]=%e\n",xs,N/2,x[N/2]);
      } else {
         printf("\n  **** save/restore test passed ****\n");
      }
/* do histogram */
      for(i=0;i<21;i++) ibin[i] = 0;
      for(i=0;i<N;i++){
         ii = xstep*x[i];
         ibin[ii]++;
      }
      printf("\n |---- Uniform test: ----| \n\n");
      printf("  Uniform histogram: compare [%e, %e]\n",
             ((float)N)/21.0-sqrt(((float)N)/21.0),
             ((float)N)/21.0+sqrt(((float)N)/21.0));
      printf("    ibin=%d\n",ibin[0]);
      for(i=1;i<21;i++) printf("         %d\n",ibin[i]);
/* now do Gaussian tests */
      printf("\n |---- Normal test: ----| \n\n");
      for(i=0;i<N;i++){
         x[i] = zuf_g();
      }
      for(i=0;i<21;i++) ibin[i] = 0;
      for(i=0;i<N;i++){
         ii = (int)(2.0*(x[i]+5.25));
         ibin[ii]++;
      }
/* compute moments */
      x1 = 0.0; x2 = 0.0; x3 = 0.0;
      x4 = 0.0; x5 = 0.0; x6 = 0.0;
      for(i=0;i<N;i++){
         x1 += (double) x[i];
         xx2 = (double) (x[i]*x[i]);
         xx4 = xx2*xx2;
         x2 += xx2;
         x3 += ((double) x[i])*xx2;
         x4 += xx4;
         x5 += ((double) x[i])*xx4;
         x6 += xx2*xx4;
      }
      fnm1 = 1.0/((double)N);
      x1 *= fnm1; x2 *= fnm1;
      x3 *= fnm1; x4 *= fnm1;
      x5 *= fnm1; x6 *= fnm1;
/* these are the theoretical deviations of these moments */
      d1 = sqrt(fnm1); 
      d2 = sqrt(2.0*fnm1); 
      d3 = sqrt(15.0*fnm1); 
      d4 = sqrt(102.0*fnm1); 
      d5 = sqrt(945.0*fnm1); 
      d6 = sqrt(10380.0*fnm1); 
      printf("  Moments:\n",
          "Compare to: (0.0 +- %7.2e)      (1.0 +- %7.2e)\n",
             "<x>    = % e  <x**2> = %e \n",
          "Compare to: (0.0 +- %7.2e)      (3.0 +- %7.2e)\n",
             "<x**3> = % e  <x**4> = %e \n"
          "Compare to: (0.0 +- %7.2e)      (15.0 +- %7.2e)\n",
             "<x**5> = % e  <x**6> = %e \n",
              d1,d2,x1,x2,d3,d4,x3,x4,d5,d6,x5,x6);
      printf("\n  Normal distr. histogram: \n");
      printf("    ibin=%d\n",ibin[0]);
      for(i=1;i<21;i++) printf("         %d\n",ibin[i]);
      printf("\n");
/* n-sphere tests */
      printf("\n |---- N-sphere test: ----| \n\n");
/* 2-D test */
      m    = 2;
      tw2m = 4;
      f2m1 = 1.0/((float)N);
      vek  = (float *) malloc(N*m*sizeof(float));
      tant = (float *) malloc(tw2m*sizeof(float));
      for(i=0;i<N;i++){
         vptr = nsphere(m,1);
         *(vek+i*m)   = *vptr;
         *(vek+i*m+1) = *(vptr+1);
      }
      tant[0] = 0.0; tant[1] = 0.0;
      tant[2] = 0.0; tant[3] = 0.0;
      for(i=0;i<N;i++){
         if(((*(vek+i*m))>=0.0)&&((*(vek+i*m+1))>=0.0)) tant[0] +=1.0;
         if(((*(vek+i*m))<0.0)&&((*(vek+i*m+1))>=0.0))  tant[1] +=1.0;
         if(((*(vek+i*m))>=0.0)&&((*(vek+i*m+1))<0.0))  tant[2] +=1.0;
         if(((*(vek+i*m))<0.0)&&((*(vek+i*m+1))<0.0))   tant[3] +=1.0;
      }
      printf("  2-D test: each quadrant should get 1/4\n");
      eps = sqrt(4.0*f2m1);
      i   = 0;
      if(fabsf(f2m1*tant[0]-0.25)>eps){
         printf("  1st quadrant frac=%e\n",f2m1*tant[0]);
         i += 1;
      }
      if(fabsf(f2m1*tant[1]-0.25)>eps){
         printf("  2nd quadrant frac=%e\n",f2m1*tant[1]);
         i += 1;
      }
      if(fabsf(f2m1*tant[2]-0.25)>eps){
         printf("  3rd quadrant frac=%e\n",f2m1*tant[2]);
         i += 1;
      }
      if(fabsf(f2m1*tant[3]-0.25)>eps){
         printf("  4th quadrant frac=%e\n",f2m1*tant[3]);
         i += 1;
      }
      if(i>0){
         printf("  %d |quadrant fraction(s) - 1/4| >< sqrt(4/N)\n",i);
      } else {
         printf("  **** 2-D test passed **** \n");
      }
      free(vek); free(tant);
/* 3-D test */
      m    = 3;
      tw2m = 8;
      f2m1 = 1.0/((float)N);
      vek  = (float *) malloc(N*m*sizeof(float));
      tant = (float *) malloc(tw2m*sizeof(float));
      for(i=0;i<N;i++){
         vptr = nsphere(m,1);
         *(vek+i*m)   = *vptr;
         *(vek+i*m+1) = *(vptr+1);
         *(vek+i*m+2) = *(vptr+2);
      }
      tant[0] = 0.0; tant[1] = 0.0;
      tant[2] = 0.0; tant[3] = 0.0;
      tant[4] = 0.0; tant[5] = 0.0;
      tant[6] = 0.0; tant[7] = 0.0;
      for(i=0;i<N;i++){
         if(((*(vek+i*m))>=0.0)&&((*(vek+i*m+1))>=0.0)&&((*(vek+i*m+2))>=0.0)) 
           tant[0] +=1.0;
         if(((*(vek+i*m))<0.0)&&((*(vek+i*m+1))>=0.0)&&((*(vek+i*m+2))>=0.0)) 
           tant[1] +=1.0;
         if(((*(vek+i*m))>=0.0)&&((*(vek+i*m+1))<0.0)&&((*(vek+i*m+2))>=0.0)) 
           tant[2] +=1.0;
         if(((*(vek+i*m))<0.0)&&((*(vek+i*m+1))<0.0)&&((*(vek+i*m+2))>=0.0)) 
           tant[3] +=1.0;
         if(((*(vek+i*m))>=0.0)&&((*(vek+i*m+1))>=0.0)&&((*(vek+i*m+2))<0.0)) 
           tant[4] +=1.0;
         if(((*(vek+i*m))<0.0)&&((*(vek+i*m+1))>=0.0)&&((*(vek+i*m+2))<0.0)) 
           tant[5] +=1.0;
         if(((*(vek+i*m))>=0.0)&&((*(vek+i*m+1))<0.0)&&((*(vek+i*m+2))<0.0)) 
           tant[6] +=1.0;
         if(((*(vek+i*m))<0.0)&&((*(vek+i*m+1))<0.0)&&((*(vek+i*m+2))<0.0))
           tant[7] +=1.0;
      }
      eps = sqrt(8.0*f2m1);
      i   = 0;
      printf("\n  3-D test: each octant should get 1/8\n");
      if(fabsf(f2m1*tant[0]-0.125)>eps){
         printf("  1st octant frac=%e\n",f2m1*tant[0]);
         i++;
      }
      if(fabsf(f2m1*tant[1]-0.125)>eps){
         printf("  2nd octant frac=%e\n",f2m1*tant[1]);
         i++;
      }
      if(fabsf(f2m1*tant[2]-0.125)>eps){
         printf("  3rd octant frac=%e\n",f2m1*tant[2]);
         i++;
      }
      if(fabsf(f2m1*tant[3]-0.125)>eps){
         printf("  4th octant frac=%e\n",f2m1*tant[3]);
         i++;
      }
      if(fabsf(f2m1*tant[4]-0.125)>eps){
         printf("  5th octant frac=%e\n",f2m1*tant[4]);
         i++;
      }
      if(fabsf(f2m1*tant[5]-0.125)>eps){
         printf("  6th octant frac=%e\n",f2m1*tant[5]);
         i++;
      }
      if(fabsf(f2m1*tant[6]-0.125)>eps){
         printf("  7th octant frac=%e\n",f2m1*tant[6]);
         i++;
      }
      if(fabsf(f2m1*tant[7]-0.125)>eps){
         printf("  8th octant frac=%e\n",f2m1*tant[7]);
         i++;
      }
      if(i>0){
         printf("  %d |octant fraction(s) - 1/8| >< sqrt(8/N)\n",i);
      } else {
         printf("  **** 3-D test passed **** \n");
      }
      free(vek); free(tant);
   } 
} 
/*                      */
/* cut here for library */
/*                      */
/* global variables */
int p,q,ptr;
int ptab[]={607,1279,2281,3217,4423,9689,19937,23209,44497};
int qtab[]={273, 418,1252, 576,2098,5502, 9842,13470,21034};
float *buf;

float zufall()
{
   void zufill();
   float rval;
   if(ptr<p-1){
     rval = buf[ptr++];
     return(rval);
   } else {
     rval = buf[ptr];
     zufill();
     ptr  = 0;
     return(rval);
   }
}
void zufill()
{
   float t;
   int i,ip,iq;
   iq = p-q;
   for(i=0;i<q;i++){
      t      = buf[i]+buf[iq++];
      if(t>1.0) t = t - 1.0;
/*    t      = t - (float)((int) t); */
      buf[i] = t;
   }
   iq = 0;
   for(i=q;i<p;i++){
      t      = buf[i]+buf[iq++];
      if(t>1.0) t = t - 1.0;
/*    t      = t - (float)((int) t); */
      buf[i] = t;
   }
}
void zufall_init(int mlev,int seed)
{
   int i;
   static int pinit=-1;
   static float ds;
   float ggl();
/* select quality level: */
   if((mlev>=0)&&(mlev<10)){
      p = ptab[mlev];
      q = qtab[mlev];
      if(p!=pinit){
         if(pinit>0){
            free(buf);
         }
         buf   = (float *) malloc(p*sizeof(float));
         pinit = p;
      }
   } else {
      printf(" zufall initialization error: 0 <= mlev=%d < 10\n",mlev);
   }
   if(seed <= 0){
      ds = 331.0;
   } else {
      ds  = (float) seed;
   }
   for(i=0;i<p;i++){
      buf[i] = ggl(&ds);
   }

   ptr = 0;
}
void zufall_spar(char *fn)
{
   int i;
   FILE *fp;
   fp = fopen(fn,"w");
   fprintf(fp,"%d %d %d\n",p,q,ptr);
   for(i=0;i<p;i++){
      fprintf(fp,"%e\n",buf[i]);
   }
   fclose(fp);
}
void zufall_rest(char *fn)
{
   int i;
   FILE *fp;
   fp = fopen(fn,"r");
   fscanf(fp,"%d %d %d\n",&p,&q,&ptr);
   for(i=0;i<p;i++){
      fscanf(fp,"%e\n",&buf[i]);
   }
   fclose(fp);
}
float ggl(float *ds)
{
/* generate u(0,1] distributed random numbers */
   double t,d2=2147483647.0;
   t   = *ds;
   t   = fmod(16807.0*t,d2);
   *ds = t;
   return(t/d2);
}
float gbuf[1024];
float zuf_g()
{
   void zuf_gfill();
   static int init=1,ptg;
   if(init){  /* buffer not initialized */
      zuf_gfill();
      init = 0;
      ptg  = 0;
   }
   if(ptg < 1024){
      ptg++;
      return(gbuf[ptg]);
   } else {
      zuf_gfill();
      ptg = 1;
      return(gbuf[0]);
   }
}
void zuf_gfill()
{
   float zufall();
   double t,u1,u2,z1,z2;
   static double twopi=6.283185307179586477;
   int i,ii;
   for(i=0;i<512;i++){
      ii = 2*i;
      u1 = (double) zufall();
      u2 = (double) zufall();
      t  = sqrt(-2.0*(log(1.0-u2)));
      z1 = cos(twopi*u1)*t; 
      z2 = sin(twopi*u1)*t;
      gbuf[ii]   = (float) z1; 
      gbuf[ii+1] = (float) z2;
   }
} 
float *nsphere(int n, int sw)
{

/* nsphere returns pointer to n-D array of elements in n-sphere
   if sw=0, points are on unit n-sphere,
   if sw>1, points are inside n-sphere and uniformly
   distributed */

   static float fnm1, *v;
   float rad,u;
   float zufall(),zuf_g();
   int i;
   static int init=-33;
   if(n!=init){
      if(init>0) free(v);
      v    = (float *)malloc(n*sizeof(double));
      init = n;
      fnm1 = 1.0/((float) n);
   }
   for(i=0;i<n;i++){
      v[i] = zuf_g();
   }
   rad = 0.0;
   for(i=0;i<n;i++){
      rad += v[i]*v[i];
   }
   rad = 1.0/((float)sqrt((double)rad));
   for(i=0;i<n;i++){
      v[i] *= rad;
   }
   if(sw==0){
      return(v);
   } else {
      u   = zufall();
      rad = (float) pow(((double)u),((double)fnm1));
      for(i=0;i<n;i++){
         v[i] *= rad;
      }
      return(v);
   }
}