//+------------------------------------------------------------------+
//|                                                   Xoshiro256.mqh |
//|                                        Copyright © 2025, Amr Ali |
//|                             https://www.mql5.com/en/users/amrali |
//+------------------------------------------------------------------+
#property copyright "Copyright © 2025, Amr Ali"
#property link      "https://www.mql5.com/en/users/amrali"
#property version   "1.09"
#property description "Random number generation using the 64-bit Xoshiro256** algorithm."
#property description "All-purpose rock-solid generator with excellent (sub-ns) speed."

// Updates:
// 2023.01.06 - v.1.01 : Updated the class to use the more robust 64-bits xoshiro256** random number generator, instead of the 32-bits xoshiro128**.
// 2025.08.01 - v.1.02 : Using a splitmix64 generator to seed the state to avoid correlation on similar seeds, as recommended by the original authors.
//                     : Added a seedable constructor with a single ulong. The four 64-bit state variables would be generated using splitmix64 function.
//                     : Increased performance of all the public methods for generating random numbers.
// 2025.08.03 - v.1.03 : Added RandomDoubleHighRes, a new method for generating high-resolution random doubles in the range of [0.0, 1.0).
// 2025.08.03 - v.1.04 : Faster generation of unbiased random integers in a range by using Lemire's fastrange method.
// 2025.08.28 - v.1.05 : Improved RNG seeding mechanism to enforce strong bit-avalanche.
// 2025.09.13 - v.1.06 : Allow richer initialization by using full entropy from all four 64-bit seeds independently.
// 2025.09.16 - v.1.07 : Implemented SHA-256 based seeding (uses CryptEncode with CRYPT_HASH_SHA256) to preserve full 256-bit entropy.
// 2025.09.19 - v.1.08 : Added sampling utilities (Shuffle, Sample, SampleWithReplacement) to the Xoshiro256 class.
// 2025.10.02 - v.1.09 : Implemented a dedicated boundedUInt32 function using the 32-bit variant of Lemire's fastrange to improve performance.

//+------------------------------------------------------------------+
//| Class Xoshiro256                                                 |
//| Usage: Provides an implementation of the Xoshiro256** algorithm. |
//+------------------------------------------------------------------+
class Xoshiro256
  {
protected:
   // RNG internal state (four 64-bit words)
   ulong             s0, s1, s2, s3;

   // For caching a 32-bit integer.
   int               has_uint32;
   uint              uinteger;

   // Internal helper functions
   ulong             mulhi(const ulong x,const ulong y,ulong &lowPart);
   ulong             rotl(const ulong value,const int count);
   uint              nextUInt32(void);
   ulong             nextUInt64(void);
   uint              boundedUInt32(const uint bound);
   ulong             boundedUInt64(const ulong bound);

public:
   // Constructor and destructor
                     Xoshiro256(void);                                   // Auto-seeded constructor
                    ~Xoshiro256(void) {   }

   // Setting the internal state for the generator
   void              Seed(ulong seed);
   void              Seed(ulong a,ulong b,ulong c,ulong d);
   string            GetName() const { return "Xoshiro 256** 64-bit"; }

   // Public methods for generating random numbers
   int               RandomInteger(void);                                // Random integer [0, INT_MAX]
   int               RandomInteger(const int max);                       // Random integer [0, max)
   int               RandomInteger(const int min,const int max);         // Random integer [min, max)

   long              RandomLong(void);                                   // Random long [0, LONG_MAX]
   long              RandomLong(const long max);                         // Random long [0, max)
   long              RandomLong(const long min,const long max);          // Random long [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)

   double            RandomDoubleHighRes(void);                          // Random double [0.0, 1.0)
   double            RandomDoubleHighRes(const double max);              // Random double [0.0, max)
   double            RandomDoubleHighRes(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 with mean 0, std dev 1)
   double            RandomNormal(const double mu,const double sigma);   // Random double (normal distribution)

   // Sampling utilities
   template<typename T>
   void              Shuffle(T &arr[]);                                  // Shuffle array arr[] in place using Fisher-Yates
   template<typename T>
   bool              Sample(const T &arr[], T &dest[], int k);           // Sample k unique items from arr[] without replacement
   template<typename T>
   bool              SampleWithReplacement(const T &arr[], T &dest[], int k);
  };
