/**************************************************************************/
/*
/* Random number generator KISS05 after a suggestion by George Marsaglia
/* in "Random numbers for C: The END?" posted on sci.crypt.random-numbers
/* in 1999
/*
/* 2005 version as in "double precision RNGs" in  sci.math.num-analysis  
/* http://sci.tech-archive.net/Archive/sci.math.num-analysis/2005-11/msg00352.html
/*
/* The  KISS (Keep It Simple Stupid) random number generator. Combines:
/* (1) The congruential generator x(n)=69069*x(n-1)+1327217885, period 2^32.
/* (2) A 3-shift shift-register generator, period 2^32-1,
/* (3) Two 16-bit multiply-with-carry generators, period 597273182964842497>2^59
/* Overall period > 2^123  
/* 
/* 
/* A call to rkiss05() gives one random real in the interval [0,1),
/* i.e., 0 <= rkiss05 < 1
/*
/* Before using rkiss05 call kissinit(seed) to initialize
/* the generator by random integers produced by Park/Millers
/* minimal standard LCG.
/* Seed should be any positive integer.
/* 
/* Implementation by Thomas Vojta, vojta@mst.edu
/* built on a Fortran module found at www.fortran.com
/* 
/* 
/* History:
/*        v0.9     Dec 11, 2010    first implementation
/*        v0.91    Dec 13, 2010    extra shuffle of seed in kissinit
/*
/*
/**************************************************************************/

#define IA 16807
#define IM 2147483647
#define IQ 127773
#define IR 2836

unsigned long x,y,z,w;

double rkiss05()
   {
   x = 69069*x+1327217885; 
   y^= (y<<13); y^=(y>>17); y^=(y<<5);
   z = 18000 * (z & 65535) + (z >> 16);
   w = 30903 * (w & 65535) + (w >> 16);
   return( (x + y + (z << 16) +w) * 2.328306436538696e-10);
   }

void kissinit(long idum)
   {
   long k;
   double d;
   idum= abs(1099087573 * idum);
   if (idum == 0) idum = 1;
   if (idum >= IM) idum = IM-1;

   k=(idum)/IQ;
   idum=IA*(idum-k*IQ)-IR*k;
   if (idum < 0) idum += IM;
   if (idum < 1) x=idum+1; else x=idum;
   k=(idum)/IQ;
   idum=IA*(idum-k*IQ)-IR*k;
   if (idum < 0) idum += IM;
   if (idum < 1) y=idum+1; else y=idum;
   k=(idum)/IQ;
   idum=IA*(idum-k*IQ)-IR*k;
   if (idum < 0) idum += IM;
   if (idum < 1) z=idum+1; else z=idum;
   k=(idum)/IQ;
   idum=IA*(idum-k*IQ)-IR*k;
   if (idum < 0) idum += IM;
   if (idum < 1) w=idum+1; else w=idum;
   
   d=rkiss05();
   }