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 */
017package org.apache.commons.rng.core.source32;
018
019import java.util.Arrays;
020import org.apache.commons.rng.core.util.NumberFactory;
021
022/**
023 * This abstract class implements the WELL class of pseudo-random number
024 * generator from François Panneton, Pierre L'Ecuyer and Makoto
025 * Matsumoto.
026 * <p>
027 * This generator is described in a paper by Fran&ccedil;ois Panneton,
028 * Pierre L'Ecuyer and Makoto Matsumoto
029 * <a href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng.pdf">
030 * Improved Long-Period Generators Based on Linear Recurrences Modulo 2</a>
031 * ACM Transactions on Mathematical Software, 32, 1 (2006).
032 * The errata for the paper are in
033 * <a href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng-errata.txt">wellrng-errata.txt</a>.
034 * </p>
035 *
036 * @see <a href="http://www.iro.umontreal.ca/~panneton/WELLRNG.html">WELL Random number generator</a>
037 *
038 * @since 1.0
039 */
040public abstract class AbstractWell extends IntProvider {
041    /** Block size. */
042    private static final int BLOCK_SIZE = 32;
043    /** Current index in the bytes pool. */
044    protected int index;
045    /** Bytes pool. */
046    protected final int[] v;
047
048    /**
049     * Creates an instance with the given {@code seed}.
050     *
051     * @param k Number of bits in the pool (not necessarily a multiple of 32).
052     * @param seed Initial seed.
053     */
054    protected AbstractWell(final int k,
055                           final int[] seed) {
056        final int r = calculateBlockCount(k);
057        v = new int[r];
058        index = 0;
059
060        // Initialize the pool content.
061        setSeedInternal(seed);
062    }
063
064    /** {@inheritDoc} */
065    @Override
066    protected byte[] getStateInternal() {
067        final int[] s = Arrays.copyOf(v, v.length + 1);
068        s[v.length] = index;
069
070        return composeStateInternal(NumberFactory.makeByteArray(s),
071                                    super.getStateInternal());
072    }
073
074    /** {@inheritDoc} */
075    @Override
076    protected void setStateInternal(byte[] s) {
077        final byte[][] c = splitStateInternal(s, (v.length + 1) * 4);
078
079        final int[] tmp = NumberFactory.makeIntArray(c[0]);
080        System.arraycopy(tmp, 0, v, 0, v.length);
081        index = tmp[v.length];
082
083        super.setStateInternal(c[1]);
084    }
085
086    /**
087     * Initializes the generator with the given {@code seed}.
088     *
089     * @param seed Seed. Cannot be null.
090     */
091    private void setSeedInternal(final int[] seed) {
092        System.arraycopy(seed, 0, v, 0, Math.min(seed.length, v.length));
093
094        if (seed.length < v.length) {
095            for (int i = seed.length; i < v.length; ++i) {
096                final long current = v[i - seed.length];
097                v[i] = (int) ((1812433253L * (current ^ (current >> 30)) + i) & 0xffffffffL);
098            }
099        }
100
101        index = 0;
102    }
103
104    /**
105     * Calculate the number of 32-bits blocks.
106     *
107     * @param k Number of bits in the pool (not necessarily a multiple of 32).
108     * @return the number of 32-bits blocks.
109     */
110    private static int calculateBlockCount(final int k) {
111        // The bits pool contains k bits, k = r w - p where r is the number
112        // of w bits blocks, w is the block size (always 32 in the original paper)
113        // and p is the number of unused bits in the last block.
114        return (k + BLOCK_SIZE - 1) / BLOCK_SIZE;
115    }
116
117    /**
118     * Inner class used to store the indirection index table which is fixed for a given
119     * type of WELL class of pseudo-random number generator.
120     */
121    protected static final class IndexTable {
122        /** Index indirection table giving for each index its predecessor taking table size into account. */
123        private final int[] iRm1;
124        /** Index indirection table giving for each index its second predecessor taking table size into account. */
125        private final int[] iRm2;
126        /** Index indirection table giving for each index the value index + m1 taking table size into account. */
127        private final int[] i1;
128        /** Index indirection table giving for each index the value index + m2 taking table size into account. */
129        private final int[] i2;
130        /** Index indirection table giving for each index the value index + m3 taking table size into account. */
131        private final int[] i3;
132
133        /** Creates a new pre-calculated indirection index table.
134         * @param k number of bits in the pool (not necessarily a multiple of 32)
135         * @param m1 first parameter of the algorithm
136         * @param m2 second parameter of the algorithm
137         * @param m3 third parameter of the algorithm
138         */
139        public IndexTable(final int k, final int m1, final int m2, final int m3) {
140
141            final int r = calculateBlockCount(k);
142
143            // precompute indirection index tables. These tables are used for optimizing access
144            // they allow saving computations like "(j + r - 2) % r" with costly modulo operations
145            iRm1 = new int[r];
146            iRm2 = new int[r];
147            i1 = new int[r];
148            i2 = new int[r];
149            i3 = new int[r];
150            for (int j = 0; j < r; ++j) {
151                iRm1[j] = (j + r - 1) % r;
152                iRm2[j] = (j + r - 2) % r;
153                i1[j] = (j + m1) % r;
154                i2[j] = (j + m2) % r;
155                i3[j] = (j + m3) % r;
156            }
157        }
158
159        /**
160         * Returns the predecessor of the given index modulo the table size.
161         * @param index the index to look at
162         * @return (index - 1) % table size
163         */
164        public int getIndexPred(final int index) {
165            return iRm1[index];
166        }
167
168        /**
169         * Returns the second predecessor of the given index modulo the table size.
170         * @param index the index to look at
171         * @return (index - 2) % table size
172         */
173        public int getIndexPred2(final int index) {
174            return iRm2[index];
175        }
176
177        /**
178         * Returns index + M1 modulo the table size.
179         * @param index the index to look at
180         * @return (index + M1) % table size
181         */
182        public int getIndexM1(final int index) {
183            return i1[index];
184        }
185
186        /**
187         * Returns index + M2 modulo the table size.
188         * @param index the index to look at
189         * @return (index + M2) % table size
190         */
191        public int getIndexM2(final int index) {
192            return i2[index];
193        }
194
195        /**
196         * Returns index + M3 modulo the table size.
197         * @param index the index to look at
198         * @return (index + M3) % table size
199         */
200        public int getIndexM3(final int index) {
201            return i3[index];
202        }
203    }
204}