Raymond | dee0849 | 2015-04-02 10:43:13 -0700 | [diff] [blame] | 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 | */ |
| 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 | */ |
| 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 | 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 |
| 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 |
| 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 |
| 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 | } |