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