/********************************************************************
/*  
/* Random number generator LFSR113
/* C++ class LFSR113
/*
/* Following a suggestion of Pierre L'Ecuyer 1997
/* "Tables of maximally equidistributed combined LFSR generators"
/* see http://www.iro.umontreal.ca/~lecuyer/papers.html
/*
/* The generator is initialized with a positive integer seed
/* when an instance is defined
/* 
/* Implementation by Thomas Vojta, vojtat@umr.edu
/* 
/* History:
/*  04 Sep 2007    v0.9    first implementation
/* 
/***********************************************************************/


#include "lfsr113_64.h"

lfsr113::lfsr113(int seed)
{
	this->lfsrinit(seed);
}


lfsr113::~lfsr113()
{

}

/** This method initializes random generator. */



int lfsr113::ilfsr113(void)
{
   unsigned int 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 ) >> 1);
}

double lfsr113::rlfsr113(void)
{
   unsigned int 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 lfsr113::lfsrinit(int idum)
{
   int k;
   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;
 
   idum=ilfsr113();
}


double lfsr113::drnd(void)
{
  return 	this->rlfsr113();
}

int lfsr113::lrnd (void)
{
	return this->ilfsr113();
}