1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.apache.commons.math.random;
18
19 import java.io.Serializable;
20
21
22 /** This class implements a powerful pseudo-random number generator
23 * developed by Makoto Matsumoto and Takuji Nishimura during
24 * 1996-1997.
25
26 * <p>This generator features an extremely long period
27 * (2<sup>19937</sup>-1) and 623-dimensional equidistribution up to 32
28 * bits accuracy. The home page for this generator is located at <a
29 * href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html">
30 * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html</a>.</p>
31
32 * <p>This generator is described in a paper by Makoto Matsumoto and
33 * Takuji Nishimura in 1998: <a
34 * href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/mt.pdf">Mersenne
35 * Twister: A 623-Dimensionally Equidistributed Uniform Pseudo-Random
36 * Number Generator</a>, ACM Transactions on Modeling and Computer
37 * Simulation, Vol. 8, No. 1, January 1998, pp 3--30</p>
38
39 * <p>This class is mainly a Java port of the 2002-01-26 version of
40 * the generator written in C by Makoto Matsumoto and Takuji
41 * Nishimura. Here is their original copyright:</p>
42
43 * <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0">
44 * <tr><td>Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
45 * All rights reserved.</td></tr>
46
47 * <tr><td>Redistribution and use in source and binary forms, with or without
48 * modification, are permitted provided that the following conditions
49 * are met:
50 * <ol>
51 * <li>Redistributions of source code must retain the above copyright
52 * notice, this list of conditions and the following disclaimer.</li>
53 * <li>Redistributions in binary form must reproduce the above copyright
54 * notice, this list of conditions and the following disclaimer in the
55 * documentation and/or other materials provided with the distribution.</li>
56 * <li>The names of its contributors may not be used to endorse or promote
57 * products derived from this software without specific prior written
58 * permission.</li>
59 * </ol></td></tr>
60
61 * <tr><td><strong>THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
62 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
63 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
64 * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
65 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS
66 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
67 * OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
68 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
69 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
70 * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
71 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
72 * USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
73 * DAMAGE.</strong></td></tr>
74 * </table>
75
76 * @version $Revision: 797246 $ $Date: 2009-07-23 18:21:46 -0400 (Thu, 23 Jul 2009) $
77 * @since 2.0
78
79 */
80 public class MersenneTwister extends BitsStreamGenerator implements Serializable {
81
82 /** Serializable version identifier. */
83 private static final long serialVersionUID = 8661194735290153518L;
84
85 /** Size of the bytes pool. */
86 private static final int N = 624;
87
88 /** Period second parameter. */
89 private static final int M = 397;
90
91 /** X * MATRIX_A for X = {0, 1}. */
92 private static final int[] MAG01 = { 0x0, 0x9908b0df };
93
94 /** Bytes pool. */
95 private int[] mt;
96
97 /** Current index in the bytes pool. */
98 private int mti;
99
100 /** Creates a new random number generator.
101 * <p>The instance is initialized using the current time as the
102 * seed.</p>
103 */
104 public MersenneTwister() {
105 mt = new int[N];
106 setSeed(System.currentTimeMillis());
107 }
108
109 /** Creates a new random number generator using a single int seed.
110 * @param seed the initial seed (32 bits integer)
111 */
112 public MersenneTwister(int seed) {
113 mt = new int[N];
114 setSeed(seed);
115 }
116
117 /** Creates a new random number generator using an int array seed.
118 * @param seed the initial seed (32 bits integers array), if null
119 * the seed of the generator will be related to the current time
120 */
121 public MersenneTwister(int[] seed) {
122 mt = new int[N];
123 setSeed(seed);
124 }
125
126 /** Creates a new random number generator using a single long seed.
127 * @param seed the initial seed (64 bits integer)
128 */
129 public MersenneTwister(long seed) {
130 mt = new int[N];
131 setSeed(seed);
132 }
133
134 /** Reinitialize the generator as if just built with the given int seed.
135 * <p>The state of the generator is exactly the same as a new
136 * generator built with the same seed.</p>
137 * @param seed the initial seed (32 bits integer)
138 */
139 public void setSeed(int seed) {
140 // we use a long masked by 0xffffffffL as a poor man unsigned int
141 long longMT = seed;
142 mt[0]= (int) longMT;
143 for (mti = 1; mti < N; ++mti) {
144 // See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier.
145 // initializer from the 2002-01-09 C version by Makoto Matsumoto
146 longMT = (1812433253l * (longMT ^ (longMT >> 30)) + mti) & 0xffffffffL;
147 mt[mti]= (int) longMT;
148 }
149 }
150
151 /** Reinitialize the generator as if just built with the given int array seed.
152 * <p>The state of the generator is exactly the same as a new
153 * generator built with the same seed.</p>
154 * @param seed the initial seed (32 bits integers array), if null
155 * the seed of the generator will be related to the current time
156 */
157 public void setSeed(int[] seed) {
158
159 if (seed == null) {
160 setSeed(System.currentTimeMillis());
161 return;
162 }
163
164 setSeed(19650218);
165 int i = 1;
166 int j = 0;
167
168 for (int k = Math.max(N, seed.length); k != 0; k--) {
169 long l0 = (mt[i] & 0x7fffffffl) | ((mt[i] < 0) ? 0x80000000l : 0x0l);
170 long l1 = (mt[i-1] & 0x7fffffffl) | ((mt[i-1] < 0) ? 0x80000000l : 0x0l);
171 long l = (l0 ^ ((l1 ^ (l1 >> 30)) * 1664525l)) + seed[j] + j; // non linear
172 mt[i] = (int) (l & 0xffffffffl);
173 i++; j++;
174 if (i >= N) {
175 mt[0] = mt[N - 1];
176 i = 1;
177 }
178 if (j >= seed.length) {
179 j = 0;
180 }
181 }
182
183 for (int k = N - 1; k != 0; k--) {
184 long l0 = (mt[i] & 0x7fffffffl) | ((mt[i] < 0) ? 0x80000000l : 0x0l);
185 long l1 = (mt[i-1] & 0x7fffffffl) | ((mt[i-1] < 0) ? 0x80000000l : 0x0l);
186 long l = (l0 ^ ((l1 ^ (l1 >> 30)) * 1566083941l)) - i; // non linear
187 mt[i] = (int) (l & 0xffffffffL);
188 i++;
189 if (i >= N) {
190 mt[0] = mt[N - 1];
191 i = 1;
192 }
193 }
194
195 mt[0] = 0x80000000; // MSB is 1; assuring non-zero initial array
196
197 }
198
199 /** Reinitialize the generator as if just built with the given long seed.
200 * <p>The state of the generator is exactly the same as a new
201 * generator built with the same seed.</p>
202 * @param seed the initial seed (64 bits integer)
203 */
204 public void setSeed(long seed) {
205 setSeed(new int[] { (int) (seed >>> 32), (int) (seed & 0xffffffffl) });
206 }
207
208 /** Generate next pseudorandom number.
209 * <p>This method is the core generation algorithm. It is used by all the
210 * public generation methods for the various primitive types {@link
211 * #nextBoolean()}, {@link #nextBytes(byte[])}, {@link #nextDouble()},
212 * {@link #nextFloat()}, {@link #nextGaussian()}, {@link #nextInt()},
213 * {@link #next(int)} and {@link #nextLong()}.</p>
214 * @param bits number of random bits to produce
215 * @return random bits generated
216 */
217 protected int next(int bits) {
218
219 int y;
220
221 if (mti >= N) { // generate N words at one time
222 int mtNext = mt[0];
223 for (int k = 0; k < N - M; ++k) {
224 int mtCurr = mtNext;
225 mtNext = mt[k + 1];
226 y = (mtCurr & 0x80000000) | (mtNext & 0x7fffffff);
227 mt[k] = mt[k + M] ^ (y >>> 1) ^ MAG01[y & 0x1];
228 }
229 for (int k = N - M; k < N - 1; ++k) {
230 int mtCurr = mtNext;
231 mtNext = mt[k + 1];
232 y = (mtCurr & 0x80000000) | (mtNext & 0x7fffffff);
233 mt[k] = mt[k + (M - N)] ^ (y >>> 1) ^ MAG01[y & 0x1];
234 }
235 y = (mtNext & 0x80000000) | (mt[0] & 0x7fffffff);
236 mt[N - 1] = mt[M - 1] ^ (y >>> 1) ^ MAG01[y & 0x1];
237
238 mti = 0;
239 }
240
241 y = mt[mti++];
242
243 // tempering
244 y ^= (y >>> 11);
245 y ^= (y << 7) & 0x9d2c5680;
246 y ^= (y << 15) & 0xefc60000;
247 y ^= (y >>> 18);
248
249 return y >>> (32 - bits);
250
251 }
252
253 }