1dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/* 2dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Licensed to the Apache Software Foundation (ASF) under one or more 3dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * contributor license agreements. See the NOTICE file distributed with 4dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * this work for additional information regarding copyright ownership. 5dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The ASF licenses this file to You under the Apache License, Version 2.0 6dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * (the "License"); you may not use this file except in compliance with 7dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the License. You may obtain a copy of the License at 8dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 9dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * http://www.apache.org/licenses/LICENSE-2.0 10dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 11dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Unless required by applicable law or agreed to in writing, software 12dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * distributed under the License is distributed on an "AS IS" BASIS, 13dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * See the License for the specific language governing permissions and 15dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * limitations under the License. 16dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 17dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpackage org.apache.commons.math.random; 18dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 19dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport java.io.Serializable; 20dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 21dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 22dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/** This abstract class implements the WELL class of pseudo-random number generator 23dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * from François Panneton, Pierre L'Ecuyer and Makoto Matsumoto. 24dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 25dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>This generator is described in a paper by François Panneton, 26dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Pierre L'Ecuyer and Makoto Matsumoto <a 27dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng.pdf">Improved 28dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Long-Period Generators Based on Linear Recurrences Modulo 2</a> ACM 29dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Transactions on Mathematical Software, 32, 1 (2006). The errata for the paper 30dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * are in <a href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng-errata.txt">wellrng-errata.txt</a>.</p> 31dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 32dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @see <a href="http://www.iro.umontreal.ca/~panneton/WELLRNG.html">WELL Random number generator</a> 33dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @version $Revision: 1003892 $ $Date: 2010-10-02 23:28:56 +0200 (sam. 02 oct. 2010) $ 34dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.2 35dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 36dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 37dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpublic abstract class AbstractWell extends BitsStreamGenerator implements Serializable { 38dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 39dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Serializable version identifier. */ 40dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static final long serialVersionUID = -817701723016583596L; 41dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 42dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Current index in the bytes pool. */ 43dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected int index; 44dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 45dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Bytes pool. */ 46dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected final int[] v; 47dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 48dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Index indirection table giving for each index its predecessor taking table size into account. */ 49dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected final int[] iRm1; 50dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 51dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Index indirection table giving for each index its second predecessor taking table size into account. */ 52dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected final int[] iRm2; 53dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 54dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Index indirection table giving for each index the value index + m1 taking table size into account. */ 55dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected final int[] i1; 56dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 57dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Index indirection table giving for each index the value index + m2 taking table size into account. */ 58dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected final int[] i2; 59dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 60dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Index indirection table giving for each index the value index + m3 taking table size into account. */ 61dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected final int[] i3; 62dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 63dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Creates a new random number generator. 64dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>The instance is initialized using the current time as the 65dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * seed.</p> 66dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param k number of bits in the pool (not necessarily a multiple of 32) 67dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param m1 first parameter of the algorithm 68dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param m2 second parameter of the algorithm 69dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param m3 third parameter of the algorithm 70dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 71dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected AbstractWell(final int k, final int m1, final int m2, final int m3) { 72dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this(k, m1, m2, m3, System.currentTimeMillis()); 73dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 74dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 75dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Creates a new random number generator using a single int seed. 76dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param k number of bits in the pool (not necessarily a multiple of 32) 77dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param m1 first parameter of the algorithm 78dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param m2 second parameter of the algorithm 79dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param m3 third parameter of the algorithm 80dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param seed the initial seed (32 bits integer) 81dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 82dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected AbstractWell(final int k, final int m1, final int m2, final int m3, final int seed) { 83dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this(k, m1, m2, m3, new int[] { seed }); 84dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 85dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 86dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Creates a new random number generator using an int array seed. 87dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param k number of bits in the pool (not necessarily a multiple of 32) 88dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param m1 first parameter of the algorithm 89dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param m2 second parameter of the algorithm 90dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param m3 third parameter of the algorithm 91dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param seed the initial seed (32 bits integers array), if null 92dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the seed of the generator will be related to the current time 93dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 94dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected AbstractWell(final int k, final int m1, final int m2, final int m3, final int[] seed) { 95dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 96dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // the bits pool contains k bits, k = r w - p where r is the number 97dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // of w bits blocks, w is the block size (always 32 in the original paper) 98dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // and p is the number of unused bits in the last block 99dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int w = 32; 100dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int r = (k + w - 1) / w; 101dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.v = new int[r]; 102dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.index = 0; 103dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 104dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // precompute indirection index tables. These tables are used for optimizing access 105dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // they allow saving computations like "(j + r - 2) % r" with costly modulo operations 106dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond iRm1 = new int[r]; 107dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond iRm2 = new int[r]; 108dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond i1 = new int[r]; 109dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond i2 = new int[r]; 110dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond i3 = new int[r]; 111dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int j = 0; j < r; ++j) { 112dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond iRm1[j] = (j + r - 1) % r; 113dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond iRm2[j] = (j + r - 2) % r; 114dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond i1[j] = (j + m1) % r; 115dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond i2[j] = (j + m2) % r; 116dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond i3[j] = (j + m3) % r; 117dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 118dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 119dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // initialize the pool content 120dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond setSeed(seed); 121dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 122dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 123dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 124dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Creates a new random number generator using a single long seed. 125dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param k number of bits in the pool (not necessarily a multiple of 32) 126dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param m1 first parameter of the algorithm 127dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param m2 second parameter of the algorithm 128dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param m3 third parameter of the algorithm 129dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param seed the initial seed (64 bits integer) 130dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 131dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected AbstractWell(final int k, final int m1, final int m2, final int m3, final long seed) { 132dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this(k, m1, m2, m3, new int[] { (int) (seed >>> 32), (int) (seed & 0xffffffffl) }); 133dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 134dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 135dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Reinitialize the generator as if just built with the given int seed. 136dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>The state of the generator is exactly the same as a new 137dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * generator built with the same seed.</p> 138dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param seed the initial seed (32 bits integer) 139dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 140dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond @Override 141dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public void setSeed(final int seed) { 142dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond setSeed(new int[] { seed }); 143dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 144dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 145dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Reinitialize the generator as if just built with the given int array seed. 146dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>The state of the generator is exactly the same as a new 147dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * generator built with the same seed.</p> 148dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param seed the initial seed (32 bits integers array), if null 149dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the seed of the generator will be related to the current time 150dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 151dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond @Override 152dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public void setSeed(final int[] seed) { 153dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 154dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (seed == null) { 155dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond setSeed(System.currentTimeMillis()); 156dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return; 157dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 158dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 159dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond System.arraycopy(seed, 0, v, 0, Math.min(seed.length, v.length)); 160dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 161dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (seed.length < v.length) { 162dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = seed.length; i < v.length; ++i) { 163dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final long l = v[i - seed.length]; 164dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond v[i] = (int) ((1812433253l * (l ^ (l >> 30)) + i) & 0xffffffffL); 165dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 166dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 167dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 168dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond index = 0; 169dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 170dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 171dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 172dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Reinitialize the generator as if just built with the given long seed. 173dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>The state of the generator is exactly the same as a new 174dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * generator built with the same seed.</p> 175dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param seed the initial seed (64 bits integer) 176dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 177dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond @Override 178dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public void setSeed(final long seed) { 179dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond setSeed(new int[] { (int) (seed >>> 32), (int) (seed & 0xffffffffl) }); 180dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 181dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 182dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 183dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond @Override 184dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected abstract int next(final int bits); 185dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 186dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond} 187