Random.hpp

#include "Common.hpp"

namespace ax
{
	namespace Random
	{
		class IRandom
		{
		public:
			virtual uint32 Next() = 0;
			float NextFloat(float min, float max)
			{
				return min + (NextFloat01() * fabs(min - max));
			}
			float NextFloat01()
			{
				constexpr float c_FMul = (1.0 / 16777216.0f);
				return float(Next() >> 8) * c_FMul;
			}
			double NextDouble01()
			{
				constexpr double c_DMul = (1.0 / 4294967296.0);
				return Next() * c_DMul;
			}
			int NextInt(int min, int max)
			{
				return min + Range(abs(min) + abs(max));
			}
			uint64 NextUlong()
			{
				return uint64(Next()) | (uint64(Next()) << 32ul);
			}
			uint64 NextUlong(uint64 min, uint64 max)
			{
				return min + (NextUlong() % (labs(min) + labs(max)));
			}
			/* returns a random number 0 <= x < bound */
			uint32 Range(uint32 bound)
			{
				uint32 threshold = -bound % bound;
				for (;;) {
					uint32 r = Next();
					if (r >= threshold) return r % bound;
				}
			}
		};
		
		class IRandom64
		{
		public:
			virtual uint64 Next() = 0;
		
			float NextFloat(float min, float max)
			{
				return min + (NextFloat01() * fabs(min - max));
			}
			float NextFloat01()
			{
				return float(NextDouble01());
			}
			double NextDouble01()
			{
				return double(Next()) / 0xffffffffu;
			}
			double NextDouble01(double min, double max)
			{
				return min + (NextDouble01() * abs(min - max));
			}
			double NextDouble()
			{
				uint64_t u64 = 0x3FF0000000000000ULL | ((Next() >> 12) | 1);
				return *(double*)&u64 - 1.0;
			}
			long NextLong(long min, long max)
			{
				return min + Range(abs(min) + abs(max));
			}
			/* returns a random number 0 <= x < bound */
			uint64 Range(uint64 bound)
			{
				return Next() % bound;
			}
		};
		
		// https://www.pcg-random.org/index.html
		// we can also add global state in a cpp file
		// compared to MT chace friendly 
		struct PCG : public IRandom
		{
			uint64 state = 0x853c49e6748fea9bULL;
			uint64 inc = 0xda3e39cb94b95bdbULL;
		public:
			PCG() {}
		
			PCG(uint64 seed)
			{
				state = 0x853c49e6748fea9bULL;
				inc = seed << 1 | 1u;
			}
		
			PCG(uint64 initstate, uint64 seed)
			{
				state = 0U;
				inc = (seed << 1u) | 1u;
				Next();
				state += initstate;
				Next();
			}
		
			// returns a random number 0 <= x < IntMax
			uint32 Next()
			{
				uint64 oldstate = state;
				// Advance internal state
				state = oldstate * 6364136223846793005ULL + (inc | 1);
				// Calculate output function (XSH RR), uses old state for max ILP
				uint32 xorshifted = ((oldstate >> 18u) ^ oldstate) >> 27u;
				uint32 rot = oldstate >> 59u;
				return (xorshifted >> rot) | (xorshifted << ((-rot) & 31));
			}
		};
		
		
		// Copyright (c) 2011, 2013 Mutsuo Saito, Makoto Matsumoto,
		// Hiroshima University and The University of Tokyo.
		// All rights reserved.
		// original implementation of mersene twister algorithm
		// generated from paper: http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/mt.pdf
		
		// important: you need to create MT classes for per thread
		// 	also I don't recommend using more than one instance in a therad
		
		class MTwister64 : public IRandom64
		{
			static constexpr int N = 624, M = 367;
			static constexpr int DIFF = N - M;
		
			uint64 m_MT[N];
			int mti = N + 1;
		
		public:
			// any non zero integer can be used as a seed
			MTwister64(uint64 seed = 4357ul)
			{
				m_MT[0] = seed & ~0ul;
				for (mti = 1; mti < N; ++mti)
					m_MT[mti] = (69069 * m_MT[mti - 1]) & ~0ul;
			}
		
			uint64 Next()
			{
				if (mti >= N) GenerateNumbers();
				uint64 x = m_MT[mti++];
				x ^= x >> 11;
				x ^= x << 7 & 0x9d2c5680ul;
				x ^= x << 15 & 0xefc60000ul;
				x ^= x >> 18;
				return x;
			}
		
		private:
		
