1 /* 2 * Scythe Statistical Library 3 * Copyright (C) 2000-2002 Andrew D. Martin and Kevin M. Quinn; 4 * 2002-present Andrew D. Martin, Kevin M. Quinn, and Daniel 5 * Pemstein. All Rights Reserved. 6 * 7 * This program is free software; you can redistribute it and/or modify 8 * under the terms of the GNU General Public License as published by 9 * Free Software Foundation; either version 2 of the License, or (at 10 * your option) any later version. See the text files COPYING 11 * and LICENSE, distributed with this source code, for further 12 * information. 13 * -------------------------------------------------------------------- 14 * scythestat/rng/mersenne.h 15 * 16 * Provides the class definition for the mersenne random number 17 * generator. This class extends the base rng class by providing an 18 * implementation of runif() based on an implementation of the 19 * mersenne twister, released under the following license: 20 * 21 * A C-program for MT19937, with initialization improved 2002/1/26. 22 * Coded by Takuji Nishimura and Makoto Matsumoto. 23 * 24 * Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, 25 * All rights reserved. 26 * 27 * Redistribution and use in source and binary forms, with or without 28 * modification, are permitted provided that the following conditions 29 * are met: 30 * 31 * 1. Redistributions of source code must retain the above copyright 32 * notice, this list of conditions and the following disclaimer. 33 * 34 * 2. Redistributions in binary form must reproduce the above 35 * copyright 36 * notice, this list of conditions and the following disclaimer 37 * in the documentation and/or other materials provided with the 38 * distribution. 39 * 40 * 3. The names of its contributors may not be used to endorse or 41 * promote products derived from this software without specific 42 * prior written permission. 43 * 44 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND 45 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, 46 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF 47 * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 48 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR 49 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT INCIDENTAL, 50 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 51 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF 52 * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED 53 * AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 54 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN 55 * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 56 * POSSIBILITY OF SUCH DAMAGE. 57 * 58 * For more information see: 59 * http://www.math.keio.ac.jp/matumoto/emt.html 60 * 61 */ 62 63 /*! \file mersenne.h 64 * \brief The Mersenne Twister random number generator. 65 * 66 * This file contains the mersenne class, a class that extends 67 * Scythe's base random number generation class (scythe::rng) by 68 * providing an implementation of scythe::rng::runif() using the 69 * Mersenne Twister algorithm. 70 */ 71 72 #ifndef SCYTHE_MERSENNE_H 73 #define SCYTHE_MERSENNE_H 74 75 #ifdef SCYTHE_COMPILE_DIRECT 76 #include "rng.h" 77 #else 78 #include "scythestat/rng.h" 79 #endif 80 81 namespace scythe { 82 83 #ifdef __MINGW32__ 84 /* constant vector a */ 85 static const unsigned long MATRIX_A = 0x9908b0dfUL; 86 87 /* most significant w-r bits */ 88 static const unsigned long UPPER_MASK = 0x80000000UL; 89 90 /* least significant r bits */ 91 static const unsigned long LOWER_MASK = 0x7fffffffUL; 92 #else 93 namespace { 94 /* constant vector a */ 95 const unsigned long MATRIX_A = 0x9908b0dfUL; 96 97 /* most significant w-r bits */ 98 const unsigned long UPPER_MASK = 0x80000000UL; 99 100 /* least significant r bits */ 101 const unsigned long LOWER_MASK = 0x7fffffffUL; 102 } 103 #endif 104 105 /*! \brief The Mersenne Twister random number generator. 106 * 107 * This class defines a random number generator, using the Mersenne 108 * Twister algorithm developed and implemented by Makoto Matsumoto 109 * and Takuji Nishimura (1997, 2002). The period of this random 110 * number generator is \f$2^{19937} - 1\f$. 111 * 112 * The mersenne class extends Scythe's basic random number 113 * generating class, scythe::rng, implementing the interface that it 114 * defines. 115 * 116 * \see rng 117 * \see lecuyer 118 * 119 */ 120 class mersenne: public rng<mersenne> 121 { 122 public: 123 124 /*! \brief Default constructor 125 * 126 * This constructor generates an unseeded and uninitialized 127 * mersenne object. It is most useful for creating arrays of 128 * random number generators. An uninitialized mersenne object 129 * will be seeded with the default seed (5489UL) automatically 130 * upon use. 131 * 132 * \see mersenne(const mersenne &m) 133 * \see initialize(unsigned long s) 134 */ mersenne()135 mersenne () 136 : rng<mersenne> (), 137 mti (N + 1) 138 {} 139 140 /*! \brief Copy constructor 141 * 142 * This constructor makes a copy of an existing mersenne 143 * object, duplicating its seed and current state exactly. 144 * 145 * \param m An existing mersenne random number generator. 146 * 147 * \see mersenne() 148 */ mersenne(const mersenne & m)149 mersenne (const mersenne &m) 150 : rng<mersenne> (), 151 mti (m.mti) 152 { 153 } 154 155 /*! \brief Sets the seed. 156 * 157 * This method sets the seed of the random number generator and 158 * readies it to begin generating random numbers. Calling this 159 * function on a mersenne object that is already in use is 160 * supported, although not suggested unless you know what you 161 * are doing. 162 * 163 * \param s A long integer seed. 164 * 165 * \see mersenne() 166 */ initialize(unsigned long s)167 void initialize (unsigned long s) 168 { 169 mt[0]= s & 0xffffffffUL; 170 for (mti=1; mti<N; mti++) { 171 mt[mti] = (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) 172 + mti); 173 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ 174 /* In the previous versions, MSBs of the seed affect */ 175 /* only MSBs of the array mt[]. */ 176 /* 2002/01/09 modified by Makoto Matsumoto */ 177 mt[mti] &= 0xffffffffUL; 178 /* for >32 bit machines */ 179 } 180 } 181 182 /*! \brief Generate a random uniform variate on (0, 1). 183 * 184 * This routine returns a random double precision floating point 185 * number from the uniform distribution on the interval (0, 186 * 1). This method overloads the pure virtual method of the 187 * same name in the rng base class. 188 * 189 * \see runif(unsigned int, unsigned int) 190 * \see genrand_int32() 191 * \see rng 192 */ runif()193 inline double runif() 194 { 195 return (((double) genrand_int32()) + 0.5) * 196 (1.0 / 4294967296.0); 197 } 198 199 /* We have to override the overloaded forms of runif because 200 * overloading the no-arg runif() hides the base class 201 * definition; C++ stops looking once it finds the above. 202 */ 203 /*! \brief Generate a Matrix of random uniform variates. 204 * 205 * This routine returns a Matrix of double precision random 206 * uniform variates. on the interval (0, 1). This method 207 * overloads the virtual method of the same name in the rng base 208 * class. 209 * 210 * This is the general template version of this method and 211 * is called through explicit template instantiation. 212 * 213 * \param rows The number of rows in the returned Matrix. 214 * \param cols The number of columns in the returned Matrix. 215 * 216 * \see runif() 217 * \see rng 218 * 219 * \note We are forced to override this overloaded method 220 * because the 1-arg version of runif() hides the base class's 221 * definition of this method from the compiler, although it 222 * probably should not. 223 */ 224 template <matrix_order O, matrix_style S> runif(unsigned int rows,unsigned int cols)225 inline Matrix<double,O,S> runif(unsigned int rows, 226 unsigned int cols) 227 { 228 return rng<mersenne>::runif<O,S>(rows, cols); 229 } 230 231 /*! \brief Generate a Matrix of random uniform variates. 232 * 233 * This routine returns a Matrix of double precision random 234 * uniform variates on the interval (0, 1). This method 235 * overloads the virtual method of the same name in the rng base 236 * class. 237 * 238 * This is the default template version of this method and 239 * is called through implicit template instantiation. 240 * 241 * \param rows The number of rows in the returned Matrix. 242 * \param cols The number of columns in the returned Matrix. 243 * 244 * \see runif() 245 * \see rng 246 * 247 * \note We are forced to override this overloaded method 248 * because the 1-arg version of runif() hides the base class's 249 * definition of this method from the compiler, although it 250 * probably should not. 251 */ runif(unsigned int rows,unsigned int cols)252 Matrix<double,Col,Concrete> runif(unsigned int rows, 253 unsigned int cols) 254 { 255 return rng<mersenne>::runif<Col,Concrete>(rows, cols); 256 } 257 258 /* generates a random number on [0,0xffffffff]-interval */ 259 /*! \brief Generate a random long integer. 260 * 261 * This method generates a random integer, drawn from the 262 * discrete uniform distribution on the interval [0,0xffffffff]. 263 * 264 * \see runif() 265 * \see initialize(unsigned long s) 266 */ genrand_int32()267 unsigned long genrand_int32() 268 { 269 unsigned long y; 270 static unsigned long mag01[2]={0x0UL, MATRIX_A}; 271 /* mag01[x] = x * MATRIX_A for x=0,1 */ 272 273 if (mti >= N) { /* generate N words at one time */ 274 int kk; 275 276 if (mti == N+1) // if init_genrand() has not been called, 277 this->initialize(5489UL); // a default initial seed is used 278 279 for (kk=0;kk<N-M;kk++) { 280 y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK); 281 mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL]; 282 } 283 for (;kk<N-1;kk++) { 284 y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK); 285 mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1UL]; 286 } 287 y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK); 288 mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL]; 289 290 mti = 0; 291 } 292 293 y = mt[mti++]; 294 295 /* Tempering */ 296 y ^= (y >> 11); 297 y ^= (y << 7) & 0x9d2c5680UL; 298 y ^= (y << 15) & 0xefc60000UL; 299 y ^= (y >> 18); 300 301 return y; 302 } 303 304 protected: 305 /* Period parameters */ 306 static const int N = 624; 307 static const int M = 398; 308 309 /* the array for the state vector */ 310 unsigned long mt[N]; 311 312 /* mti==N+1 means mt[N] is not initialized */ 313 int mti; 314 }; 315 316 } 317 318 #endif /* SCYTHE_MERSENNE_H */ 319