AliasMethodDiscreteSampler.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.sampling.distribution;

  18. import org.apache.commons.rng.UniformRandomProvider;

  19. import java.util.Arrays;

  20. /**
  21.  * Distribution sampler that uses the <a
  22.  * href="https://en.wikipedia.org/wiki/Alias_method">Alias method</a>. It can be used to
  23.  * sample from {@code n} values each with an associated probability. If all unique items
  24.  * are assigned the same probability it is more efficient to use the {@link DiscreteUniformSampler}.
  25.  *
  26.  * <p>This implementation is based on the detailed explanation of the alias method by
  27.  * Keith Schartz and implements Vose's algorithm.</p>
  28.  *
  29.  * <ul>
  30.  *  <li>
  31.  *   <blockquote>
  32.  *    Vose, M.D.,
  33.  *    <i>A linear algorithm for generating random numbers with a given distribution,</i>
  34.  *     IEEE Transactions on Software Engineering, 17, 972-975, 1991.
  35.  *   </blockquote>
  36.  *  </li>
  37.  * </ul>
  38.  *
  39.  * <p>The algorithm will sample values in {@code O(1)} time after a pre-processing step of
  40.  * {@code O(n)} time.</p>
  41.  *
  42.  * <p>The alias tables are constructed using fraction probabilities with an assumed denominator
  43.  * of 2<sup>53</sup>. In the generic case sampling uses {@link UniformRandomProvider#nextInt(int)}
  44.  * and the upper 53-bits from {@link UniformRandomProvider#nextLong()}.</p>
  45.  *
  46.  * <p>Zero padding the input probabilities can be used to make more sampling more efficient.
  47.  * Any zero entry will always be aliased removing the requirement to compute a {@code long}.
  48.  * Increased sampling speed comes at the cost of increased storage space. The algorithm requires
  49.  * approximately 12 bytes of storage per input probability, that is {@code n * 12} for size
  50.  * {@code n}. Zero-padding only requires 4 bytes of storage per padded value as the probability is
  51.  * known to be zero. A table can be padded to a power of 2 using the utility function
  52.  * {@link #of(UniformRandomProvider, double[], int)} to construct the sampler.</p>
  53.  *
  54.  * <p>An optimisation is performed for small table sizes that are a power of 2. In this case the
  55.  * sampling uses 1 or 2 calls from {@link UniformRandomProvider#nextInt()} to generate up to
  56.  * 64-bits for creation of an 11-bit index and 53-bits for the {@code long}. This optimisation
  57.  * requires a generator with a high cycle length for the lower order bits.</p>
  58.  *
  59.  * <p>Larger table sizes that are a power of 2 will benefit from fast algorithms for
  60.  * {@link UniformRandomProvider#nextInt(int)} that exploit the power of 2.</p>
  61.  *
  62.  * @see <a href="https://en.wikipedia.org/wiki/Alias_method">Alias Method</a>
  63.  * @see <a href="http://www.keithschwarz.com/darts-dice-coins/">Darts, Dice, and Coins:
  64.  * Sampling from a Discrete Distribution by Keith Schwartz</a>
  65.  * @see <a href="https://ieeexplore.ieee.org/document/92917">Vose (1991) IEEE Transactions
  66.  * on Software Engineering 17, 972-975.</a>
  67.  * @since 1.3
  68.  */
  69. public class AliasMethodDiscreteSampler
  70.     implements SharedStateDiscreteSampler {
  71.     /**
  72.      * The default alpha factor for zero-padding an input probability table. The default
  73.      * value will pad the probabilities by to the next power-of-2.
  74.      */
  75.     private static final int DEFAULT_ALPHA = 0;
  76.     /** The value zero for a {@code double}. */
  77.     private static final double ZERO = 0.0;
  78.     /** The value 1.0 represented as the numerator of a fraction with denominator 2<sup>53</sup>. */
  79.     private static final long ONE_AS_NUMERATOR = 1L << 53;
  80.     /**
  81.      * The multiplier to convert a {@code double} probability in the range {@code [0, 1]}
  82.      * to the numerator of a fraction with denominator 2<sup>53</sup>.
  83.      */
  84.     private static final double CONVERT_TO_NUMERATOR = ONE_AS_NUMERATOR;
  85.     /**
  86.      * The maximum size of the small alias table. This is 2<sup>11</sup>.
  87.      */
  88.     private static final int MAX_SMALL_POWER_2_SIZE = 1 << 11;

  89.     /** Underlying source of randomness. */
  90.     protected final UniformRandomProvider rng;

  91.     /**
  92.      * The probability table. During sampling a random index into this table is selected.
  93.      * A random probability is compared to the value at this index: if lower then the sample is the
  94.      * index; if higher then the sample uses the corresponding entry in the alias table.
  95.      *
  96.      * <p>This has entries up to the last non-zero element since there is no need to store
  97.      * probabilities of zero. This is an optimisation for zero-padded input. Any zero value will
  98.      * always be aliased so any look-up index outside this table always uses the alias.</p>
  99.      *
  100.      * <p>Note that a uniform double in the range [0,1) can be generated using 53-bits from a long
  101.      * to sample all the dyadic rationals with a denominator of 2<sup>53</sup>
  102.      * (e.g. see org.apache.commons.rng.core.utils.NumberFactory.makeDouble(long)). To avoid
  103.      * computation of a double and comparison to the probability as a double the probabilities are
  104.      * stored as 53-bit longs to use integer arithmetic. This is the equivalent of storing the
  105.      * numerator of a fraction with the denominator of 2<sup>53</sup>.</p>
  106.      *
  107.      * <p>During conversion of the probability to a double it is rounded up to the next integer
  108.      * value. This ensures the functionality of comparing a uniform deviate distributed evenly on
  109.      * the interval 1/2^53 to the unevenly distributed probability is equivalent, i.e. a uniform
  110.      * deviate is either below the probability or above it:
  111.      *
  112.      * <pre>
  113.      * Uniform deviate
  114.      *  1/2^53    2/2^53    3/2^53    4/2^53
  115.      * --|---------|---------|---------|---
  116.      *      ^
  117.      *      |
  118.      *  probability
  119.      *             ^
  120.      *             |
  121.      *         rounded up
  122.      * </pre>
  123.      *
  124.      * <p>Round-up ensures a non-zero probability is always non-zero and zero probability remains
  125.      * zero. Thus any item with a non-zero input probability can always be sampled, and a zero
  126.      * input probability cannot be sampled.</p>
  127.      *
  128.      * @see <a href="https://en.wikipedia.org/wiki/Dyadic_rational">Dyadic rational</a>
  129.      */
  130.     protected final long[] probability;

  131.     /**
  132.      * The alias table. During sampling if the random probability is not below the entry in the
  133.      * probability table then the sample is the alias.
  134.      */
  135.     protected final int[] alias;

  136.     /**
  137.      * Sample from the computed tables exploiting the small power-of-two table size.
  138.      * This implements a variant of the optimised algorithm as per Vose (1991):
  139.      *
  140.      * <pre>
  141.      * bits = obtained required number of random bits
  142.      * v = (some of the bits) * constant1
  143.      * j = (rest of the bits) * constant2
  144.      * if v &lt; prob[j] then
  145.      *   return j
  146.      * else
  147.      *   return alias[j]
  148.      * </pre>
  149.      *
  150.      * <p>This is a variant because the bits are not multiplied by constants. In the case of
  151.      * {@code v} the constant is a scale that is pre-applied to the probability table. In the
  152.      * case of {@code j} the constant is not used to scale a deviate to an index; the index is
  153.      * from a power-of-2 range and so the bits are used directly.</p>
  154.      *
  155.      * <p>This is implemented using up to 64 bits from the random generator.
  156.      * The index for the table is computed using a mask to extract up to 11 of the lower bits
  157.      * from an integer. The probability is computed using a second integer combined with the
  158.      * remaining bits to create 53-bits for the numerator of a fraction with denominator
  159.      * 2<sup>53</sup>. This is only computed on demand.</p>
  160.      *
  161.      * <p>Note: This supports a table size of up to 2^11, or 2048, exclusive. Any larger requires
  162.      * consuming more than 64-bits and the algorithm is not more efficient than the
  163.      * {@link AliasMethodDiscreteSampler}.</p>
  164.      *
  165.      * <p>Sampling uses 1 or 2 calls to {@link UniformRandomProvider#nextInt()}.</p>
  166.      */
  167.     private static final class SmallTableAliasMethodDiscreteSampler extends AliasMethodDiscreteSampler {
  168.         /** The mask to isolate the lower bits. */
  169.         private final int mask;

  170.         /**
  171.          * Create a new instance.
  172.          *
  173.          * @param rng Generator of uniformly distributed random numbers.
  174.          * @param probability Probability table.
  175.          * @param alias Alias table.
  176.          */
  177.         SmallTableAliasMethodDiscreteSampler(final UniformRandomProvider rng,
  178.                                              final long[] probability,
  179.                                              final int[] alias) {
  180.             super(rng, probability, alias);
  181.             // Assume the table size is a power of 2 and create the mask
  182.             mask = alias.length - 1;
  183.         }

  184.         @Override
  185.         public int sample() {
  186.             final int bits = rng.nextInt();
  187.             // Isolate lower bits
  188.             final int j = bits & mask;

  189.             // Optimisation for zero-padded input tables
  190.             if (j >= probability.length) {
  191.                 // No probability must use the alias
  192.                 return alias[j];
  193.             }

  194.             // Create a uniform random deviate as a long.
  195.             // This replicates functionality from the o.a.c.rng.core.utils.NumberFactory.makeLong
  196.             final long longBits = (((long) rng.nextInt()) << 32) | (bits & 0xffffffffL);

  197.             // Choose between the two. Use a 53-bit long for the probability.
  198.             return (longBits >>> 11) < probability[j] ? j : alias[j];
  199.         }

  200.         /** {@inheritDoc} */
  201.         @Override
  202.         public SharedStateDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) {
  203.             return new SmallTableAliasMethodDiscreteSampler(rng, probability, alias);
  204.         }
  205.     }

  206.     /**
  207.      * Creates a sampler.
  208.      *
  209.      * <p>The input parameters are not validated and must be correctly computed alias tables.</p>
  210.      *
  211.      * @param rng Generator of uniformly distributed random numbers.
  212.      * @param probability Probability table.
  213.      * @param alias Alias table.
  214.      */
  215.     AliasMethodDiscreteSampler(final UniformRandomProvider rng,
  216.                                final long[] probability,
  217.                                final int[] alias) {
  218.         this.rng = rng;
  219.         // Deliberate direct storage of input arrays
  220.         this.probability = probability;
  221.         this.alias = alias;
  222.     }

  223.     /** {@inheritDoc} */
  224.     @Override
  225.     public int sample() {
  226.         // This implements the algorithm as per Vose (1991):
  227.         // v = uniform()  in [0, 1)
  228.         // j = uniform(n) in [0, n)
  229.         // if v < prob[j] then
  230.         //   return j
  231.         // else
  232.         //   return alias[j]

  233.         final int j = rng.nextInt(alias.length);

  234.         // Optimisation for zero-padded input tables
  235.         if (j >= probability.length) {
  236.             // No probability must use the alias
  237.             return alias[j];
  238.         }

  239.         // Note: We could check the probability before computing a deviate.
  240.         // p(j) == 0  => alias[j]
  241.         // p(j) == 1  => j
  242.         // However it is assumed these edge cases are rare:
  243.         //
  244.         // The probability table will be 1 for approximately 1/n samples, i.e. only the
  245.         // last unpaired probability. This is only worth checking for when the table size (n)
  246.         // is small. But in that case the user should zero-pad the table for performance.
  247.         //
  248.         // The probability table will be 0 when an input probability was zero. We
  249.         // will assume this is also rare if modelling a discrete distribution where
  250.         // all samples are possible. The edge case for zero-padded tables is handled above.

  251.         // Choose between the two. Use a 53-bit long for the probability.
  252.         return (rng.nextLong() >>> 11) < probability[j] ? j : alias[j];
  253.     }

  254.     /** {@inheritDoc} */
  255.     @Override
  256.     public String toString() {
  257.         return "Alias method [" + rng.toString() + "]";
  258.     }

  259.     /** {@inheritDoc} */
  260.     @Override
  261.     public SharedStateDiscreteSampler withUniformRandomProvider(UniformRandomProvider rng) {
  262.         return new AliasMethodDiscreteSampler(rng, probability, alias);
  263.     }

  264.     /**
  265.      * Creates a sampler.
  266.      *
  267.      * <p>The probabilities will be normalised using their sum. The only requirement
  268.      * is the sum is strictly positive.</p>
  269.      *
  270.      * <p>Where possible this method zero-pads the probabilities so the length is the next
  271.      * power-of-two. Padding is bounded by the upper limit on the size of an array.</p>
  272.      *
  273.      * <p>To avoid zero-padding use the
  274.      * {@link #of(UniformRandomProvider, double[], int)} method with a negative
  275.      * {@code alpha} factor.</p>
  276.      *
  277.      * @param rng Generator of uniformly distributed random numbers.
  278.      * @param probabilities The list of probabilities.
  279.      * @return the sampler
  280.      * @throws IllegalArgumentException if {@code probabilities} is null or empty, a
  281.      * probability is negative, infinite or {@code NaN}, or the sum of all
  282.      * probabilities is not strictly positive.
  283.      * @see #of(UniformRandomProvider, double[], int)
  284.      */
  285.     public static SharedStateDiscreteSampler of(final UniformRandomProvider rng,
  286.                                                 final double[] probabilities) {
  287.         return of(rng, probabilities, DEFAULT_ALPHA);
  288.     }

  289.     /**
  290.      * Creates a sampler.
  291.      *
  292.      * <p>The probabilities will be normalised using their sum. The only requirement
  293.      * is the sum is strictly positive.</p>
  294.      *
  295.      * <p>Where possible this method zero-pads the probabilities to improve sampling
  296.      * efficiency. Padding is bounded by the upper limit on the size of an array and
  297.      * controlled by the {@code alpha} argument. Set to negative to disable
  298.      * padding.</p>
  299.      *
  300.      * <p>For each zero padded value an entry is added to the tables which is always
  301.      * aliased. This can be sampled with fewer bits required from the
  302.      * {@link UniformRandomProvider}. Increasing the padding of zeros increases the
  303.      * chance of using this fast path to selecting a sample. The penalty is
  304.      * two-fold: initialisation is bounded by {@code O(n)} time with {@code n} the
  305.      * size <strong>after</strong> padding; an additional memory cost of 4 bytes per
  306.      * padded value.</p>
  307.      *
  308.      * <p>Zero padding to any length improves performance; using a power of 2 allows
  309.      * the index into the tables to be more efficiently generated. The argument
  310.      * {@code alpha} controls the level of padding. Positive values of {@code alpha}
  311.      * represent a scale factor in powers of 2. The size of the input array will be
  312.      * increased by a factor of 2<sup>alpha</sup> and then rounded-up to the next
  313.      * power of 2. Padding is bounded by the upper limit on the size of an
  314.      * array.</p>
  315.      *
  316.      * <p>The chance of executing the slow path is upper bounded at
  317.      * 2<sup>-alpha</sup> when padding is enabled. Each successive doubling of
  318.      * padding will have diminishing performance gains.</p>
  319.      *
  320.      * @param rng Generator of uniformly distributed random numbers.
  321.      * @param probabilities The list of probabilities.
  322.      * @param alpha The alpha factor controlling the zero padding.
  323.      * @return the sampler
  324.      * @throws IllegalArgumentException if {@code probabilities} is null or empty, a
  325.      * probability is negative, infinite or {@code NaN}, or the sum of all
  326.      * probabilities is not strictly positive.
  327.      */
  328.     public static SharedStateDiscreteSampler of(final UniformRandomProvider rng,
  329.                                                 final double[] probabilities,
  330.                                                 int alpha) {
  331.         // The Alias method balances N categories with counts around the mean into N sections,
  332.         // each allocated 'mean' observations.
  333.         //
  334.         // Consider 4 categories with counts 6,3,2,1. The histogram can be balanced into a
  335.         // 2D array as 4 sections with a height of the mean:
  336.         //
  337.         // 6
  338.         // 6
  339.         // 6
  340.         // 63   => 6366   --
  341.         // 632     6326    |-- mean
  342.         // 6321    6321   --
  343.         //
  344.         // section abcd
  345.         //
  346.         // Each section is divided as:
  347.         // a: 6=1/1
  348.         // b: 3=1/1
  349.         // c: 2=2/3; 6=1/3   (6 is the alias)
  350.         // d: 1=1/3; 6=2/3   (6 is the alias)
  351.         //
  352.         // The sample is obtained by randomly selecting a section, then choosing which category
  353.         // from the pair based on a uniform random deviate.

  354.         final double sumProb = InternalUtils.validateProbabilities(probabilities);

  355.         // Allow zero-padding
  356.         final int n = computeSize(probabilities.length, alpha);

  357.         // Partition into small and large by splitting on the average.
  358.         final double mean = sumProb / n;
  359.         // The cardinality of smallSize + largeSize = n.
  360.         // So fill the same array from either end.
  361.         final int[] indices = new int[n];
  362.         int large = n;
  363.         int small = 0;
  364.         for (int i = 0; i < probabilities.length; i++) {
  365.             if (probabilities[i] >= mean) {
  366.                 indices[--large] = i;
  367.             } else {
  368.                 indices[small++] = i;
  369.             }
  370.         }

  371.         small = fillRemainingIndices(probabilities.length, indices, small);

  372.         // This may be smaller than the input length if the probabilities were already padded.
  373.         final int nonZeroIndex = findLastNonZeroIndex(probabilities);

  374.         // The probabilities are modified so use a copy.
  375.         // Note: probabilities are required only up to last nonZeroIndex
  376.         final double[] remainingProbabilities = Arrays.copyOf(probabilities, nonZeroIndex + 1);

  377.         // Allocate the final tables.
  378.         // Probability table may be truncated (when zero padded).
  379.         // The alias table is full length.
  380.         final long[] probability = new long[remainingProbabilities.length];
  381.         final int[] alias = new int[n];

  382.         // This loop uses each large in turn to fill the alias table for small probabilities that
  383.         // do not reach the requirement to fill an entire section alone (i.e. p < mean).
  384.         // Since the sum of the small should be less than the sum of the large it should use up
  385.         // all the small first. However floating point round-off can result in
  386.         // misclassification of items as small or large. The Vose algorithm handles this using
  387.         // a while loop conditioned on the size of both sets and a subsequent loop to use
  388.         // unpaired items.
  389.         while (large != n && small != 0) {
  390.             // Index of the small and the large probabilities.
  391.             final int j = indices[--small];
  392.             final int k = indices[large++];

  393.             // Optimisation for zero-padded input:
  394.             // p(j) = 0 above the last nonZeroIndex
  395.             if (j > nonZeroIndex) {
  396.                 // The entire amount for the section is taken from the alias.
  397.                 remainingProbabilities[k] -= mean;
  398.             } else {
  399.                 final double pj = remainingProbabilities[j];

  400.                 // Item j is a small probability that is below the mean.
  401.                 // Compute the weight of the section for item j: pj / mean.
  402.                 // This is scaled by 2^53 and the ceiling function used to round-up
  403.                 // the probability to a numerator of a fraction in the range [1,2^53].
  404.                 // Ceiling ensures non-zero values.
  405.                 probability[j] = (long) Math.ceil(CONVERT_TO_NUMERATOR * (pj / mean));

  406.                 // The remaining amount for the section is taken from the alias.
  407.                 // Effectively: probabilities[k] -= (mean - pj)
  408.                 remainingProbabilities[k] += pj - mean;
  409.             }

  410.             // If not j then the alias is k
  411.             alias[j] = k;

  412.             // Add the remaining probability from large to the appropriate list.
  413.             if (remainingProbabilities[k] >= mean) {
  414.                 indices[--large] = k;
  415.             } else {
  416.                 indices[small++] = k;
  417.             }
  418.         }

  419.         // Final loop conditions to consume unpaired items.
  420.         // Note: The large set should never be non-empty but this can occur due to round-off
  421.         // error so consume from both.
  422.         fillTable(probability, alias, indices, 0, small);
  423.         fillTable(probability, alias, indices, large, n);

  424.         // Change the algorithm for small power of 2 sized tables
  425.         return isSmallPowerOf2(n) ?
  426.             new SmallTableAliasMethodDiscreteSampler(rng, probability, alias) :
  427.             new AliasMethodDiscreteSampler(rng, probability, alias);
  428.     }

  429.     /**
  430.      * Allocate the remaining indices from zero padding as small probabilities. The
  431.      * number to add is from the length of the probability array to the length of
  432.      * the padded probability array (which is the same length as the indices array).
  433.      *
  434.      * @param length Length of probability array.
  435.      * @param indices Indices.
  436.      * @param small Number of small indices.
  437.      * @return the updated number of small indices
  438.      */
  439.     private static int fillRemainingIndices(final int length, final int[] indices, int small) {
  440.         int updatedSmall = small;
  441.         for (int i = length; i < indices.length; i++) {
  442.             indices[updatedSmall++] = i;
  443.         }
  444.         return updatedSmall;
  445.     }

  446.     /**
  447.      * Find the last non-zero index in the probabilities. This may be smaller than
  448.      * the input length if the probabilities were already padded.
  449.      *
  450.      * @param probabilities The list of probabilities.
  451.      * @return the index
  452.      */
  453.     private static int findLastNonZeroIndex(final double[] probabilities) {
  454.         // No bounds check is performed when decrementing as the array contains at least one
  455.         // value above zero.
  456.         int nonZeroIndex = probabilities.length - 1;
  457.         while (probabilities[nonZeroIndex] == ZERO) {
  458.             nonZeroIndex--;
  459.         }
  460.         return nonZeroIndex;
  461.     }

  462.     /**
  463.      * Compute the size after padding. A value of {@code alpha < 0} disables
  464.      * padding. Otherwise the length will be increased by 2<sup>alpha</sup>
  465.      * rounded-up to the next power of 2.
  466.      *
  467.      * @param length Length of probability array.
  468.      * @param alpha The alpha factor controlling the zero padding.
  469.      * @return the padded size
  470.      */
  471.     private static int computeSize(int length, int alpha) {
  472.         if (alpha < 0) {
  473.             // No padding
  474.             return length;
  475.         }
  476.         // Use the number of leading zeros function to find the next power of 2,
  477.         // i.e. ceil(log2(x))
  478.         int pow2 = 32 - Integer.numberOfLeadingZeros(length - 1);
  479.         // Increase by the alpha. Clip this to limit to a positive integer (2^30)
  480.         pow2 = Math.min(30, pow2 + alpha);
  481.         // Use max to handle a length above the highest possible power of 2
  482.         return Math.max(length, 1 << pow2);
  483.     }

  484.     /**
  485.      * Fill the tables using unpaired items that are in the range between {@code start} inclusive
  486.      * and {@code end} exclusive.
  487.      *
  488.      * <p>Anything left must fill the entire section so the probability table is set
  489.      * to 1 and there is no alias. This will occur for 1/n samples, i.e. the last
  490.      * remaining unpaired probability. Note: When the tables are zero-padded the
  491.      * remaining indices are from an input probability that is above zero so the
  492.      * index will be allowed in the truncated probability array and no
  493.      * index-out-of-bounds exception will occur.
  494.      *
  495.      * @param probability Probability table.
  496.      * @param alias Alias table.
  497.      * @param indices Unpaired indices.
  498.      * @param start Start position.
  499.      * @param end End position.
  500.      */
  501.     private static void fillTable(long[] probability, int[] alias, int[] indices, int start, int end) {
  502.         for (int i = start; i < end; i++) {
  503.             final int index = indices[i];
  504.             probability[index] = ONE_AS_NUMERATOR;
  505.             alias[index] = index;
  506.         }
  507.     }

  508.     /**
  509.      * Checks if the size is a small power of 2 so can be supported by the
  510.      * {@link SmallTableAliasMethodDiscreteSampler}.
  511.      *
  512.      * @param n Size of the alias table.
  513.      * @return true if supported by {@link SmallTableAliasMethodDiscreteSampler}
  514.      */
  515.     private static boolean isSmallPowerOf2(int n) {
  516.         return n <= MAX_SMALL_POWER_2_SIZE && (n & (n - 1)) == 0;
  517.     }
  518. }