View Javadoc
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;
18  
19  import java.util.Arrays;
20  import org.apache.commons.rng.core.util.NumberFactory;
21  
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 border="0" width="80%" cellpadding="10" style="background-color: #E0E0E0" summary="Mersenne Twister licence">
39   * <tr><td>Copyright (C) 2004, Makoto Matsumoto and Takuji Nishimura,
40   *     All rights reserved.</td></tr>
41   *
42   * <tr><td>Redistribution and use in source and binary forms, with or without
43   * modification, are permitted provided that the following conditions
44   * are met:
45   * <ol>
46   *   <li>Redistributions of source code must retain the above copyright
47   *       notice, this list of conditions and the following disclaimer.</li>
48   *   <li>Redistributions in binary form must reproduce the above copyright
49   *       notice, this list of conditions and the following disclaimer in the
50   *       documentation and/or other materials provided with the distribution.</li>
51   *   <li>The names of its contributors may not be used to endorse or promote
52   *       products derived from this software without specific prior written
53   *       permission.</li>
54   * </ol></td></tr>
55   *
56   * <tr><td><strong>THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
57   * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
58   * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
59   * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
60   * DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS
61   * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
62   * OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
63   * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
64   * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
65   * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
66   * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
67   * USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
68   * DAMAGE.</strong></td></tr>
69   * </table>
70   *
71   * @see <a href="https://en.wikipedia.org/wiki/Mersenne_Twister">Mersenne Twister (Wikipedia)</a>
72   * @since 1.0
73   */
74  public class MersenneTwister64 extends LongProvider {
75      /** Size of the bytes pool. */
76      private static final int NN = 312;
77      /** Period second parameter. */
78      private static final int MM = 156;
79      /** X * MATRIX_A for X = {0, 1}. */
80      private static final long[] MAG01 = { 0x0, 0xb5026f5aa96619e9L };
81      /** Most significant 33 bits. */
82      private static final long UM = 0xffffffff80000000L;
83      /** Least significant 31 bits. */
84      private static final long LM = 0x7fffffffL;
85      /** Bytes pool. */
86      private long[] mt = new long[NN];
87      /** Current index in the bytes pool. */
88      private int mti;
89  
90      /**
91       * Creates a new random number generator.
92       *
93       * @param seed Initial seed.
94       */
95      public MersenneTwister64(long[] seed) {
96          setSeedInternal(seed);
97      }
98  
99      /** {@inheritDoc} */
100     @Override
101     protected byte[] getStateInternal() {
102         final long[] s = Arrays.copyOf(mt, NN + 1);
103         s[NN] = mti;
104 
105         return composeStateInternal(NumberFactory.makeByteArray(s),
106                                     super.getStateInternal());
107     }
108 
109     /** {@inheritDoc} */
110     @Override
111     protected void setStateInternal(byte[] s) {
112         final byte[][] c = splitStateInternal(s, (NN + 1) * 8);
113 
114         final long[] tmp = NumberFactory.makeLongArray(c[0]);
115         System.arraycopy(tmp, 0, mt, 0, NN);
116         mti = (int) tmp[NN];
117 
118         super.setStateInternal(c[1]);
119     }
120 
121     /**
122      * Initializes the generator with the given seed.
123      *
124      * @param seed Initial seed.
125      */
126     private void setSeedInternal(long[] seed) {
127         if (seed.length == 0) {
128             // Accept empty seed.
129             seed = new long[1];
130         }
131 
132         initState(19650218L);
133         int i = 1;
134         int j = 0;
135 
136         for (int k = Math.max(NN, seed.length); k != 0; k--) {
137             final long mm1 = mt[i - 1];
138             mt[i] = (mt[i] ^ ((mm1 ^ (mm1 >>> 62)) * 0x369dea0f31a53f85L)) + seed[j] + j; // non linear
139             i++;
140             j++;
141             if (i >= NN) {
142                 mt[0] = mt[NN - 1];
143                 i = 1;
144             }
145             if (j >= seed.length) {
146                 j = 0;
147             }
148         }
149         for (int k = NN - 1; k != 0; k--) {
150             final long mm1 = mt[i - 1];
151             mt[i] = (mt[i] ^ ((mm1 ^ (mm1 >>> 62)) * 0x27bb2ee687b0b0fdL)) - i; // non linear
152             i++;
153             if (i >= NN) {
154                 mt[0] = mt[NN - 1];
155                 i = 1;
156             }
157         }
158 
159         mt[0] = 0x8000000000000000L; // MSB is 1; assuring non-zero initial array
160     }
161 
162     /**
163      * Initialize the internal state of this instance.
164      *
165      * @param seed Seed.
166      */
167     private void initState(long seed) {
168         mt[0] = seed;
169         for (mti = 1; mti < NN; mti++) {
170             final long mm1 = mt[mti - 1];
171             mt[mti] = 0x5851f42d4c957f2dL * (mm1 ^ (mm1 >>> 62)) + mti;
172         }
173     }
174 
175     /** {@inheritDoc} */
176     @Override
177     public long next() {
178         long x;
179 
180         if (mti >= NN) { // generate NN words at one time
181             for (int i = 0; i < NN - MM; i++) {
182                 x = (mt[i] & UM) | (mt[i + 1] & LM);
183                 mt[i] = mt[i + MM] ^ (x >>> 1) ^ MAG01[(int)(x & 0x1L)];
184             }
185             for (int i = NN - MM; i < NN - 1; i++) {
186                 x = (mt[i] & UM) | (mt[i + 1] & LM);
187                 mt[i] = mt[ i + (MM - NN)] ^ (x >>> 1) ^ MAG01[(int)(x & 0x1L)];
188             }
189 
190             x = (mt[NN - 1] & UM) | (mt[0] & LM);
191             mt[NN - 1] = mt[MM - 1] ^ (x >>> 1) ^ MAG01[(int)(x & 0x1L)];
192 
193             mti = 0;
194         }
195 
196         x = mt[mti++];
197 
198         x ^= (x >>> 29) & 0x5555555555555555L;
199         x ^= (x << 17) & 0x71d67fffeda60000L;
200         x ^= (x << 37) & 0xfff7eee000000000L;
201         x ^= x >>> 43;
202 
203         return x;
204     }
205 }