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.  *      https://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.  */

  18. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */
  22. package org.hipparchus.random;

  24. import java.io.Serializable;

  26. import org.hipparchus.util.FastMath;


  29. /**
  30.  * This class implements a powerful pseudo-random number generator
  31.  * developed by Makoto Matsumoto and Takuji Nishimura during
  32.  * 1996-1997.
  33.  * <p>
  34.  * <b>Caveat:</b> It is recommended to use one of WELL generators rather
  35.  * than the MersenneTwister generator (see
  36.  * <a href="http://www.iro.umontreal.ca/~panneton/WELLRNG.html">
  37.  * this paper</a> for more information).
  38.  * </p>
  39.  * <p>
  40.  * This generator features an extremely long period
  41.  * (2<sup>19937</sup>-1) and 623-dimensional equidistribution up to 32
  42.  * bits accuracy. The home page for this generator is located at <a
  43.  * href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html">
  44.  * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html</a>.
  45.  * </p>
  46.  * <p>
  47.  * This generator is described in a paper by Makoto Matsumoto and
  48.  * Takuji Nishimura in 1998: <a
  49.  * href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/mt.pdf">Mersenne
  50.  * Twister: A 623-Dimensionally Equidistributed Uniform Pseudo-Random
  51.  * Number Generator</a>, ACM Transactions on Modeling and Computer
  52.  * Simulation, Vol. 8, No. 1, January 1998, pp 3--30.
  53.  * </p>
  54.  * <p>
  55.  * This class is mainly a Java port of the 2002-01-26 version of
  56.  * the generator written in C by Makoto Matsumoto and Takuji
  57.  * Nishimura. Here is their original copyright:
  58.  * </p>
  59.  * <blockquote>
  60.  * <p>Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
  61.  *     All rights reserved.</p>
  62.  * <p>Redistribution and use in source and binary forms, with or without
  63.  * modification, are permitted provided that the following conditions
  64.  * are met:</p>
  65.  * <ol>
  66.  *   <li>Redistributions of source code must retain the above copyright
  67.  *       notice, this list of conditions and the following disclaimer.</li>
  68.  *   <li>Redistributions in binary form must reproduce the above copyright
  69.  *       notice, this list of conditions and the following disclaimer in the
  70.  *       documentation and/or other materials provided with the distribution.</li>
  71.  *   <li>The names of its contributors may not be used to endorse or promote
  72.  *       products derived from this software without specific prior written
  73.  *       permission.</li>
  74.  * </ol>
  75.  *
  76.  * <p><strong>THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
  77.  * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
  78.  * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
  79.  * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
  80.  * DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS
  81.  * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
  82.  * OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
  83.  * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
  84.  * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
  85.  * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  86.  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
  87.  * USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
  88.  * DAMAGE.</strong></p>
  89.  * </blockquote>
  90.  */
  91. public class MersenneTwister extends IntRandomGenerator implements Serializable {

  93.     /** Serializable version identifier. */
  94.     private static final long serialVersionUID = 20160529L;

  96.     /** Size of the bytes pool. */
  97.     private static final int   N     = 624;

  99.     /** Period second parameter. */
  100.     private static final int   M     = 397;

  102.     /** X * MATRIX_A for X = {0, 1}. */
  103.     private static final int[] MAG01 = { 0x0, 0x9908b0df };

  105.     /** Bytes pool. */
  106.     private int[] mt;

  108.     /** Current index in the bytes pool. */
  109.     private int   mti;

  111.     /**
  112.      * Creates a new random number generator.
  113.      * <p>
  114.      * The instance is initialized using the current time plus the
  115.      * system identity hash code of this instance as the seed.
  116.      */
  117.     public MersenneTwister() {
  118.         mt = new int[N];
  119.         setSeed(System.currentTimeMillis() + System.identityHashCode(this));
  120.     }

  122.     /**
  123.      * Creates a new random number generator using a single int seed.
  124.      * @param seed the initial seed (32 bits integer)
  125.      */
  126.     public MersenneTwister(int seed) {
  127.         mt = new int[N];
  128.         setSeed(seed);
  129.     }

  131.     /**
  132.      * Creates a new random number generator using an int array seed.
  133.      * @param seed the initial seed (32 bits integers array), if null
  134.      * the seed of the generator will be related to the current time
  135.      */
  136.     public MersenneTwister(int[] seed) {
  137.         mt = new int[N];
  138.         setSeed(seed);
  139.     }

  141.     /**
  142.      * Creates a new random number generator using a single long seed.
  143.      * @param seed the initial seed (64 bits integer)
  144.      */
  145.     public MersenneTwister(long seed) {
  146.         mt = new int[N];
  147.         setSeed(seed);
  148.     }

  150.     /**
  151.      * Reinitialize the generator as if just built with the given int seed.
  152.      * <p>
  153.      * The state of the generator is exactly the same as a new
  154.      * generator built with the same seed.
  155.      *
  156.      * @param seed the initial seed (32 bits integer)
  157.      */
  158.     @Override
  159.     public void setSeed(int seed) {
  160.         // we use a long masked by 0xffffffffL as a poor man unsigned int
  161.         long longMT = seed & 0xffffffffL;
  162.         // NB: unlike original C code, we are working with java longs,
  163.         // the cast below makes masking unnecessary
  164.         mt[0]= (int) longMT;
  165.         for (mti = 1; mti < N; ++mti) {
  166.             // See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier.
  167.             // initializer from the 2002-01-09 C version by Makoto Matsumoto
  168.             longMT = (1812433253l * (longMT ^ (longMT >> 30)) + mti) & 0xffffffffL;
  169.             mt[mti]= (int) longMT;
  170.         }

  172.         clearCache(); // Clear normal deviate cache
  173.     }

  175.     /**
  176.      * Reinitialize the generator as if just built with the given int array seed.
  177.      * <p>
  178.      * The state of the generator is exactly the same as a new
  179.      * generator built with the same seed.
  180.      *
  181.      * @param seed the initial seed (32 bits integers array), if null
  182.      * the seed of the generator will be the current system time plus the
  183.      * system identity hash code of this instance
  184.      */
  185.     @Override
  186.     public void setSeed(int[] seed) {

  188.         if (seed == null) {
  189.             setSeed(System.currentTimeMillis() + System.identityHashCode(this));
  190.             return;
  191.         }

  193.         setSeed(19650218);
  194.         int i = 1;
  195.         int j = 0;

  197.         for (int k = FastMath.max(N, seed.length); k != 0; k--) {
  198.             long l0 = (mt[i] & 0x7fffffffl)   | ((mt[i]   < 0) ? 0x80000000l : 0x0l);
  199.             long l1 = (mt[i-1] & 0x7fffffffl) | ((mt[i-1] < 0) ? 0x80000000l : 0x0l);
  200.             long l  = (l0 ^ ((l1 ^ (l1 >> 30)) * 1664525l)) + seed[j] + j; // non linear
  201.             mt[i]   = (int) (l & 0xffffffffl);
  202.             i++; j++;
  203.             if (i >= N) {
  204.                 mt[0] = mt[N - 1];
  205.                 i = 1;
  206.             }
  207.             if (j >= seed.length) {
  208.                 j = 0;
  209.             }
  210.         }

  212.         for (int k = N - 1; k != 0; k--) {
  213.             long l0 = (mt[i] & 0x7fffffffl)   | ((mt[i]   < 0) ? 0x80000000l : 0x0l);
  214.             long l1 = (mt[i-1] & 0x7fffffffl) | ((mt[i-1] < 0) ? 0x80000000l : 0x0l);
  215.             long l  = (l0 ^ ((l1 ^ (l1 >> 30)) * 1566083941l)) - i; // non linear
  216.             mt[i]   = (int) (l & 0xffffffffL);
  217.             i++;
  218.             if (i >= N) {
  219.                 mt[0] = mt[N - 1];
  220.                 i = 1;
  221.             }
  222.         }

  224.         mt[0] = 0x80000000; // MSB is 1; assuring non-zero initial array

  226.         clearCache(); // Clear normal deviate cache
  227.     }

  229.     /** {@inheritDoc} */
  230.     @Override
  231.     public int nextInt() {

  233.         int y;

  235.         if (mti >= N) { // generate N words at one time
  236.             int mtNext = mt[0];
  237.             for (int k = 0; k < N - M; ++k) {
  238.                 int mtCurr = mtNext;
  239.                 mtNext = mt[k + 1];
  240.                 y = (mtCurr & 0x80000000) | (mtNext & 0x7fffffff);
  241.                 mt[k] = mt[k + M] ^ (y >>> 1) ^ MAG01[y & 0x1];
  242.             }
  243.             for (int k = N - M; k < N - 1; ++k) {
  244.                 int mtCurr = mtNext;
  245.                 mtNext = mt[k + 1];
  246.                 y = (mtCurr & 0x80000000) | (mtNext & 0x7fffffff);
  247.                 mt[k] = mt[k + (M - N)] ^ (y >>> 1) ^ MAG01[y & 0x1];
  248.             }
  249.             y = (mtNext & 0x80000000) | (mt[0] & 0x7fffffff);
  250.             mt[N - 1] = mt[M - 1] ^ (y >>> 1) ^ MAG01[y & 0x1];

  252.             mti = 0;
  253.         }

  255.         y = mt[mti++];

  257.         // tempering
  258.         y ^=  y >>> 11;
  259.         y ^= (y <<   7) & 0x9d2c5680;
  260.         y ^= (y <<  15) & 0xefc60000;
  261.         y ^=  y >>> 18;

  263.         return y;
  264.     }

  266. }
Created with JaCoCo 0.8.12.202403310830