    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     */
    017    package org.apache.commons.math3.random;
    018
    019    import java.io.Serializable;
    020
    021    import org.apache.commons.math3.util.FastMath;
    022
    023
    024    /** This class implements a powerful pseudo-random number generator
    025     * developed by Makoto Matsumoto and Takuji Nishimura during
    026     * 1996-1997.
    027
    028     * <p>This generator features an extremely long period
    029     * (2<sup>19937</sup>-1) and 623-dimensional equidistribution up to 32
    030     * bits accuracy. The home page for this generator is located at <a
    031     * href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html">
    032     * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html</a>.</p>
    033
    034     * <p>This generator is described in a paper by Makoto Matsumoto and
    035     * Takuji Nishimura in 1998: <a
    036     * href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/mt.pdf">Mersenne
    037     * Twister: A 623-Dimensionally Equidistributed Uniform Pseudo-Random
    038     * Number Generator</a>, ACM Transactions on Modeling and Computer
    039     * Simulation, Vol. 8, No. 1, January 1998, pp 3--30</p>
    040
    041     * <p>This class is mainly a Java port of the 2002-01-26 version of
    042     * the generator written in C by Makoto Matsumoto and Takuji
    043     * Nishimura. Here is their original copyright:</p>
    044
    045     * <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0">
    046     * <tr><td>Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
    047     *     All rights reserved.</td></tr>
    048
    049     * <tr><td>Redistribution and use in source and binary forms, with or without
    050     * modification, are permitted provided that the following conditions
    051     * are met:
    052     * <ol>
    053     *   <li>Redistributions of source code must retain the above copyright
    054     *       notice, this list of conditions and the following disclaimer.</li>
    055     *   <li>Redistributions in binary form must reproduce the above copyright
    056     *       notice, this list of conditions and the following disclaimer in the
    057     *       documentation and/or other materials provided with the distribution.</li>
    058     *   <li>The names of its contributors may not be used to endorse or promote
    059     *       products derived from this software without specific prior written
    060     *       permission.</li>
    061     * </ol></td></tr>
    062
    063     * <tr><td><strong>THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
    064     * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
    065     * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
    066     * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
    067     * DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS
    068     * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
    069     * OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
    070     * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
    071     * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
    072     * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
    073     * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
    074     * USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
    075     * DAMAGE.</strong></td></tr>
    076     * </table>
    077
    078     * @version $Id: MersenneTwister.java 1416643 2012-12-03 19:37:14Z tn $
    079     * @since 2.0
    080
    081     */
    082    public class MersenneTwister extends BitsStreamGenerator implements Serializable {
    083
    084        /** Serializable version identifier. */
    085        private static final long serialVersionUID = 8661194735290153518L;
    086
    087        /** Size of the bytes pool. */
    088        private static final int   N     = 624;
    089
    090        /** Period second parameter. */
    091        private static final int   M     = 397;
    092
    093        /** X * MATRIX_A for X = {0, 1}. */
    094        private static final int[] MAG01 = { 0x0, 0x9908b0df };
    095
    096        /** Bytes pool. */
    097        private int[] mt;
    098
    099        /** Current index in the bytes pool. */
    100        private int   mti;
    101
    102        /** Creates a new random number generator.
    103         * <p>The instance is initialized using the current time plus the
    104         * system identity hash code of this instance as the seed.</p>
    105         */
    106        public MersenneTwister() {
    107            mt = new int[N];
    108            setSeed(System.currentTimeMillis() + System.identityHashCode(this));
    109        }
    110
    111        /** Creates a new random number generator using a single int seed.
    112         * @param seed the initial seed (32 bits integer)
    113         */
    114        public MersenneTwister(int seed) {
    115            mt = new int[N];
    116            setSeed(seed);
    117        }
    118
    119        /** Creates a new random number generator using an int array seed.
    120         * @param seed the initial seed (32 bits integers array), if null
    121         * the seed of the generator will be related to the current time
    122         */
    123        public MersenneTwister(int[] seed) {
    124            mt = new int[N];
    125            setSeed(seed);
    126        }
    127
    128        /** Creates a new random number generator using a single long seed.
    129         * @param seed the initial seed (64 bits integer)
    130         */
    131        public MersenneTwister(long seed) {
    132            mt = new int[N];
    133            setSeed(seed);
    134        }
    135
    136        /** Reinitialize the generator as if just built with the given int seed.
    137         * <p>The state of the generator is exactly the same as a new
    138         * generator built with the same seed.</p>
    139         * @param seed the initial seed (32 bits integer)
    140         */
    141        @Override
    142        public void setSeed(int seed) {
    143            // we use a long masked by 0xffffffffL as a poor man unsigned int
    144            long longMT = seed;
    145            // NB: unlike original C code, we are working with java longs, the cast below makes masking unnecessary
    146            mt[0]= (int) longMT;
    147            for (mti = 1; mti < N; ++mti) {
    148                // See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier.
    149                // initializer from the 2002-01-09 C version by Makoto Matsumoto
    150                longMT = (1812433253l * (longMT ^ (longMT >> 30)) + mti) & 0xffffffffL;
    151                mt[mti]= (int) longMT;
    152            }
    153
    154            clear(); // Clear normal deviate cache
    155        }
    156
    157        /** Reinitialize the generator as if just built with the given int array seed.
    158         * <p>The state of the generator is exactly the same as a new
    159         * generator built with the same seed.</p>
    160         * @param seed the initial seed (32 bits integers array), if null
    161         * the seed of the generator will be the current system time plus the
    162         * system identity hash code of this instance
    163         */
    164        @Override
    165        public void setSeed(int[] seed) {
    166
    167            if (seed == null) {
    168                setSeed(System.currentTimeMillis() + System.identityHashCode(this));
    169                return;
    170            }
    171
    172            setSeed(19650218);
    173            int i = 1;
    174            int j = 0;
    175
    176            for (int k = FastMath.max(N, seed.length); k != 0; k--) {
    177                long l0 = (mt[i] & 0x7fffffffl)   | ((mt[i]   < 0) ? 0x80000000l : 0x0l);
    178                long l1 = (mt[i-1] & 0x7fffffffl) | ((mt[i-1] < 0) ? 0x80000000l : 0x0l);
    179                long l  = (l0 ^ ((l1 ^ (l1 >> 30)) * 1664525l)) + seed[j] + j; // non linear
    180                mt[i]   = (int) (l & 0xffffffffl);
    181                i++; j++;
    182                if (i >= N) {
    183                    mt[0] = mt[N - 1];
    184                    i = 1;
    185                }
    186                if (j >= seed.length) {
    187                    j = 0;
    188                }
    189            }
    190
    191            for (int k = N - 1; k != 0; k--) {
    192                long l0 = (mt[i] & 0x7fffffffl)   | ((mt[i]   < 0) ? 0x80000000l : 0x0l);
    193                long l1 = (mt[i-1] & 0x7fffffffl) | ((mt[i-1] < 0) ? 0x80000000l : 0x0l);
    194                long l  = (l0 ^ ((l1 ^ (l1 >> 30)) * 1566083941l)) - i; // non linear
    195                mt[i]   = (int) (l & 0xffffffffL);
    196                i++;
    197                if (i >= N) {
    198                    mt[0] = mt[N - 1];
    199                    i = 1;
    200                }
    201            }
    202
    203            mt[0] = 0x80000000; // MSB is 1; assuring non-zero initial array
    204
    205            clear(); // Clear normal deviate cache
    206
    207        }
    208
    209        /** Reinitialize the generator as if just built with the given long seed.
    210         * <p>The state of the generator is exactly the same as a new
    211         * generator built with the same seed.</p>
    212         * @param seed the initial seed (64 bits integer)
    213         */
    214        @Override
    215        public void setSeed(long seed) {
    216            setSeed(new int[] { (int) (seed >>> 32), (int) (seed & 0xffffffffl) });
    217        }
    218
    219        /** Generate next pseudorandom number.
    220         * <p>This method is the core generation algorithm. It is used by all the
    221         * public generation methods for the various primitive types {@link
    222         * #nextBoolean()}, {@link #nextBytes(byte[])}, {@link #nextDouble()},
    223         * {@link #nextFloat()}, {@link #nextGaussian()}, {@link #nextInt()},
    224         * {@link #next(int)} and {@link #nextLong()}.</p>
    225         * @param bits number of random bits to produce
    226         * @return random bits generated
    227         */
    228        @Override
    229        protected int next(int bits) {
    230
    231            int y;
    232
    233            if (mti >= N) { // generate N words at one time
    234                int mtNext = mt[0];
    235                for (int k = 0; k < N - M; ++k) {
    236                    int mtCurr = mtNext;
    237                    mtNext = mt[k + 1];
    238                    y = (mtCurr & 0x80000000) | (mtNext & 0x7fffffff);
    239                    mt[k] = mt[k + M] ^ (y >>> 1) ^ MAG01[y & 0x1];
    240                }
    241                for (int k = N - M; k < N - 1; ++k) {
    242                    int mtCurr = mtNext;
    243                    mtNext = mt[k + 1];
    244                    y = (mtCurr & 0x80000000) | (mtNext & 0x7fffffff);
    245                    mt[k] = mt[k + (M - N)] ^ (y >>> 1) ^ MAG01[y & 0x1];
    246                }
    247                y = (mtNext & 0x80000000) | (mt[0] & 0x7fffffff);
    248                mt[N - 1] = mt[M - 1] ^ (y >>> 1) ^ MAG01[y & 0x1];
    249
    250                mti = 0;
    251            }
    252
    253            y = mt[mti++];
    254
    255            // tempering
    256            y ^=  y >>> 11;
    257            y ^= (y <<   7) & 0x9d2c5680;
    258            y ^= (y <<  15) & 0xefc60000;
    259            y ^=  y >>> 18;
    260
    261            return y >>> (32 - bits);
    262
    263        }
    264
    265    }