/************************** MOTHER.CPP ****************** AgF 1999-03-03 *
*  'Mother-of-All' random number generator                               *
*                                                                        *
*  This is a multiply-with-carry type of random number generator         *
*  invented by George Marsaglia.  The algorithm is:                      *
*  S = 2111111111*X[n-4] + 1492*X[n-3] + 1776*X[n-2] + 5115*X[n-1] + C   *
*  X[n] = S modulo 2^32                                                  *
*  C = floor(S / 2^32)                                                   *
*                                                                        *
*  This implementation uses a long double for C. Note that not all C++   *
*  compilers support the long double (80-bit) floating point format.     *
*  You will get an error message if your compiler doesn't support this.  *
*************************************************************************/

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

typedef unsigned long uint32;            // 32 bit unsigned integer

class TMotherRandomGenerator {           // 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
  TMotherRandomGenerator(uint32 seed=-1);// constructor
  protected:
  double x[5];                           // history buffer
  int imin, iinterval;                   // interval for iRandom
  };


// constructor:
TMotherRandomGenerator::TMotherRandomGenerator(uint32 seed) {
  // check that compiler supports 80-bit long double
  assert(sizeof(long double)>9);
  // initialize
  RandomInit(seed);  SetInterval(0, 99);}


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


// returns a random number between 0 and 1:
double TMotherRandomGenerator::Random() {
  long double c;
  c = 2111111111.0 * x[3] +
      1492.0 * (x[3] = x[2]) +
      1776.0 * (x[2] = x[1]) +
      5115.0 * (x[1] = x[0]) +
      x[4];
  x[4] = floorl(c);
  x[0] = c - x[4];
  x[4] = x[4] * (1./(65536.*65536.));
  return x[0];}


// returns integer random number in desired interval:
int TMotherRandomGenerator::iRandom() {
  int i = iinterval * Random();
  if (i >= iinterval) i = iinterval;
  return imin + i;}


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

  // make random numbers and put them into the buffer
  for (i=0; i<5; i++) {
    seed ^= seed << 13; seed ^= seed >> 17; seed ^= seed << 5;
    x[i] = seed * (1./(65536.*65536.));}
  // randomize some more
  for (i=0; i<19; i++) Random();}


//////////////////////////////////////////////////////////////////////////
// Example showing how to use the random number generator:
void main() {
  TMotherRandomGenerator 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<100; i++) {
    // generate 100 random floating point numbers and 100 random integers
    printf("\n%14.10f   %2i", (double)RG.Random(), RG.iRandom());}

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