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