Random.h

// Random.h: Definition and Implementation of Random Number Distribution Class
//           Ref: Richard Saucier, "Computer Generation of Statistical
//                Distributions," ARL-TR-2168, US Army Research Laboratory,
//                Aberdeen Proving Ground, MD, 21005-5068, March 2000.
 
#ifndef RANDOM_H
#define RANDOM_H
 
#include <array>
#include <cstdint>    // for int32_t and its limits
using namespace std;  // NOTE: to be removed!
 
// Constants for Tausworthe random bit generator
// Ref: Tausworthe, Robert C., "Random Numbers Generated by Linear Recurrence
//      Modulo Two," Mathematics of Computation, vol. 19, pp. 201-209, 1965.
 
static const unsigned DEGREE_MAX = 32;   // max degree (bits per word)
 
static const unsigned BIT[ 1 + DEGREE_MAX ] = {
 
// Hexadecimal       Value      Degree
// -----------       -----      ------
   0x00000000,    // 0          0
   0x00000001,    // 2^0        1
   0x00000002,    // 2^1        2
   0x00000004,    // 2^2        3
   0x00000008,    // 2^3        4
   0x00000010,    // 2^4        5
   0x00000020,    // 2^5        6
   0x00000040,    // 2^6        7
   0x00000080,    // 2^7        8
   0x00000100,    // 2^8        9
   0x00000200,    // 2^9        10
   0x00000400,    // 2^10       11
   0x00000800,    // 2^11       12
   0x00001000,    // 2^12       13
   0x00002000,    // 2^13       14
   0x00004000,    // 2^14       15
   0x00008000,    // 2^15       16
   0x00010000,    // 2^16       17
   0x00020000,    // 2^17       18
   0x00040000,    // 2^18       19
   0x00080000,    // 2^19       20
   0x00100000,    // 2^20       21
   0x00200000,    // 2^21       22
   0x00400000,    // 2^22       23
   0x00800000,    // 2^23       24
   0x01000000,    // 2^24       25
   0x02000000,    // 2^25       26
   0x04000000,    // 2^26       27
   0x08000000,    // 2^27       28
   0x10000000,    // 2^28       29
   0x20000000,    // 2^29       30
   0x40000000,    // 2^30       31
   0x80000000     // 2^31       32
};
 
// Coefficients that define a primitive polynomial (mod 2)
// Ref: Watson, E. J., "Primitive Polynomials (Mod 2),"
//      Mathematics of Computation, vol. 16, pp. 368-369, 1962.
 
static const unsigned MASK[ 1 + DEGREE_MAX ] = {
 
   BIT[0],                                      // 0
   BIT[0],                                      // 1
   BIT[1],                                      // 2
   BIT[1],                                      // 3
   BIT[1],                                      // 4
   BIT[2],                                      // 5
   BIT[1],                                      // 6
   BIT[1],                                      // 7
   BIT[4] + BIT[3] + BIT[2],                    // 8
   BIT[4],                                      // 9
   BIT[3],                                      // 10
   BIT[2],                                      // 11
   BIT[6] + BIT[4] + BIT[1],                    // 12
   BIT[4] + BIT[3] + BIT[1],                    // 13
   BIT[5] + BIT[3] + BIT[1],                    // 14
   BIT[1],                                      // 15
   BIT[5] + BIT[3] + BIT[2],                    // 16
   BIT[3],                                      // 17
   BIT[5] + BIT[2] + BIT[1],                    // 18
   BIT[5] + BIT[2] + BIT[1],                    // 19
   BIT[3],                                      // 20
   BIT[2],                                      // 21
   BIT[1],                                      // 22
   BIT[5],                                      // 23
   BIT[4] + BIT[3] + BIT[1],                    // 24
   BIT[3],                                      // 25
   BIT[6] + BIT[2] + BIT[1],                    // 26
   BIT[5] + BIT[2] + BIT[1],                    // 27
   BIT[3],                                      // 28
   BIT[2],                                      // 29
   BIT[6] + BIT[4] + BIT[1],                    // 30
   BIT[3],                                      // 31
   BIT[7] + BIT[5] + BIT[3] + BIT[2] + BIT[1]   // 32
};
 
 
 
class Random {
 
public:
   // overloaded relational operators
   friend bool operator==(const Random& a, const Random& b);
   inline friend bool operator!=(const Random& a, const Random& b) { return !(a == b); }
 
   Random( int32_t seed );                 
   Random( );                              
   Random( int32_t seed, bool f );         
   Random( const Random& r ) = default;            
   Random& operator=( const Random& r ) = default; 
 
   void reset( int32_t seed );   
   void reset( );            
   double beta( double v, double w,
                double xMin = 0., double xMax = 1. );      
   double gamma( double a, double b, double c );           
   double exponential( double a = 0., double c = 1. );     
   double normal( double mu = 0., double sigma = 1. );     
 
   double uniform( double xMin = 0., double xMax = 1. );   
   int poisson( double mu );   
   int uniformDiscrete( int i, int j );   
 
   int32_t getSeed(); 
 
private:
 
/* Variables used in normalModified function */
   struct {
      bool f; 
      double pR2; 
      std::array<double,2> pXY;
   } _gauss;
 
   static constexpr int32_t   _M    = 0x7fffffff; // 2147483647 (Mersenne prime 2^31-1)
   static constexpr int32_t   A_    = 0x10ff5;    // 69621
   static constexpr int32_t   Q_    = 0x787d;     // 30845
   static constexpr int32_t   R_    = 0x5d5e;     // 23902
   static constexpr double    _F    = 1. / _M;
   static constexpr short     _NTAB = 32;         // arbitrary length of shuffle table
   static constexpr int32_t   _DIV  = 1+(_M-1)/_NTAB;
   int32_t        _table[ _NTAB ];            // shuffle table of seeds
   int32_t        _next;                      // seed to be used as index into table
   int32_t        _seed;                      // current random number seed
   unsigned       _seed2;                     // seed for tausworthe random bit
 
 
   void _seedTable( );          // seeds the shuffle table
   double _u( );                // uniform rng
 
   double normalModified( double mu = 0.,                  // Normal
                              double sigma = 1. );
};
 
#endif