Implement RNG class
This commit is contained in:
99
include/J3ML/Algorithm/RNG.h
Normal file
99
include/J3ML/Algorithm/RNG.h
Normal file
@@ -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:
|
||||
<ul>
|
||||
<li> Easy to implement.
|
||||
<li> Fast.
|
||||
</ul>
|
||||
|
||||
Cons:
|
||||
<ul>
|
||||
<li> NOT safe for cryptography because of the easily calculatable sequential
|
||||
correlation between successive calls. A case study:
|
||||
http://www.cigital.com/papers/download/developer_gambling.php
|
||||
|
||||
<li> Tends to have less random low-order bits (compared to the high-order bits)
|
||||
Thus, NEVER do something like this:
|
||||
|
||||
u32 numBetween1And10 = 1 + LCGRand.Int() % 10;
|
||||
|
||||
Instead, take into account EVERY bit of the generated number, like this:
|
||||
|
||||
u32 numBetween1And10 = 1 + (int)(10.0 * (double)LCGRand.Int()
|
||||
/(LCGRand.Max()+1.0));
|
||||
or simply
|
||||
|
||||
u32 numBetween1And10 = LCGRand.Float(1.f, 10.f);
|
||||
</ul> */
|
||||
|
||||
|
||||
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;
|
||||
};
|
||||
}
|
||||
154
src/J3ML/Algorithm/RNG.cpp
Normal file
154
src/J3ML/Algorithm/RNG.cpp
Normal file
@@ -0,0 +1,154 @@
|
||||
#include <J3ML/Algorithm/RNG.h>
|
||||
#include <stdexcept>
|
||||
#include <cassert>
|
||||
|
||||
|
||||
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 <typename To, typename From>
|
||||
union ReinterpretOp {
|
||||
To to;
|
||||
From from;
|
||||
};
|
||||
|
||||
template <typename To, typename From>
|
||||
To ReinterpretAs(From input)
|
||||
{
|
||||
ReinterpretOp<To, From> fi {};
|
||||
fi.to = input;
|
||||
return fi.from;
|
||||
}
|
||||
|
||||
float RNG::Float() {
|
||||
u32 i = ((u32)Int() & 0x007FFFFF /* random mantissa */) | 0x3F800000 /* fixed exponent */;
|
||||
auto f = ReinterpretAs<float, u32>(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<float, u32>(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<float, u32>(i); // f is now in range ]-2, -1[ union [1, 2].
|
||||
float fone = ReinterpretAs<float, u32>(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;
|
||||
}
|
||||
}
|
||||
Reference in New Issue
Block a user