MersenneTwister64.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.source64;

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

  22. /**
  23.  * This class provides the 64-bits version of the originally 32-bits
  24.  * {@link org.apache.commons.rng.core.source32.MersenneTwister
  25.  * Mersenne Twister}.
  26.  *
  27.  * <p>
  28.  * This class is mainly a Java port of
  29.  * <a href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt64.html">
  30.  *  the 2014/2/23 version of the generator
  31.  * </a> written in C by Takuji Nishimura and Makoto Matsumoto.
  32.  * </p>
  33.  *
  34.  * <p>
  35.  * Here is their original copyright:
  36.  * </p>
  37.  *
  38.  * <table style="background-color: #E0E0E0; width: 80%">
  39.  * <caption>Mersenne Twister licence</caption>
  40.  * <tr><td style="padding: 10px">Copyright (C) 2004, Makoto Matsumoto and Takuji Nishimura,
  41.  *     All rights reserved.</td></tr>
  42.  *
  43.  * <tr><td style="padding: 10px">Redistribution and use in source and binary forms, with or without
  44.  * modification, are permitted provided that the following conditions
  45.  * are met:
  46.  * <ol>
  47.  *   <li>Redistributions of source code must retain the above copyright
  48.  *       notice, this list of conditions and the following disclaimer.</li>
  49.  *   <li>Redistributions in binary form must reproduce the above copyright
  50.  *       notice, this list of conditions and the following disclaimer in the
  51.  *       documentation and/or other materials provided with the distribution.</li>
  52.  *   <li>The names of its contributors may not be used to endorse or promote
  53.  *       products derived from this software without specific prior written
  54.  *       permission.</li>
  55.  * </ol></td></tr>
  56.  *
  57.  * <tr><td style="padding: 10px"><strong>THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
  58.  * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
  59.  * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
  60.  * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
  61.  * DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS
  62.  * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
  63.  * OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
  64.  * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
  65.  * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
  66.  * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  67.  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
  68.  * USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
  69.  * DAMAGE.</strong></td></tr>
  70.  * </table>
  71.  *
  72.  * @see <a href="https://en.wikipedia.org/wiki/Mersenne_Twister">Mersenne Twister (Wikipedia)</a>
  73.  * @since 1.0
  74.  */
  75. public class MersenneTwister64 extends LongProvider {
  76.     /** Size of the bytes pool. */
  77.     private static final int NN = 312;
  78.     /** Period second parameter. */
  79.     private static final int MM = 156;
  80.     /** X * MATRIX_A for X = {0, 1}. */
  81.     private static final long[] MAG01 = {0x0L, 0xb5026f5aa96619e9L};
  82.     /** Most significant 33 bits. */
  83.     private static final long UM = 0xffffffff80000000L;
  84.     /** Least significant 31 bits. */
  85.     private static final long LM = 0x7fffffffL;
  86.     /** Bytes pool. */
  87.     private final long[] mt = new long[NN];
  88.     /** Current index in the bytes pool. */
  89.     private int mti;

  91.     /**
  92.      * Creates a new random number generator.
  93.      *
  94.      * @param seed Initial seed.
  95.      */
  96.     public MersenneTwister64(long[] seed) {
  97.         setSeedInternal(seed);
  98.     }

  100.     /** {@inheritDoc} */
  101.     @Override
  102.     protected byte[] getStateInternal() {
  103.         final long[] s = Arrays.copyOf(mt, NN + 1);
  104.         s[NN] = mti;

  106.         return composeStateInternal(NumberFactory.makeByteArray(s),
  107.                                     super.getStateInternal());
  108.     }

  110.     /** {@inheritDoc} */
  111.     @Override
  112.     protected void setStateInternal(byte[] s) {
  113.         final byte[][] c = splitStateInternal(s, (NN + 1) * 8);

  115.         final long[] tmp = NumberFactory.makeLongArray(c[0]);
  116.         System.arraycopy(tmp, 0, mt, 0, NN);
  117.         mti = (int) tmp[NN];

  119.         super.setStateInternal(c[1]);
  120.     }

  122.     /**
  123.      * Initializes the generator with the given seed.
  124.      *
  125.      * @param inputSeed Initial seed.
  126.      */
  127.     private void setSeedInternal(long[] inputSeed) {
  128.         // Accept empty seed.
  129.         final long[] seed = (inputSeed.length == 0) ? new long[1] : inputSeed;

  131.         initState(19650218L);
  132.         int i = 1;
  133.         int j = 0;

  135.         for (int k = Math.max(NN, seed.length); k != 0; k--) {
  136.             final long mm1 = mt[i - 1];
  137.             mt[i] = (mt[i] ^ ((mm1 ^ (mm1 >>> 62)) * 0x369dea0f31a53f85L)) + seed[j] + j; // non linear
  138.             i++;
  139.             j++;
  140.             if (i >= NN) {
  141.                 mt[0] = mt[NN - 1];
  142.                 i = 1;
  143.             }
  144.             if (j >= seed.length) {
  145.                 j = 0;
  146.             }
  147.         }
  148.         for (int k = NN - 1; k != 0; k--) {
  149.             final long mm1 = mt[i - 1];
  150.             mt[i] = (mt[i] ^ ((mm1 ^ (mm1 >>> 62)) * 0x27bb2ee687b0b0fdL)) - i; // non linear
  151.             i++;
  152.             if (i >= NN) {
  153.                 mt[0] = mt[NN - 1];
  154.                 i = 1;
  155.             }
  156.         }

  158.         mt[0] = 0x8000000000000000L; // MSB is 1; assuring non-zero initial array
  159.     }

  161.     /**
  162.      * Initialize the internal state of this instance.
  163.      *
  164.      * @param seed Seed.
  165.      */
  166.     private void initState(long seed) {
  167.         mt[0] = seed;
  168.         for (mti = 1; mti < NN; mti++) {
  169.             final long mm1 = mt[mti - 1];
  170.             mt[mti] = 0x5851f42d4c957f2dL * (mm1 ^ (mm1 >>> 62)) + mti;
  171.         }
  172.     }

  174.     /** {@inheritDoc} */
  175.     @Override
  176.     public long next() {
  177.         long x;

  179.         if (mti >= NN) { // generate NN words at one time
  180.             for (int i = 0; i < NN - MM; i++) {
  181.                 x = (mt[i] & UM) | (mt[i + 1] & LM);
  182.                 mt[i] = mt[i + MM] ^ (x >>> 1) ^ MAG01[(int)(x & 0x1L)];
  183.             }
  184.             for (int i = NN - MM; i < NN - 1; i++) {
  185.                 x = (mt[i] & UM) | (mt[i + 1] & LM);
  186.                 mt[i] = mt[ i + (MM - NN)] ^ (x >>> 1) ^ MAG01[(int)(x & 0x1L)];
  187.             }

  189.             x = (mt[NN - 1] & UM) | (mt[0] & LM);
  190.             mt[NN - 1] = mt[MM - 1] ^ (x >>> 1) ^ MAG01[(int)(x & 0x1L)];

  192.             mti = 0;
  193.         }

  195.         x = mt[mti++];

  197.         x ^= (x >>> 29) & 0x5555555555555555L;
  198.         x ^= (x << 17) & 0x71d67fffeda60000L;
  199.         x ^= (x << 37) & 0xfff7eee000000000L;
  200.         x ^= x >>> 43;

  202.         return x;
  203.     }
  204. }
Created with JaCoCo 0.8.12.202403310830