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}