001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.rng.core.source64;
018
019import java.util.Arrays;
020import org.apache.commons.rng.core.util.NumberFactory;
021
022/**
023 * This class provides the 64-bits version of the originally 32-bits
024 * {@link org.apache.commons.rng.core.source32.MersenneTwister
025 * Mersenne Twister}.
026 *
027 * <p>
028 * This class is mainly a Java port of
029 * <a href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt64.html">
030 *  the 2014/2/23 version of the generator
031 * </a> written in C by Takuji Nishimura and Makoto Matsumoto.
032 * </p>
033 *
034 * <p>
035 * Here is their original copyright:
036 * </p>
037 *
038 * <table style="background-color: #E0E0E0; width: 80%">
039 * <caption>Mersenne Twister licence</caption>
040 * <tr><td style="padding: 10px">Copyright (C) 2004, Makoto Matsumoto and Takuji Nishimura,
041 *     All rights reserved.</td></tr>
042 *
043 * <tr><td style="padding: 10px">Redistribution and use in source and binary forms, with or without
044 * modification, are permitted provided that the following conditions
045 * are met:
046 * <ol>
047 *   <li>Redistributions of source code must retain the above copyright
048 *       notice, this list of conditions and the following disclaimer.</li>
049 *   <li>Redistributions in binary form must reproduce the above copyright
050 *       notice, this list of conditions and the following disclaimer in the
051 *       documentation and/or other materials provided with the distribution.</li>
052 *   <li>The names of its contributors may not be used to endorse or promote
053 *       products derived from this software without specific prior written
054 *       permission.</li>
055 * </ol></td></tr>
056 *
057 * <tr><td style="padding: 10px"><strong>THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
058 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
059 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
060 * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
061 * DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS
062 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
063 * OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
064 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
065 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
066 * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
067 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
068 * USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
069 * DAMAGE.</strong></td></tr>
070 * </table>
071 *
072 * @see <a href="https://en.wikipedia.org/wiki/Mersenne_Twister">Mersenne Twister (Wikipedia)</a>
073 * @since 1.0
074 */
075public class MersenneTwister64 extends LongProvider {
076    /** Size of the bytes pool. */
077    private static final int NN = 312;
078    /** Period second parameter. */
079    private static final int MM = 156;
080    /** X * MATRIX_A for X = {0, 1}. */
081    private static final long[] MAG01 = {0x0L, 0xb5026f5aa96619e9L};
082    /** Most significant 33 bits. */
083    private static final long UM = 0xffffffff80000000L;
084    /** Least significant 31 bits. */
085    private static final long LM = 0x7fffffffL;
086    /** Bytes pool. */
087    private final long[] mt = new long[NN];
088    /** Current index in the bytes pool. */
089    private int mti;
090
091    /**
092     * Creates a new random number generator.
093     *
094     * @param seed Initial seed.
095     */
096    public MersenneTwister64(long[] seed) {
097        setSeedInternal(seed);
098    }
099
100    /** {@inheritDoc} */
101    @Override
102    protected byte[] getStateInternal() {
103        final long[] s = Arrays.copyOf(mt, NN + 1);
104        s[NN] = mti;
105
106        return composeStateInternal(NumberFactory.makeByteArray(s),
107                                    super.getStateInternal());
108    }
109
110    /** {@inheritDoc} */
111    @Override
112    protected void setStateInternal(byte[] s) {
113        final byte[][] c = splitStateInternal(s, (NN + 1) * 8);
114
115        final long[] tmp = NumberFactory.makeLongArray(c[0]);
116        System.arraycopy(tmp, 0, mt, 0, NN);
117        mti = (int) tmp[NN];
118
119        super.setStateInternal(c[1]);
120    }
121
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;
130
131        initState(19650218L);
132        int i = 1;
133        int j = 0;
134
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        }
157
158        mt[0] = 0x8000000000000000L; // MSB is 1; assuring non-zero initial array
159    }
160
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    }
173
174    /** {@inheritDoc} */
175    @Override
176    public long next() {
177        long x;
178
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            }
188
189            x = (mt[NN - 1] & UM) | (mt[0] & LM);
190            mt[NN - 1] = mt[MM - 1] ^ (x >>> 1) ^ MAG01[(int)(x & 0x1L)];
191
192            mti = 0;
193        }
194
195        x = mt[mti++];
196
197        x ^= (x >>> 29) & 0x5555555555555555L;
198        x ^= (x << 17) & 0x71d67fffeda60000L;
199        x ^= (x << 37) & 0xfff7eee000000000L;
200        x ^= x >>> 43;
201
202        return x;
203    }
204}