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