//+------------------------------------------------------------------+
//|                                                       Random.mqh |
//|                                        Copyright © 2018, Amr Ali |
//|                             https://www.mql5.com/en/users/amrali |
//+------------------------------------------------------------------+
#property description "Random number generation using the 32-bit PCG generator."
//+------------------------------------------------------------------+
//| Class CRandom                                                    |
//| Usage: Random number generation using the 32-bit PCG generator.  |
//+------------------------------------------------------------------+
class CRandom
  {
protected:
   //--- The generator uses two 64-bit integers for its internal state
   ulong             m_state;
   ulong             m_inc;
   //---
   void              pcg32_srandom_r(const ulong initstate,const ulong initseq);
   uint              pcg32_random_r(void);
   uint              pcg32_boundedrand_r(const uint bound);
   ulong             pcg32x2_random_r(void);

public:
   //--- Constructor and destructor
                     CRandom(void);                                      // Constructor (auto-seed)
                     CRandom(const ulong initstate,const ulong initseq); // Constructor (custom seed)
                    ~CRandom(void) {   }

   //--- Methods for generation of random numbers
   int               RandomInteger(void);                                // Random integer [0,2147483648)
   int               RandomInteger(const int max);                       // Random integer [0,max)
   int               RandomInteger(const int min,const int max);         // Random integer [min,max)

   double            RandomDouble(void);                                 // Random double  [0.0,1.0)
   double            RandomDouble(const double max);                     // Random double  [0.0,max)
   double            RandomDouble(const double min,const double max);    // Random double  [min,max)

   bool              RandomBoolean(void);                                // Random true/false (equal probability)
   bool              RandomBoolean(const double prob_true);              // Random true/false (with prob_true)

   double            RandomNormal(void);                                 // Random double (normal distribution)
   double            RandomNormal(const double mu,const double sigma);   // Random double (normal distribution)
  };
//+------------------------------------------------------------------+
//| Initializes a new instance of the CRandom class with auto-seed.  |
//+------------------------------------------------------------------+
void CRandom::CRandom(void)
  {
   pcg32_srandom_r(GetTickCount(), ChartID() * MathRand());
  }
//+------------------------------------------------------------------+
//| Initializes a new instance of the CRandom class and set the      |
//| initial values for the generator of random numbers.              |
//| Seed the rng. Specified in two parts, state initializer and a    |
//| sequence selection constant (a.k.a. stream id)                   |
//+------------------------------------------------------------------+
void CRandom::CRandom(const ulong initstate,const ulong initseq)
  {
   pcg32_srandom_r(initstate, initseq);
  }
//+------------------------------------------------------------------+
//| Uniformly distributed random integer in the range [0,2147483648) |
//| The maximum possible random integer is 2147483647 (INT_MAX).     |
//+------------------------------------------------------------------+
int CRandom::RandomInteger(void)
  {
//--- return result
   return (int)(pcg32_random_r() >> 1);
  }
//+------------------------------------------------------------------+
//| Uniformly distributed random integer in the range [0,max)        |
//| max : exclusive and must be greater than zero.                   |
//+------------------------------------------------------------------+
int CRandom::RandomInteger(const int max)
  {
   if(max <= 0)
      return(-1);
//--- return result
   return (int)pcg32_boundedrand_r(max);
  }
//+------------------------------------------------------------------+
//| Uniformly distributed random integer in the range [min,max)      |
//| min : inclusive                                                  |
//| max : exclusive and must be greater than min.                    |
//+------------------------------------------------------------------+
int CRandom::RandomInteger(const int min,const int max)
  {
   if(max <= min)
      return(-1);
//--- return result
   return (int)pcg32_boundedrand_r(max - min) + min;
  }
//+------------------------------------------------------------------+
//| Uniformly distributed random double in the range [0.0, 1.0)      |
//+------------------------------------------------------------------+
//| https://docs.python.org/3/library/random.html#recipes            |
//| http://mumble.net/~campbell/tmp/random_real.c                    |
//| The mantissa comes from a uniform distribution of integers in    |
//| the range 0 <= mantissa < 2^52. The exponent comes from a        |
//| geometric distribution where exponents smaller than -53 occur    |
//| half as often as the next larger exponent. The recipe is         |
//| conceptually equivalent to an algorithm that chooses from all    |
//| the multiples of 2^-1074 in the range 0.0 = x < 1.0. All floats  |
//| in the interval [0, 1) are possible selections.                  |
//+------------------------------------------------------------------+
double CRandom::RandomDouble(void)
  {
   ulong mantissa = pcg32x2_random_r() >> 12;
   double significand = 1.0 + mantissa * DBL_EPSILON;

//--- sample p from geometric(2) by finding the first 1-bit
//--- in a stream of bits generated uniformly at random and
//--- interpret it as a power of two, 00101000000000000000.
   double power2 = 2.0;
   uint x;
   while((x = pcg32_random_r()) == 0)
     {
      power2 *= 4294967296.0; // 2^32
     }

//--- x & -x creates a bitmask for the rightmost set bit in x.
//--- power2 is 2 with probability 1/2, 4 with probability 1/4,
//--- 8 with probability 1/8, and so on.
   power2 *= (x & (0-x));

//--- scale 2^52 numbers from [1, 2) into a smaller binade,
//--- this constructs fp number in [0, 1) as dyadic rational.
   return significand / power2;
  }