			void GenerateNumbers()
			{
				static const uint64 mag01[2] = { 0x0, 0x9908b0dful };
		
				int kk = 0;
				uint64 y;
		
				while (kk < N - M) // unrolled for chace line optimizations
				{
					y = (m_MT[kk] & 0x80000000ul) | (m_MT[kk + 1] & 0x7ffffffful);
					m_MT[kk] = m_MT[kk + M] ^ (y >> 1) ^ mag01[y & 0x1];
					kk++;
					y = (m_MT[kk] & 0x80000000ul) | (m_MT[kk + 1] & 0x7ffffffful);
					m_MT[kk] = m_MT[kk + M] ^ (y >> 1) ^ mag01[y & 0x1];
					kk++;
					y = (m_MT[kk] & 0x80000000ul) | (m_MT[kk + 1] & 0x7ffffffful);
					m_MT[kk] = m_MT[kk + M] ^ (y >> 1) ^ mag01[y & 0x1];
					kk++;
				}
		
				while (kk < N - 1)
				{
					y = (m_MT[kk] & 0x80000000ul) | (m_MT[kk + 1] & 0x7ffffffful);
					m_MT[kk] = m_MT[kk + (M - N)] ^ (y >> 1) ^ mag01[y & 0x1];
					++kk;
					y = (m_MT[kk] & 0x80000000ul) | (m_MT[kk + 1] & 0x7ffffffful);
					m_MT[kk] = m_MT[kk + (M - N)] ^ (y >> 1) ^ mag01[y & 0x1];
					++kk;
					y = (m_MT[kk] & 0x80000000ul) | (m_MT[kk + 1] & 0x7ffffffful);
					m_MT[kk] = m_MT[kk + (M - N)] ^ (y >> 1) ^ mag01[y & 0x1];
					++kk;
				}
		
				y = (m_MT[N - 1] & 0x80000000ul) | (m_MT[0] & 0x7ffffffful);
				m_MT[N - 1] = m_MT[kk + (M - 1)] ^ (y >> 1) ^ mag01[y & 0x1];
		
				mti = 0;
			}
		};
		
		class MTwister : public IRandom
		{
			static constexpr int SIZE = 624, PERIOD = 397;
			static constexpr int DIFF = SIZE - PERIOD;
		
			static constexpr uint32 MAGIC = 0x9908b0df;
		
			uint32 MT[SIZE];
			int index = SIZE;
		
		public:
		
			MTwister(uint32 value = 4586u)
			{
				MT[0] = value;
				index = SIZE;
		
				for (uint32 i = 1; i < SIZE; ++i)
					MT[i] = 0x6c078965 * (MT[i - 1] ^ MT[i - 1] >> 30) + i;
			}
		
			uint32 Next()
			{
				if (index == SIZE) {
					GenerateNumbers();
					index = 0;
				}
				uint32 y = MT[index++];
				y ^= y >> 11;
				y ^= y << 7 & 0x9d2c5680;
				y ^= y << 15 & 0xefc60000;
				y ^= y >> 18;
				return y;
			}
		
		private:
			void GenerateNumbers()
			{
				size_t i = 0;
				uint32_t y;
		
				while (i < DIFF) {
					y = (0x80000000 & MT[i]) | (0x7FFFFFFF & MT[i + 1]);
					MT[i] = MT[i + PERIOD] ^ (y >> 1) ^ (((int(y) << 31) >> 31) & MAGIC);
					++i;
					y = (0x80000000 & MT[i]) | (0x7FFFFFFF & MT[i + 1]);
					MT[i] = MT[i + PERIOD] ^ (y >> 1) ^ (((int(y) << 31) >> 31) & MAGIC);
					++i;
				}
		
				while (i < SIZE - 1) 
				{
					y = (0x80000000 & MT[i]) | (0x7FFFFFFF & MT[i + 1]);
					MT[i] = MT[i - DIFF] ^ (y >> 1) ^ (((int(y) << 31) >> 31) & MAGIC); 
					++i;
					y = (0x80000000 & MT[i]) | (0x7FFFFFFF & MT[i + 1]);
					MT[i] = MT[i - DIFF] ^ (y >> 1) ^ (((int(y) << 31) >> 31) & MAGIC);
					++i;
					y = (0x80000000 & MT[i]) | (0x7FFFFFFF & MT[i + 1]);
					MT[i] = MT[i - DIFF] ^ (y >> 1) ^ (((int(y) << 31) >> 31) & MAGIC);
					++i;
				}
		
				// i = 623, last step rolls over
				y = (0x80000000 & MT[SIZE - 1]) | (0x7FFFFFFF & MT[0]);
				MT[SIZE - 1] = MT[PERIOD - 1] ^ (y >> 1) ^ (((int32_t(y) << 31) >> 31) & MAGIC);
		
				index = 0;
			}
		};

	} // namespace Random end
}