MarsagliaTsangWangDiscreteSampler.java

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

import org.apache.commons.rng.UniformRandomProvider;

/**
 * Sampler for a discrete distribution using an optimised look-up table.
 *
 * <ul>
 *  <li>
 *   The method requires 30-bit integer probabilities that sum to 2<sup>30</sup> as described
 *   in George Marsaglia, Wai Wan Tsang, Jingbo Wang (2004) Fast Generation of Discrete
 *   Random Variables. Journal of Statistical Software. Vol. 11, Issue. 3, pp. 1-11.
 *  </li>
 * </ul>
 *
 * <p>Sampling uses 1 call to {@link UniformRandomProvider#nextInt()}.</p>
 *
 * <p>Memory requirements depend on the maximum number of possible sample values, {@code n},
 * and the values for the probabilities. Storage is optimised for {@code n}. The worst case
 * scenario is a uniform distribution of the maximum sample size. This is capped at 0.06MB for
 * {@code n <= } 2<sup>8</sup>, 17.0MB for {@code n <= } 2<sup>16</sup>, and 4.3GB for
 * {@code n <=} 2<sup>30</sup>. Realistic requirements will be in the kB range.</p>
 *
 * <p>The sampler supports the following distributions:</p>
 *
 * <ul>
 *  <li>Enumerated distribution (probabilities must be provided for each sample)
 *  <li>Poisson distribution up to {@code mean = 1024}
 *  <li>Binomial distribution up to {@code trials = 65535}
 * </ul>
 *
 * @see <a href="http://dx.doi.org/10.18637/jss.v011.i03">Margsglia, et al (2004) JSS Vol.
 * 11, Issue 3</a>
 * @since 1.3
 */
public final class MarsagliaTsangWangDiscreteSampler {
    /** The value 2<sup>8</sup> as an {@code int}. */
    private static final int INT_8 = 1 << 8;
    /** The value 2<sup>16</sup> as an {@code int}. */
    private static final int INT_16 = 1 << 16;
    /** The value 2<sup>30</sup> as an {@code int}. */
    private static final int INT_30 = 1 << 30;
    /** The value 2<sup>31</sup> as a {@code double}. */
    private static final double DOUBLE_31 = 1L << 31;

    // =========================================================================
    // Implementation note:
    //
    // This sampler uses prepared look-up tables that are searched using a single
    // random int variate. The look-up tables contain the sample value. The tables
    // are constructed using probabilities that sum to 2^30. The original paper
    // by Marsaglia, et al (2004) describes the use of 5, 3, or 2 look-up tables
    // indexed using digits of base 2^6, 2^10 or 2^15. Currently only base 64 (2^6)
    // is supported using 5 look-up tables.
    //
    // The implementations use 8, 16 or 32 bit storage tables to support different
    // distribution sizes with optimal storage. Separate class implementations of
    // the same algorithm allow array storage to be accessed directly from 1D tables.
    // This provides a performance gain over using: abstracted storage accessed via
    // an interface; or a single 2D table.
    //
    // To allow the optimal implementation to be chosen the sampler is created
    // using factory methods. The sampler supports any probability distribution
    // when provided via an array of probabilities and the Poisson and Binomial
    // distributions for a restricted set of parameters. The restrictions are
    // imposed by the requirement to compute the entire probability distribution
    // from the controlling parameter(s) using a recursive method. Factory
    // constructors return a SharedStateDiscreteSampler instance. Each distribution
    // type is contained in an inner class.
    // =========================================================================

    /**
     * The base class for Marsaglia-Tsang-Wang samplers.
     */
    private abstract static class AbstractMarsagliaTsangWangDiscreteSampler
            implements SharedStateDiscreteSampler {
        /** Underlying source of randomness. */
        protected final UniformRandomProvider rng;

        /** The name of the distribution. */
        private final String distributionName;

        /**
         * @param rng Generator of uniformly distributed random numbers.
         * @param distributionName Distribution name.
         */
        AbstractMarsagliaTsangWangDiscreteSampler(UniformRandomProvider rng,
                                                  String distributionName) {
            this.rng = rng;
            this.distributionName = distributionName;
        }

        /**
         * @param rng Generator of uniformly distributed random numbers.
         * @param source Source to copy.
         */
        AbstractMarsagliaTsangWangDiscreteSampler(UniformRandomProvider rng,
                                                  AbstractMarsagliaTsangWangDiscreteSampler source) {
            this.rng = rng;
            this.distributionName = source.distributionName;
        }

        /** {@inheritDoc} */
        @Override
        public String toString() {
            return "Marsaglia Tsang Wang " + distributionName + " deviate [" + rng.toString() + "]";
        }
    }