//+------------------------------------------------------------------+
//| Initializes a new instance of the Xoshiro class with auto-seed.  |
//+------------------------------------------------------------------+
void Xoshiro256::Xoshiro256(void) : has_uint32(0), uinteger(0)
  {
   // collect multiple entropy sources (non-cryptographic)
   ulong a = (ulong)MathRand() ^ (ulong)_RandomSeed;
   ulong b = (ulong)GetTickCount64();
   ulong c = (ulong)ChartID();
   ulong d = (ulong)MQLInfoInteger(MQL_GLOBAL_COUNTER);

   // initialize RNG state
   Seed(a, b, c, d);
  }
//+------------------------------------------------------------------+
//| Sets the internal state for the generator with multiple          |
//| 64-bit seeds.                                                    |
//+------------------------------------------------------------------+
void Xoshiro256::Seed(ulong a, ulong b, ulong c, ulong d)
  {
   uchar nullKey[]={};
   union ULongWords { ulong words[4]; uchar bytes[32]; };

   ULongWords entropy;
   entropy.words[0] = a;
   entropy.words[1] = b;
   entropy.words[2] = c;
   entropy.words[3] = d;

   // compute SHA-256 (32-byte) digest of entropy bytes
   ULongWords digest;
   int result = CryptEncode(CRYPT_HASH_SHA256, entropy.bytes, nullKey, digest.bytes);
   if(result != 32)
     {
      Print("Hashing failed. Error code: ", GetLastError());
      s0 = a; s1 = b; s2 = c; s3 = 0x9E3779B97F4A7C15;  // fallback
      return;
     }

   // fill RNG state from the hash digest
   s0 = digest.words[0];
   s1 = digest.words[1];
   s2 = digest.words[2];
   s3 = digest.words[3];
  }
//+------------------------------------------------------------------+
//| Sets the internal state for the generator with a single          |
//| 64-bit seed (e.g., for reproducibility).                         |
//+------------------------------------------------------------------+
void Xoshiro256::Seed(ulong seed)
  {
   Seed(seed, 0UL, 0UL, 0UL);
  }
//+------------------------------------------------------------------+
//| Uniformly distributed random 32-bit INT, in range [0, INT_MAX]   |
//+------------------------------------------------------------------+
int Xoshiro256::RandomInteger(void)
  {
   return (int)(nextUInt32() >> 1);
  }
//+------------------------------------------------------------------+
//| Uniformly distributed random 32-bit INT, in the range [0, max)   |
//| i.e., inclusive of 0 and exclusive of max.                       |
//+------------------------------------------------------------------+
int Xoshiro256::RandomInteger(const int max)
  {
   return (max > 0) ? (int)boundedUInt32(max) : 0;
  }
//+------------------------------------------------------------------+
//| Uniformly distributed random 32-bit INT, in the range [min, max) |
//| i.e., inclusive of min and exclusive of max.                     |
//+------------------------------------------------------------------+
int Xoshiro256::RandomInteger(const int min,const int max)
  {
   return (max > min) ? (int)boundedUInt32(max - min) + min : min;
  }
//+------------------------------------------------------------------+
//| Uniformly distributed random 64-bit LONG, in range [0, LONG_MAX] |
//+------------------------------------------------------------------+
long Xoshiro256::RandomLong(void)
  {
   return (long)(nextUInt64() >> 1);
  }
