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 */ 017package org.apache.commons.rng.core.source32; 018 019import java.util.Arrays; 020import org.apache.commons.rng.core.util.NumberFactory; 021 022/** 023 * This class implements a powerful pseudo-random number generator 024 * developed by Makoto Matsumoto and Takuji Nishimura during 025 * 1996-1997. 026 * 027 * <p> 028 * This generator features an extremely long period 029 * (2<sup>19937</sup>-1) and 623-dimensional equidistribution up to 030 * 32 bits accuracy. The home page for this generator is located at 031 * <a 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>. 033 * </p> 034 * 035 * <p> 036 * This generator is described in a paper by Makoto Matsumoto and 037 * Takuji Nishimura in 1998: 038 * <a href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/mt.pdf"> 039 * Mersenne Twister: A 623-Dimensionally Equidistributed Uniform Pseudo-Random 040 * Number Generator</a>, 041 * ACM Transactions on Modeling and Computer Simulation, Vol. 8, No. 1, 042 * January 1998, pp 3--30 043 * </p> 044 * 045 * <p> 046 * This class is mainly a Java port of the 047 * <a href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html"> 048 * 2002-01-26 version of the generator</a> written in C by Makoto Matsumoto 049 * and Takuji Nishimura. Here is their original copyright: 050 * </p> 051 * 052 * <table border="0" width="80%" cellpadding="10" style="background-color: #E0E0E0" summary="Mersenne Twister licence"> 053 * <tr><td>Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, 054 * All rights reserved.</td></tr> 055 * 056 * <tr><td>Redistribution and use in source and binary forms, with or without 057 * modification, are permitted provided that the following conditions 058 * are met: 059 * <ol> 060 * <li>Redistributions of source code must retain the above copyright 061 * notice, this list of conditions and the following disclaimer.</li> 062 * <li>Redistributions in binary form must reproduce the above copyright 063 * notice, this list of conditions and the following disclaimer in the 064 * documentation and/or other materials provided with the distribution.</li> 065 * <li>The names of its contributors may not be used to endorse or promote 066 * products derived from this software without specific prior written 067 * permission.</li> 068 * </ol></td></tr> 069 * 070 * <tr><td><strong>THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND 071 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, 072 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF 073 * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 074 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS 075 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, 076 * OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 077 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 078 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY 079 * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 080 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE 081 * USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH 082 * DAMAGE.</strong></td></tr> 083 * </table> 084 * 085 * @see <a href="https://en.wikipedia.org/wiki/Mersenne_Twister">Mersenne Twister (Wikipedia)</a> 086 * @since 1.0 087 */ 088public class MersenneTwister extends IntProvider { 089 /** Mask 32 most significant bits. */ 090 private static final long INT_MASK_LONG = 0xffffffffL; 091 /** Most significant w-r bits. */ 092 private static final long UPPER_MASK_LONG = 0x80000000L; 093 /** Least significant r bits. */ 094 private static final long LOWER_MASK_LONG = 0x7fffffffL; 095 /** Most significant w-r bits. */ 096 private static final int UPPER_MASK = 0x80000000; 097 /** Least significant r bits. */ 098 private static final int LOWER_MASK = 0x7fffffff; 099 /** 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 NumberFactory.makeByteArray(s); 126 } 127 128 /** {@inheritDoc} */ 129 @Override 130 protected void setStateInternal(byte[] s) { 131 checkStateSize(s, (N + 1) * 4); 132 133 final int[] tmp = NumberFactory.makeIntArray(s); 134 System.arraycopy(tmp, 0, mt, 0, N); 135 mti = tmp[N]; 136 } 137 138 /** 139 * Initializes the generator with the given seed. 140 * 141 * @param seed Initial seed. 142 */ 143 private void setSeedInternal(int[] seed) { 144 fillStateMersenneTwister(mt, seed); 145 146 // Initial index. 147 mti = N; 148 } 149 150 /** 151 * Utility for wholly filling a {@code state} array with non-zero 152 * bytes, even if the {@code seed} has a smaller size. 153 * The procedure is the one defined by the standard implementation 154 * of the algorithm. 155 * 156 * @param state State to be filled (must be allocated). 157 * @param seed Seed (cannot be {@code null}). 158 */ 159 private static void fillStateMersenneTwister(int[] state, 160 int[] seed) { 161 if (seed.length == 0) { 162 // Accept empty seed. 163 seed = new int[1]; 164 } 165 166 final int stateSize = state.length; 167 168 long mt = 19650218 & INT_MASK_LONG; 169 state[0] = (int) mt; 170 for (int i = 1; i < stateSize; i++) { 171 mt = (1812433253L * (mt ^ (mt >> 30)) + i) & INT_MASK_LONG; 172 state[i] = (int) mt; 173 } 174 175 int i = 1; 176 int j = 0; 177 178 for (int k = Math.max(stateSize, seed.length); k > 0; k--) { 179 final long a = (state[i] & LOWER_MASK_LONG) | ((state[i] < 0) ? UPPER_MASK_LONG : 0); 180 final long b = (state[i - 1] & LOWER_MASK_LONG) | ((state[i - 1] < 0) ? UPPER_MASK_LONG : 0); 181 final long c = (a ^ ((b ^ (b >> 30)) * 1664525L)) + seed[j] + j; // Non linear. 182 state[i] = (int) (c & INT_MASK_LONG); 183 i++; 184 j++; 185 if (i >= stateSize) { 186 state[0] = state[stateSize - 1]; 187 i = 1; 188 } 189 if (j >= seed.length) { 190 j = 0; 191 } 192 } 193 194 for (int k = stateSize - 1; k > 0; k--) { 195 final long a = (state[i] & LOWER_MASK_LONG) | ((state[i] < 0) ? UPPER_MASK_LONG : 0); 196 final long b = (state[i - 1] & LOWER_MASK_LONG) | ((state[i - 1] < 0) ? UPPER_MASK_LONG : 0); 197 final long c = (a ^ ((b ^ (b >> 30)) * 1566083941L)) - i; // Non linear. 198 state[i] = (int) (c & INT_MASK_LONG); 199 i++; 200 if (i >= stateSize) { 201 state[0] = state[stateSize - 1]; 202 i = 1; 203 } 204 } 205 206 state[0] = (int) UPPER_MASK_LONG; // MSB is 1, ensuring non-zero initial array. 207 } 208 209 /** {@inheritDoc} */ 210 @Override 211 public int next() { 212 int y; 213 214 if (mti >= N) { // Generate N words at one time. 215 int mtNext = mt[0]; 216 for (int k = 0; k < N - M; ++k) { 217 int mtCurr = mtNext; 218 mtNext = mt[k + 1]; 219 y = (mtCurr & UPPER_MASK) | (mtNext & LOWER_MASK); 220 mt[k] = mt[k + M] ^ (y >>> 1) ^ MAG01[y & 1]; 221 } 222 for (int k = N - M; k < N - 1; ++k) { 223 int mtCurr = mtNext; 224 mtNext = mt[k + 1]; 225 y = (mtCurr & UPPER_MASK) | (mtNext & LOWER_MASK); 226 mt[k] = mt[k + (M - N)] ^ (y >>> 1) ^ MAG01[y & 1]; 227 } 228 y = (mtNext & UPPER_MASK) | (mt[0] & LOWER_MASK); 229 mt[N - 1] = mt[M - 1] ^ (y >>> 1) ^ MAG01[y & 1]; 230 231 mti = 0; 232 } 233 234 y = mt[mti++]; 235 236 // Tempering. 237 y ^= y >>> 11; 238 y ^= (y << 7) & 0x9d2c5680; 239 y ^= (y << 15) & 0xefc60000; 240 y ^= y >>> 18; 241 242 return y; 243 } 244}