Java tutorial
//package org.aitools.util.math; import java.io.IOException; import java.io.ObjectInputStream; import java.io.ObjectOutputStream; import java.io.Serializable; /** * <h3>Mersenne Twister and MersenneTwisterFast</h3> * <p> * <b>Version 3 </b>, based on version MT199937(99/10/29) of the Mersenne * Twister algorithm found at <a * href="http://www.math.keio.ac.jp/matumoto/emt.html"> The Mersenne Twister * Home Page </a>. By Sean Luke, June 2000. * <p> * <b>MersenneTwister </b> is a drop-in subclass replacement for * java.util.Random. It is properly synchronized and can be used in a * multithreaded environment. * <p> * <b>MersenneTwisterFast </b> is not a subclass of java.util.Random. It has the * same public methods as Random does, however, and it is algorithmically * identical to MersenneTwister. MersenneTwisterFast has hard-code inlined all * of its methods directly, and made all of them final (well, the ones of * consequence anyway). Further, these methods are <i>not </i> synchronized, so * the same MersenneTwisterFast instance cannot be shared by multiple threads. * But all this helps MersenneTwisterFast achieve over twice the speed of * MersenneTwister. * <h3>About the Mersenne Twister</h3> * <p> * This is a Java version of the C-program for MT19937: Integer version. The * MT19937 algorithm was created by Makoto Matsumoto and Takuji Nishimura, who * ask: "When you use this, send an email to: matumoto@math.keio.ac.jp with an * appropriate reference to your work". Indicate that this is a translation of * their algorithm into Java. * <p> * <b>Reference. </b> Makato Matsumoto and Takuji Nishimura, "Mersenne Twister: * A 623-Dimensionally Equidistributed Uniform Pseudo-Random Number Generator", * <i>ACM Transactions on Modeling and Computer Simulation, </i> Vol. 8, No. 1, * January 1998, pp 3--30. * <h3>About this Version</h3> * This version of the code implements the MT19937 Mersenne Twister algorithm, * with the 99/10/29 seeding mechanism. The original mechanism did not permit 0 * as a seed, and odd numbers were not good seed choices. The new version * permits any 32-bit signed integer. This algorithm is identical to the MT19937 * integer algorithm; real values conform to Sun's float and double random * number generator standards rather than attempting to implement the half-open * or full-open MT19937-1 and MT199937-2 algorithms. * <p> * This code is based on standard MT19937 C/C++ code by Takuji Nishimura, with * suggestions from Topher Cooper and Marc Rieffel, July 1997. The code was * originally translated into Java by Michael Lecuyer, January 1999, and is * Copyright (c) 1999 by Michael Lecuyer. The included license is as follows: * <blockquote><font size="-1"> The basic algorithmic work of this library * (appearing in nextInt() and setSeed()) is free software; you can redistribute * it and or modify it under the terms of the GNU Library General Public License * as published by the Free Software Foundation; either version 2 of the * License, or (at your option) any later version. * <p> * This library is distributed in the hope that it will be useful, but WITHOUT * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * FOR A PARTICULAR PURPOSE. See the GNU Library General Public License for more * details. You should have received a copy of the GNU Library General Public * License along with this library; if not, write to the Free Foundation, Inc., * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA </font> </blockquote> * <h3>Bug Fixes</h3> * <p> * This implementation implements the bug fixes made in Java 1.2's version of * Random, which means it can be used with earlier versions of Java. See <a * href="http://www.javasoft.com/products/jdk/1.2/docs/api/java/util/Random.html"> * the JDK 1.2 java.util.Random documentation </a> for further documentation on * the random-number generation contracts made. Additionally, there's an * undocumented bug in the JDK java.util.Random.nextBytes() method, which this * code fixes. * <h3>Important Note on Seeds</h3> * <p> * Just like java.util.Random, this generator accepts a long seed but doesn't * use all of it. java.util.Random uses 48 bits. The Mersenne Twister instead * uses 32 bits (int size). So it's best if your seed does not exceed the int * range. * <h3>Timings On Different Java Versions</h3> * <p> * MersenneTwister can be used reliably on JDK version 1.1.5 or above. Earlier * Java versions have serious bugs in java.util.Random; only MersenneTwisterFast * (and not MersenneTwister nor java.util.Random) should be used with them. And * why would you use 'em anyway? They're very slow, as you'll see. Here are some * timings in milliseconds on a Sun Creator3D/Ultra 60 running SunOS 5.6. * <dl> * <dt><i>Standard C Version (gcc -O2) </i> * <dd>1070 * <dt><i>Standard C Version (Solaris cc -O) </i> * <dd>1210 * <dt><i>JDK 1.2.2 w/Hotspot Compiler (-O) </i> * <dd>MTF: 1785, MT: 3699, java.util.Random: 4849 * <dt><i>JDK 1.2.1/1.2.2 (-O) </i> * <dd>MTF: 1827, MT: 3868, java.util.Random: 4194 * <dt><i>JDK 1.1.8 (-O) </i> * <dd>MTF: 40509, MT: 45853, java.util.Random: 24604 <br> * Beats me why it's so slow... * <dt><i>JDK 1.1.5 (-O) </i> * <dd>MTF: 4056, MT: 20478, java.util.Random: 19692 * <dt><i>JDK 1.0.2 (-O) </i> * <dd>MTF: 71640, MT: 66176, java.util.Random: 67269 <br> * <i>Important note: </i> Do not MersenneTwister.java or java.util.Random on a * Java version this early! Random number generation in versions less than 1.1.5 * has serious bugs. * </dl> * * @version 3 */ public class MersenneTwister extends java.util.Random implements Serializable { // Period parameters private static final int N = 624; private static final int M = 397; private static final int MATRIX_A = 0x9908b0df; // private static final * constant vector a private static final int UPPER_MASK = 0x80000000; // most significant w-r bits private static final int LOWER_MASK = 0x7fffffff; // least significant r bits // Tempering parameters private static final int TEMPERING_MASK_B = 0x9d2c5680; private static final int TEMPERING_MASK_C = 0xefc60000; // #define TEMPERING_SHIFT_U(y) (y >>> 11) // #define TEMPERING_SHIFT_S(y) (y << 7) // #define TEMPERING_SHIFT_T(y) (y << 15) // #define TEMPERING_SHIFT_L(y) (y >>> 18) private int mt[]; // the array for the state vector private int mti; // mti==N+1 means mt[N] is not initialized private int mag01[]; // a good initial seed (of int size, though stored in a long) private static final long GOOD_SEED = 4357; /* * implemented here because there's a bug in Random's implementation of the * Gaussian code (divide by zero, and log(0), ugh!), yet its gaussian * variables are private so we can't access them here. :-( */ private double __nextNextGaussian; private boolean __haveNextNextGaussian; /** * Constructor using the default seed. */ public MersenneTwister() { super(GOOD_SEED); /* just in case */ setSeed(GOOD_SEED); } /** * Constructor using a given seed. Though you pass this seed in as a long, * it's best to make sure it's actually an integer. * * @param seed the seed to use */ public MersenneTwister(final long seed) { super(seed); /* just in case */ setSeed(seed); } /** * Initalize the pseudo random number generator. This is the old * seed-setting mechanism for the original Mersenne Twister algorithm. You * must not use 0 as your seed, and don't pass in a long that's bigger than * an int (Mersenne Twister only uses the first 32 bits for its seed). Also * it's suggested that for you avoid even-numbered seeds in this older * seed-generation procedure. * * @param seed the seed to use */ synchronized public void setSeedOld(final long seed) { // it's always good style to call super super.setSeed(seed); // Due to a bug in java.util.Random clear up to 1.2, we're // doing our own Gaussian variable. this.__haveNextNextGaussian = false; this.mt = new int[N]; // setting initial seeds to mt[N] using // the generator Line 25 of Table 1 in // [KNUTH 1981, The Art of Computer Programming // Vol. 2 (2nd Ed.), pp102] // the 0xffffffff is commented out because in Java // ints are always 32 bits; hence i & 0xffffffff == i this.mt[0] = ((int) seed); // & 0xffffffff; for (this.mti = 1; this.mti < N; this.mti++) this.mt[this.mti] = (69069 * this.mt[this.mti - 1]); // & // 0xffffffff; // mag01[x] = x * MATRIX_A for x=0,1 this.mag01 = new int[2]; this.mag01[0] = 0x0; this.mag01[1] = MATRIX_A; } /** * An alternative, more complete, method of seeding the pseudo random number * generator. array must be an array of 624 ints, and they can be any value * as long as they're not *all* zero. * * @param array an array of 624 ints */ synchronized public void setSeed(final int[] array) { // it's always good style to call super -- // we'll use MT's canonical random number, but it doesn't // really matter. super.setSeed(4357); // Due to a bug in java.util.Random clear up to 1.2, we're // doing our own Gaussian variable. this.__haveNextNextGaussian = false; this.mt = new int[N]; System.arraycopy(array, 0, this.mt, 0, N); this.mti = N; // mag01[x] = x * MATRIX_A for x=0,1 this.mag01 = new int[2]; this.mag01[0] = 0x0; this.mag01[1] = MATRIX_A; } /** * Initalize the pseudo random number generator. Don't pass in a long that's * bigger than an int (Mersenne Twister only uses the first 32 bits for its * seed). * * @param seed the seed to use */ @Override synchronized public void setSeed(final long seed) { // it's always good style to call super super.setSeed(seed); // seed needs to be casted into an int first for this to work int _seed = (int) seed; // Due to a bug in java.util.Random clear up to 1.2, we're // doing our own Gaussian variable. this.__haveNextNextGaussian = false; this.mt = new int[N]; for (int i = 0; i < N; i++) { this.mt[i] = _seed & 0xffff0000; _seed = 69069 * _seed + 1; this.mt[i] |= (_seed & 0xffff0000) >>> 16; _seed = 69069 * _seed + 1; } this.mti = N; // mag01[x] = x * MATRIX_A for x=0,1 this.mag01 = new int[2]; this.mag01[0] = 0x0; this.mag01[1] = MATRIX_A; } /** * @param bits the number of bits to use * @return an integer with <i>bits</i> bits filled with a random number */ @Override synchronized protected int next(final int bits) { int y; if (this.mti >= N) // generate N words at one time { int kk; for (kk = 0; kk < N - M; kk++) { y = (this.mt[kk] & UPPER_MASK) | (this.mt[kk + 1] & LOWER_MASK); this.mt[kk] = this.mt[kk + M] ^ (y >>> 1) ^ this.mag01[y & 0x1]; } for (; kk < N - 1; kk++) { y = (this.mt[kk] & UPPER_MASK) | (this.mt[kk + 1] & LOWER_MASK); this.mt[kk] = this.mt[kk + (M - N)] ^ (y >>> 1) ^ this.mag01[y & 0x1]; } y = (this.mt[N - 1] & UPPER_MASK) | (this.mt[0] & LOWER_MASK); this.mt[N - 1] = this.mt[M - 1] ^ (y >>> 1) ^ this.mag01[y & 0x1]; this.mti = 0; } y = this.mt[this.mti++]; y ^= y >>> 11; // TEMPERING_SHIFT_U(y) y ^= (y << 7) & TEMPERING_MASK_B; // TEMPERING_SHIFT_S(y) y ^= (y << 15) & TEMPERING_MASK_C; // TEMPERING_SHIFT_T(y) y ^= (y >>> 18); // TEMPERING_SHIFT_L(y) return y >>> (32 - bits); // hope that's right! } /* * If you've got a truly old version of Java, you can omit these two next * methods. */ private synchronized void writeObject(final ObjectOutputStream out) throws IOException { // just so we're synchronized. out.defaultWriteObject(); } private synchronized void readObject(final ObjectInputStream in) throws IOException, ClassNotFoundException { // just so we're synchronized. in.defaultReadObject(); } /** * This generates a coin flip with a probability <tt>probability</tt> of * returning true, else returning false. <tt>probability</tt> must be * between 0.0 and 1.0, inclusive. Not as precise a random real event as * nextBoolean(double), but twice as fast. To explicitly use this, remember * you may need to cast to float first. * * @param probability the probability to use (between 0.0 and 1.0) * @return the coin flip result */ public boolean nextBoolean(final float probability) { if (probability < 0.0f || probability > 1.0f) throw new IllegalArgumentException("probability must be between 0.0 and 1.0 inclusive."); return nextFloat() < probability; } /** * This generates a coin flip with a probability <tt>probability</tt> of * returning true, else returning false. <tt>probability</tt> must be * between 0.0 and 1.0, inclusive. * * @param probability must be between 0.0 and 1.0 * @return the result of the coin flip */ public boolean nextBoolean(final double probability) { if (probability < 0.0 || probability > 1.0) throw new IllegalArgumentException("probability must be between 0.0 and 1.0 inclusive."); return nextDouble() < probability; } /** * A bug fix for all versions of the JDK. The JDK appears to use all four * bytes in an integer as independent byte values! Totally wrong. I've * submitted a bug report. * * @param bytes the bytes for which to get the next bytes (?) */ @Override public void nextBytes(final byte[] bytes) { for (int x = 0; x < bytes.length; x++) bytes[x] = (byte) next(8); } /** * For completeness' sake, though it's not in java.util.Random. * * @return the next char */ public char nextChar() { // chars are 16-bit UniCode values return (char) (next(16)); } /** * For completeness' sake, though it's not in java.util.Random. * * @return the next short */ public short nextShort() { return (short) (next(16)); } /** * For completeness' sake, though it's not in java.util.Random. * * @return the next byte */ public byte nextByte() { return (byte) (next(8)); } /** * A bug fix for all JDK code including 1.2. nextGaussian can theoretically * ask for the log of 0 and divide it by 0! See Java bug <a * href="http://developer.java.sun.com/developer/bugParade/bugs/4254501.html"> * http://developer.java.sun.com/developer/bugParade/bugs/4254501.html </a> * * @return the next Gaussian */ @Override synchronized public double nextGaussian() { if (this.__haveNextNextGaussian) { this.__haveNextNextGaussian = false; return this.__nextNextGaussian; } // (otherwise...) double v1, v2, s; do { v1 = 2 * nextDouble() - 1; // between -1.0 and 1.0 v2 = 2 * nextDouble() - 1; // between -1.0 and 1.0 s = v1 * v1 + v2 * v2; } while (s >= 1 || s == 0); double multiplier = Math.sqrt(-2 * Math.log(s) / s); this.__nextNextGaussian = v2 * multiplier; this.__haveNextNextGaussian = true; return v1 * multiplier; } }