Skip to content
This repository
Newer
Older
100644 178 lines (156 sloc) 5.986 kb
65814005 » Laurent Sansonetti
2010-02-13 integrate Mersenne Twister based random number generator from 1.9 ups…
1 /*
2 * This is based on trimmed version of MT19937. To get the original version,
3 * contact <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html>.
4 *
5 * The original copyright notice follows.
6 *
7 * A C-program for MT19937, with initialization improved 2002/2/10.
8 * Coded by Takuji Nishimura and Makoto Matsumoto.
9 * This is a faster version by taking Shawn Cokus's optimization,
10 * Matthe Bellew's simplification, Isaku Wada's real version.
11 *
12 * Before using, initialize the state by using init_genrand(mt, seed)
13 * or init_by_array(mt, init_key, key_length).
14 *
15 * Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
16 * All rights reserved.
17 *
18 * Redistribution and use in source and binary forms, with or without
19 * modification, are permitted provided that the following conditions
20 * are met:
21 *
22 * 1. Redistributions of source code must retain the above copyright
23 * notice, this list of conditions and the following disclaimer.
24 *
25 * 2. Redistributions in binary form must reproduce the above copyright
26 * notice, this list of conditions and the following disclaimer in the
27 * documentation and/or other materials provided with the distribution.
28 *
29 * 3. The names of its contributors may not be used to endorse or promote
30 * products derived from this software without specific prior written
31 * permission.
32 *
33 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
34 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
35 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
36 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
37 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
38 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
39 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
40 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
41 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
42 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
43 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
44 *
45 *
46 * Any feedback is very welcome.
47 * http://www.math.keio.ac.jp/matumoto/emt.html
48 * email: matumoto@math.keio.ac.jp
49 */
50
51 #include <limits.h>
52 typedef int int_must_be_32bit_at_least[sizeof(int) * CHAR_BIT < 32 ? -1 : 1];
53
54 /* Period parameters */
55 #define N 624
56 #define M 397
57 #define MATRIX_A 0x9908b0dfU /* constant vector a */
58 #define UMASK 0x80000000U /* most significant w-r bits */
59 #define LMASK 0x7fffffffU /* least significant r bits */
60 #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) )
61 #define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1U ? MATRIX_A : 0U))
62
63 enum {MT_MAX_STATE = N};
64
65 struct MT {
66 /* assume int is enough to store 32bits */
67 unsigned int state[N]; /* the array for the state vector */
68 unsigned int *next;
69 int left;
70 };
71
72 #define genrand_initialized(mt) ((mt)->next != 0)
73 #define uninit_genrand(mt) ((mt)->next = 0)
74
75 /* initializes state[N] with a seed */
76 static void
77 init_genrand(struct MT *mt, unsigned int s)
78 {
79 int j;
80 mt->state[0] = s & 0xffffffffU;
81 for (j=1; j<N; j++) {
82 mt->state[j] = (1812433253U * (mt->state[j-1] ^ (mt->state[j-1] >> 30)) + j);
83 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
84 /* In the previous versions, MSBs of the seed affect */
85 /* only MSBs of the array state[]. */
86 /* 2002/01/09 modified by Makoto Matsumoto */
87 mt->state[j] &= 0xffffffff; /* for >32 bit machines */
88 }
89 mt->left = 1;
90 mt->next = mt->state + N;
91 }
92
93 /* initialize by an array with array-length */
94 /* init_key is the array for initializing keys */
95 /* key_length is its length */
96 /* slight change for C++, 2004/2/26 */
97 static void
98 init_by_array(struct MT *mt, unsigned int init_key[], int key_length)
99 {
100 int i, j, k;
101 init_genrand(mt, 19650218U);
102 i=1; j=0;
103 k = (N>key_length ? N : key_length);
104 for (; k; k--) {
105 mt->state[i] = (mt->state[i] ^ ((mt->state[i-1] ^ (mt->state[i-1] >> 30)) * 1664525U))
106 + init_key[j] + j; /* non linear */
107 mt->state[i] &= 0xffffffffU; /* for WORDSIZE > 32 machines */
108 i++; j++;
109 if (i>=N) { mt->state[0] = mt->state[N-1]; i=1; }
110 if (j>=key_length) j=0;
111 }
112 for (k=N-1; k; k--) {
113 mt->state[i] = (mt->state[i] ^ ((mt->state[i-1] ^ (mt->state[i-1] >> 30)) * 1566083941U))
114 - i; /* non linear */
115 mt->state[i] &= 0xffffffffU; /* for WORDSIZE > 32 machines */
116 i++;
117 if (i>=N) { mt->state[0] = mt->state[N-1]; i=1; }
118 }
119
120 mt->state[0] = 0x80000000U; /* MSB is 1; assuring non-zero initial array */
121 }
122
123 static void
124 next_state(struct MT *mt)
125 {
126 unsigned int *p = mt->state;
127 int j;
128
129 /* if init_genrand() has not been called, */
130 /* a default initial seed is used */
131 if (!genrand_initialized(mt)) init_genrand(mt, 5489U);
132
133 mt->left = N;
134 mt->next = mt->state;
135
136 for (j=N-M+1; --j; p++)
137 *p = p[M] ^ TWIST(p[0], p[1]);
138
139 for (j=M; --j; p++)
140 *p = p[M-N] ^ TWIST(p[0], p[1]);
141
142 *p = p[M-N] ^ TWIST(p[0], mt->state[0]);
143 }
144
145 /* generates a random number on [0,0xffffffff]-interval */
146 static unsigned int
147 genrand_int32(struct MT *mt)
148 {
149 unsigned int y;
150
151 if (--mt->left <= 0) next_state(mt);
152 y = *mt->next++;
153
154 /* Tempering */
155 y ^= (y >> 11);
156 y ^= (y << 7) & 0x9d2c5680;
157 y ^= (y << 15) & 0xefc60000;
158 y ^= (y >> 18);
159
160 return y;
161 }
162
163 /* generates a random number on [0,1) with 53-bit resolution*/
164 static double
165 genrand_real(struct MT *mt)
166 {
167 unsigned int a = genrand_int32(mt)>>5, b = genrand_int32(mt)>>6;
168 return(a*67108864.0+b)*(1.0/9007199254740992.0);
169 }
170
171 /* generates a random number on [0,1] with 53-bit resolution*/
172 // TODO
173 #define genrand_real2 genrand_real
174
175 /* These real versions are due to Isaku Wada, 2002/01/09 added */
176
177 #undef N
178 #undef M
Something went wrong with that request. Please try again.