/********************************************************************
/* Random number generator RLFSR113  - real version
/*
/* Following a suggestion of Pierre L'Ecuyer 1997
/* "Tables of maximally equidistributed combined LFSR generators"
/* see http://www.iro.umontreal.ca/~lecuyer/papers.html
/*
/* A call to rlfsr113() gives one random double in the 
/* open interval (0,1).
/*
/* Before using rlfsr113 call lfsrinit(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@physik.tu-chemnitz.de
/* 
/* History:
/*  04 Feb 1998    v0.9    first implementation
/* 
/***********************************************************************/

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

unsigned long z1, z2, z3, z4;

double rlfsr113 ()
   {
   unsigned long b;
   b  = (((z1 << 6) ^ z1) >> 13);
   z1 = (((z1 & 4294967294) << 18) ^b);
   b  = (((z2 << 2) ^ z2) >> 27);
   z2 = (((z2 & 4294967288) << 2) ^b);
   b  = (((z3 << 13) ^ z3) >> 21);
   z3 = (((z3 & 4294967280) << 7) ^b);
   b  = (((z4 << 3) ^ z4) >> 12);
   z4 = (((z4 & 4294967168) << 13) ^b);
   return((z1 ^ z2 ^ z3 ^ z4 ) * 2.3283064365e-10);
   }

void lfsrinit(long idum)
   {
   long k;
   double d;
   if (idum <= 0) idum = 1;
   k=(idum)/IQ;
   idum=IA*(idum-k*IQ)-IR*k;
   if (idum < 0) idum += IM;
   if (idum < 2) z1=idum+2; else z1=idum;
   k=(idum)/IQ;
   idum=IA*(idum-k*IQ)-IR*k;
   if (idum < 0) idum += IM;
   if (idum < 8) z2=idum+8; else z2=idum;
   k=(idum)/IQ;
   idum=IA*(idum-k*IQ)-IR*k;
   if (idum < 0) idum += IM;
   if (idum < 16) z3=idum+16; else z3=idum;
   k=(idum)/IQ;
   idum=IA*(idum-k*IQ)-IR*k;
   if (idum < 0) idum += IM;
   if (idum < 128) z4=idum+128; else z4=idum;
   
   d=rlfsr113();
   }