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 NumberFactory.makeByteArray(s);
106     }
107 
108     /** {@inheritDoc} */
109     @Override
110     protected void setStateInternal(byte[] s) {
111         checkStateSize(s, (NN + 1) * 8);
112 
113         final long[] tmp = NumberFactory.makeLongArray(s);
114         System.arraycopy(tmp, 0, mt, 0, NN);
115         mti = (int) tmp[NN];
116     }
117 
118     /**
119      * Initializes the generator with the given seed.
120      *
121      * @param seed Initial seed.
122      */
123     private void setSeedInternal(long[] seed) {
124         if (seed.length == 0) {
125             // Accept empty seed.
126             seed = new long[1];
127         }
128 
129         initState(19650218L);
130         int i = 1;
131         int j = 0;
132 
133         for (int k = Math.max(NN, seed.length); k != 0; k--) {
134             final long mm1 = mt[i - 1];
135             mt[i] = (mt[i] ^ ((mm1 ^ (mm1 >>> 62)) * 0x369dea0f31a53f85L)) + seed[j] + j; // non linear
136             i++;
137             j++;
138             if (i >= NN) {
139                 mt[0] = mt[NN - 1];
140                 i = 1;
141             }
142             if (j >= seed.length) {
143                 j = 0;
144             }
145         }
146         for (int k = NN - 1; k != 0; k--) {
147             final long mm1 = mt[i - 1];
148             mt[i] = (mt[i] ^ ((mm1 ^ (mm1 >>> 62)) * 0x27bb2ee687b0b0fdL)) - i; // non linear
149             i++;
150             if (i >= NN) {
151                 mt[0] = mt[NN - 1];
152                 i = 1;
153             }
154         }
155 
156         mt[0] = 0x8000000000000000L; // MSB is 1; assuring non-zero initial array
157     }
158 
159     /**
160      * Initialize the internal state of this instance.
161      *
162      * @param seed Seed.
163      */
164     private void initState(long seed) {
165         mt[0] = seed;
166         for (mti = 1; mti < NN; mti++) {
167             final long mm1 = mt[mti - 1];
168             mt[mti] = 0x5851f42d4c957f2dL * (mm1 ^ (mm1 >>> 62)) + mti;
169         }
170     }
171 
172     /** {@inheritDoc} */
173     @Override
174     public long next() {
175         long x;
176 
177         if (mti >= NN) { // generate NN words at one time
178             for (int i = 0; i < NN - MM; i++) {
179                 x = (mt[i] & UM) | (mt[i + 1] & LM);
180                 mt[i] = mt[i + MM] ^ (x >>> 1) ^ MAG01[(int)(x & 0x1L)];
181             }
182             for (int i = NN - MM; i < NN - 1; i++) {
183                 x = (mt[i] & UM) | (mt[i + 1] & LM);
184                 mt[i] = mt[ i + (MM - NN)] ^ (x >>> 1) ^ MAG01[(int)(x & 0x1L)];
185             }
186 
187             x = (mt[NN - 1] & UM) | (mt[0] & LM);
188             mt[NN - 1] = mt[MM - 1] ^ (x >>> 1) ^ MAG01[(int)(x & 0x1L)];
189 
190             mti = 0;
191         }
192 
193         x = mt[mti++];
194 
195         x ^= (x >>> 29) & 0x5555555555555555L;
196         x ^= (x << 17) & 0x71d67fffeda60000L;
197         x ^= (x << 37) & 0xfff7eee000000000L;
198         x ^= x >>> 43;
199 
200         return x;
201     }
202 }