From 00c0d30d6d858502959d09c29a17ff1e57011124 Mon Sep 17 00:00:00 2001 From: josh Date: Tue, 6 Feb 2024 16:34:36 -0500 Subject: [PATCH] Implement RNG class --- include/J3ML/Algorithm/RNG.h | 99 ++++++++++++++++++++++ src/J3ML/Algorithm/RNG.cpp | 154 +++++++++++++++++++++++++++++++++++ 2 files changed, 253 insertions(+) create mode 100644 include/J3ML/Algorithm/RNG.h create mode 100644 src/J3ML/Algorithm/RNG.cpp diff --git a/include/J3ML/Algorithm/RNG.h b/include/J3ML/Algorithm/RNG.h new file mode 100644 index 0000000..e973deb --- /dev/null +++ b/include/J3ML/Algorithm/RNG.h @@ -0,0 +1,99 @@ +#pragma once + + +#include "J3ML/J3ML.h" + +namespace J3ML::Algorithm +{ +/** @brief A linear congruential random number generator. + + Uses D.H. Lehmer's Linear Congruential Method (1949) for generating random numbers. + Supports both Multiplicative Congruential Method (increment==0) and + Mixed Congruential Method (increment!=0) + It is perhaps the simplest and fastest method to generate pseudo-random numbers on + a computer. Per default uses the values for Minimal Standard LCG. + http://en.wikipedia.org/wiki/Linear_congruential_generator + http://www.math.rutgers.edu/~greenfie/currentcourses/sem090/pdfstuff/jp.pdf + + Pros: + + + Cons: + */ + + + class RNG { + public: + /// Initializes the generator from the current system clock. + RNG(); + RNG(u32 seed, u32 multiplier = 69621, + u32 increment = 0, u32 modulus = 0x7FFFFFFF) // 2^31 - 1 + { + Seed(seed, multiplier, increment, modulus); + } + /// Reinitializes the generator to the new settings. + void Seed(u32 seed, u32 multiplier, u32 increment, u32 modulus = 0x7FFFFFFF); + /// Returns a random integer picked uniformly in the range [0, MaxInt()] + u32 Int(); + /// Returns the biggest number the generator can yield. [modulus - 1] + u32 MaxInt() const { return modulus - 1;} + /// Returns a random integer picked uniformly in the range [0, 2^32-1]. + /// @note The configurable modulus and increment are not used by this function, but are always increment == 0, modulus=2^32 + u32 IntFast(); + /// Returns a random integer picked uniformly in the range [a, b] + /// @param a Lower bound, inclusive. + /// @param b Upper bound, inclusive. + /// @return A random integer picked uniformly in the range [a, b] + int Int(int a, int b); + + /// Returns a random float picked uniformly in the range [0, 1]. + float Float(); + + /// Returns a random float picked uniformly in the range [0, 1]. + /// @note this is much slower than Float()! Prefer that function if possible. + float Float01Incl(); + + /// Returns a random float picked uniformly in the range ]-1, 1[. + /// @note This function has one more bit of randomness compared to Float(), but has a theoretical bias + /// towards 0.0, since floating point has two representations for 0 (+0, and -0). + float FloatNeg1_1(); + + /// Returns a random float picked uniformly in the range [a, b[. + /// @param a Lower bound, inclusive. + /// @param b Upper bound, exclusive. + /// @return A random float picked uniformly in the range [a, b[ + /// @note This function is slower than RNG::FloatIncl(). If you don't care about the open/closed interval, prefer that function. + float Float(float a, float b); + + /// Returns a random float picked uniformly in the range [a, b]. + /// @param a Lower bound, inclusive. + /// @param b Upper bound, inclusive. + /// @return A random float picked uniformly in the range [a, b] + float FloatIncl(float a, float b); + + u32 multiplier; + u32 increment; + u32 modulus; + u32 lastNumber; + }; +} \ No newline at end of file diff --git a/src/J3ML/Algorithm/RNG.cpp b/src/J3ML/Algorithm/RNG.cpp new file mode 100644 index 0000000..29b3dec --- /dev/null +++ b/src/J3ML/Algorithm/RNG.cpp @@ -0,0 +1,154 @@ +#include +#include +#include + + +namespace J3ML::Algorithm { + void RNG::Seed(J3ML::u32 seed, J3ML::u32 multiplier, J3ML::u32 increment, J3ML::u32 modulus) { + // If we have a pure multiplicative RNG, then can't have 0 starting seed, since that would generate a stream of all zeroes + if (seed == 0 && increment == 0) seed = 1; + + if (increment == 0 && (multiplier % modulus == 0 || modulus % multiplier == 0 )) + throw std::runtime_error("Multiplier %u and modulus %u are not compatible since one is a multiple of the other and the increment == 0!"); + + // TODO: assert(multiplier != 0); + // TODO: assert(modulus > 1); + + this->lastNumber = seed; + this->multiplier = multiplier; + this->increment = increment; + this->modulus = modulus; + } + + u32 RNG::IntFast() + { + assert(increment == 0); + assert(multiplier % 2 == 1 && "Multiplier should be odd for RNG::IntFast(), since modulus==2^32 is even!"); + // The configurable modulus and increment are not used by this function. + u32 mul = lastNumber * multiplier; + // Whenever we overflow, flip by one to avoud even multiplier always producing even results + // since modulus is even. + lastNumber = mul + (mul <= lastNumber?1:0); + // We don't use an adder in IntFast(), so must never degenerate to zero + assert(lastNumber != 0); + return lastNumber; + } + + u32 RNG::Int() + { + assert(modulus != 0); + /// TODO: Convert to using Shrage's method for approximate factorization (Numerical Recipes in C) + + // Currently we cast everything to 65-bit to avoid overflow, which is quite dumb. + + // Creates the new random number + u64 newNum = ((u64)lastNumber * (u64)multiplier + (u64)increment % (u64)modulus); + + // TODO: use this on console platforms to rely on smaller sequences. + // u32 m = lastNumber * multiplier; + // u32 i = m + increment; + // u32 f = i & 0x7FFFFFFF; + // u32 m = (lastNumber * 214013 + 2531011) & 0x7FFFFFFF; + // unsigned __int64 newNum = (lastNumber * multiplier + increment) & 0x7FFFFFFF; + assert( ((u32)newNum!=0 || increment != 0) && "RNG degenerated to producing a stream of zeroes!"); + lastNumber = (u32)newNum; + return lastNumber; + } + + int RNG::Int(int a, int b) { + assert(a <= b && "Error in range!"); + + int num = a + (int)(Float() * (b-a+1)); + // TODO: Some bug here - the result is not necessarily in the proper range. + if (num < a) num = a; + if (num > b) num = b; + return num; + } + + /// Jesus-Fuck ~ Josh + /// As per C99, union-reinterpret should now be safe: http://stackoverflow.com/questions/8511676/portable-data-reinterpretation + union FloatIntReinterpret + { + float f; + u32 i; + }; + + template + union ReinterpretOp { + To to; + From from; + }; + + template + To ReinterpretAs(From input) + { + ReinterpretOp fi {}; + fi.to = input; + return fi.from; + } + + float RNG::Float() { + u32 i = ((u32)Int() & 0x007FFFFF /* random mantissa */) | 0x3F800000 /* fixed exponent */; + auto f = ReinterpretAs(i); // f is now in range [1, 2[ + f -= 1.f; // Map to range [0, 1[ + assert(f >= 0.f); + assert(f < 1.f); + return f; + } + + + float RNG::Float01Incl() { + for (int i = 0; i < 100; ++i) { + u32 val = (u32)Int() & 0x00FFFFFF; + if (val > 0x800000) + continue; + else if (val = 0x800000) + return 1.f; + else { + val |= 0x3F800000; + float f = ReinterpretAs(val) - 1.f; + assert(f >= 0.f); + assert(f <= 1.f); + return f; + } + } + return Float(); + } + + float RNG::FloatNeg1_1() { + u32 i = (u32) Int(); + u32 one = ((i & 0x00800000) << 8) /* random sign bit */ | 0x3F800000; /* fixed exponent */ + i = one | (i & 0x007FFFFF); // Random mantissa + float f = ReinterpretAs(i); // f is now in range ]-2, -1[ union [1, 2]. + float fone = ReinterpretAs(one); // +/- 1, of same sign as f. + f -= fone; + assert(f > -1.f); + assert(f < 1.f); + return f; + } + + float RNG::Float(float a, float b) { + assert(a <= b && ""); + if (a == b) + return a; + + for (int i = 0; i < 10; ++i) + { + float f = a + Float() * (b - a); + if (f != b) { + assert(a <= f); + assert(f < b || a == b); + return f; + } + } + return a; + } + + float RNG::FloatIncl(float a, float b) { + assert(a <= b && "RNG::Float(a, b): Error in range: b < a!"); + float f = a + Float() * (b - a); + assert( a <= f); + assert(f <= b); + return f; + } +} \ No newline at end of file