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.rng.core.source32; 18 19 import java.util.Arrays; 20 import org.apache.commons.rng.core.util.NumberFactory; 21 22 /** 23 * This class implements a powerful pseudo-random number generator 24 * developed by Makoto Matsumoto and Takuji Nishimura during 25 * 1996-1997. 26 * 27 * <p> 28 * This generator features an extremely long period 29 * (2<sup>19937</sup>-1) and 623-dimensional equidistribution up to 30 * 32 bits accuracy. The home page for this generator is located at 31 * <a 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>. 33 * </p> 34 * 35 * <p> 36 * This generator is described in a paper by Makoto Matsumoto and 37 * Takuji Nishimura in 1998: 38 * <a href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/mt.pdf"> 39 * Mersenne Twister: A 623-Dimensionally Equidistributed Uniform Pseudo-Random 40 * Number Generator</a>, 41 * ACM Transactions on Modeling and Computer Simulation, Vol. 8, No. 1, 42 * January 1998, pp 3--30 43 * </p> 44 * 45 * <p> 46 * This class is mainly a Java port of the 47 * <a href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html"> 48 * 2002-01-26 version of the generator</a> written in C by Makoto Matsumoto 49 * and Takuji Nishimura. Here is their original copyright: 50 * </p> 51 * 52 * <table style="background-color: #E0E0E0; width: 80%"> 53 * <caption>Mersenne Twister licence</caption> 54 * <tr><td style="padding: 10px">Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, 55 * All rights reserved.</td></tr> 56 * 57 * <tr><td style="padding: 10px">Redistribution and use in source and binary forms, with or without 58 * modification, are permitted provided that the following conditions 59 * are met: 60 * <ol> 61 * <li>Redistributions of source code must retain the above copyright 62 * notice, this list of conditions and the following disclaimer.</li> 63 * <li>Redistributions in binary form must reproduce the above copyright 64 * notice, this list of conditions and the following disclaimer in the 65 * documentation and/or other materials provided with the distribution.</li> 66 * <li>The names of its contributors may not be used to endorse or promote 67 * products derived from this software without specific prior written 68 * permission.</li> 69 * </ol></td></tr> 70 * 71 * <tr><td style="padding: 10px"><strong>THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND 72 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, 73 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF 74 * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 75 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS 76 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, 77 * OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 78 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 79 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY 80 * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 81 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE 82 * USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH 83 * DAMAGE.</strong></td></tr> 84 * </table> 85 * 86 * @see <a href="https://en.wikipedia.org/wiki/Mersenne_Twister">Mersenne Twister (Wikipedia)</a> 87 * @since 1.0 88 */ 89 public class MersenneTwister extends IntProvider { 90 /** Mask 32 most significant bits. */ 91 private static final long INT_MASK_LONG = 0xffffffffL; 92 /** Most significant w-r bits. */ 93 private static final long UPPER_MASK_LONG = 0x80000000L; 94 /** Least significant r bits. */ 95 private static final long LOWER_MASK_LONG = 0x7fffffffL; 96 /** Most significant w-r bits. */ 97 private static final int UPPER_MASK = 0x80000000; 98 /** Least significant r bits. */ 99 private static final int LOWER_MASK = 0x7fffffff; 100 /** Size of the bytes pool. */ 101 private static final int N = 624; 102 /** Period second parameter. */ 103 private static final int M = 397; 104 /** X * MATRIX_A for X = {0, 1}. */ 105 private static final int[] MAG01 = {0x0, 0x9908b0df}; 106 /** Bytes pool. */ 107 private final int[] mt = new int[N]; 108 /** Current index in the bytes pool. */ 109 private int mti; 110 111 /** 112 * Creates a new random number generator. 113 * 114 * @param seed Initial seed. 115 */ 116 public MersenneTwister(int[] seed) { 117 setSeedInternal(seed); 118 } 119 120 /** {@inheritDoc} */ 121 @Override 122 protected byte[] getStateInternal() { 123 final int[] s = Arrays.copyOf(mt, N + 1); 124 s[N] = mti; 125 126 return composeStateInternal(NumberFactory.makeByteArray(s), 127 super.getStateInternal()); 128 } 129 130 /** {@inheritDoc} */ 131 @Override 132 protected void setStateInternal(byte[] s) { 133 final byte[][] c = splitStateInternal(s, (N + 1) * 4); 134 135 final int[] tmp = NumberFactory.makeIntArray(c[0]); 136 System.arraycopy(tmp, 0, mt, 0, N); 137 mti = tmp[N]; 138 139 super.setStateInternal(c[1]); 140 } 141 142 /** 143 * Initializes the generator with the given seed. 144 * 145 * @param seed Initial seed. 146 */ 147 private void setSeedInternal(int[] seed) { 148 fillStateMersenneTwister(mt, seed); 149 150 // Initial index. 151 mti = N; 152 } 153 154 /** 155 * Utility for wholly filling a {@code state} array with non-zero 156 * bytes, even if the {@code seed} has a smaller size. 157 * The procedure is the one defined by the standard implementation 158 * of the algorithm. 159 * 160 * @param state State to be filled (must be allocated). 161 * @param inputSeed Seed (cannot be {@code null}). 162 */ 163 private static void fillStateMersenneTwister(int[] state, 164 int[] inputSeed) { 165 // Accept empty seed. 166 final int[] seed = (inputSeed.length == 0) ? new int[1] : inputSeed; 167 168 initializeState(state); 169 170 final int nextIndex = mixSeedAndState(state, seed); 171 172 mixState(state, nextIndex); 173 174 state[0] = (int) UPPER_MASK_LONG; // MSB is 1, ensuring non-zero initial array. 175 } 176 177 /** 178 * Fill the state using a defined pseudo-random sequence. 179 * 180 * @param state State to be filled (must be allocated). 181 */ 182 private static void initializeState(int[] state) { 183 long mt = 19650218 & INT_MASK_LONG; 184 state[0] = (int) mt; 185 for (int i = 1; i < state.length; i++) { 186 mt = (1812433253L * (mt ^ (mt >> 30)) + i) & INT_MASK_LONG; 187 state[i] = (int) mt; 188 } 189 } 190 191 /** 192 * Mix the seed into the state using a non-linear combination. The procedure 193 * uses {@code k} steps where {@code k = max(state.length, seed.length)}. If 194 * the seed is smaller than the state it is wrapped to obtain enough values. 195 * If the seed is larger than the state then the procedure visits entries in 196 * the state multiple times. 197 * 198 * <p>Returns the index immediately after the most recently visited position 199 * in the state array.</p> 200 * 201 * @param state State to be filled (must be allocated). 202 * @param seed Seed (must be at least length 1). 203 * @return the next index 204 */ 205 private static int mixSeedAndState(int[] state, final int[] seed) { 206 final int stateSize = state.length; 207 208 int i = 1; 209 int j = 0; 210 211 for (int k = Math.max(stateSize, seed.length); k > 0; k--) { 212 final long a = (state[i] & LOWER_MASK_LONG) | ((state[i] < 0) ? UPPER_MASK_LONG : 0); 213 final long b = (state[i - 1] & LOWER_MASK_LONG) | ((state[i - 1] < 0) ? UPPER_MASK_LONG : 0); 214 final long c = (a ^ ((b ^ (b >> 30)) * 1664525L)) + seed[j] + j; // Non linear. 215 state[i] = (int) (c & INT_MASK_LONG); 216 i++; 217 j++; 218 if (i >= stateSize) { 219 state[0] = state[stateSize - 1]; 220 i = 1; 221 } 222 if (j >= seed.length) { 223 j = 0; 224 } 225 } 226 return i; 227 } 228 229 /** 230 * Mix each position of the state using a non-linear combination. The 231 * procedure starts from the specified index in the state array and wraps 232 * iteration through the array if required. 233 * 234 * @param state State to be filled (must be allocated). 235 * @param startIndex The index to begin within the state array. 236 */ 237 private static void mixState(int[] state, int startIndex) { 238 final int stateSize = state.length; 239 240 int i = startIndex; 241 for (int k = stateSize - 1; k > 0; k--) { 242 final long a = (state[i] & LOWER_MASK_LONG) | ((state[i] < 0) ? UPPER_MASK_LONG : 0); 243 final long b = (state[i - 1] & LOWER_MASK_LONG) | ((state[i - 1] < 0) ? UPPER_MASK_LONG : 0); 244 final long c = (a ^ ((b ^ (b >> 30)) * 1566083941L)) - i; // Non linear. 245 state[i] = (int) (c & INT_MASK_LONG); 246 i++; 247 if (i >= stateSize) { 248 state[0] = state[stateSize - 1]; 249 i = 1; 250 } 251 } 252 } 253 254 /** {@inheritDoc} */ 255 @Override 256 public int next() { 257 int y; 258 259 if (mti >= N) { // Generate N words at one time. 260 int mtNext = mt[0]; 261 for (int k = 0; k < N - M; ++k) { 262 final int mtCurr = mtNext; 263 mtNext = mt[k + 1]; 264 y = (mtCurr & UPPER_MASK) | (mtNext & LOWER_MASK); 265 mt[k] = mt[k + M] ^ (y >>> 1) ^ MAG01[y & 1]; 266 } 267 for (int k = N - M; k < N - 1; ++k) { 268 final int mtCurr = mtNext; 269 mtNext = mt[k + 1]; 270 y = (mtCurr & UPPER_MASK) | (mtNext & LOWER_MASK); 271 mt[k] = mt[k + (M - N)] ^ (y >>> 1) ^ MAG01[y & 1]; 272 } 273 y = (mtNext & UPPER_MASK) | (mt[0] & LOWER_MASK); 274 mt[N - 1] = mt[M - 1] ^ (y >>> 1) ^ MAG01[y & 1]; 275 276 mti = 0; 277 } 278 279 y = mt[mti++]; 280 281 // Tempering. 282 y ^= y >>> 11; 283 y ^= (y << 7) & 0x9d2c5680; 284 y ^= (y << 15) & 0xefc60000; 285 y ^= y >>> 18; 286 287 return y; 288 } 289 }