001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017 package org.apache.commons.math3.random;
018
019 import java.io.Serializable;
020
021 import org.apache.commons.math3.util.FastMath;
022
023
024 /** This class implements a powerful pseudo-random number generator
025 * developed by Makoto Matsumoto and Takuji Nishimura during
026 * 1996-1997.
027
028 * <p>This generator features an extremely long period
029 * (2<sup>19937</sup>-1) and 623-dimensional equidistribution up to 32
030 * bits accuracy. The home page for this generator is located at <a
031 * href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html">
032 * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html</a>.</p>
033
034 * <p>This generator is described in a paper by Makoto Matsumoto and
035 * Takuji Nishimura in 1998: <a
036 * href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/mt.pdf">Mersenne
037 * Twister: A 623-Dimensionally Equidistributed Uniform Pseudo-Random
038 * Number Generator</a>, ACM Transactions on Modeling and Computer
039 * Simulation, Vol. 8, No. 1, January 1998, pp 3--30</p>
040
041 * <p>This class is mainly a Java port of the 2002-01-26 version of
042 * the generator written in C by Makoto Matsumoto and Takuji
043 * Nishimura. Here is their original copyright:</p>
044
045 * <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0">
046 * <tr><td>Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
047 * All rights reserved.</td></tr>
048
049 * <tr><td>Redistribution and use in source and binary forms, with or without
050 * modification, are permitted provided that the following conditions
051 * are met:
052 * <ol>
053 * <li>Redistributions of source code must retain the above copyright
054 * notice, this list of conditions and the following disclaimer.</li>
055 * <li>Redistributions in binary form must reproduce the above copyright
056 * notice, this list of conditions and the following disclaimer in the
057 * documentation and/or other materials provided with the distribution.</li>
058 * <li>The names of its contributors may not be used to endorse or promote
059 * products derived from this software without specific prior written
060 * permission.</li>
061 * </ol></td></tr>
062
063 * <tr><td><strong>THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
064 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
065 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
066 * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
067 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS
068 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
069 * OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
070 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
071 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
072 * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
073 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
074 * USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
075 * DAMAGE.</strong></td></tr>
076 * </table>
077
078 * @version $Id: MersenneTwister.java 1416643 2012-12-03 19:37:14Z tn $
079 * @since 2.0
080
081 */
082 public class MersenneTwister extends BitsStreamGenerator implements Serializable {
083
084 /** Serializable version identifier. */
085 private static final long serialVersionUID = 8661194735290153518L;
086
087 /** Size of the bytes pool. */
088 private static final int N = 624;
089
090 /** Period second parameter. */
091 private static final int M = 397;
092
093 /** X * MATRIX_A for X = {0, 1}. */
094 private static final int[] MAG01 = { 0x0, 0x9908b0df };
095
096 /** Bytes pool. */
097 private int[] mt;
098
099 /** Current index in the bytes pool. */
100 private int mti;
101
102 /** Creates a new random number generator.
103 * <p>The instance is initialized using the current time plus the
104 * system identity hash code of this instance as the seed.</p>
105 */
106 public MersenneTwister() {
107 mt = new int[N];
108 setSeed(System.currentTimeMillis() + System.identityHashCode(this));
109 }
110
111 /** Creates a new random number generator using a single int seed.
112 * @param seed the initial seed (32 bits integer)
113 */
114 public MersenneTwister(int seed) {
115 mt = new int[N];
116 setSeed(seed);
117 }
118
119 /** Creates a new random number generator using an int array seed.
120 * @param seed the initial seed (32 bits integers array), if null
121 * the seed of the generator will be related to the current time
122 */
123 public MersenneTwister(int[] seed) {
124 mt = new int[N];
125 setSeed(seed);
126 }
127
128 /** Creates a new random number generator using a single long seed.
129 * @param seed the initial seed (64 bits integer)
130 */
131 public MersenneTwister(long seed) {
132 mt = new int[N];
133 setSeed(seed);
134 }
135
136 /** Reinitialize the generator as if just built with the given int seed.
137 * <p>The state of the generator is exactly the same as a new
138 * generator built with the same seed.</p>
139 * @param seed the initial seed (32 bits integer)
140 */
141 @Override
142 public void setSeed(int seed) {
143 // we use a long masked by 0xffffffffL as a poor man unsigned int
144 long longMT = seed;
145 // NB: unlike original C code, we are working with java longs, the cast below makes masking unnecessary
146 mt[0]= (int) longMT;
147 for (mti = 1; mti < N; ++mti) {
148 // See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier.
149 // initializer from the 2002-01-09 C version by Makoto Matsumoto
150 longMT = (1812433253l * (longMT ^ (longMT >> 30)) + mti) & 0xffffffffL;
151 mt[mti]= (int) longMT;
152 }
153
154 clear(); // Clear normal deviate cache
155 }
156
157 /** Reinitialize the generator as if just built with the given int array seed.
158 * <p>The state of the generator is exactly the same as a new
159 * generator built with the same seed.</p>
160 * @param seed the initial seed (32 bits integers array), if null
161 * the seed of the generator will be the current system time plus the
162 * system identity hash code of this instance
163 */
164 @Override
165 public void setSeed(int[] seed) {
166
167 if (seed == null) {
168 setSeed(System.currentTimeMillis() + System.identityHashCode(this));
169 return;
170 }
171
172 setSeed(19650218);
173 int i = 1;
174 int j = 0;
175
176 for (int k = FastMath.max(N, seed.length); k != 0; k--) {
177 long l0 = (mt[i] & 0x7fffffffl) | ((mt[i] < 0) ? 0x80000000l : 0x0l);
178 long l1 = (mt[i-1] & 0x7fffffffl) | ((mt[i-1] < 0) ? 0x80000000l : 0x0l);
179 long l = (l0 ^ ((l1 ^ (l1 >> 30)) * 1664525l)) + seed[j] + j; // non linear
180 mt[i] = (int) (l & 0xffffffffl);
181 i++; j++;
182 if (i >= N) {
183 mt[0] = mt[N - 1];
184 i = 1;
185 }
186 if (j >= seed.length) {
187 j = 0;
188 }
189 }
190
191 for (int k = N - 1; k != 0; k--) {
192 long l0 = (mt[i] & 0x7fffffffl) | ((mt[i] < 0) ? 0x80000000l : 0x0l);
193 long l1 = (mt[i-1] & 0x7fffffffl) | ((mt[i-1] < 0) ? 0x80000000l : 0x0l);
194 long l = (l0 ^ ((l1 ^ (l1 >> 30)) * 1566083941l)) - i; // non linear
195 mt[i] = (int) (l & 0xffffffffL);
196 i++;
197 if (i >= N) {
198 mt[0] = mt[N - 1];
199 i = 1;
200 }
201 }
202
203 mt[0] = 0x80000000; // MSB is 1; assuring non-zero initial array
204
205 clear(); // Clear normal deviate cache
206
207 }
208
209 /** Reinitialize the generator as if just built with the given long seed.
210 * <p>The state of the generator is exactly the same as a new
211 * generator built with the same seed.</p>
212 * @param seed the initial seed (64 bits integer)
213 */
214 @Override
215 public void setSeed(long seed) {
216 setSeed(new int[] { (int) (seed >>> 32), (int) (seed & 0xffffffffl) });
217 }
218
219 /** Generate next pseudorandom number.
220 * <p>This method is the core generation algorithm. It is used by all the
221 * public generation methods for the various primitive types {@link
222 * #nextBoolean()}, {@link #nextBytes(byte[])}, {@link #nextDouble()},
223 * {@link #nextFloat()}, {@link #nextGaussian()}, {@link #nextInt()},
224 * {@link #next(int)} and {@link #nextLong()}.</p>
225 * @param bits number of random bits to produce
226 * @return random bits generated
227 */
228 @Override
229 protected int next(int bits) {
230
231 int y;
232
233 if (mti >= N) { // generate N words at one time
234 int mtNext = mt[0];
235 for (int k = 0; k < N - M; ++k) {
236 int mtCurr = mtNext;
237 mtNext = mt[k + 1];
238 y = (mtCurr & 0x80000000) | (mtNext & 0x7fffffff);
239 mt[k] = mt[k + M] ^ (y >>> 1) ^ MAG01[y & 0x1];
240 }
241 for (int k = N - M; k < N - 1; ++k) {
242 int mtCurr = mtNext;
243 mtNext = mt[k + 1];
244 y = (mtCurr & 0x80000000) | (mtNext & 0x7fffffff);
245 mt[k] = mt[k + (M - N)] ^ (y >>> 1) ^ MAG01[y & 0x1];
246 }
247 y = (mtNext & 0x80000000) | (mt[0] & 0x7fffffff);
248 mt[N - 1] = mt[M - 1] ^ (y >>> 1) ^ MAG01[y & 0x1];
249
250 mti = 0;
251 }
252
253 y = mt[mti++];
254
255 // tempering
256 y ^= y >>> 11;
257 y ^= (y << 7) & 0x9d2c5680;
258 y ^= (y << 15) & 0xefc60000;
259 y ^= y >>> 18;
260
261 return y >>> (32 - bits);
262
263 }
264
265 }