Source code

001/*
002 * Copyright 2005, Nick Galbreath -- nickg [at] modp [dot] com
003 * All rights reserved.
004 *
005 * Redistribution and use in source and binary forms, with or without
006 * modification, are permitted provided that the following conditions are
007 * met:
008 *
009 *   Redistributions of source code must retain the above copyright
010 *   notice, this list of conditions and the following disclaimer.
011 *
012 *   Redistributions in binary form must reproduce the above copyright
013 *   notice, this list of conditions and the following disclaimer in the
014 *   documentation and/or other materials provided with the distribution.
015 *
016 *   Neither the name of the modp.com nor the names of its
017 *   contributors may be used to endorse or promote products derived from
018 *   this software without specific prior written permission.
019 *
020 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
021 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
022 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
023 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
024 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
025 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
026 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
027 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
028 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
029 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
030 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
031 *
032 * This is the standard "new" BSD license:
033 * http://www.opensource.org/licenses/bsd-license.php
034 *
035 * Portions may also be
036 * Copyright (C) 2004, Makoto Matsumoto and Takuji Nishimura,
037 * All rights reserved.
038 * (and covered under the BSD license)
039 * See http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/VERSIONS/C-LANG/mt19937-64.c
040 */
041package squidpony.squidmath;
042
043import squidpony.StringKit;
044
045import java.io.Serializable;
046import java.util.Arrays;
047
048/**
049 * Mersenne Twister, 64-bit version as a RandomnessSource.
050 * <br>
051 * Similar to the regular 32-bit Mersenne Twister but implemented with 64-bit
052 * values (Java <code>long</code>), and with different output. This generator is probably
053 * not the best to use because of known statistical problems and low speed, but its period
054 * is absurdly high, {@code pow(2, 19937) - 1}. {@link LongPeriodRNG} has significantly
055 * better speed and statistical quality, and also has a large period, {@code pow(2, 1024) - 1}.
056 * {@link IsaacRNG} is slower, but offers impeccable quality, and from its webpage, "Cycles
057 * are guaranteed to be at least {@code pow(2, 40)} values long, and they are
058 * {@code pow(2, 8295)} values long on average." IsaacRNG should be your choice if security is a
059 * concern, LongPeriodRNG if quality and speed are important, and MersenneTwister should be used
060 * if period is the only criterion to judge an RNG on. Keep in mind that extremely long periods
061 * are not always a good thing; there are known states for the Mersenne Twister that produce dozens
062 * of {@code 0} outputs in a row, which is fundamentally impossible for {@link LightRNG} or
063 * {@link DiverRNG}. It also would take longer to exhaust the period of a 128-bit-state generator
064 * (generating 100 gigabytes per second) than the amount of time humans have walked the Earth.
065 * This has 19968 bits of state... so if 128 is more than is possible to exhaust, why do you need
066 * 19968 bits, again? Consider {@link GoatRNG} if you want a high-quality 128-bit generator.
067 * <br>
068 * This is mostly a straight port of the
069 * <a href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/VERSIONS/C-LANG/mt19937-64.c">
070 * C-code (mt19937-64.c)</a>, but made more "java-like". This version was originally found
071 * at <a href="https://github.com/javabeanz/javarng/blob/master/com/modp/random/MersenneTwister64.java">
072 * an archived version of a Google Code repo</a>, and it is licensed under the 3-clause BSD license, like
073 * the Mersenne Twister.
074 * <br>
075 * References:
076 * <br>
077 * <ul>
078 * <li>
079 * <a href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html">
080 * The Mersenne Twister Home Page </a>
081 * </li>
082 * <li>  Makato Matsumoto and Takuji Nishimura,
083 * <a href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/mt.pdf">"Mersenne Twister:A 623-Dimensionally Equidistributed Uniform Pseudo-Random Number Generator"</a>,
084 * <i>ACM Transactions on Modeling and Computer Simulation, </i> Vol. 8, No. 1,
085 * January 1998, pp 3--30.</li>
086 * </ul>
087 *
088 * @author Nick Galbreath -- nickg [at] modp [dot] com
089 * @author Tommy Ettinger
090 * @version 1.1 -- 07-Oct-2017
091 */
092public class MersenneTwister implements Serializable, RandomnessSource {
093
094    private static final long serialVersionUID = 1001000L;
095
096    private static final int NN = 312;
097
098    private static final int MM = 156;
099
100    private static final long MATRIX_A = 0xB5026F5AA96619E9L;
101
102    /**
103     * Mask: Most significant 33 bits
104     */
105    private static final long UM = 0xFFFFFFFF80000000L;
106
107    /**
108     * Mask: Least significant 31 bits
109     */
110    private static final long LM = 0x7FFFFFFFL;
111
112    /**
113     * Mersenne Twister data.
114     */
115    private final long[] mt = new long[NN];
116
117    /**
118     * Mersenne Twister Index.
119     */
120    private int mti = NN + 1;
121
122    /**
123     * Internal to hold 64 bits, that might
124     * used to generate two 32 bit values.
125     */
126    private long extra;
127
128    /**
129     * Set to true if we need to generate another 64 bits, false if we have enough bits available for an int.
130     */
131    private boolean bitState = true;
132
133    /**
134     * Seeds this using two calls to Math.random().
135     */
136    public MersenneTwister() {
137        setSeed((long) ((Math.random() * 2.0 - 1.0) * 0x8000000000000L)
138                ^ (long) ((Math.random() * 2.0 - 1.0) * 0x8000000000000000L));
139    }
140
141    /**
142     * Seeds this with the given long, which will be used to affect the large state.
143     *
144     * @param seed any long
145     */
146    public MersenneTwister(final long seed) {
147        setSeed(seed);
148    }
149
150    /**
151     * Seeds this with the given long array, which will be used to affect the large state, and not used directly.
152     *
153     * @param seed a long array; generally should be non-null
154     */
155    public MersenneTwister(final long[] seed) {
156        setSeed(seed);
157    }
158
159    /**
160     * Initalize the pseudo random number generator with 64 bits.
161     * Not the same as setState() in StatefulRandomness; this changes the seed quite a bit.
162     *
163     * @param seed any long
164     */
165    public void setSeed(final long seed) {
166        mt[0] = seed;
167        for (mti = 1; mti < NN; mti++) {
168            mt[mti] = (6364136223846793005L * (mt[mti - 1] ^ (mt[mti - 1] >>> 62)) + mti);
169        }
170    }
171
172    /**
173     * Initalize the pseudo random number generator with a long array of any size, which should not be null but can be.
174     * Not the same as setState() in StatefulRandomness; this changes the seed quite a bit.
175     *
176     * @param array any long array
177     */
178    public void setSeed(final long[] array) {
179        setSeed(19650218L);
180        if (array == null)
181            return;
182        int i = 1;
183        int j = 0;
184        int k = (Math.max(NN, array.length));
185        for (; k != 0; k--) {
186            mt[i] = (mt[i] ^ ((mt[i - 1] ^ (mt[i - 1] >>> 62)) * 3935559000370003845L))
187                    + array[j] + j;
188            i++;
189            j++;
190            if (i >= NN) {
191                mt[0] = mt[NN - 1];
192                i = 1;
193            }
194            if (j >= array.length)
195                j = 0;
196        }
197        for (k = NN - 1; k != 0; k--) {
198            mt[i] = (mt[i] ^ ((mt[i - 1] ^ (mt[i - 1] >>> 62)) * 2862933555777941757L))
199                    - i;
200            i++;
201            if (i >= NN) {
202                mt[0] = mt[NN - 1];
203                i = 1;
204            }
205        }
206
207        mt[0] = 0x8000000000000000L; /* MSB is 1; assuring non-zero initial array */
208    }
209
210    /**
211     * Returns up to 32 random bits.
212     * <br>
213     * The implementation splits a 64-bit long into two 32-bit chunks.
214     * @param bits the number of bits to output, between 1 and 32 (both inclusive)
215     * @return an int with the specified number of pseudo-random bits
216     */
217    @Override
218    public int next(final int bits) {
219        //return ((int)nextLong()) >>> (32 - numbits);
220
221        if (bitState) {
222            extra = nextLong();
223            bitState = false;
224            return (int) (extra >>> (64 - bits));
225        } else {
226            bitState = true;
227            return ((int) extra) >>> (32 - bits);
228        }
229    }
230
231    /**
232     * Returns 64 random bits.
233     * @return a pseudo-random long, which can have any 64-bit value, positive or negative
234     */
235    @Override
236    public long nextLong() {
237        int i;
238        long x;
239        if (mti >= NN) { /* generate NN words at one time */
240
241            for (i = 0; i < NN - MM; i++) {
242                x = (mt[i] & UM) | (mt[i + 1] & LM);
243                mt[i] = mt[i + MM] ^ (x >>> 1) ^ (x & 1L) * MATRIX_A;
244            }
245            for (; i < NN - 1; i++) {
246                x = (mt[i] & UM) | (mt[i + 1] & LM);
247                mt[i] = mt[i + (MM - NN)] ^ (x >>> 1) ^ (x & 1L) * MATRIX_A;
248            }
249            x = (mt[NN - 1] & UM) | (mt[0] & LM);
250            mt[NN - 1] = mt[MM - 1] ^ (x >>> 1) ^ (x & 1L) * MATRIX_A;
251
252            mti = 0;
253        }
254
255        x = mt[mti++];
256
257        x ^= (x >>> 29) & 0x5555555555555555L;
258        x ^= (x << 17) & 0x71D67FFFEDA60000L;
259        x ^= (x << 37) & 0xFFF7EEE000000000L;
260        x ^= (x >>> 43);
261
262        return x;
263    }
264
265    @Override
266    public final MersenneTwister copy() {
267        MersenneTwister f = new MersenneTwister(MATRIX_A);
268        f.mti = mti;
269        f.extra = extra;
270        f.bitState = bitState;
271        System.arraycopy(mt, 0, f.mt, 0, mt.length);
272        return f;
273    }
274
275    @Override
276    public String toString() {
277        return "MersenneTwister with state hashed as " + StringKit.hexHash(mt) +
278                " and index " + mti;
279    }
280
281    @Override
282    public boolean equals(Object o) {
283        if (this == o) return true;
284        if (o == null || getClass() != o.getClass()) return false;
285
286        MersenneTwister mt64RNG = (MersenneTwister) o;
287
288        return mti == mt64RNG.mti && extra == mt64RNG.extra && bitState == mt64RNG.bitState && Arrays.equals(mt, mt64RNG.mt);
289    }
290
291    @Override
292    public int hashCode() {
293        int result = CrossHash.hash(mt);
294        result = 31 * result + mti;
295        result = 31 * result + (int) (extra ^ (extra >>> 32));
296        result = 31 * result + (bitState ? 1 : 0);
297        return result;
298    }
299
300}