#include <stdio.h>
#include <math.h>
#include <time.h>
#include "rngheader.h"
#define N 100000
main()
{
/* 
   Example of G. H. Gonnet's RNG, this variant does 16/time:
   that is, 16 independent streams.
   
   Header file rngheader.h is required: "rng" is defined
   in this header by a comma operator procedure.
                                             wpp 18/11/2003 
*/ 
   double fnm1,xv,x1,x2,x3,x4,x[N];
   float t1,t2;
   int i,j,initseed=331;
   zufall_init(initseed,&rng_buf_ptr,rng_buf);
   x1 = 0.0; x2 = 0.0; x3 = 0.0; x4 = 0.0; 
   fnm1 = 1.0/((double) N);

   t1 = ((float)clock())/((float) CLOCKS_PER_SEC);
   for(i=0;i<N;i++){
      x[i] = rng;
   }
   t2 = ((float)clock())/((float) CLOCKS_PER_SEC);
   t2 -= t1;
   printf("Time for 100000 ghg = %e seconds\n",t2);
   for(i=0;i<N;i++){
      xv   = x[i] - 0.5;
      x1  += xv;
      x2  += xv*xv;
      x3  += xv*xv*xv;
      x4  += xv*xv*xv*xv;
   }

/* x1 := <(x-1/2)>,   x2 := <(x-1/2)^2> 
   x3 := <(x-1/2)^3>, x4 := <(x-1/2)^4> */
   x1 *= fnm1; x2 *= fnm1; x3 *= fnm1, x4 *= fnm1; 
   x4  = x4 - 4.0*x1*x3 + 6.0*x1*x1*x2 - 4.0*x1*x1*x1*x1;
   x3  = x3 - 3.0*x1*x2 + 2.0*x1*x1*x1;
   x2 -= x1*x1;
   printf("Sample size N = %d, initial seed = %d \n",N,initseed);
   printf("      <x-1/2> = %13.6e (compare to 0.0)\n",x1);
   printf("12<(x-1/2)^2> = %13.6e (compare to 1.0)\n",12.0*x2);
   printf("  <(x-1/2)^3> = %13.6e (compare to 0.0)\n",x3);
   printf("80<(x-1/2)^4> = %13.6e (compare to 1.0)\n",80.0*x4);
}
#define ADD_ERR(aa,bb,cc) { double tmpa,tmpb,tmpe; \
   tmpa = aa; tmpb = bb; \
   aa += tmpb; \
   if(abs(tmpa) > abs(tmpb)) { \
      tmpe = aa-tmpa; cc += tmpb-tmpe; \
   } else { \
      tmpe = aa-tmpb; cc += tmpa-tmpe; \
   } }

#define bij(i,k) buf[i+16*k]
void zufall(int *p,double *buf)
{

/* 
   Gaston H. Gonnet's (30 Apr., 2003) uniform random number generator,
   a variant on the idea of Green-Smith-Klem: see D. Knuth, vol. 2, pages 
   26-27 (1969 edition). GHG verified this polynomial for irreducibility 
   on the field of integers modulo 2^53 using Paul Zimmerman's version
   of the LLL algorithm of A. K. Lenstra, H. W. Lenstra, and L. Lovasz, 
   Math. Ann. 261, 1982. This variant with 16/time: wpp, Nov. 2003.  
*/

   double t,r0,r1;
   int i,j0,j1,j2,j3,j;
   int off[4]  = {-5, -43, -97, -128};
   double pwr2[4] = {8.0e0,1.28e2, 1.5625e-2, 1.0};
   for(j=0;j<128;j++){
      j0 = (j-off[0])&(127);
      j1 = (j-off[1])&(127);
      j2 = (j-off[2])&(127);
      j3 = (j-off[3])&(127);
      for(i=0;i<16;i++){
         r0    = bij(i,j);
         r1    = 0;
         ADD_ERR(r0,bij(i,j0)*pwr2[0],r1); 
         ADD_ERR(r0,bij(i,j1)*pwr2[1],r1); 
         ADD_ERR(r0,bij(i,j2)*pwr2[2],r1); 
         ADD_ERR(r0,bij(i,j3)*pwr2[3],r1);
         r0 -= floor(r0+r1);
         t   = r0 + r1;
         if(t > 1){
            t -= 1;
         } else if(t < 0) {
            t += 1;
         }
         bij(i,j) = t;
      }
   }
   *p = 0;
}
void zufall_init(int initseed, int *p,double *buf)
{
    int i,j;
    double ggl(double *);
    static double seed;
    seed = (double) initseed; 
/* initialize buffer */
    for(j=0;j<128;j++){
       for(i=0;i<16;i++){
          bij(i,j) = ggl(&seed);
       }
    }
    *p = 2048;
}
#undef bij
double ggl(double *ds)
{
/* Generator for u(0.,1.] distributed random numbers. 
   This version of ggl (same as RNUM in IMSL) by
   W. Petersen and M. Troyer, ETHZ, 25 Oct. 2002. */
   double t,d2=0.2147483647e10;
   t   = *ds;
   t   = fmod(0.16807e5*t,d2);
   *ds = t;
   return((t-1.0e0)/(d2-1.0e0));
}