//+------------------------------------------------------------------+
//| Uniformly distributed random 64-bit LONG, in the range [0, max)  |
//| i.e., inclusive of 0 and exclusive of max.                       |
//+------------------------------------------------------------------+
long Xoshiro256::RandomLong(const long max)
  {
   return (max > 0) ? (long)boundedUInt64(max) : 0;
  }
//+------------------------------------------------------------------+
//| Uniformly distributed random 64-bit LONG, in the range [min, max)|
//| i.e., inclusive of min and exclusive of max.                     |
//+------------------------------------------------------------------+
long Xoshiro256::RandomLong(const long min,const long max)
  {
   return (max > min) ? (long)boundedUInt64(max - min) + min : min;
  }
//+------------------------------------------------------------------+
//| Uniformly distributed random double in the range [0.0, 1.0)      |
//| i.e., inclusive of 0.0 and exclusive of 1.0.                     |
//| Doubles are rounded down to the nearest multiple of 1/2^53.      |
//+------------------------------------------------------------------+
double Xoshiro256::RandomDouble()
  {
   return (nextUInt64() >> 11) * (1.0 / (1ull << 53));
  }
//+------------------------------------------------------------------+
//| Uniformly distributed random double in the range [0.0, max)      |
//| i.e., inclusive of 0.0 and exclusive of max.                     |
//+------------------------------------------------------------------+
double Xoshiro256::RandomDouble(const double max)
  {
   if(max <= 0) return 0;
   double r;
   do {
      r = RandomDouble() * max;
   } while(r >= max); // Check for rounding error
   return r;
  }
//+------------------------------------------------------------------+
//| Uniformly distributed random double in the range [min, max)      |
//| i.e., inclusive of min and exclusive of max.                     |
//+------------------------------------------------------------------+
double Xoshiro256::RandomDouble(const double min,const double max)
  {
   if(max <= min) return min;
   double r;
   double range = max - min;
   do {
      r = RandomDouble() * range + min;
   } while(r >= max);
   return r;
  }
//+------------------------------------------------------------------+
//| Uniformly distributed random double in the range [0.0, 1.0)      |
//| i.e., inclusive of 0.0 and exclusive of 1.0.                     |
//|                                                                  |
//| https://docs.python.org/3/library/random.html#recipe             |
//| http://mumble.net/~campbell/tmp/random_real.c                    |
//|                                                                  |
//| Uses an alternative precise sampling method that is capable of   |
//| generating all possible values in the interval [0, 1) that can   |
//| be represented by a double precision float, a feature important  |
//| for scientific simulation and other high-precision applications. |
//| Note however that this method is slower than RandomDouble().     |
//+------------------------------------------------------------------+
double Xoshiro256::RandomDoubleHighRes()
  {
   ulong mantissa = 0x10000000000000 | nextUInt64() >> 12; // 2^52...2^53-1
   double power2  = 0x20000000000000; // 2^53
   uint x;
   while((x = nextUInt32()) == 0)
     {
      power2 *= 0x100000000; // 2^32
     }
   // bitmask for rightmost 1-bit (geometric distribution)
   x &= (0u-x);
   power2 *= x;
   return mantissa / power2;
  }
//+------------------------------------------------------------------+
//| Uniformly distributed random double in the range [0.0, max)      |
//| i.e., inclusive of 0.0 and exclusive of max.                     |
//+------------------------------------------------------------------+
double Xoshiro256::RandomDoubleHighRes(const double max)
  {
   if(max <= 0) return 0;
   double r;
   do {
      r = RandomDoubleHighRes() * max;
   } while(r >= max);
   return r;
  }
//+------------------------------------------------------------------+
//| Uniformly distributed random double in the range [min, max)      |
//| i.e., inclusive of min and exclusive of max.                     |
//+------------------------------------------------------------------+
double Xoshiro256::RandomDoubleHighRes(const double min,const double max)
  {
   if(max <= min) return min;
   double r;
   double range = max - min;
   do {
      r = RandomDoubleHighRes() * range + min;
   } while(r >= max);
   return r;
  }
