MersenneTwister.java

  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;

  19. import java.util.Arrays;
  20. import org.apache.commons.rng.core.util.NumberFactory;

  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 style="background-color: #E0E0E0; width: 80%">
  53.  * <caption>Mersenne Twister licence</caption>
  54.  * <tr><td style="padding: 10px">Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
  55.  *     All rights reserved.</td></tr>
  56.  *
  57.  * <tr><td style="padding: 10px">Redistribution and use in source and binary forms, with or without
  58.  * modification, are permitted provided that the following conditions
  59.  * are met:
  60.  * <ol>
  61.  *   <li>Redistributions of source code must retain the above copyright
  62.  *       notice, this list of conditions and the following disclaimer.</li>
  63.  *   <li>Redistributions in binary form must reproduce the above copyright
  64.  *       notice, this list of conditions and the following disclaimer in the
  65.  *       documentation and/or other materials provided with the distribution.</li>
  66.  *   <li>The names of its contributors may not be used to endorse or promote
  67.  *       products derived from this software without specific prior written
  68.  *       permission.</li>
  69.  * </ol></td></tr>
  70.  *
  71.  * <tr><td style="padding: 10px"><strong>THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
  72.  * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
  73.  * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
  74.  * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
  75.  * DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS
  76.  * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
  77.  * OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
  78.  * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
  79.  * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
  80.  * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  81.  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
  82.  * USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
  83.  * DAMAGE.</strong></td></tr>
  84.  * </table>
  85.  *
  86.  * @see <a href="https://en.wikipedia.org/wiki/Mersenne_Twister">Mersenne Twister (Wikipedia)</a>
  87.  * @since 1.0
  88.  */
  89. public class MersenneTwister extends IntProvider {
  90.     /** Mask 32 most significant bits. */
  91.     private static final long INT_MASK_LONG = 0xffffffffL;
  92.     /** Most significant w-r bits. */
  93.     private static final long UPPER_MASK_LONG = 0x80000000L;
  94.     /** Least significant r bits. */
  95.     private static final long LOWER_MASK_LONG = 0x7fffffffL;
  96.     /** Most significant w-r bits. */
  97.     private static final int UPPER_MASK = 0x80000000;
  98.     /** Least significant r bits. */
  99.     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;

  111.     /**
  112.      * Creates a new random number generator.
  113.      *
  114.      * @param seed Initial seed.
  115.      */
  116.     public MersenneTwister(int[] seed) {
  117.         setSeedInternal(seed);
  118.     }

  120.     /** {@inheritDoc} */
  121.     @Override
  122.     protected byte[] getStateInternal() {
  123.         final int[] s = Arrays.copyOf(mt, N + 1);
  124.         s[N] = mti;

  126.         return composeStateInternal(NumberFactory.makeByteArray(s),
  127.                                     super.getStateInternal());
  128.     }

  130.     /** {@inheritDoc} */
  131.     @Override
  132.     protected void setStateInternal(byte[] s) {
  133.         final byte[][] c = splitStateInternal(s, (N + 1) * 4);

  135.         final int[] tmp = NumberFactory.makeIntArray(c[0]);
  136.         System.arraycopy(tmp, 0, mt, 0, N);
  137.         mti = tmp[N];

  139.         super.setStateInternal(c[1]);
  140.     }

  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);

  150.         // Initial index.
  151.         mti = N;
  152.     }

  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;

  168.         initializeState(state);

  170.         final int nextIndex = mixSeedAndState(state, seed);

  172.         mixState(state, nextIndex);

  174.         state[0] = (int) UPPER_MASK_LONG; // MSB is 1, ensuring non-zero initial array.
  175.     }

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

  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;

  208.         int i = 1;
  209.         int j = 0;

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

  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;

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

  254.     /** {@inheritDoc} */
  255.     @Override
  256.     public int next() {
  257.         int y;

  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];

  276.             mti = 0;
  277.         }

  279.         y = mt[mti++];

  281.         // Tempering.
  282.         y ^=  y >>> 11;
  283.         y ^= (y << 7) & 0x9d2c5680;
  284.         y ^= (y << 15) & 0xefc60000;
  285.         y ^=  y >>> 18;

  287.         return y;
  288.     }
  289. }
Created with JaCoCo 0.8.12.202403310830Code Coverage Report for Apache Commons RNG 1.6