View Javadoc
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 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 }