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