//+------------------------------------------------------------------+
//| Uniformly distributed random true/false with equal probability   |
//+------------------------------------------------------------------+
bool Xoshiro256::RandomBoolean(void)
  {
   // Use a high bit since the low bits are less random.
   return (nextUInt32() & 0x80000000) != 0;
  }
//+------------------------------------------------------------------+
//| Random true/false where 'true' has the probability of prob_true  |
//+------------------------------------------------------------------+
bool Xoshiro256::RandomBoolean(const double prob_true)
  {
   if (prob_true < 0.0) return false;
   if (prob_true >= 1.0) return true;
   return RandomDouble() < prob_true;
  }
//+------------------------------------------------------------------+
//| Random normal sample with specified mean and standard deviation  |
//+------------------------------------------------------------------+
double Xoshiro256::RandomNormal(const double mu,const double sigma)
  {
   return mu + sigma * RandomNormal();
  }
//+------------------------------------------------------------------+
//| Random normal sample with mean 0 and standard deviation 1        |
//+------------------------------------------------------------------+
double Xoshiro256::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);
  }
//+------------------------------------------------------------------+
//| Shuffle the order of array elements in-place.                    |
//| To shuffle to a new array, use Sample(a, dest, ArraySize(a)).    |
//+------------------------------------------------------------------+
template<typename T>
void Xoshiro256::Shuffle(T &arr[])
  {
   // Fisher-Yates shuffle
   int size = ArraySize(arr);
   for(int i = 0; i < size - 1; i++)
     {
      // Select random index and swap items [j] and [i]
      int j = RandomInteger(i, size);
      T tmp = arr[i];
      arr[i] = arr[j];
      arr[j] = tmp;
     }
  }
//+------------------------------------------------------------------+
//| Sample k unique items from array (without replacement)           |
//| k cannot be larger than the size of arr[].                       |
//+------------------------------------------------------------------+
template<typename T>
bool Xoshiro256::Sample(const T &arr[], T &dest[], int k)
  {
   int size = ArraySize(arr);
   if(size == 0 || k <= 0) return false;
   if(k > size) k = size; // cannot sample more than available
   if(ArrayResize(dest, k) != k) return false;

   // Make a temporary index array
   int idx[];
   if(ArrayResize(idx, size) != size) return false;
   for(int i=0; i<size; i++) idx[i] = i;

   // Select the first k random unique indices
   for(int i=0; i<k; i++)
     {
      // Fisher-Yates shuffle on indices
      int j = RandomInteger(i, size);
      int tmp = idx[i];
      idx[i] = idx[j];
      idx[j] = tmp;
      dest[i] = arr[idx[i]];
     }
   return true;
  }
//+------------------------------------------------------------------+
//| Sample k items from array with replacement.                      |
//| A sampled item is replaced back into arr[], so that the same     |
//| item can be selected more than once (duplicate selections).      |
//+------------------------------------------------------------------+
template<typename T>
bool Xoshiro256::SampleWithReplacement(const T &arr[], T &dest[], int k)
  {
   int size = ArraySize(arr);
   if(size == 0 || k <= 0) return false;
   if(ArrayResize(dest, k) != k) return false;

   // Pick a random index from arr[]
   for(int i=0; i<k; i++)
     {
      int ind = RandomInteger(0, size);
      dest[i] = arr[ind];
     }
   return true;
  }
//+------------------------------------------------------------------+
//| Copyright 2018 by David Blackman and Sebastiano Vigna            |
//|                                                                  |
//|     https://prng.di.unimi.it/xoshiro256starstar.c                |
//|                                                                  |
//| This is xoshiro256** 1.0, one of our all-purpose, rock-solid     |
//| generators. It has excellent (sub-ns) speed, a state (256 bits)  |
//| that is large enough for any parallel application, and it        |
//| passes all tests we are aware of.                                |
//|                                                                  |
//| The state must be seeded so that it is not everywhere zero.      |
//| If you have a 64-bit seed, we suggest to seed a splitmix64       |
//| generator and use its output to fill s.                          |
//+------------------------------------------------------------------+

