/************************* RANROTA.CPP ****************** AgF 1999-03-03 *
*  Random Number generator 'RANROT' type A                               *
*                                                                        *
*  This is a lagged-Fibonacci type of random number generator with       *
*  rotation of bits.  The algorithm is:                                  *
*  X[n] = ((X[n-j] + X[n-k]) modulo 2^b) rotl r                          *
*                                                                        *
*  The last k values of X are stored in a circular buffer named          *
*  randbuffer.                                                           *
*  The code includes a self-test facility which will detect any          *
*  repetition of previous states.                                        *
*  The function uses a fast method for conversion to floating point.     *
*  This method relies on floating point numbers being stored in the      *
*  standard 64-bit IEEE format.                                          *
*                                                                        *
*  The theory of the RANROT type of generators and the reason for the    *
*  self-test are described at www.agner.org/random/ranrot.htm            *
*                                                                        *
*************************************************************************/

// names of header files may differ
#include <stdio.h>
#include <conio.h>
#include <math.h>
#include <mem.h>
#include <stdlib.h>
#include <assert.h>

// define if you want self-test:
#define SELF_TEST

// define 32-bits unsigned integer
typedef unsigned long uint32;

// If your system doesn't have a rotate function for 32 bits integers,
// then define it thus:
// uint32 _lrotl (uint32 x, int r) {
//   return (x << r) | (x >> (sizeof(x)*8-r));}

// define parameters
const int KK = 11;
const int JJ =  7;
const int RR = 13;


class TRanrotAGenerator {             // encapsulate random number generator
  public:
  void RandomInit(uint32 seed);       // initialization
  void SetInterval(int min, int max); // set interval for iRandom
  int iRandom();                      // get integer random number
  double Random();                    // get floating point random number
  TRanrotAGenerator(uint32 seed=-1);  // constructor
  protected:
  int p1, p2;                         // indexes into buffer
  int imin, iinterval;                // interval for iRandom
  uint32 randbuffer[KK];              // history buffer
#ifdef SELF_TEST
  uint32 randbufcopy[KK*2];           // used for self-test
#endif
};


// constructor:
TRanrotAGenerator::TRanrotAGenerator(uint32 seed) {
  RandomInit(seed);  SetInterval(0, 99);}


// returns a random number between 0 and 1:
double TRanrotAGenerator::Random() {
  uint32 x;
  // generate next random number
  x = randbuffer[p1] = _lrotl(randbuffer[p1] + randbuffer[p2], RR);
  // rotate list pointers
  if (--p1 < 0) p1 = KK - 1;
  if (--p2 < 0) p2 = KK - 1;
#ifdef SELF_TEST
  // perform self-test
  if (randbuffer[p1] == randbufcopy[0] &&
    memcmp(randbuffer, randbufcopy+KK-p1, KK*sizeof(uint32)) == 0) {
      // self-test failed
      if ((p2 + KK - p1) % KK != JJ) {
        // note: the way of printing error messages depends on system
        // In Windows you may use FatalAppExit
        printf("Random number generator not initialized");}
      else {
        printf("Random number generator returned to initial state");}
      exit(1);}
#endif
  // fast conversion to float:
  union {
    double randp1;
    uint32 randbits[2];};
  randbits[0] = x << 20;
  randbits[1] = (x >> 12) | 0x3FF00000;
  return randp1 - 1.0;}


int TRanrotAGenerator::iRandom() {
  // get integer random number
  int i = iinterval * Random();
  if (i >= iinterval) i = iinterval;
  return imin + i;}


// set interval for iRandom
void TRanrotAGenerator::SetInterval(int min, int max) {
  imin = min; iinterval = max - min + 1;}


void TRanrotAGenerator::RandomInit (uint32 seed) {
  // this function initializes the random number generator.
  int i;
  // make sure seed != 0
  if (seed==0) seed--;

  // make random numbers and put them into the buffer
  for (i=0; i<KK; i++) {
    seed ^= seed << 13; seed ^= seed >> 17; seed ^= seed << 5;
    randbuffer[i] = seed;}

  // check that the right data formats are used by compiler:
  union {
    double randp1;
    uint32 randbits[2];};
  randp1 = 1.5; assert(randbits[1]==0x3FF80000);

  // initialize pointers to circular buffer
  p1 = 0;  p2 = JJ;
#ifdef SELF_TEST
  // store state for self-test
  memcpy (randbufcopy, randbuffer, KK*sizeof(uint32));
  memcpy (randbufcopy+KK, randbuffer, KK*sizeof(uint32));
#endif
  // randomize some more
  for (i=0; i<9; i++) Random();
}


//////////////////////////////////////////////////////////////////////////
// Example showing how to use the random number generator:
void main() {
  TRanrotAGenerator RG;    // define random number generator

  RG.RandomInit(time(0)); // use time as random seed

  RG.SetInterval(1, 10);  // set interval for integer output

  for (int i=0; i<20; i++) {
    // generate 10 random floating point numbers and 10 random integers
    printf("\n%14.10f   %2i", RG.Random(), RG.iRandom());}

  getch();  // wait for user to press any key
  }