15821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)/* 25821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) A C-program for MT19937, with initialization improved 2002/1/26. 35821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) Coded by Takuji Nishimura and Makoto Matsumoto. 45821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 55821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) Before using, initialize the state by using init_genrand(seed) 65821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) or init_by_array(init_key, key_length). 75821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 85821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, 95821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) All rights reserved. 105821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 115821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) Redistribution and use in source and binary forms, with or without 125821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) modification, are permitted provided that the following conditions 135821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) are met: 145821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 155821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 1. Redistributions of source code must retain the above copyright 165821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) notice, this list of conditions and the following disclaimer. 175821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 185821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 2. Redistributions in binary form must reproduce the above copyright 195821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) notice, this list of conditions and the following disclaimer in the 205821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) documentation and/or other materials provided with the distribution. 215821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 225821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 3. The names of its contributors may not be used to endorse or promote 235821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) products derived from this software without specific prior written 245821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) permission. 255821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 265821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 275821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 285821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 295821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR 305821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 315821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 325821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 335821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 345821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 355821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 365821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 375821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 385821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 395821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) Any feedback is very welcome. 405821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html 415821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) 425821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)*/ 435821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 445821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)#include "mt19937ar.h" 455821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 465821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)/* Period parameters */ 475821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)#define N 624 485821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)#define M 397 495821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)#define MATRIX_A 0x9908b0dfUL /* constant vector a */ 505821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)#define UPPER_MASK 0x80000000UL /* most significant w-r bits */ 515821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)#define LOWER_MASK 0x7fffffffUL /* least significant r bits */ 525821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 535821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)MersenneTwister::MersenneTwister() : mt(N), mti(N+1) { 545821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)} 555821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 565821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)MersenneTwister::~MersenneTwister() { 575821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)} 585821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 595821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)/* initializes mt[N] with a seed */ 605821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)void MersenneTwister::init_genrand(uint32 s) 615821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles){ 625821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) mt[0]= s & 0xffffffffUL; 635821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) for (mti=1; mti<N; mti++) { 645821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) mt[mti] = 655821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti); 665821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ 675821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) /* In the previous versions, MSBs of the seed affect */ 685821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) /* only MSBs of the array mt[]. */ 695821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) /* 2002/01/09 modified by Makoto Matsumoto */ 705821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) mt[mti] &= 0xffffffffUL; 715821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) /* for >32 bit machines */ 725821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) } 735821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)} 745821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 755821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)/* initialize by an array with array-length */ 765821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)/* init_key is the array for initializing keys */ 775821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)/* key_length is its length */ 785821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)/* slight change for C++, 2004/2/26 */ 795821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)void MersenneTwister::init_by_array(uint32 init_key[], int key_length) 805821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles){ 815821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) int i, j, k; 825821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) init_genrand(19650218UL); 835821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) i=1; j=0; 845821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) k = (N>key_length ? N : key_length); 855821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) for (; k; k--) { 865821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL)) 875821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) + init_key[j] + j; /* non linear */ 885821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ 895821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) i++; j++; 905821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) if (i>=N) { mt[0] = mt[N-1]; i=1; } 915821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) if (j>=key_length) j=0; 925821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) } 935821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) for (k=N-1; k; k--) { 945821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL)) 955821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) - i; /* non linear */ 965821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ 975821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) i++; 985821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) if (i>=N) { mt[0] = mt[N-1]; i=1; } 995821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) } 1005821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 1015821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */ 1025821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)} 1035821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 1045821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)/* generates a random number on [0,0xffffffff]-interval */ 1055821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)uint32 MersenneTwister::genrand_int32(void) 1065821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles){ 1075821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) uint32 y; 1085821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) static uint32 mag01[2]={0x0UL, MATRIX_A}; 1095821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) /* mag01[x] = x * MATRIX_A for x=0,1 */ 1105821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 1115821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) if (mti >= N) { /* generate N words at one time */ 1125821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) int kk; 1135821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 1145821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) if (mti == N+1) /* if init_genrand() has not been called, */ 1155821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) init_genrand(5489UL); /* a default initial seed is used */ 1165821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 1175821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) for (kk=0;kk<N-M;kk++) { 1185821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK); 1195821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL]; 1205821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) } 1215821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) for (;kk<N-1;kk++) { 1225821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK); 1235821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1UL]; 1245821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) } 1255821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK); 1265821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL]; 1275821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 1285821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) mti = 0; 1295821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) } 1305821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 1315821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) y = mt[mti++]; 1325821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 1335821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) /* Tempering */ 1345821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) y ^= (y >> 11); 1355821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) y ^= (y << 7) & 0x9d2c5680UL; 1365821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) y ^= (y << 15) & 0xefc60000UL; 1375821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) y ^= (y >> 18); 1385821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) 1395821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles) return y; 1405821806d5e7f356e8fa4b058a389a808ea183019Torne (Richard Coles)} 141