libmove3d
3.13.0
|
00001 // MersenneTwister.h 00002 // Mersenne Twister random number generator -- a C++ class MTRand 00003 // Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus 00004 // Richard J. Wagner v1.0 15 May 2003 rjwagner@writeme.com 00005 00006 // The Mersenne Twister is an algorithm for generating random numbers. It 00007 // was designed with consideration of the flaws in various other generators. 00008 // The period, 2^19937-1, and the order of equidistribution, 623 dimensions, 00009 // are far greater. The generator is also fast; it avoids multiplication and 00010 // division, and it benefits from caches and pipelines. For more information 00011 // see the inventors' web page at http://www.math.keio.ac.jp/~matumoto/emt.html 00012 00013 // Reference 00014 // M. Matsumoto and T. Nishimura, "Mersenne Twister: A 623-Dimensionally 00015 // Equidistributed Uniform Pseudo-Random Number Generator", ACM Transactions on 00016 // Modeling and Computer Simulation, Vol. 8, No. 1, January 1998, pp 3-30. 00017 00018 // Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, 00019 // Copyright (C) 2000 - 2003, Richard J. Wagner 00020 // All rights reserved. 00021 // 00022 // Redistribution and use in source and binary forms, with or without 00023 // modification, are permitted provided that the following conditions 00024 // are met: 00025 // 00026 // 1. Redistributions of source code must retain the above copyright 00027 // notice, this list of conditions and the following disclaimer. 00028 // 00029 // 2. Redistributions in binary form must reproduce the above copyright 00030 // notice, this list of conditions and the following disclaimer in the 00031 // documentation and/or other materials provided with the distribution. 00032 // 00033 // 3. The names of its contributors may not be used to endorse or promote 00034 // products derived from this software without specific prior written 00035 // permission. 00036 // 00037 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 00038 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 00039 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 00040 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR 00041 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 00042 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 00043 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 00044 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 00045 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 00046 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 00047 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00048 00049 // The original code included the following notice: 00050 // 00051 // When you use this, send an email to: matumoto@math.keio.ac.jp 00052 // with an appropriate reference to your work. 00053 // 00054 // It would be nice to CC: rjwagner@writeme.com and Cokus@math.washington.edu 00055 // when you write. 00056 00057 #ifndef MERSENNETWISTER_H 00058 #define MERSENNETWISTER_H 00059 00060 //Note: for inclusion in MSL, the original file, MersenneTwister.h, was 00061 //split into MTRand.h and MTRand.C and the "inline" declarations were 00062 //removed. The original code can be found at: 00063 //http://www-personal.engin.umich.edu/~wagnerr/MersenneTwister.html 00064 00065 // Not thread safe (unless auto-initialization is avoided and each thread has 00066 // its own MTRand object) 00067 00068 #include <iostream> 00069 #include <limits.h> 00070 #include <cstdio> 00071 #include <ctime> 00072 #include <math.h> 00073 #include <limits.h> 00074 00075 //using namespace std; 00076 class MTRand { 00077 // Data 00078 public: 00079 typedef unsigned long uint32; // unsigned integer type, at least 32 bits 00080 00081 enum { N = 624 }; // length of state vector 00082 enum { SAVE = N + 1 }; // length of array for save() 00083 00084 protected: 00085 enum { M = 397 }; // period parameter 00086 00087 uint32 state[N]; // internal state 00088 uint32 *pNext; // next value to get from state 00089 int left; // number of values left before reload needed 00090 00091 00092 //Methods 00093 public: 00094 MTRand( const uint32& oneSeed ); // initialize with a simple uint32 00095 MTRand( uint32 *const bigSeed, uint32 const seedLength = N ); // or an array 00096 MTRand(); // auto-initialize with /dev/urandom or time() and clock() 00097 00098 // Do NOT use for CRYPTOGRAPHY without securely hashing several returned 00099 // values together, otherwise the generator state can be learned after 00100 // reading 624 consecutive values. 00101 00102 // Access to 32-bit random numbers 00103 double rand(); // real number in [0,1] 00104 double rand( const double& n ); // real number in [0,n] 00105 double randExc(); // real number in [0,1) 00106 double randExc( const double& n ); // real number in [0,n) 00107 double randDblExc(); // real number in (0,1) 00108 double randDblExc( const double& n ); // real number in (0,n) 00109 uint32 randInt(); // integer in [0,2^32-1] 00110 uint32 randInt( const uint32& n ); // integer in [0,n] for n < 2^32 00111 double operator()() { return rand(); } // same as rand() 00112 00113 // Access to 53-bit random numbers (capacity of IEEE double precision) 00114 double rand53(); // real number in [0,1) 00115 00116 // Access to nonuniform random number distributions 00117 double randNorm( const double& mean = 0.0, const double& variance = 0.0 ); 00118 00119 // Re-seeding functions with same behavior as initializers 00120 void seed( const uint32 oneSeed ); 00121 void seed( uint32 *const bigSeed, const uint32 seedLength = N ); 00122 void seed(); 00123 00124 // Saving and loading generator state 00125 void save( uint32* saveArray ) const; // to array of size SAVE 00126 void load( uint32 *const loadArray ); // from such array 00127 friend std::ostream& operator<<( std::ostream& os, const MTRand& mtrand ); 00128 friend std::istream& operator>>( std::istream& is, MTRand& mtrand ); 00129 00130 protected: 00131 void initialize( const uint32 oneSeed ); 00132 void reload(); 00133 uint32 hiBit( const uint32& u ) const { return u & 0x80000000UL; } 00134 uint32 loBit( const uint32& u ) const { return u & 0x00000001UL; } 00135 uint32 loBits( const uint32& u ) const { return u & 0x7fffffffUL; } 00136 uint32 mixBits( const uint32& u, const uint32& v ) const 00137 { return hiBit(u) | loBits(v); } 00138 uint32 twist( const uint32& m, const uint32& s0, const uint32& s1 ) const 00139 { return m ^ (mixBits(s0,s1)>>1) ^ (-loBit(s1) & 0x9908b0dfUL); } 00140 static uint32 hash( time_t t, clock_t c ); 00141 }; 00142 00143 // Change log: 00144 // 00145 // v0.1 - First release on 15 May 2000 00146 // - Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus 00147 // - Translated from C to C++ 00148 // - Made completely ANSI compliant 00149 // - Designed convenient interface for initialization, seeding, and 00150 // obtaining numbers in default or user-defined ranges 00151 // - Added automatic seeding from /dev/urandom or time() and clock() 00152 // - Provided functions for saving and loading generator state 00153 // 00154 // v0.2 - Fixed bug which reloaded generator one step too late 00155 // 00156 // v0.3 - Switched to clearer, faster reload() code from Matthew Bellew 00157 // 00158 // v0.4 - Removed trailing newline in saved generator format to be consistent 00159 // with output format of built-in types 00160 // 00161 // v0.5 - Improved portability by replacing static const int's with enum's and 00162 // clarifying return values in seed(); suggested by Eric Heimburg 00163 // - Removed MAXINT constant; use 0xffffffffUL instead 00164 // 00165 // v0.6 - Eliminated seed overflow when uint32 is larger than 32 bits 00166 // - Changed integer [0,n] generator to give better uniformity 00167 // 00168 // v0.7 - Fixed operator precedence ambiguity in reload() 00169 // - Added access for real numbers in (0,1) and (0,n) 00170 // 00171 // v0.8 - Included time.h header to properly support time_t and clock_t 00172 // 00173 // v1.0 - Revised seeding to match 26 Jan 2002 update of Nishimura and Matsumoto 00174 // - Allowed for seeding with arrays of any length 00175 // - Added access for real numbers in [0,1) with 53-bit resolution 00176 // - Added access for real numbers from normal (Gaussian) distributions 00177 // - Increased overall speed by optimizing twist() 00178 // - Doubled speed of integer [0,n] generation 00179 // - Fixed out-of-range number generation on 64-bit machines 00180 // - Improved portability by substituting literal constants for long enum's 00181 // - Changed license from GNU LGPL to BSD 00182 00183 #endif // MERSENNETWISTER_H