From 7dbcb078e5b9cac2c01c80fdaad479eae93d3901 Mon Sep 17 00:00:00 2001 From: cignoni Date: Thu, 17 Apr 2014 08:19:06 +0000 Subject: [PATCH] Standardized the generate method of the marsenne twister random generator in order to get also a unsigned capped random generation (like all the other generate() of the other random generators) --- vcg/math/random_generator.h | 426 ++++++++++++++++++------------------ 1 file changed, 215 insertions(+), 211 deletions(-) diff --git a/vcg/math/random_generator.h b/vcg/math/random_generator.h index 442cdd38..49fb804d 100644 --- a/vcg/math/random_generator.h +++ b/vcg/math/random_generator.h @@ -40,29 +40,29 @@ class RandomGenerator // construction public: - RandomGenerator(){} + RandomGenerator(){} - virtual ~RandomGenerator() - {} + virtual ~RandomGenerator() + {} // public methods public: - /// (Re-)initialize with a given seed. - virtual void initialize(unsigned int seed)=0; + /// (Re-)initialize with a given seed. + virtual void initialize(unsigned int seed)=0; - /// Return a random number in the given range (note that not all the RNG can handle a given limit). - virtual unsigned int generate(unsigned int limit)=0; + /// Return a random number in the given range (note that not all the RNG can handle a given limit). + virtual unsigned int generate(unsigned int limit)=0; - /// Return a random number in the [0,1) real interval. - virtual double generate01()=0; + /// Return a random number in the [0,1) real interval. + virtual double generate01()=0; - /// Returns a random number in the [0,1] real interval. - virtual double generate01closed()=0; + /// Returns a random number in the [0,1] real interval. + virtual double generate01closed()=0; - /// Generates a random number in the (0,1) real interval. - virtual double generate01open()=0; - virtual double generateRange(double minV, double maxV) { return minV+(maxV-minV)*generate01(); } + /// Generates a random number in the (0,1) real interval. + virtual double generate01open()=0; + virtual double generateRange(double minV, double maxV) { return minV+(maxV-minV)*generate01(); } }; @@ -81,9 +81,9 @@ vcg::Point3 GenerateBarycentricUniform(GeneratorType &rnd) interp[2] = 1.0 - interp[2]; } - assert(interp[1] + interp[2] <= 1.0); - interp[0]=1.0-(interp[1] + interp[2]); - return interp; + assert(interp[1] + interp[2] <= 1.0); + interp[0]=1.0-(interp[1] + interp[2]); + return interp; } /// \brief Generate a random point insidie a box with uniform distribution @@ -159,80 +159,80 @@ class SubtractiveRingRNG : public RandomGenerator // private data member private: - // Subtractive Ring RNG status variables - unsigned int _M_table[55]; - size_t _M_index1; - size_t _M_index2; + // Subtractive Ring RNG status variables + unsigned int _M_table[55]; + size_t _M_index1; + size_t _M_index2; // construction public: - // ctor - SubtractiveRingRNG(int default_seed=161803398u) - { - initialize(default_seed); - } + // ctor + SubtractiveRingRNG(int default_seed=161803398u) + { + initialize(default_seed); + } - virtual ~SubtractiveRingRNG() - {} + virtual ~SubtractiveRingRNG() + {} // public methods public: - /// (Re-)initialize with a given seed. - void initialize(unsigned int seed) - { - unsigned int __k = 1; - _M_table[54] = seed; - size_t __i; - for (__i = 0; __i < 54; __i++) - { - size_t __ii = (21 * (__i + 1) % 55) - 1; - _M_table[__ii] = __k; - __k = seed - __k; - seed = _M_table[__ii]; - } - for (int __loop = 0; __loop < 4; __loop++) - { - for (__i = 0; __i < 55; __i++) - _M_table[__i] = _M_table[__i] - _M_table[(1 + __i + 30) % 55]; - } - _M_index1 = 0; - _M_index2 = 31; - } + /// (Re-)initialize with a given seed. + void initialize(unsigned int seed) + { + unsigned int __k = 1; + _M_table[54] = seed; + size_t __i; + for (__i = 0; __i < 54; __i++) + { + size_t __ii = (21 * (__i + 1) % 55) - 1; + _M_table[__ii] = __k; + __k = seed - __k; + seed = _M_table[__ii]; + } + for (int __loop = 0; __loop < 4; __loop++) + { + for (__i = 0; __i < 55; __i++) + _M_table[__i] = _M_table[__i] - _M_table[(1 + __i + 30) % 55]; + } + _M_index1 = 0; + _M_index2 = 31; + } - /// Return a random number in the given range (limit) using the Subtractive Ring method. - unsigned int generate(unsigned int limit= 0xffffffffu) - { - _M_index1 = (_M_index1 + 1) % 55; - _M_index2 = (_M_index2 + 1) % 55; - _M_table[_M_index1] = _M_table[_M_index1] - _M_table[_M_index2]; - return _M_table[_M_index1] % limit; - } + /// Return a random number in the given range (limit) using the Subtractive Ring method. + unsigned int generate(unsigned int limit= 0xffffffffu) + { + _M_index1 = (_M_index1 + 1) % 55; + _M_index2 = (_M_index2 + 1) % 55; + _M_table[_M_index1] = _M_table[_M_index1] - _M_table[_M_index2]; + return _M_table[_M_index1] % limit; + } - /// Return a random number in the [0,1) real interval using the Subtractive Ring method. - double generate01() - { - const unsigned int lmt = 0xffffffffu; - unsigned int number = generate(lmt); - return static_cast(number) / static_cast(lmt); - } + /// Return a random number in the [0,1) real interval using the Subtractive Ring method. + double generate01() + { + const unsigned int lmt = 0xffffffffu; + unsigned int number = generate(lmt); + return static_cast(number) / static_cast(lmt); + } - /// Returns a random number in the [0,1] real interval using the Subtractive Ring method. - double generate01closed() - { - const unsigned int lmt = 0xffffffffu; - unsigned int number = generate(lmt); - return static_cast(number) / static_cast(0xfffffffEu); - } + /// Returns a random number in the [0,1] real interval using the Subtractive Ring method. + double generate01closed() + { + const unsigned int lmt = 0xffffffffu; + unsigned int number = generate(lmt); + return static_cast(number) / static_cast(0xfffffffEu); + } - /// Generates a random number in the (0,1) real interval using the Subtractive Ring method. - double generate01open() - { - const unsigned int lmt = 0xffffffffu; - unsigned int number = generate(lmt); - return (static_cast(number) + 0.5) * (1.0/static_cast(lmt)); - } + /// Generates a random number in the (0,1) real interval using the Subtractive Ring method. + double generate01open() + { + const unsigned int lmt = 0xffffffffu; + unsigned int number = generate(lmt); + return (static_cast(number) + 0.5) * (1.0/static_cast(lmt)); + } }; @@ -253,174 +253,178 @@ class MarsenneTwisterRNG : public RandomGenerator // definitions private: - static const int N = 624; - static const int M = 397; - static const unsigned int MATRIX_A = 0x9908b0dfu; // constant vector a - static const unsigned int UPPER_MASK = 0x80000000u; // most significant w-r bits - static const unsigned int LOWER_MASK = 0x7fffffffu; // least significant r bits + static const int N = 624; + static const int M = 397; + static const unsigned int MATRIX_A = 0x9908b0dfu; // constant vector a + static const unsigned int UPPER_MASK = 0x80000000u; // most significant w-r bits + static const unsigned int LOWER_MASK = 0x7fffffffu; // least significant r bits // private data member private: - // Improved Marsenne-Twister RNG status variables - unsigned int mt[N]; // the array for the state vector - int mti; + // Improved Marsenne-Twister RNG status variables + unsigned int mt[N]; // the array for the state vector + int mti; // construction public: - // ctor - MarsenneTwisterRNG() - { - initialize(5489u); - } + // ctor + MarsenneTwisterRNG() + { + initialize(5489u); + } - MarsenneTwisterRNG(unsigned int seed) - { - initialize(seed); - } + MarsenneTwisterRNG(unsigned int seed) + { + initialize(seed); + } - virtual ~MarsenneTwisterRNG() - {} + virtual ~MarsenneTwisterRNG() + {} // public methods public: - /// (Re-)initialize with the given seed. - void initialize(unsigned int seed) - { - mt[0]= seed & 0xffffffffu; - for (mti=1; mti> 30)) + mti); - /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ - /* In the previous versions, MSBs of the seed affect */ - /* only MSBs of the array mt[]. */ - /* 2002/01/09 modified by Makoto Matsumoto */ - mt[mti] &= 0xffffffffu; - /* for >32 bit machines */ - } - } + /// (Re-)initialize with the given seed. + void initialize(unsigned int seed) + { + mt[0]= seed & 0xffffffffu; + for (mti=1; mti> 30)) + mti); + /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ + /* In the previous versions, MSBs of the seed affect */ + /* only MSBs of the array mt[]. */ + /* 2002/01/09 modified by Makoto Matsumoto */ + mt[mti] &= 0xffffffffu; + /* for >32 bit machines */ + } + } - /** - * Initialize by an array with array-length. - * - * init_key is the array for initializing keys - * key_length is its length - */ - void initializeByArray(unsigned int init_key[], int key_length) - { - int i, j, k; - initialize(19650218u); - i=1; j=0; - k = (N>key_length ? N : key_length); - for (; k; k--) - { - mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525u)) + init_key[j] + j; /* non linear */ - mt[i] &= 0xffffffffu; /* for WORDSIZE > 32 machines */ - i++; j++; + /** + * Initialize by an array with array-length. + * + * init_key is the array for initializing keys + * key_length is its length + */ + void initializeByArray(unsigned int init_key[], int key_length) + { + int i, j, k; + initialize(19650218u); + i=1; j=0; + k = (N>key_length ? N : key_length); + for (; k; k--) + { + mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525u)) + init_key[j] + j; /* non linear */ + mt[i] &= 0xffffffffu; /* for WORDSIZE > 32 machines */ + i++; j++; - if (i>=N) - { - mt[0] = mt[N-1]; - i=1; - } + if (i>=N) + { + mt[0] = mt[N-1]; + i=1; + } - if (j>=key_length) j=0; - } + if (j>=key_length) j=0; + } - for (k=N-1; k; k--) - { - mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941u)) - i; /* non linear */ - mt[i] &= 0xffffffffu; /* for WORDSIZE > 32 machines */ - i++; - if (i>=N) - { - mt[0] = mt[N-1]; - i=1; - } - } + for (k=N-1; k; k--) + { + mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941u)) - i; /* non linear */ + mt[i] &= 0xffffffffu; /* for WORDSIZE > 32 machines */ + i++; + if (i>=N) + { + mt[0] = mt[N-1]; + i=1; + } + } - mt[0] = 0x80000000u; /* MSB is 1; assuring non-zero initial array */ - } + mt[0] = 0x80000000u; /* MSB is 1; assuring non-zero initial array */ + } - /** - * Return a random number in the [0,0xffffffff] interval using the improved Marsenne Twister algorithm. - * - * NOTE: Limit is not considered, the interval is fixed. - */ - unsigned int generate(unsigned int /*limit*/) - { - unsigned int y; - static unsigned int mag01[2]={0x0u, MATRIX_A}; - /* mag01[x] = x * MATRIX_A for x=0,1 */ + unsigned int generate(unsigned int limit) + { + return generate()%limit; + } - if (mti >= N) // generate N words at one time - { - int kk; + /** + * Return a random number in the [0,0xffffffff] interval using the improved Marsenne Twister algorithm. + * + * NOTE: Limit is not considered, the interval is fixed. + */ + unsigned int generate() + { + unsigned int y; + static unsigned int mag01[2]={0x0u, MATRIX_A}; + /* mag01[x] = x * MATRIX_A for x=0,1 */ - for (kk=0;kk> 1) ^ mag01[y & 0x1u]; - } + if (mti >= N) // generate N words at one time + { + int kk; - for (;kk> 1) ^ mag01[y & 0x1u]; - } + for (kk=0;kk> 1) ^ mag01[y & 0x1u]; + } - y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK); - mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1u]; + for (;kk> 1) ^ mag01[y & 0x1u]; + } - mti = 0; - } + y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK); + mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1u]; - y = mt[mti++]; + mti = 0; + } - /* Tempering */ - y ^= (y >> 11); - y ^= (y << 7) & 0x9d2c5680u; - y ^= (y << 15) & 0xefc60000u; - y ^= (y >> 18); + y = mt[mti++]; - return y; - } + /* Tempering */ + y ^= (y >> 11); + y ^= (y << 7) & 0x9d2c5680u; + y ^= (y << 15) & 0xefc60000u; + y ^= (y >> 18); - /// Returns a random number in the [0,1] real interval using the improved Marsenne-Twister. - double generate01closed() - { - return generate(0)*(1.0/4294967295.0); - } + return y; + } - /// Returns a random number in the [0,1) real interval using the improved Marsenne-Twister. - double generate01() - { - return generate(0)*(1.0/4294967296.0); - } + /// Returns a random number in the [0,1] real interval using the improved Marsenne-Twister. + double generate01closed() + { + return generate()*(1.0/4294967295.0); + } - /// Generates a random number in the (0,1) real interval using the improved Marsenne-Twister. - double generate01open() - { - return (((double)generate(0)) + 0.5)*(1.0/4294967296.0); - } + /// Returns a random number in the [0,1) real interval using the improved Marsenne-Twister. + double generate01() + { + return generate()*(1.0/4294967296.0); + } - /// Generate a random triple of baricentric coords - template - void generateBarycentric(PointType &p){ - p[1] = this->generate01(); - p[2] = this->generate01(); + /// Generates a random number in the (0,1) real interval using the improved Marsenne-Twister. + double generate01open() + { + return (((double)generate()) + 0.5)*(1.0/4294967296.0); + } - if(p[1] + p[2] > 1.0){ - p[1] = 1.0 - p[1]; - p[2] = 1.0 - p[2]; - } - p[0]=1.0-(p[1] + p[2]); - } -}; + /// Generate a random triple of baricentric coords + template + void generateBarycentric(PointType &p){ + p[1] = this->generate01(); + p[2] = this->generate01(); + if(p[1] + p[2] > 1.0){ + p[1] = 1.0 - p[1]; + p[2] = 1.0 - p[2]; + } + p[0]=1.0-(p[1] + p[2]); + } +}; // end class MarsenneTwisterRNG /* Returns a value with normal distribution with mean m, standard deviation s *