Source code

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 style="background-color: #E0E0E0; width: 80%">
053 * <caption>Mersenne Twister licence</caption>
054 * <tr><td style="padding: 10px">Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
055 *     All rights reserved.</td></tr>
056 *
057 * <tr><td style="padding: 10px">Redistribution and use in source and binary forms, with or without
058 * modification, are permitted provided that the following conditions
059 * are met:
060 * <ol>
061 *   <li>Redistributions of source code must retain the above copyright
062 *       notice, this list of conditions and the following disclaimer.</li>
063 *   <li>Redistributions in binary form must reproduce the above copyright
064 *       notice, this list of conditions and the following disclaimer in the
065 *       documentation and/or other materials provided with the distribution.</li>
066 *   <li>The names of its contributors may not be used to endorse or promote
067 *       products derived from this software without specific prior written
068 *       permission.</li>
069 * </ol></td></tr>
070 *
071 * <tr><td style="padding: 10px"><strong>THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
072 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
073 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
074 * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
075 * DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS
076 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
077 * OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
078 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
079 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
080 * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
081 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
082 * USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
083 * DAMAGE.</strong></td></tr>
084 * </table>
085 *
086 * @see <a href="https://en.wikipedia.org/wiki/Mersenne_Twister">Mersenne Twister (Wikipedia)</a>
087 * @since 1.0
088 */
089public class MersenneTwister extends IntProvider {
090    /** Mask 32 most significant bits. */
091    private static final long INT_MASK_LONG = 0xffffffffL;
092    /** Most significant w-r bits. */
093    private static final long UPPER_MASK_LONG = 0x80000000L;
094    /** Least significant r bits. */
095    private static final long LOWER_MASK_LONG = 0x7fffffffL;
096    /** Most significant w-r bits. */
097    private static final int UPPER_MASK = 0x80000000;
098    /** Least significant r bits. */
099    private static final int LOWER_MASK = 0x7fffffff;
100    /** Size of the bytes pool. */
101    private static final int N = 624;
102    /** Period second parameter. */
103    private static final int M = 397;
104    /** X * MATRIX_A for X = {0, 1}. */
105    private static final int[] MAG01 = {0x0, 0x9908b0df};
106    /** Bytes pool. */
107    private final int[] mt = new int[N];
108    /** Current index in the bytes pool. */
109    private int mti;
110
111    /**
112     * Creates a new random number generator.
113     *
114     * @param seed Initial seed.
115     */
116    public MersenneTwister(int[] seed) {
117        setSeedInternal(seed);
118    }
119
120    /** {@inheritDoc} */
121    @Override
122    protected byte[] getStateInternal() {
123        final int[] s = Arrays.copyOf(mt, N + 1);
124        s[N] = mti;
125
126        return composeStateInternal(NumberFactory.makeByteArray(s),
127                                    super.getStateInternal());
128    }
129
130    /** {@inheritDoc} */
131    @Override
132    protected void setStateInternal(byte[] s) {
133        final byte[][] c = splitStateInternal(s, (N + 1) * 4);
134
135        final int[] tmp = NumberFactory.makeIntArray(c[0]);
136        System.arraycopy(tmp, 0, mt, 0, N);
137        mti = tmp[N];
138
139        super.setStateInternal(c[1]);
140    }
141
142    /**
143     * Initializes the generator with the given seed.
144     *
145     * @param seed Initial seed.
146     */
147    private void setSeedInternal(int[] seed) {
148        fillStateMersenneTwister(mt, seed);
149
150        // Initial index.
151        mti = N;
152    }
153
154    /**
155     * Utility for wholly filling a {@code state} array with non-zero
156     * bytes, even if the {@code seed} has a smaller size.
157     * The procedure is the one defined by the standard implementation
158     * of the algorithm.
159     *
160     * @param state State to be filled (must be allocated).
161     * @param inputSeed Seed (cannot be {@code null}).
162     */
163    private static void fillStateMersenneTwister(int[] state,
164                                                 int[] inputSeed) {
165        // Accept empty seed.
166        final int[] seed = (inputSeed.length == 0) ? new int[1] : inputSeed;
167
168        initializeState(state);
169
170        final int nextIndex = mixSeedAndState(state, seed);
171
172        mixState(state, nextIndex);
173
174        state[0] = (int) UPPER_MASK_LONG; // MSB is 1, ensuring non-zero initial array.
175    }
176
177    /**
178     * Fill the state using a defined pseudo-random sequence.
179     *
180     * @param state State to be filled (must be allocated).
181     */
182    private static void initializeState(int[] state) {
183        long mt = 19650218 & INT_MASK_LONG;
184        state[0] = (int) mt;
185        for (int i = 1; i < state.length; i++) {
186            mt = (1812433253L * (mt ^ (mt >> 30)) + i) & INT_MASK_LONG;
187            state[i] = (int) mt;
188        }
189    }
190
191    /**
192     * Mix the seed into the state using a non-linear combination. The procedure
193     * uses {@code k} steps where {@code k = max(state.length, seed.length)}. If
194     * the seed is smaller than the state it is wrapped to obtain enough values.
195     * If the seed is larger than the state then the procedure visits entries in
196     * the state multiple times.
197     *
198     * <p>Returns the index immediately after the most recently visited position
199     * in the state array.</p>
200     *
201     * @param state State to be filled (must be allocated).
202     * @param seed Seed (must be at least length 1).
203     * @return the next index
204     */
205    private static int mixSeedAndState(int[] state, final int[] seed) {
206        final int stateSize = state.length;
207
208        int i = 1;
209        int j = 0;
210
211        for (int k = Math.max(stateSize, seed.length); k > 0; k--) {
212            final long a = (state[i] & LOWER_MASK_LONG) | ((state[i] < 0) ? UPPER_MASK_LONG : 0);
213            final long b = (state[i - 1] & LOWER_MASK_LONG) | ((state[i - 1] < 0) ? UPPER_MASK_LONG : 0);
214            final long c = (a ^ ((b ^ (b >> 30)) * 1664525L)) + seed[j] + j; // Non linear.
215            state[i] = (int) (c & INT_MASK_LONG);
216            i++;
217            j++;
218            if (i >= stateSize) {
219                state[0] = state[stateSize - 1];
220                i = 1;
221            }
222            if (j >= seed.length) {
223                j = 0;
224            }
225        }
226        return i;
227    }
228
229    /**
230     * Mix each position of the state using a non-linear combination. The
231     * procedure starts from the specified index in the state array and wraps
232     * iteration through the array if required.
233     *
234     * @param state State to be filled (must be allocated).
235     * @param startIndex The index to begin within the state array.
236     */
237    private static void mixState(int[] state, int startIndex) {
238        final int stateSize = state.length;
239
240        int i = startIndex;
241        for (int k = stateSize - 1; k > 0; k--) {
242            final long a = (state[i] & LOWER_MASK_LONG) | ((state[i] < 0) ? UPPER_MASK_LONG : 0);
243            final long b = (state[i - 1] & LOWER_MASK_LONG) | ((state[i - 1] < 0) ? UPPER_MASK_LONG : 0);
244            final long c = (a ^ ((b ^ (b >> 30)) * 1566083941L)) - i; // Non linear.
245            state[i] = (int) (c & INT_MASK_LONG);
246            i++;
247            if (i >= stateSize) {
248                state[0] = state[stateSize - 1];
249                i = 1;
250            }
251        }
252    }
253
254    /** {@inheritDoc} */
255    @Override
256    public int next() {
257        int y;
258
259        if (mti >= N) { // Generate N words at one time.
260            int mtNext = mt[0];
261            for (int k = 0; k < N - M; ++k) {
262                final int mtCurr = mtNext;
263                mtNext = mt[k + 1];
264                y = (mtCurr & UPPER_MASK) | (mtNext & LOWER_MASK);
265                mt[k] = mt[k + M] ^ (y >>> 1) ^ MAG01[y & 1];
266            }
267            for (int k = N - M; k < N - 1; ++k) {
268                final int mtCurr = mtNext;
269                mtNext = mt[k + 1];
270                y = (mtCurr & UPPER_MASK) | (mtNext & LOWER_MASK);
271                mt[k] = mt[k + (M - N)] ^ (y >>> 1) ^ MAG01[y & 1];
272            }
273            y = (mtNext & UPPER_MASK) | (mt[0] & LOWER_MASK);
274            mt[N - 1] = mt[M - 1] ^ (y >>> 1) ^ MAG01[y & 1];
275
276            mti = 0;
277        }
278
279        y = mt[mti++];
280
281        // Tempering.
282        y ^=  y >>> 11;
283        y ^= (y << 7) & 0x9d2c5680;
284        y ^= (y << 15) & 0xefc60000;
285        y ^=  y >>> 18;
286
287        return y;
288    }
289}