godot/thirdparty/thekla_atlas/nvmath/Random.h

377 lines
9.4 KiB
C++
Raw Normal View History

// This code is in the public domain -- castanyo@yahoo.es
#pragma once
#ifndef NV_MATH_RANDOM_H
#define NV_MATH_RANDOM_H
#include "nvmath.h"
#include "nvcore/Utils.h" // nextPowerOfTwo
namespace nv
{
/// Interface of the random number generators.
class Rand
{
public:
virtual ~Rand() {}
enum time_e { Time };
/// Provide a new seed.
virtual void seed( uint s ) { /* empty */ };
/// Get an integer random number.
virtual uint get() = 0;
/// Get a random number on [0, max] interval.
uint getRange( uint max )
{
if (max == 0) return 0;
if (max == NV_UINT32_MAX) return get();
const uint np2 = nextPowerOfTwo( max+1 ); // @@ This fails if max == NV_UINT32_MAX
const uint mask = np2 - 1;
uint n;
do { n = get() & mask; } while( n > max );
return n;
}
/// Random number on [0.0, 1.0] interval.
float getFloat()
{
union
{
uint32 i;
float f;
} pun;
pun.i = 0x3f800000UL | (get() & 0x007fffffUL);
return pun.f - 1.0f;
}
float getFloatRange(float min, float max) {
return getFloat() * (max - min) + min;
}
/*
/// Random number on [0.0, 1.0] interval.
double getReal()
{
return double(get()) * (1.0/4294967295.0); // 2^32-1
}
/// Random number on [0.0, 1.0) interval.
double getRealExclusive()
{
return double(get()) * (1.0/4294967296.0); // 2^32
}
*/
/// Get the max value of the random number.
uint max() const { return NV_UINT32_MAX; }
// Get a random seed.
static uint randomSeed();
};
/// Very simple random number generator with low storage requirements.
class SimpleRand : public Rand
{
public:
/// Constructor that uses the current time as the seed.
SimpleRand( time_e )
{
seed(randomSeed());
}
/// Constructor that uses the given seed.
SimpleRand( uint s = 0 )
{
seed(s);
}
/// Set the given seed.
virtual void seed( uint s )
{
current = s;
}
/// Get a random number.
virtual uint get()
{
return current = current * 1103515245 + 12345;
}
private:
uint current;
};
/// Mersenne twister random number generator.
class MTRand : public Rand
{
public:
enum { N = 624 }; // length of state vector
enum { M = 397 };
/// Constructor that uses the current time as the seed.
MTRand( time_e )
{
seed(randomSeed());
}
/// Constructor that uses the given seed.
MTRand( uint s = 0 )
{
seed(s);
}
/// Constructor that uses the given seeds.
NVMATH_API MTRand( const uint * seed_array, uint length );
/// Provide a new seed.
virtual void seed( uint s )
{
initialize(s);
reload();
}
/// Get a random number between 0 - 65536.
virtual uint get()
{
// Pull a 32-bit integer from the generator state
// Every other access function simply transforms the numbers extracted here
if( left == 0 ) {
reload();
}
left--;
uint s1;
s1 = *next++;
s1 ^= (s1 >> 11);
s1 ^= (s1 << 7) & 0x9d2c5680U;
s1 ^= (s1 << 15) & 0xefc60000U;
return ( s1 ^ (s1 >> 18) );
};
private:
NVMATH_API void initialize( uint32 seed );
NVMATH_API void reload();
uint hiBit( uint u ) const { return u & 0x80000000U; }
uint loBit( uint u ) const { return u & 0x00000001U; }
uint loBits( uint u ) const { return u & 0x7fffffffU; }
uint mixBits( uint u, uint v ) const { return hiBit(u) | loBits(v); }
uint twist( uint m, uint s0, uint s1 ) const { return m ^ (mixBits(s0,s1)>>1) ^ ((~loBit(s1)+1) & 0x9908b0dfU); }
private:
uint state[N]; // internal state
uint * next; // next value to get from state
int left; // number of values left before reload needed
};
/** George Marsaglia's random number generator.
* Code based on Thatcher Ulrich public domain source code:
* http://cvs.sourceforge.net/viewcvs.py/tu-testbed/tu-testbed/base/tu_random.cpp?rev=1.7&view=auto
*
* PRNG code adapted from the complimentary-multiply-with-carry
* code in the article: George Marsaglia, "Seeds for Random Number
* Generators", Communications of the ACM, May 2003, Vol 46 No 5,
* pp90-93.
*
* The article says:
*
* "Any one of the choices for seed table size and multiplier will
* provide a RNG that has passed extensive tests of randomness,
* particularly those in [3], yet is simple and fast --
* approximately 30 million random 32-bit integers per second on a
* 850MHz PC. The period is a*b^n, where a is the multiplier, n
* the size of the seed table and b=2^32-1. (a is chosen so that
* b is a primitive root of the prime a*b^n + 1.)"
*
* [3] Marsaglia, G., Zaman, A., and Tsang, W. Toward a universal
* random number generator. _Statistics and Probability Letters
* 8_ (1990), 35-39.
*/
class GMRand : public Rand
{
public:
enum { SEED_COUNT = 8 };
// const uint64 a = 123471786; // for SEED_COUNT=1024
// const uint64 a = 123554632; // for SEED_COUNT=512
// const uint64 a = 8001634; // for SEED_COUNT=255
// const uint64 a = 8007626; // for SEED_COUNT=128
// const uint64 a = 647535442; // for SEED_COUNT=64
// const uint64 a = 547416522; // for SEED_COUNT=32
// const uint64 a = 487198574; // for SEED_COUNT=16
// const uint64 a = 716514398U; // for SEED_COUNT=8
enum { a = 716514398U };
GMRand( time_e )
{
seed(randomSeed());
}
GMRand(uint s = 987654321)
{
seed(s);
}
/// Provide a new seed.
virtual void seed( uint s )
{
c = 362436;
i = SEED_COUNT - 1;
for(int i = 0; i < SEED_COUNT; i++) {
s = s ^ (s << 13);
s = s ^ (s >> 17);
s = s ^ (s << 5);
Q[i] = s;
}
}
/// Get a random number between 0 - 65536.
virtual uint get()
{
const uint32 r = 0xFFFFFFFE;
uint64 t;
uint32 x;
i = (i + 1) & (SEED_COUNT - 1);
t = a * Q[i] + c;
c = uint32(t >> 32);
x = uint32(t + c);
if( x < c ) {
x++;
c++;
}
uint32 val = r - x;
Q[i] = val;
return val;
};
private:
uint32 c;
uint32 i;
uint32 Q[8];
};
/** Random number implementation from the GNU Sci. Lib. (GSL).
* Adapted from Nicholas Chapman version:
*
* Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
* This is the Unix rand48() generator. The generator returns the
* upper 32 bits from each term of the sequence,
*
* x_{n+1} = (a x_n + c) mod m
*
* using 48-bit unsigned arithmetic, with a = 0x5DEECE66D , c = 0xB
* and m = 2^48. The seed specifies the upper 32 bits of the initial
* value, x_1, with the lower 16 bits set to 0x330E.
*
* The theoretical value of x_{10001} is 244131582646046.
*
* The period of this generator is ? FIXME (probably around 2^48).
*/
class Rand48 : public Rand
{
public:
Rand48( time_e )
{
seed(randomSeed());
}
Rand48( uint s = 0x1234ABCD )
{
seed(s);
}
/** Set the given seed. */
virtual void seed( uint s ) {
vstate.x0 = 0x330E;
vstate.x1 = uint16(s & 0xFFFF);
vstate.x2 = uint16((s >> 16) & 0xFFFF);
}
/** Get a random number. */
virtual uint get() {
advance();
uint x1 = vstate.x1;
uint x2 = vstate.x2;
return (x2 << 16) + x1;
}
private:
void advance()
{
/* work with unsigned long ints throughout to get correct integer
promotions of any unsigned short ints */
const uint32 x0 = vstate.x0;
const uint32 x1 = vstate.x1;
const uint32 x2 = vstate.x2;
uint32 a;
a = a0 * x0 + c0;
vstate.x0 = uint16(a & 0xFFFF);
a >>= 16;
/* although the next line may overflow we only need the top 16 bits
in the following stage, so it does not matter */
a += a0 * x1 + a1 * x0;
vstate.x1 = uint16(a & 0xFFFF);
a >>= 16;
a += a0 * x2 + a1 * x1 + a2 * x0;
vstate.x2 = uint16(a & 0xFFFF);
}
private:
NVMATH_API static const uint16 a0, a1, a2, c0;
struct rand48_state_t {
uint16 x0, x1, x2;
} vstate;
};
} // nv namespace
#endif // NV_MATH_RANDOM_H