    /**
     * An implementation for the sample algorithm based on the decomposition of the
     * index in the range {@code [0,2^30)} into 5 base-64 digits with 8-bit backing storage.
     */
    private static final class MarsagliaTsangWangBase64Int8DiscreteSampler
        extends AbstractMarsagliaTsangWangDiscreteSampler {
        /** The mask to convert a {@code byte} to an unsigned 8-bit integer. */
        private static final int MASK = 0xff;

        /** Limit for look-up table 1. */
        private final int t1;
        /** Limit for look-up table 2. */
        private final int t2;
        /** Limit for look-up table 3. */
        private final int t3;
        /** Limit for look-up table 4. */
        private final int t4;

        /** Look-up table table1. */
        private final byte[] table1;
        /** Look-up table table2. */
        private final byte[] table2;
        /** Look-up table table3. */
        private final byte[] table3;
        /** Look-up table table4. */
        private final byte[] table4;
        /** Look-up table table5. */
        private final byte[] table5;

        /**
         * @param rng Generator of uniformly distributed random numbers.
         * @param distributionName Distribution name.
         * @param prob The probabilities.
         * @param offset The offset (must be positive).
         */
        MarsagliaTsangWangBase64Int8DiscreteSampler(UniformRandomProvider rng,
                                                    String distributionName,
                                                    int[] prob,
                                                    int offset) {
            super(rng, distributionName);

            // Get table sizes for each base-64 digit
            int n1 = 0;
            int n2 = 0;
            int n3 = 0;
            int n4 = 0;
            int n5 = 0;
            for (final int m : prob) {
                n1 += getBase64Digit(m, 1);
                n2 += getBase64Digit(m, 2);
                n3 += getBase64Digit(m, 3);
                n4 += getBase64Digit(m, 4);
                n5 += getBase64Digit(m, 5);
            }

            table1 = new byte[n1];
            table2 = new byte[n2];
            table3 = new byte[n3];
            table4 = new byte[n4];
            table5 = new byte[n5];

            // Compute offsets
            t1 = n1 << 24;
            t2 = t1 + (n2 << 18);
            t3 = t2 + (n3 << 12);
            t4 = t3 + (n4 << 6);
            n1 = n2 = n3 = n4 = n5 = 0;

            // Fill tables
            for (int i = 0; i < prob.length; i++) {
                final int m = prob[i];
                // Primitive type conversion will extract lower 8 bits
                final byte k = (byte) (i + offset);
                n1 = fill(table1, n1, n1 + getBase64Digit(m, 1), k);
                n2 = fill(table2, n2, n2 + getBase64Digit(m, 2), k);
                n3 = fill(table3, n3, n3 + getBase64Digit(m, 3), k);
                n4 = fill(table4, n4, n4 + getBase64Digit(m, 4), k);
                n5 = fill(table5, n5, n5 + getBase64Digit(m, 5), k);
            }
        }

        /**
         * @param rng Generator of uniformly distributed random numbers.
         * @param source Source to copy.
         */
        private MarsagliaTsangWangBase64Int8DiscreteSampler(UniformRandomProvider rng,
                MarsagliaTsangWangBase64Int8DiscreteSampler source) {
            super(rng, source);
            t1 = source.t1;
            t2 = source.t2;
            t3 = source.t3;
            t4 = source.t4;
            table1 = source.table1;
            table2 = source.table2;
            table3 = source.table3;
            table4 = source.table4;
            table5 = source.table5;
        }

        /**
         * Fill the table with the value.
         *
         * @param table Table.
         * @param from Lower bound index (inclusive)
         * @param to Upper bound index (exclusive)
         * @param value Value.
         * @return the upper bound index
         */
        private static int fill(byte[] table, int from, int to, byte value) {
            for (int i = from; i < to; i++) {
                table[i] = value;
            }
            return to;
        }

        @Override
        public int sample() {
            final int j = rng.nextInt() >>> 2;
            if (j < t1) {
                return table1[j >>> 24] & MASK;
            }
            if (j < t2) {
                return table2[(j - t1) >>> 18] & MASK;
            }
            if (j < t3) {
                return table3[(j - t2) >>> 12] & MASK;
            }
            if (j < t4) {
                return table4[(j - t3) >>> 6] & MASK;
            }
            // Note the tables are filled on the assumption that the sum of the probabilities.
            // is >=2^30. If this is not true then the final table table5 will be smaller by the
            // difference. So the tables *must* be constructed correctly.
            return table5[j - t4] & MASK;
        }

        @Override
        public SharedStateDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) {
            return new MarsagliaTsangWangBase64Int8DiscreteSampler(rng, this);
        }
    }

    /**
     * An implementation for the sample algorithm based on the decomposition of the
     * index in the range {@code [0,2^30)} into 5 base-64 digits with 16-bit backing storage.
     */
    private static final class MarsagliaTsangWangBase64Int16DiscreteSampler
        extends AbstractMarsagliaTsangWangDiscreteSampler {
        /** The mask to convert a {@code byte} to an unsigned 16-bit integer. */
        private static final int MASK = 0xffff;

        /** Limit for look-up table 1. */
        private final int t1;
        /** Limit for look-up table 2. */
        private final int t2;
        /** Limit for look-up table 3. */
        private final int t3;
        /** Limit for look-up table 4. */
        private final int t4;

        /** Look-up table table1. */
        private final short[] table1;
        /** Look-up table table2. */
        private final short[] table2;
        /** Look-up table table3. */
        private final short[] table3;
        /** Look-up table table4. */
        private final short[] table4;
        /** Look-up table table5. */
        private final short[] table5;

        /**
         * @param rng Generator of uniformly distributed random numbers.
         * @param distributionName Distribution name.
         * @param prob The probabilities.
         * @param offset The offset (must be positive).
         */
        MarsagliaTsangWangBase64Int16DiscreteSampler(UniformRandomProvider rng,
                                                     String distributionName,
                                                     int[] prob,
                                                     int offset) {
            super(rng, distributionName);

            // Get table sizes for each base-64 digit
            int n1 = 0;
            int n2 = 0;
            int n3 = 0;
            int n4 = 0;
            int n5 = 0;
            for (final int m : prob) {
                n1 += getBase64Digit(m, 1);
                n2 += getBase64Digit(m, 2);
                n3 += getBase64Digit(m, 3);
                n4 += getBase64Digit(m, 4);
                n5 += getBase64Digit(m, 5);
            }

            table1 = new short[n1];
            table2 = new short[n2];
            table3 = new short[n3];
            table4 = new short[n4];
            table5 = new short[n5];

            // Compute offsets
            t1 = n1 << 24;
            t2 = t1 + (n2 << 18);
            t3 = t2 + (n3 << 12);
            t4 = t3 + (n4 << 6);
            n1 = n2 = n3 = n4 = n5 = 0;

            // Fill tables
            for (int i = 0; i < prob.length; i++) {
                final int m = prob[i];
                // Primitive type conversion will extract lower 16 bits
                final short k = (short) (i + offset);
                n1 = fill(table1, n1, n1 + getBase64Digit(m, 1), k);
                n2 = fill(table2, n2, n2 + getBase64Digit(m, 2), k);
                n3 = fill(table3, n3, n3 + getBase64Digit(m, 3), k);
                n4 = fill(table4, n4, n4 + getBase64Digit(m, 4), k);
                n5 = fill(table5, n5, n5 + getBase64Digit(m, 5), k);
            }
        }

        /**
         * @param rng Generator of uniformly distributed random numbers.
         * @param source Source to copy.
         */
        private MarsagliaTsangWangBase64Int16DiscreteSampler(UniformRandomProvider rng,
                MarsagliaTsangWangBase64Int16DiscreteSampler source) {
            super(rng, source);
            t1 = source.t1;
            t2 = source.t2;
            t3 = source.t3;
            t4 = source.t4;
            table1 = source.table1;
            table2 = source.table2;
            table3 = source.table3;
            table4 = source.table4;
            table5 = source.table5;
        }

        /**
         * Fill the table with the value.
         *
         * @param table Table.
         * @param from Lower bound index (inclusive)
         * @param to Upper bound index (exclusive)
         * @param value Value.
         * @return the upper bound index
         */
        private static int fill(short[] table, int from, int to, short value) {
            for (int i = from; i < to; i++) {
                table[i] = value;
            }
            return to;
        }

        @Override
        public int sample() {
            final int j = rng.nextInt() >>> 2;
            if (j < t1) {
                return table1[j >>> 24] & MASK;
            }
            if (j < t2) {
                return table2[(j - t1) >>> 18] & MASK;
            }
            if (j < t3) {
                return table3[(j - t2) >>> 12] & MASK;
            }
            if (j < t4) {
                return table4[(j - t3) >>> 6] & MASK;
            }
            // Note the tables are filled on the assumption that the sum of the probabilities.
            // is >=2^30. If this is not true then the final table table5 will be smaller by the
            // difference. So the tables *must* be constructed correctly.
            return table5[j - t4] & MASK;
        }

        @Override
        public SharedStateDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) {
            return new MarsagliaTsangWangBase64Int16DiscreteSampler(rng, this);
        }
    }

    /**
     * An implementation for the sample algorithm based on the decomposition of the
     * index in the range {@code [0,2^30)} into 5 base-64 digits with 32-bit backing storage.
     */
    private static final class MarsagliaTsangWangBase64Int32DiscreteSampler
        extends AbstractMarsagliaTsangWangDiscreteSampler {
        /** Limit for look-up table 1. */
        private final int t1;
        /** Limit for look-up table 2. */
        private final int t2;
        /** Limit for look-up table 3. */
        private final int t3;
        /** Limit for look-up table 4. */
        private final int t4;

        /** Look-up table table1. */
        private final int[] table1;
        /** Look-up table table2. */
        private final int[] table2;
        /** Look-up table table3. */
        private final int[] table3;
        /** Look-up table table4. */
        private final int[] table4;
        /** Look-up table table5. */
        private final int[] table5;

        /**
         * @param rng Generator of uniformly distributed random numbers.
         * @param distributionName Distribution name.
         * @param prob The probabilities.
         * @param offset The offset (must be positive).
         */
        MarsagliaTsangWangBase64Int32DiscreteSampler(UniformRandomProvider rng,
                                                     String distributionName,
                                                     int[] prob,
                                                     int offset) {
            super(rng, distributionName);

            // Get table sizes for each base-64 digit
            int n1 = 0;
            int n2 = 0;
            int n3 = 0;
            int n4 = 0;
            int n5 = 0;
            for (final int m : prob) {
                n1 += getBase64Digit(m, 1);
                n2 += getBase64Digit(m, 2);
                n3 += getBase64Digit(m, 3);
                n4 += getBase64Digit(m, 4);
                n5 += getBase64Digit(m, 5);
            }

            table1 = new int[n1];
            table2 = new int[n2];
            table3 = new int[n3];
            table4 = new int[n4];
            table5 = new int[n5];

            // Compute offsets
            t1 = n1 << 24;
            t2 = t1 + (n2 << 18);
            t3 = t2 + (n3 << 12);
            t4 = t3 + (n4 << 6);
            n1 = n2 = n3 = n4 = n5 = 0;

            // Fill tables
            for (int i = 0; i < prob.length; i++) {
                final int m = prob[i];
                final int k = i + offset;
                n1 = fill(table1, n1, n1 + getBase64Digit(m, 1), k);
                n2 = fill(table2, n2, n2 + getBase64Digit(m, 2), k);
                n3 = fill(table3, n3, n3 + getBase64Digit(m, 3), k);
                n4 = fill(table4, n4, n4 + getBase64Digit(m, 4), k);
                n5 = fill(table5, n5, n5 + getBase64Digit(m, 5), k);
            }
        }

        /**
         * @param rng Generator of uniformly distributed random numbers.
         * @param source Source to copy.
         */
        private MarsagliaTsangWangBase64Int32DiscreteSampler(UniformRandomProvider rng,
                MarsagliaTsangWangBase64Int32DiscreteSampler source) {
            super(rng, source);
            t1 = source.t1;
            t2 = source.t2;
            t3 = source.t3;
            t4 = source.t4;
            table1 = source.table1;
            table2 = source.table2;
            table3 = source.table3;
            table4 = source.table4;
            table5 = source.table5;
        }

        /**
         * Fill the table with the value.
         *
         * @param table Table.
         * @param from Lower bound index (inclusive)
         * @param to Upper bound index (exclusive)
         * @param value Value.
         * @return the upper bound index
         */
        private static int fill(int[] table, int from, int to, int value) {
            for (int i = from; i < to; i++) {
                table[i] = value;
            }
            return to;
        }

        @Override
        public int sample() {
            final int j = rng.nextInt() >>> 2;
            if (j < t1) {
                return table1[j >>> 24];
            }
            if (j < t2) {
                return table2[(j - t1) >>> 18];
            }
            if (j < t3) {
                return table3[(j - t2) >>> 12];
            }
            if (j < t4) {
                return table4[(j - t3) >>> 6];
            }
            // Note the tables are filled on the assumption that the sum of the probabilities.
            // is >=2^30. If this is not true then the final table table5 will be smaller by the
            // difference. So the tables *must* be constructed correctly.
            return table5[j - t4];
        }

        @Override
        public SharedStateDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) {
            return new MarsagliaTsangWangBase64Int32DiscreteSampler(rng, this);
        }
    }



    /** Class contains only static methods. */
    private MarsagliaTsangWangDiscreteSampler() {}

    /**
     * Gets the k<sup>th</sup> base 64 digit of {@code m}.
     *
     * @param m the value m.
     * @param k the digit.
     * @return the base 64 digit
     */
    private static int getBase64Digit(int m, int k) {
        return (m >>> (30 - 6 * k)) & 63;
    }

    /**
     * Convert the probability to an integer in the range [0,2^30]. This is the numerator of
     * a fraction with assumed denominator 2<sup>30</sup>.
     *
     * @param p Probability.
     * @return the fraction numerator
     */
    private static int toUnsignedInt30(double p) {
        return (int) (p * INT_30 + 0.5);
    }

    /**
     * Create a new instance for probabilities {@code p(i)} where the sample value {@code x} is
     * {@code i + offset}.
     *
     * <p>The sum of the probabilities must be {@code >=} 2<sup>30</sup>. Only the
     * values for cumulative probability up to 2<sup>30</sup> will be sampled.</p>
     *
     * @param rng Generator of uniformly distributed random numbers.
     * @param distributionName Distribution name.
     * @param prob The probabilities.
     * @param offset The offset (must be positive).
     * @return Sampler.
     */
    private static SharedStateDiscreteSampler createSampler(UniformRandomProvider rng,
                                                            String distributionName,
                                                            int[] prob,
                                                            int offset) {
        // Note: No argument checks for private method.

        // Choose implementation based on the maximum index
        final int maxIndex = prob.length + offset - 1;
        if (maxIndex < INT_8) {
            return new MarsagliaTsangWangBase64Int8DiscreteSampler(rng, distributionName, prob, offset);
        }
        if (maxIndex < INT_16) {
            return new MarsagliaTsangWangBase64Int16DiscreteSampler(rng, distributionName, prob, offset);
        }
        return new MarsagliaTsangWangBase64Int32DiscreteSampler(rng, distributionName, prob, offset);
    }

    // =========================================================================
    // The following public classes provide factory methods to construct a sampler for:
    // - Enumerated probability distribution (from provided double[] probabilities)
    // - Poisson distribution for mean <= 1024
    // - Binomial distribution for trials <= 65535
    // =========================================================================

    /**
     * Create a sampler for an enumerated distribution of {@code n} values each with an
     * associated probability.
     * The samples corresponding to each probability are assumed to be a natural sequence
     * starting at zero.
     */
    public static final class Enumerated {
        /** The name of the enumerated probability distribution. */
        private static final String ENUMERATED_NAME = "Enumerated";

        /** Class contains only static methods. */
        private Enumerated() {}

        /**
         * Creates a sampler for an enumerated distribution of {@code n} values each with an
         * associated probability.
         *
         * <p>The probabilities will be normalised using their sum. The only requirement
         * is the sum is positive.</p>
         *
         * <p>The sum of the probabilities is normalised to 2<sup>30</sup>. Note that
         * probabilities are adjusted to the nearest 2<sup>-30</sup> due to round-off during
         * the normalisation conversion. Consequently any probability less than 2<sup>-31</sup>
         * will not be observed in samples.</p>
         *
         * @param rng Generator of uniformly distributed random numbers.
         * @param probabilities The list of probabilities.
         * @return Sampler.
         * @throws IllegalArgumentException if {@code probabilities} is null or empty, a
         * probability is negative, infinite or {@code NaN}, or the sum of all
         * probabilities is not strictly positive.
         */
        public static SharedStateDiscreteSampler of(UniformRandomProvider rng,
                                                    double[] probabilities) {
            return createSampler(rng, ENUMERATED_NAME, normaliseProbabilities(probabilities), 0);
        }

        /**
         * Normalise the probabilities to integers that sum to 2<sup>30</sup>.
         *
         * @param probabilities The list of probabilities.
         * @return the normalised probabilities.
         * @throws IllegalArgumentException if {@code probabilities} is null or empty, a
         * probability is negative, infinite or {@code NaN}, or the sum of all
         * probabilities is not strictly positive.
         */
        private static int[] normaliseProbabilities(double[] probabilities) {
            final double sumProb = InternalUtils.validateProbabilities(probabilities);

            // Compute the normalisation: 2^30 / sum
            final double normalisation = INT_30 / sumProb;
            final int[] prob = new int[probabilities.length];
            int sum = 0;
            int max = 0;
            int mode = 0;
            for (int i = 0; i < prob.length; i++) {
                // Add 0.5 for rounding
                final int p = (int) (probabilities[i] * normalisation + 0.5);
                sum += p;
                // Find the mode (maximum probability)
                if (max < p) {
                    max = p;
                    mode = i;
                }
                prob[i] = p;
            }

            // The sum must be >= 2^30.
            // Here just compensate the difference onto the highest probability.
            prob[mode] += INT_30 - sum;

            return prob;
        }
    }

    /**
     * Create a sampler for the Poisson distribution.
     */
    public static final class Poisson {
        /** The name of the Poisson distribution. */
        private static final String POISSON_NAME = "Poisson";

        /**
         * Upper bound on the mean for the Poisson distribution.
         *
         * <p>The original source code provided in Marsaglia, et al (2004) has no explicit
         * limit but the code fails at mean {@code >= 1941} as the transform to compute p(x=mode)
         * produces infinity. Use a conservative limit of 1024.</p>
         */

        private static final double MAX_MEAN = 1024;
        /**
         * The threshold for the mean of the Poisson distribution to switch the method used
         * to compute the probabilities. This is taken from the example software provided by
         * Marsaglia, et al (2004).
         */
        private static final double MEAN_THRESHOLD = 21.4;

        /** Class contains only static methods. */
        private Poisson() {}

        /**
         * Creates a sampler for the Poisson distribution.
         *
         * <p>Any probability less than 2<sup>-31</sup> will not be observed in samples.</p>
         *
         * <p>Storage requirements depend on the tabulated probability values. Example storage
         * requirements are listed below.</p>
         *
         * <pre>
         * mean      table size     kB
         * 0.25      882            0.88
         * 0.5       1135           1.14
         * 1         1200           1.20
         * 2         1451           1.45
         * 4         1955           1.96
         * 8         2961           2.96
         * 16        4410           4.41
         * 32        6115           6.11
         * 64        8499           8.50
         * 128       11528          11.53
         * 256       15935          31.87
         * 512       20912          41.82
         * 1024      30614          61.23
         * </pre>
         *
         * <p>Note: Storage changes to 2 bytes per index between {@code mean=128} and {@code mean=256}.</p>
         *
         * @param rng Generator of uniformly distributed random numbers.
         * @param mean Mean.
         * @return Sampler.
         * @throws IllegalArgumentException if {@code mean <= 0} or {@code mean > 1024}.
         */
        public static SharedStateDiscreteSampler of(UniformRandomProvider rng,
                                                    double mean) {
            validatePoissonDistributionParameters(mean);

            // Create the distribution either from X=0 or from X=mode when the mean is high.
            return mean < MEAN_THRESHOLD ?
                createPoissonDistributionFromX0(rng, mean) :
                createPoissonDistributionFromXMode(rng, mean);
        }

        /**
         * Validate the Poisson distribution parameters.
         *
         * @param mean Mean.
         * @throws IllegalArgumentException if {@code mean <= 0} or {@code mean > 1024}.
         */
        private static void validatePoissonDistributionParameters(double mean) {
            InternalUtils.requireStrictlyPositive(mean, "mean");
            if (mean > MAX_MEAN) {
                throw new IllegalArgumentException("mean " + mean + " > " + MAX_MEAN);
            }
        }

        /**
         * Creates the Poisson distribution by computing probabilities recursively from {@code X=0}.
         *
         * @param rng Generator of uniformly distributed random numbers.
         * @param mean Mean.
         * @return Sampler.
         */
        private static SharedStateDiscreteSampler createPoissonDistributionFromX0(
                UniformRandomProvider rng, double mean) {
            final double p0 = Math.exp(-mean);

            // Recursive update of Poisson probability until the value is too small
            // p(x + 1) = p(x) * mean / (x + 1)
            double p = p0;
            int i = 1;
            while (p * DOUBLE_31 >= 1) {
                p *= mean / i++;
            }

            // Probabilities are 30-bit integers, assumed denominator 2^30
            final int size = i - 1;
            final int[] prob = new int[size];

            p = p0;
            prob[0] = toUnsignedInt30(p);
            // The sum must exceed 2^30. In edges cases this is false due to round-off.
            int sum = prob[0];
            for (i = 1; i < prob.length; i++) {
                p *= mean / i;
                prob[i] = toUnsignedInt30(p);
                sum += prob[i];
            }

            // If the sum is < 2^30 add the remaining sum to the mode (floor(mean)).
            prob[(int) mean] += Math.max(0, INT_30 - sum);

            // Note: offset = 0
            return createSampler(rng, POISSON_NAME, prob, 0);
        }

        /**
         * Creates the Poisson distribution by computing probabilities recursively upward and downward
         * from {@code X=mode}, the location of the largest p-value.
         *
         * @param rng Generator of uniformly distributed random numbers.
         * @param mean Mean.
         * @return Sampler.
         */
        private static SharedStateDiscreteSampler createPoissonDistributionFromXMode(
                UniformRandomProvider rng, double mean) {
            // If mean >= 21.4, generate from largest p-value up, then largest down.
            // The largest p-value will be at the mode (floor(mean)).

            // Find p(x=mode)
            final int mode = (int) mean;
            // This transform is stable until mean >= 1941 where p will result in Infinity
            // before the divisor i is large enough to start reducing the product (i.e. i > c).
            final double c = mean * Math.exp(-mean / mode);
            double p = 1.0;
            for (int i = 1; i <= mode; i++) {
                p *= c / i;
            }
            final double pMode = p;

            // Find the upper limit using recursive computation of the p-value.
            // Note this will exit when i overflows to negative so no check on the range
            int i = mode + 1;
            while (p * DOUBLE_31 >= 1) {
                p *= mean / i++;
            }
            final int last = i - 2;

            // Find the lower limit using recursive computation of the p-value.
            p = pMode;
            int j = -1;
            for (i = mode - 1; i >= 0; i--) {
                p *= (i + 1) / mean;
                if (p * DOUBLE_31 < 1) {
                    j = i;
                    break;
                }
            }

            // Probabilities are 30-bit integers, assumed denominator 2^30.
            // This is the minimum sample value: prob[x - offset] = p(x)
            final int offset = j + 1;
            final int size = last - offset + 1;
            final int[] prob = new int[size];

            p = pMode;
            prob[mode - offset] = toUnsignedInt30(p);
            // The sum must exceed 2^30. In edges cases this is false due to round-off.
            int sum = prob[mode - offset];
            // From mode to upper limit
            for (i = mode + 1; i <= last; i++) {
                p *= mean / i;
                prob[i - offset] = toUnsignedInt30(p);
                sum += prob[i - offset];
            }
            // From mode to lower limit
            p = pMode;
            for (i = mode - 1; i >= offset; i--) {
                p *= (i + 1) / mean;
                prob[i - offset] = toUnsignedInt30(p);
                sum += prob[i - offset];
            }

            // If the sum is < 2^30 add the remaining sum to the mode.
            // If above 2^30 then the effect is truncation of the long tail of the distribution.
            prob[mode - offset] += Math.max(0, INT_30 - sum);

            return createSampler(rng, POISSON_NAME, prob, offset);
        }
    }

    /**
     * Create a sampler for the Binomial distribution.
     */
    public static final class Binomial {
        /** The name of the Binomial distribution. */
        private static final String BINOMIAL_NAME = "Binomial";

        /**
         * Return a fixed result for the Binomial distribution. This is a special class to handle
         * an edge case of probability of success equal to 0 or 1.
         */
        private static final class MarsagliaTsangWangFixedResultBinomialSampler
            extends AbstractMarsagliaTsangWangDiscreteSampler {
            /** The result. */
            private final int result;

            /**
             * @param result Result.
             */
            MarsagliaTsangWangFixedResultBinomialSampler(int result) {
                super(null, BINOMIAL_NAME);
                this.result = result;
            }

            @Override
            public int sample() {
                return result;
            }

            @Override
            public String toString() {
                return BINOMIAL_NAME + " deviate";
            }

            @Override
            public SharedStateDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) {
                // No shared state
                return this;
            }
        }

        /**
         * Return an inversion result for the Binomial distribution. This assumes the
         * following:
         *
         * <pre>
         * Binomial(n, p) = 1 - Binomial(n, 1 - p)
         * </pre>
         */
        private static final class MarsagliaTsangWangInversionBinomialSampler
            extends AbstractMarsagliaTsangWangDiscreteSampler {
            /** The number of trials. */
            private final int trials;
            /** The Binomial distribution sampler. */
            private final SharedStateDiscreteSampler sampler;

            /**
             * @param trials Number of trials.
             * @param sampler Binomial distribution sampler.
             */
            MarsagliaTsangWangInversionBinomialSampler(int trials,
                                                       SharedStateDiscreteSampler sampler) {
                super(null, BINOMIAL_NAME);
                this.trials = trials;
                this.sampler = sampler;
            }

            @Override
            public int sample() {
                return trials - sampler.sample();
            }

            @Override
            public String toString() {
                return sampler.toString();
            }

            @Override
            public SharedStateDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) {
                return new MarsagliaTsangWangInversionBinomialSampler(this.trials,
                    this.sampler.withUniformRandomProvider(rng));
            }
        }

        /** Class contains only static methods. */
        private Binomial() {}

        /**
         * Creates a sampler for the Binomial distribution.
         *
         * <p>Any probability less than 2<sup>-31</sup> will not be observed in samples.</p>
         *
         * <p>Storage requirements depend on the tabulated probability values. Example storage
         * requirements are listed below (in kB).</p>
         *
         * <pre>
         *          p
         * trials   0.5    0.1   0.01  0.001
         *    4    0.06   0.63   0.44   0.44
         *   16    0.69   1.14   0.76   0.44
         *   64    4.73   2.40   1.14   0.51
         *  256    8.63   5.17   1.89   0.82
         * 1024   31.12   9.45   3.34   0.89
         * </pre>
         *
         * <p>The method requires that the Binomial distribution probability at {@code x=0} can be computed.
         * This will fail when {@code (1 - p)^trials == 0} which requires {@code trials} to be large
         * and/or {@code p} to be small. In this case an exception is raised.</p>
         *
         * @param rng Generator of uniformly distributed random numbers.
         * @param trials Number of trials.
         * @param probabilityOfSuccess Probability of success (p).
         * @return Sampler.
         * @throws IllegalArgumentException if {@code trials < 0} or {@code trials >= 2^16},
         * {@code p} is not in the range {@code [0-1]}, or the probability distribution cannot
         * be computed.
         */
        public static SharedStateDiscreteSampler of(UniformRandomProvider rng,
                                                    int trials,
                                                    double probabilityOfSuccess) {
            validateBinomialDistributionParameters(trials, probabilityOfSuccess);

            // Handle edge cases
            if (probabilityOfSuccess == 0) {
                return new MarsagliaTsangWangFixedResultBinomialSampler(0);
            }
            if (probabilityOfSuccess == 1) {
                return new MarsagliaTsangWangFixedResultBinomialSampler(trials);
            }

            // Check the supported size.
            if (trials >= INT_16) {
                throw new IllegalArgumentException("Unsupported number of trials: " + trials);
            }

            return createBinomialDistributionSampler(rng, trials, probabilityOfSuccess);
        }

        /**
         * Validate the Binomial distribution parameters.
         *
         * @param trials Number of trials.
         * @param probabilityOfSuccess Probability of success (p).
         * @throws IllegalArgumentException if {@code trials < 0} or
         * {@code p} is not in the range {@code [0-1]}
         */
        private static void validateBinomialDistributionParameters(int trials, double probabilityOfSuccess) {
            if (trials < 0) {
                throw new IllegalArgumentException("Trials is not positive: " + trials);
            }
            InternalUtils.requireRangeClosed(0, 1, probabilityOfSuccess, "probability of success");
        }

        /**
         * Creates the Binomial distribution sampler.
         *
         * <p>This assumes the parameters for the distribution are valid. The method
         * will only fail if the initial probability for {@code X=0} is zero.</p>
         *
         * @param rng Generator of uniformly distributed random numbers.
         * @param trials Number of trials.
         * @param probabilityOfSuccess Probability of success (p).
         * @return Sampler.
         * @throws IllegalArgumentException if the probability distribution cannot be
         * computed.
         */
        private static SharedStateDiscreteSampler createBinomialDistributionSampler(
                UniformRandomProvider rng, int trials, double probabilityOfSuccess) {

            // The maximum supported value for Math.exp is approximately -744.
            // This occurs when trials is large and p is close to 1.
            // Handle this by using an inversion: generate j=Binomial(n,1-p), return n-j
            final boolean useInversion = probabilityOfSuccess > 0.5;
            final double p = useInversion ? 1 - probabilityOfSuccess : probabilityOfSuccess;

            // Check if the distribution can be computed
            final double p0 = Math.exp(trials * Math.log(1 - p));
            if (p0 < Double.MIN_VALUE) {
                throw new IllegalArgumentException("Unable to compute distribution");
            }

            // First find size of probability array
            double t = p0;
            final double h = p / (1 - p);
            // Find first probability above the threshold of 2^-31
            int begin = 0;
            if (t * DOUBLE_31 < 1) {
                // Somewhere after p(0)
                // Note:
                // If this loop is entered p(0) is < 2^-31.
                // This has been tested at the extreme for p(0)=Double.MIN_VALUE and either
                // p=0.5 or trials=2^16-1 and does not fail to find the beginning.
                for (int i = 1; i <= trials; i++) {
                    t *= (trials + 1 - i) * h / i;
                    if (t * DOUBLE_31 >= 1) {
                        begin = i;
                        break;
                    }
                }
            }
            // Find last probability
            int end = trials;
            for (int i = begin + 1; i <= trials; i++) {
                t *= (trials + 1 - i) * h / i;
                if (t * DOUBLE_31 < 1) {
                    end = i - 1;
                    break;
                }
            }

            return createBinomialDistributionSamplerFromRange(rng, trials, p, useInversion,
                    p0, begin, end);
        }

        /**
         * Creates the Binomial distribution sampler using only the probability values for {@code X}
         * between the begin and the end (inclusive).
         *
         * @param rng Generator of uniformly distributed random numbers.
         * @param trials Number of trials.
         * @param p Probability of success (p).
         * @param useInversion Set to {@code true} if the probability was inverted.
         * @param p0 Probability at {@code X=0}
         * @param begin Begin value {@code X} for the distribution.
         * @param end End value {@code X} for the distribution.
         * @return Sampler.
         */
        private static SharedStateDiscreteSampler createBinomialDistributionSamplerFromRange(
                UniformRandomProvider rng, int trials, double p,
                boolean useInversion, double p0, int begin, int end) {

            // Assign probability values as 30-bit integers
            final int size = end - begin + 1;
            final int[] prob = new int[size];
            double t = p0;
            final double h = p / (1 - p);
            for (int i = 1; i <= begin; i++) {
                t *= (trials + 1 - i) * h / i;
            }
            int sum = toUnsignedInt30(t);
            prob[0] = sum;
            for (int i = begin + 1; i <= end; i++) {
                t *= (trials + 1 - i) * h / i;
                prob[i - begin] = toUnsignedInt30(t);
                sum += prob[i - begin];
            }

            // If the sum is < 2^30 add the remaining sum to the mode (floor((n+1)p))).
            // If above 2^30 then the effect is truncation of the long tail of the distribution.
            final int mode = (int) ((trials + 1) * p) - begin;
            prob[mode] += Math.max(0, INT_30 - sum);

            final SharedStateDiscreteSampler sampler = createSampler(rng, BINOMIAL_NAME, prob, begin);

            // Check if an inversion was made
            return useInversion ?
                   new MarsagliaTsangWangInversionBinomialSampler(trials, sampler) :
                   sampler;
        }
    }
}