00001 // mtrand.h 00002 // C++ include file for MT19937, with initialization improved 2002/1/26. 00003 // Coded by Takuji Nishimura and Makoto Matsumoto. 00004 // Ported to C++ by Jasper Bedaux 2003/1/1 (see http://www.bedaux.net/mtrand/). 00005 // The generators returning floating point numbers are based on 00006 // a version by Isaku Wada, 2002/01/09 00007 // 00008 // Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, 00009 // All rights reserved. 00010 // 00011 // Redistribution and use in source and binary forms, with or without 00012 // modification, are permitted provided that the following conditions 00013 // are met: 00014 // 00015 // 1. Redistributions of source code must retain the above copyright 00016 // notice, this list of conditions and the following disclaimer. 00017 // 00018 // 2. Redistributions in binary form must reproduce the above copyright 00019 // notice, this list of conditions and the following disclaimer in the 00020 // documentation and/or other materials provided with the distribution. 00021 // 00022 // 3. The names of its contributors may not be used to endorse or promote 00023 // products derived from this software without specific prior written 00024 // permission. 00025 // 00026 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 00027 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 00028 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 00029 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR 00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00037 // 00038 // Any feedback is very welcome. 00039 // http://www.math.keio.ac.jp/matumoto/emt.html 00040 // email: matumoto@math.keio.ac.jp 00041 // 00042 // Feedback about the C++ port should be sent to Jasper Bedaux, 00043 // see http://www.bedaux.net/mtrand/ for e-mail address and info. 00044 00045 #ifndef MTRAND_H 00046 #define MTRAND_H 00047 00048 class MTRand_int32 { // Mersenne Twister random number generator 00049 public: 00050 // default constructor: uses default seed only if this is the first instance 00051 MTRand_int32() { if (!init) seed(5489UL); init = true; } 00052 // constructor with 32 bit int as seed 00053 MTRand_int32(unsigned long s) { seed(s); init = true; } 00054 // constructor with array of size 32 bit ints as seed 00055 MTRand_int32(const unsigned long* array, int size) { seed(array, size); init = true; } 00056 // the two seed functions 00057 void seed(unsigned long); // seed with 32 bit integer 00058 void seed(const unsigned long*, int size); // seed with array 00059 // overload operator() to make this a generator (functor) 00060 unsigned long operator()() { return rand_int32(); } 00061 // 2007-02-11: made the destructor virtual; thanks "double more" for pointing this out 00062 virtual ~MTRand_int32() {} // destructor 00063 protected: // used by derived classes, otherwise not accessible; use the ()-operator 00064 unsigned long rand_int32(); // generate 32 bit random integer 00065 private: 00066 static const int n = 624, m = 397; // compile time constants 00067 // the variables below are static (no duplicates can exist) 00068 static unsigned long state[n]; // state vector array 00069 static int p; // position in state array 00070 static bool init; // true if init function is called 00071 // private functions used to generate the pseudo random numbers 00072 unsigned long twiddle(unsigned long, unsigned long); // used by gen_state() 00073 void gen_state(); // generate new state 00074 // make copy constructor and assignment operator unavailable, they don't make sense 00075 MTRand_int32(const MTRand_int32&); // copy constructor not defined 00076 void operator=(const MTRand_int32&); // assignment operator not defined 00077 }; 00078 00079 // inline for speed, must therefore reside in header file 00080 inline unsigned long MTRand_int32::twiddle(unsigned long u, unsigned long v) { 00081 return (((u & 0x80000000UL) | (v & 0x7FFFFFFFUL)) >> 1) 00082 ^ ((v & 1UL) ? 0x9908B0DFUL : 0x0UL); 00083 } 00084 00085 inline unsigned long MTRand_int32::rand_int32() { // generate 32 bit random int 00086 if (p == n) gen_state(); // new state vector needed 00087 // gen_state() is split off to be non-inline, because it is only called once 00088 // in every 624 calls and otherwise irand() would become too big to get inlined 00089 unsigned long x = state[p++]; 00090 x ^= (x >> 11); 00091 x ^= (x << 7) & 0x9D2C5680UL; 00092 x ^= (x << 15) & 0xEFC60000UL; 00093 return x ^ (x >> 18); 00094 } 00095 00096 // generates double floating point numbers in the half-open interval [0, 1) 00097 class MTRand : public MTRand_int32 { 00098 public: 00099 MTRand() : MTRand_int32() {} 00100 MTRand(unsigned long seed) : MTRand_int32(seed) {} 00101 MTRand(const unsigned long* seed, int size) : MTRand_int32(seed, size) {} 00102 ~MTRand() {} 00103 double operator()() { 00104 return static_cast<double>(rand_int32()) * (1. / 4294967296.); } // divided by 2^32 00105 private: 00106 MTRand(const MTRand&); // copy constructor not defined 00107 void operator=(const MTRand&); // assignment operator not defined 00108 }; 00109 00110 // generates double floating point numbers in the closed interval [0, 1] 00111 class MTRand_closed : public MTRand_int32 { 00112 public: 00113 MTRand_closed() : MTRand_int32() {} 00114 MTRand_closed(unsigned long seed) : MTRand_int32(seed) {} 00115 MTRand_closed(const unsigned long* seed, int size) : MTRand_int32(seed, size) {} 00116 ~MTRand_closed() {} 00117 double operator()() { 00118 return static_cast<double>(rand_int32()) * (1. / 4294967295.); } // divided by 2^32 - 1 00119 private: 00120 MTRand_closed(const MTRand_closed&); // copy constructor not defined 00121 void operator=(const MTRand_closed&); // assignment operator not defined 00122 }; 00123 00124 // generates double floating point numbers in the open interval (0, 1) 00125 class MTRand_open : public MTRand_int32 { 00126 public: 00127 MTRand_open() : MTRand_int32() {} 00128 MTRand_open(unsigned long seed) : MTRand_int32(seed) {} 00129 MTRand_open(const unsigned long* seed, int size) : MTRand_int32(seed, size) {} 00130 ~MTRand_open() {} 00131 double operator()() { 00132 return (static_cast<double>(rand_int32()) + .5) * (1. / 4294967296.); } // divided by 2^32 00133 private: 00134 MTRand_open(const MTRand_open&); // copy constructor not defined 00135 void operator=(const MTRand_open&); // assignment operator not defined 00136 }; 00137 00138 // generates 53 bit resolution doubles in the half-open interval [0, 1) 00139 class MTRand53 : public MTRand_int32 { 00140 public: 00141 MTRand53() : MTRand_int32() {} 00142 MTRand53(unsigned long seed) : MTRand_int32(seed) {} 00143 MTRand53(const unsigned long* seed, int size) : MTRand_int32(seed, size) {} 00144 ~MTRand53() {} 00145 double operator()() { 00146 return (static_cast<double>(rand_int32() >> 5) * 67108864. + 00147 static_cast<double>(rand_int32() >> 6)) * (1. / 9007199254740992.); } 00148 private: 00149 MTRand53(const MTRand53&); // copy constructor not defined 00150 void operator=(const MTRand53&); // assignment operator not defined 00151 }; 00152 00153 #endif // MTRAND_H