#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <stdlib.h>
#define NPTS 200000
#define DIM_RAND 256
double rand_a[DIM_RAND];

main()
/* 
   Simple timing test program for G. H. Gonnet's "ghg" 
   random number generator. ggl is used to initialize the 
   buffer for ghg. 
                                           wpp 30/4/2003.
*/
{
    double a[NPTS];
    int seed,ia[20];
    int i,ii;
    double ggl(double*),ghg(void);
    double t1;
    static double fseed;
/* initialize buffer */
    fseed = (double) 331;
    for(i=0;i<DIM_RAND;i++){
        rand_a[i] = ggl(&fseed);
    }
    t1 = (double) clock();
    for(i=0;i<NPTS;i++){
        a[i] = ghg();
    }
    t1 = (double) clock() - t1;
    t1 = t1/CLOCKS_PER_SEC;
    printf(" time for %d sample = %e\n",NPTS,t1);
    for(i=0;i<20;i++) ia[i] = 0;
    for(i=0;i<NPTS;i++){
       ii = (int) (a[i]*20.0);
       ia[ii]++;
    }
    for(i=0;i<20;i++){
       printf("ia[%d] = %d\n",i,ia[i]);
    }
}
void ghginit()
{
    int seed;
    int i;
    double ggl(double *);
    static double fseed;
/* initialize buffer */
    fseed = (double) 66790;
    for(i=0;i<DIM_RAND;i++){
        rand_a[i] = ggl(&fseed);
    }
}
#define DIM_RAND 128
#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; \
   } }

int off[3]  = {19, 67, 128};
double pwr2[3] = {3.2e1, 7.8125e-3, 1.0};
int randi=0;
double ghg()
{
/* 
   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 Maple and the 
   LLL algorithm of A. K. Lenstra, H. W. Lenstra, and L. Lovasz, 
   Math. Ann. 261, 1982. This variant: wpp, Oct. 2002.  
*/
   double t,r0,r1;
   int j;
/* randi = (randi+1) & (DIM_RAND-1); wpp fix: increment at end 5/9/2003 */ 
   r0    = rand_a[randi];
   r1    = 0;
   ADD_ERR(r0,rand_a[(randi-off[0])&(DIM_RAND-1)]*pwr2[0],r1); 
   ADD_ERR(r0,rand_a[(randi-off[1])&(DIM_RAND-1)]*pwr2[1],r1); 
   ADD_ERR(r0,rand_a[(randi-off[2])&(DIM_RAND-1)]*pwr2[2],r1); 
   r0 -= floor(r0+r1);
   t   = r0 + r1;
   if(t > 1){
      t -= 1;
   } else if(t < 0) {
      t += 1;
   }
   rand_a[randi] = t;
   randi = (randi+1) & (DIM_RAND-1);
   return(t);
}
#include <math.h>
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));
}