Branch data Line data Source code
1 : : /*
2 : : A C-program for MT19937, with initialization improved 2002/2/10.
3 : : Coded by Takuji Nishimura and Makoto Matsumoto.
4 : : This is a faster version by taking Shawn Cokus's optimization,
5 : : Matthe Bellew's simplification, Isaku Wada's real version.
6 : :
7 : : Before using, initialize the state by using init_genrand(seed)
8 : : or init_by_array(init_key, key_length).
9 : :
10 : : Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
11 : : All rights reserved.
12 : :
13 : : Redistribution and use in source and binary forms, with or without
14 : : modification, are permitted provided that the following conditions
15 : : are met:
16 : :
17 : : 1. Redistributions of source code must retain the above copyright
18 : : notice, this list of conditions and the following disclaimer.
19 : :
20 : : 2. Redistributions in binary form must reproduce the above copyright
21 : : notice, this list of conditions and the following disclaimer in the
22 : : documentation and/or other materials provided with the distribution.
23 : :
24 : : 3. The names of its contributors may not be used to endorse or promote
25 : : products derived from this software without specific prior written
26 : : permission.
27 : :
28 : : THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
29 : : "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
30 : : LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
31 : : A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
32 : : CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
33 : : EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
34 : : PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
35 : : PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
36 : : LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
37 : : NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
38 : : SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
39 : :
40 : :
41 : : Any feedback is very welcome.
42 : : http://www.math.keio.ac.jp/matumoto/emt.html
43 : : email: [email protected]
44 : : */
45 : :
46 : : #include "RandomMersenne.hpp"
47 : :
48 : : /* Period parameters */
49 : : #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) )
50 : : #define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1UL ? MATRIX_A : 0UL))
51 : :
52 : : /* initializes state[N] with a seed */
53 : 43 : void RandomMersenne::init_genrand(unsigned long s)
54 : : {
55 : 43 : mt[0]= s & 0xffffffffUL;
56 [ + + ]: 26832 : for (mti=1; mti<N; mti++) {
57 : 26789 : mt[mti] =
58 : 26789 : (1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
59 : : /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
60 : : /* In the previous versions, MSBs of the seed affect */
61 : : /* only MSBs of the array mt[]. */
62 : : /* 2002/01/09 modified by Makoto Matsumoto */
63 : 26789 : mt[mti] &= 0xffffffffUL;
64 : : /* for >32 bit machines */
65 : : }
66 : 43 : }
67 : :
68 : : /* initialize by an array with array-length */
69 : : /* init_key is the array for initializing keys */
70 : : /* key_length is its length */
71 : 0 : void RandomMersenne::init_by_array(unsigned long init_key[], unsigned long key_length)
72 : : {
73 : : unsigned long i, j, k;
74 : 0 : const unsigned long n = (unsigned long)N;
75 : 0 : init_genrand(19650218UL);
76 : 0 : i=1UL; j=0UL;
77 : 0 : k = (n>key_length ? n : key_length);
78 [ # # ]: 0 : for (; k; k--) {
79 : 0 : mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL))
80 : 0 : + init_key[j] + j; /* non linear */
81 : 0 : mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
82 : 0 : i++; j++;
83 [ # # ]: 0 : if (i>=n) { mt[0] = mt[n-1]; i=1; }
84 [ # # ]: 0 : if (j>=key_length) j=0;
85 : : }
86 [ # # ]: 0 : for (k=n-1; k; k--) {
87 : 0 : mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL))
88 : 0 : - i; /* non linear */
89 : 0 : mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
90 : 0 : i++;
91 [ # # ]: 0 : if (i>=n) { mt[0] = mt[n-1]; i=1; }
92 : : }
93 : :
94 : 0 : mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
95 : 0 : }
96 : :
97 : 954 : unsigned long RandomMersenne::genrand_int32()
98 : : {
99 : : unsigned long y;
100 : : static unsigned long mag01[2]={0x0UL, MATRIX_A};
101 : : /* mag01[x] = x * MATRIX_A for x=0,1 */
102 : :
103 [ + + ]: 954 : if (mti >= N) { /* generate N words at one time */
104 : : unsigned kk;
105 : :
106 : : //if (mti == N+1) /* if init_genrand() has not been called, */
107 : : //init_genrand(5489UL); /* a default initial seed is used */
108 : :
109 [ + + ]: 9804 : for (kk=0;kk<N-M;kk++) {
110 : 9761 : y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
111 : 9761 : mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL];
112 : : }
113 [ + + ]: 17071 : for (;kk<N-1;kk++) {
114 : 17028 : y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
115 : 17028 : mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
116 : : }
117 : 43 : y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
118 : 43 : mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
119 : :
120 : 43 : mti = 0;
121 : : }
122 : :
123 : 954 : y = mt[mti++];
124 : :
125 : : /* Tempering */
126 : 954 : y ^= (y >> 11);
127 : 954 : y ^= (y << 7) & 0x9d2c5680UL;
128 : 954 : y ^= (y << 15) & 0xefc60000UL;
129 : 954 : y ^= (y >> 18);
130 : :
131 : 954 : return y;
132 : : }
133 : :
134 : : #ifdef TEST_RANDOM
135 : :
136 : : #include <cstdio>
137 : :
138 : : int main(void)
139 : : {
140 : : int i;
141 : : unsigned long init[4]={0x123, 0x234, 0x345, 0x456}, length=4;
142 : : RandomMersenne rm(init, length);
143 : : printf("1000 outputs of genrand_int32()\n");
144 : : for (i=0; i<1000; i++) {
145 : : printf("%10lu ", rm.genrand_int32());
146 : : if (i%5==4) printf("\n");
147 : : }
148 : : printf("\n1000 outputs of genrand_real2()\n");
149 : : for (i=0; i<1000; i++) {
150 : : printf("%10.8f ", rm.genrand_real2());
151 : : if (i%5==4) printf("\n");
152 : : }
153 : : return 0;
154 : : }
155 : : #endif
156 : :
|
