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&ccedil;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