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 border="0" width="80%" cellpadding="10" style="background-color: #E0E0E0" summary="Mersenne Twister licence"> 53 * <tr><td>Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, 54 * All rights reserved.</td></tr> 55 * 56 * <tr><td>Redistribution and use in source and binary forms, with or without 57 * modification, are permitted provided that the following conditions 58 * are met: 59 * <ol> 60 * <li>Redistributions of source code must retain the above copyright 61 * notice, this list of conditions and the following disclaimer.</li> 62 * <li>Redistributions in binary form must reproduce the above copyright 63 * notice, this list of conditions and the following disclaimer in the 64 * documentation and/or other materials provided with the distribution.</li> 65 * <li>The names of its contributors may not be used to endorse or promote 66 * products derived from this software without specific prior written 67 * permission.</li> 68 * </ol></td></tr> 69 * 70 * <tr><td><strong>THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND 71 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, 72 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF 73 * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 74 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS 75 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, 76 * OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 77 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 78 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY 79 * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 80 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE 81 * USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH 82 * DAMAGE.</strong></td></tr> 83 * </table> 84 * 85 * @see <a href="https://en.wikipedia.org/wiki/Mersenne_Twister">Mersenne Twister (Wikipedia)</a> 86 * @since 1.0 87 */ 88 public class MersenneTwister extends IntProvider { 89 /** Mask 32 most significant bits. */ 90 private static final long INT_MASK_LONG = 0xffffffffL; 91 /** Most significant w-r bits. */ 92 private static final long UPPER_MASK_LONG = 0x80000000L; 93 /** Least significant r bits. */ 94 private static final long LOWER_MASK_LONG = 0x7fffffffL; 95 /** Most significant w-r bits. */ 96 private static final int UPPER_MASK = 0x80000000; 97 /** Least significant r bits. */ 98 private static final int LOWER_MASK = 0x7fffffff; 99 /** Size of the bytes pool. */ 100 private static final int N = 624; 101 /** Period second parameter. */ 102 private static final int M = 397; 103 /** X * MATRIX_A for X = {0, 1}. */ 104 private static final int[] MAG01 = { 0x0, 0x9908b0df }; 105 /** Bytes pool. */ 106 private int[] mt = new int[N]; 107 /** Current index in the bytes pool. */ 108 private int mti; 109 110 /** 111 * Creates a new random number generator. 112 * 113 * @param seed Initial seed. 114 */ 115 public MersenneTwister(int[] seed) { 116 setSeedInternal(seed); 117 } 118 119 /** {@inheritDoc} */ 120 @Override 121 protected byte[] getStateInternal() { 122 final int[] s = Arrays.copyOf(mt, N + 1); 123 s[N] = mti; 124 125 return composeStateInternal(NumberFactory.makeByteArray(s), 126 super.getStateInternal()); 127 } 128 129 /** {@inheritDoc} */ 130 @Override 131 protected void setStateInternal(byte[] s) { 132 final byte[][] c = splitStateInternal(s, (N + 1) * 4); 133 134 final int[] tmp = NumberFactory.makeIntArray(c[0]); 135 System.arraycopy(tmp, 0, mt, 0, N); 136 mti = tmp[N]; 137 138 super.setStateInternal(c[1]); 139 } 140 141 /** 142 * Initializes the generator with the given seed. 143 * 144 * @param seed Initial seed. 145 */ 146 private void setSeedInternal(int[] seed) { 147 fillStateMersenneTwister(mt, seed); 148 149 // Initial index. 150 mti = N; 151 } 152 153 /** 154 * Utility for wholly filling a {@code state} array with non-zero 155 * bytes, even if the {@code seed} has a smaller size. 156 * The procedure is the one defined by the standard implementation 157 * of the algorithm. 158 * 159 * @param state State to be filled (must be allocated). 160 * @param seed Seed (cannot be {@code null}). 161 */ 162 private static void fillStateMersenneTwister(int[] state, 163 int[] seed) { 164 if (seed.length == 0) { 165 // Accept empty seed. 166 seed = new int[1]; 167 } 168 169 final int stateSize = state.length; 170 171 long mt = 19650218 & INT_MASK_LONG; 172 state[0] = (int) mt; 173 for (int i = 1; i < stateSize; i++) { 174 mt = (1812433253L * (mt ^ (mt >> 30)) + i) & INT_MASK_LONG; 175 state[i] = (int) mt; 176 } 177 178 int i = 1; 179 int j = 0; 180 181 for (int k = Math.max(stateSize, seed.length); k > 0; k--) { 182 final long a = (state[i] & LOWER_MASK_LONG) | ((state[i] < 0) ? UPPER_MASK_LONG : 0); 183 final long b = (state[i - 1] & LOWER_MASK_LONG) | ((state[i - 1] < 0) ? UPPER_MASK_LONG : 0); 184 final long c = (a ^ ((b ^ (b >> 30)) * 1664525L)) + seed[j] + j; // Non linear. 185 state[i] = (int) (c & INT_MASK_LONG); 186 i++; 187 j++; 188 if (i >= stateSize) { 189 state[0] = state[stateSize - 1]; 190 i = 1; 191 } 192 if (j >= seed.length) { 193 j = 0; 194 } 195 } 196 197 for (int k = stateSize - 1; k > 0; k--) { 198 final long a = (state[i] & LOWER_MASK_LONG) | ((state[i] < 0) ? UPPER_MASK_LONG : 0); 199 final long b = (state[i - 1] & LOWER_MASK_LONG) | ((state[i - 1] < 0) ? UPPER_MASK_LONG : 0); 200 final long c = (a ^ ((b ^ (b >> 30)) * 1566083941L)) - i; // Non linear. 201 state[i] = (int) (c & INT_MASK_LONG); 202 i++; 203 if (i >= stateSize) { 204 state[0] = state[stateSize - 1]; 205 i = 1; 206 } 207 } 208 209 state[0] = (int) UPPER_MASK_LONG; // MSB is 1, ensuring non-zero initial array. 210 } 211 212 /** {@inheritDoc} */ 213 @Override 214 public int next() { 215 int y; 216 217 if (mti >= N) { // Generate N words at one time. 218 int mtNext = mt[0]; 219 for (int k = 0; k < N - M; ++k) { 220 int mtCurr = mtNext; 221 mtNext = mt[k + 1]; 222 y = (mtCurr & UPPER_MASK) | (mtNext & LOWER_MASK); 223 mt[k] = mt[k + M] ^ (y >>> 1) ^ MAG01[y & 1]; 224 } 225 for (int k = N - M; k < N - 1; ++k) { 226 int mtCurr = mtNext; 227 mtNext = mt[k + 1]; 228 y = (mtCurr & UPPER_MASK) | (mtNext & LOWER_MASK); 229 mt[k] = mt[k + (M - N)] ^ (y >>> 1) ^ MAG01[y & 1]; 230 } 231 y = (mtNext & UPPER_MASK) | (mt[0] & LOWER_MASK); 232 mt[N - 1] = mt[M - 1] ^ (y >>> 1) ^ MAG01[y & 1]; 233 234 mti = 0; 235 } 236 237 y = mt[mti++]; 238 239 // Tempering. 240 y ^= y >>> 11; 241 y ^= (y << 7) & 0x9d2c5680; 242 y ^= (y << 15) & 0xefc60000; 243 y ^= y >>> 18; 244 245 return y; 246 } 247 }