//+------------------------------------------------------------------+
//| Uniformly distributed random double in the range [0.0,max)       |
//| max : exclusive and must be greater than zero.                   |
//+------------------------------------------------------------------+
double CRandom::RandomDouble(const double max)
  {
   if(max <= 0)
      return(-1);
//--- return result
   for(;;)
     {
      double r = RandomDouble() * max;
      //--- check for rounding error
      if(r < max)
         return r;
     }
  }
//+------------------------------------------------------------------+
//| Uniformly distributed random double in the range [min,max)       |
//| min : inclusive                                                  |
//| max : exclusive and must be greater than min.                    |
//+------------------------------------------------------------------+
double CRandom::RandomDouble(const double min,const double max)
  {
   if(max <= min)
      return(-1);
//--- return result
   for(;;)
     {
      double r = RandomDouble() * (max - min) + min;
      //--- check for rounding error
      if(r < max)
         return r;
     }
  }
//+------------------------------------------------------------------+
//| Uniformly distributed random true/false with equal probability   |
//+------------------------------------------------------------------+
bool CRandom::RandomBoolean(void)
  {
   return pcg32_random_r() <= INT_MAX;
  }
//+------------------------------------------------------------------+
//| Random true/false where 'true' has the probability of prob_true  |
//+------------------------------------------------------------------+
bool CRandom::RandomBoolean(const double prob_true)
  {
   return pcg32_random_r() <= prob_true * UINT_MAX;
  }
//+------------------------------------------------------------------+
//| Random normal sample with specified mean and standard deviation  |
//+------------------------------------------------------------------+
double CRandom::RandomNormal(const double mu,const double sigma)
  {
   return mu + sigma * RandomNormal();
  }
//+------------------------------------------------------------------+
//| Random normal sample with mean 0 and standard deviation 1        |
//+------------------------------------------------------------------+
double CRandom::RandomNormal(void)
  {
// Use Box-Muller algorithm
   double u1 = RandomDouble();
   double u2 = RandomDouble();
   double r=::MathSqrt(-2.0*::MathLog(u1));
   double theta=2.0*M_PI*u2;
   return r * ::MathSin(theta);
  }
//+------------------------------------------------------------------+
//| PCG Random Number Generation, Minimal C Edition                  |
//|                                                                  |
//|     http://www.pcg-random.org                                    |
//|                                                                  |
//| Copyright 2014 Melissa O'Neill <oneill@pcg-random.org>           |
//+------------------------------------------------------------------+
//+------------------------------------------------------------------+
//| Set the initial values for the generator of random numbers.      |
//| Seed the rng. Specified in two parts, state initializer and a    |
//| sequence selection constant (a.k.a. stream id)                   |
//+------------------------------------------------------------------+
void CRandom::pcg32_srandom_r(const ulong initstate,const ulong initseq)
  {
   m_state = 0;
   m_inc = (initseq << 1) | 1; // Must *always* be odd.
   pcg32_random_r();
   m_state += initstate;
   pcg32_random_r();
  }
//+------------------------------------------------------------------+
//| Uniformly distributed 32-bit unsigned integer [0,UINT_MAX]       |
//| This is the heart of the generator.                              |
//+------------------------------------------------------------------+
uint CRandom::pcg32_random_r(void)
  {
   ulong oldstate = m_state;
   m_state = oldstate * (ulong)6364136223846793005 + m_inc;
   uint xorshifted = (uint)(((oldstate >> 18) ^ oldstate) >> 27);
   uint rot = (uint)(oldstate >> 59);
   return (xorshifted >> rot) | (xorshifted << ((0-rot) & 31));
  }
//+------------------------------------------------------------------+
//| Uniformly distributed random integer, r, where 0 <= r < bound    |
//| Generates random numbers less than a given bound.                |
//+------------------------------------------------------------------+
uint CRandom::pcg32_boundedrand_r(const uint bound)
  {
   // To avoid bias, we need to make the range of the RNG a multiple of
   // bound, which we do by dropping output less than a threshold.
   // A naive scheme to calculate the threshold would be to do
   //
   //     uint32_t threshold = 0x100000000ull % bound;
   //
   // but 64-bit div/mod is slower than 32-bit div/mod (especially on
   // 32-bit platforms).  In essence, we do
   //
   //     uint32_t threshold = (0x100000000ull-bound) % bound;
   //
   // because this version will calculate the same modulus, but the LHS
   // value is less than 2^32.

   uint threshold = (0-bound) % bound;

   // Uniformity guarantees that this loop will terminate.  In practice, it
   // should usually terminate quickly; on average (assuming all bounds are
   // equally likely), 82.25% of the time, we can expect it to require just
   // one iteration.  In the worst case, someone passes a bound of 2^31 + 1
   // (i.e., 2147483649), which invalidates almost 50% of the range.  In
   // practice, bounds are typically small and only a tiny amount of the range
   // is eliminated.
   for(;;)
     {
      uint r = pcg32_random_r();
      if(r >= threshold)
         return r % bound;
     }
  }
//+------------------------------------------------------------------+
//| Uniformly distributed 64-bit unsigned long integer [0,ULONG_MAX] |
//| We build a 64-bit number by tying together two 32-bit numbers.   |
//+------------------------------------------------------------------+
ulong CRandom::pcg32x2_random_r(void)
  {
   return ((ulong)(pcg32_random_r()) << 32)
          | pcg32_random_r();
  }
//+------------------------------------------------------------------+