AbstractWell.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.apache.commons.rng.core.source32;

  18. import java.util.Arrays;
  19. import org.apache.commons.rng.core.util.NumberFactory;

  20. /**
  21.  * This abstract class implements the WELL class of pseudo-random number
  22.  * generator from François Panneton, Pierre L'Ecuyer and Makoto
  23.  * Matsumoto.
  24.  * <p>
  25.  * This generator is described in a paper by Fran&ccedil;ois Panneton,
  26.  * Pierre L'Ecuyer and Makoto Matsumoto
  27.  * <a href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng.pdf">
  28.  * Improved Long-Period Generators Based on Linear Recurrences Modulo 2</a>
  29.  * ACM Transactions on Mathematical Software, 32, 1 (2006).
  30.  * The errata for the paper are in
  31.  * <a href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng-errata.txt">wellrng-errata.txt</a>.
  32.  * </p>
  33.  *
  34.  * @see <a href="http://www.iro.umontreal.ca/~panneton/WELLRNG.html">WELL Random number generator</a>
  35.  *
  36.  * @since 1.0
  37.  */
  38. public abstract class AbstractWell extends IntProvider {
  39.     /** Block size. */
  40.     private static final int BLOCK_SIZE = 32;
  41.     /** Current index in the bytes pool. */
  42.     protected int index;
  43.     /** Bytes pool. */
  44.     protected final int[] v;

  45.     /**
  46.      * Creates an instance with the given {@code seed}.
  47.      *
  48.      * @param k Number of bits in the pool (not necessarily a multiple of 32).
  49.      * @param seed Initial seed.
  50.      */
  51.     protected AbstractWell(final int k,
  52.                            final int[] seed) {
  53.         final int r = calculateBlockCount(k);
  54.         v = new int[r];
  55.         index = 0;

  56.         // Initialize the pool content.
  57.         setSeedInternal(seed);
  58.     }

  59.     /** {@inheritDoc} */
  60.     @Override
  61.     protected byte[] getStateInternal() {
  62.         final int[] s = Arrays.copyOf(v, v.length + 1);
  63.         s[v.length] = index;

  64.         return composeStateInternal(NumberFactory.makeByteArray(s),
  65.                                     super.getStateInternal());
  66.     }

  67.     /** {@inheritDoc} */
  68.     @Override
  69.     protected void setStateInternal(byte[] s) {
  70.         final byte[][] c = splitStateInternal(s, (v.length + 1) * 4);

  71.         final int[] tmp = NumberFactory.makeIntArray(c[0]);
  72.         System.arraycopy(tmp, 0, v, 0, v.length);
  73.         index = tmp[v.length];

  74.         super.setStateInternal(c[1]);
  75.     }

  76.     /**
  77.      * Initializes the generator with the given {@code seed}.
  78.      *
  79.      * @param seed Seed. Cannot be null.
  80.      */
  81.     private void setSeedInternal(final int[] seed) {
  82.         System.arraycopy(seed, 0, v, 0, Math.min(seed.length, v.length));

  83.         if (seed.length < v.length) {
  84.             for (int i = seed.length; i < v.length; ++i) {
  85.                 final long current = v[i - seed.length];
  86.                 v[i] = (int) ((1812433253L * (current ^ (current >> 30)) + i) & 0xffffffffL);
  87.             }
  88.         }

  89.         index = 0;
  90.     }

  91.     /**
  92.      * Calculate the number of 32-bits blocks.
  93.      *
  94.      * @param k Number of bits in the pool (not necessarily a multiple of 32).
  95.      * @return the number of 32-bits blocks.
  96.      */
  97.     private static int calculateBlockCount(final int k) {
  98.         // The bits pool contains k bits, k = r w - p where r is the number
  99.         // of w bits blocks, w is the block size (always 32 in the original paper)
  100.         // and p is the number of unused bits in the last block.
  101.         return (k + BLOCK_SIZE - 1) / BLOCK_SIZE;
  102.     }

  103.     /**
  104.      * Inner class used to store the indirection index table which is fixed for a given
  105.      * type of WELL class of pseudo-random number generator.
  106.      */
  107.     protected static final class IndexTable {
  108.         /** Index indirection table giving for each index its predecessor taking table size into account. */
  109.         private final int[] iRm1;
  110.         /** Index indirection table giving for each index its second predecessor taking table size into account. */
  111.         private final int[] iRm2;
  112.         /** Index indirection table giving for each index the value index + m1 taking table size into account. */
  113.         private final int[] i1;
  114.         /** Index indirection table giving for each index the value index + m2 taking table size into account. */
  115.         private final int[] i2;
  116.         /** Index indirection table giving for each index the value index + m3 taking table size into account. */
  117.         private final int[] i3;

  118.         /** Creates a new pre-calculated indirection index table.
  119.          * @param k number of bits in the pool (not necessarily a multiple of 32)
  120.          * @param m1 first parameter of the algorithm
  121.          * @param m2 second parameter of the algorithm
  122.          * @param m3 third parameter of the algorithm
  123.          */
  124.         public IndexTable(final int k, final int m1, final int m2, final int m3) {

  125.             final int r = calculateBlockCount(k);

  126.             // precompute indirection index tables. These tables are used for optimizing access
  127.             // they allow saving computations like "(j + r - 2) % r" with costly modulo operations
  128.             iRm1 = new int[r];
  129.             iRm2 = new int[r];
  130.             i1 = new int[r];
  131.             i2 = new int[r];
  132.             i3 = new int[r];
  133.             for (int j = 0; j < r; ++j) {
  134.                 iRm1[j] = (j + r - 1) % r;
  135.                 iRm2[j] = (j + r - 2) % r;
  136.                 i1[j] = (j + m1) % r;
  137.                 i2[j] = (j + m2) % r;
  138.                 i3[j] = (j + m3) % r;
  139.             }
  140.         }

  141.         /**
  142.          * Returns the predecessor of the given index modulo the table size.
  143.          * @param index the index to look at
  144.          * @return (index - 1) % table size
  145.          */
  146.         public int getIndexPred(final int index) {
  147.             return iRm1[index];
  148.         }

  149.         /**
  150.          * Returns the second predecessor of the given index modulo the table size.
  151.          * @param index the index to look at
  152.          * @return (index - 2) % table size
  153.          */
  154.         public int getIndexPred2(final int index) {
  155.             return iRm2[index];
  156.         }

  157.         /**
  158.          * Returns index + M1 modulo the table size.
  159.          * @param index the index to look at
  160.          * @return (index + M1) % table size
  161.          */
  162.         public int getIndexM1(final int index) {
  163.             return i1[index];
  164.         }

  165.         /**
  166.          * Returns index + M2 modulo the table size.
  167.          * @param index the index to look at
  168.          * @return (index + M2) % table size
  169.          */
  170.         public int getIndexM2(final int index) {
  171.             return i2[index];
  172.         }

  173.         /**
  174.          * Returns index + M3 modulo the table size.
  175.          * @param index the index to look at
  176.          * @return (index + M3) % table size
  177.          */
  178.         public int getIndexM3(final int index) {
  179.             return i3[index];
  180.         }
  181.     }
  182. }