//+------------------------------------------------------------------+
//| Emulate 64x64 bit multiplication to 128 bit result.              |
//| Returns high part of product, and low part is set in `lowPart`.  |
//| github.com/numpy/numpy/blob/main/numpy/random/src/pcg64/pcg64.h  |
//+------------------------------------------------------------------+
ulong Xoshiro256::mulhi(const ulong x,const ulong y,ulong &lowPart)
  {
   lowPart = x * y; // Lower 64 bits are straightforward clock-arithmetic.

   ulong x0 = x & 0xFFFFFFFF;
   ulong x1 = x >> 32;
   ulong y0 = y & 0xFFFFFFFF;
   ulong y1 = y >> 32;
   ulong w0 = x0 * y0;
   ulong t = x1 * y0 + (w0 >> 32);
   ulong w1 = t & 0xFFFFFFFF;
   ulong w2 = t >> 32;
   w1 += x0 * y1;
   ulong hi = x1 * y1 + w2 + (w1 >> 32);

   return hi;
  }
//+------------------------------------------------------------------+
//| Rotates bits of a 64-bit value to the left by count positions.   |
//+------------------------------------------------------------------+
ulong Xoshiro256::rotl(const ulong value,const int count)
  {
   return (value << count) | (value >> (64 - count));
  }
//+------------------------------------------------------------------+
//| Uniformly distributed 64-bit unsigned long integer [0,ULONG_MAX] |
//| This is the heart of the generator.                              |
//+------------------------------------------------------------------+
ulong Xoshiro256::nextUInt64(void)
  {
   const ulong result = rotl(s1 * 5, 7) * 9;

   // Update the generator state
   const ulong t = s1 << 17;
   s2 ^= s0;
   s3 ^= s1;
   s1 ^= s2;
   s0 ^= s3;
   s2 ^= t;
   s3 = rotl(s3, 45);

   return result;
  }
//+------------------------------------------------------------------+
//| Uniformly distributed 32-bit unsigned integer [0, UINT_MAX]      |
//+------------------------------------------------------------------+
uint Xoshiro256::nextUInt32(void)
  {
   ulong next;
   if(has_uint32)
     {
      has_uint32 = 0;
      return uinteger;
     }
   next = nextUInt64();
   has_uint32 = 1;
   uinteger = (uint)(next >> 32);
   return (uint)(next & 0xffffffff);
  }
//+------------------------------------------------------------------------+
//| Uniformly distributed random 64-bit ULONG, r, where 0 <= r < bound     |
//| This implementation uses a fast rejection method (Lemire's fastrange). |
//| Source: https://github.com/lemire/fastrange                            |
//+------------------------------------------------------------------------+
ulong Xoshiro256::boundedUInt64(const ulong bound)
  {
   ulong l;
   ulong m = mulhi(nextUInt64(), bound, l);
   if(l < bound)
     {
      ulong t = (0ull-bound) % bound;
      while(l < t)
        {
         m = mulhi(nextUInt64(), bound, l);
        }
     }
   return m;
  }
//+------------------------------------------------------------------------+
//| Uniformly distributed random 32-bit UINT, r, where 0 <= r < bound      |
//| 32-bit variant of Lemire's fastrange.                                  |
//+------------------------------------------------------------------------+
uint Xoshiro256::boundedUInt32(const uint bound)
  {
   ulong m = nextUInt32() * (ulong)bound;
   uint l = (uint) m;
   if(l < bound)
     {
      uint t = (0u-bound) % bound;
      while(l < t)
        {
         m = nextUInt32() * (ulong)bound;
         l = (uint) m;
        }
     }
   return (uint) (m >> 32);
  }