    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     */
    017
    018    package org.apache.commons.math.random;
    019
    020    import java.io.Serializable;
    021    import java.security.MessageDigest;
    022    import java.security.NoSuchAlgorithmException;
    023    import java.security.NoSuchProviderException;
    024    import java.security.SecureRandom;
    025    import java.util.Collection;
    026
    027    import org.apache.commons.math.MathException;
    028    import org.apache.commons.math.distribution.BetaDistributionImpl;
    029    import org.apache.commons.math.distribution.BinomialDistributionImpl;
    030    import org.apache.commons.math.distribution.CauchyDistributionImpl;
    031    import org.apache.commons.math.distribution.ChiSquaredDistributionImpl;
    032    import org.apache.commons.math.distribution.ContinuousDistribution;
    033    import org.apache.commons.math.distribution.FDistributionImpl;
    034    import org.apache.commons.math.distribution.GammaDistributionImpl;
    035    import org.apache.commons.math.distribution.HypergeometricDistributionImpl;
    036    import org.apache.commons.math.distribution.IntegerDistribution;
    037    import org.apache.commons.math.distribution.PascalDistributionImpl;
    038    import org.apache.commons.math.distribution.TDistributionImpl;
    039    import org.apache.commons.math.distribution.WeibullDistributionImpl;
    040    import org.apache.commons.math.distribution.ZipfDistributionImpl;
    041    import org.apache.commons.math.exception.MathInternalError;
    042    import org.apache.commons.math.exception.NotStrictlyPositiveException;
    043    import org.apache.commons.math.exception.NumberIsTooLargeException;
    044    import org.apache.commons.math.exception.util.LocalizedFormats;
    045    import org.apache.commons.math.util.FastMath;
    046    import org.apache.commons.math.util.MathUtils;
    047
    048    /**
    049     * Implements the {@link RandomData} interface using a {@link RandomGenerator}
    050     * instance to generate non-secure data and a {@link java.security.SecureRandom}
    051     * instance to provide data for the <code>nextSecureXxx</code> methods. If no
    052     * <code>RandomGenerator</code> is provided in the constructor, the default is
    053     * to use a generator based on {@link java.util.Random}. To plug in a different
    054     * implementation, either implement <code>RandomGenerator</code> directly or
    055     * extend {@link AbstractRandomGenerator}.
    056     * <p>
    057     * Supports reseeding the underlying pseudo-random number generator (PRNG). The
    058     * <code>SecurityProvider</code> and <code>Algorithm</code> used by the
    059     * <code>SecureRandom</code> instance can also be reset.
    060     * </p>
    061     * <p>
    062     * For details on the default PRNGs, see {@link java.util.Random} and
    063     * {@link java.security.SecureRandom}.
    064     * </p>
    065     * <p>
    066     * <strong>Usage Notes</strong>:
    067     * <ul>
    068     * <li>
    069     * Instance variables are used to maintain <code>RandomGenerator</code> and
    070     * <code>SecureRandom</code> instances used in data generation. Therefore, to
    071     * generate a random sequence of values or strings, you should use just
    072     * <strong>one</strong> <code>RandomDataImpl</code> instance repeatedly.</li>
    073     * <li>
    074     * The "secure" methods are *much* slower. These should be used only when a
    075     * cryptographically secure random sequence is required. A secure random
    076     * sequence is a sequence of pseudo-random values which, in addition to being
    077     * well-dispersed (so no subsequence of values is an any more likely than other
    078     * subsequence of the the same length), also has the additional property that
    079     * knowledge of values generated up to any point in the sequence does not make
    080     * it any easier to predict subsequent values.</li>
    081     * <li>
    082     * When a new <code>RandomDataImpl</code> is created, the underlying random
    083     * number generators are <strong>not</strong> initialized. If you do not
    084     * explicitly seed the default non-secure generator, it is seeded with the
    085     * current time in milliseconds on first use. The same holds for the secure
    086     * generator. If you provide a <code>RandomGenerator</code> to the constructor,
    087     * however, this generator is not reseeded by the constructor nor is it reseeded
    088     * on first use.</li>
    089     * <li>
    090     * The <code>reSeed</code> and <code>reSeedSecure</code> methods delegate to the
    091     * corresponding methods on the underlying <code>RandomGenerator</code> and
    092     * <code>SecureRandom</code> instances. Therefore, <code>reSeed(long)</code>
    093     * fully resets the initial state of the non-secure random number generator (so
    094     * that reseeding with a specific value always results in the same subsequent
    095     * random sequence); whereas reSeedSecure(long) does <strong>not</strong>
    096     * reinitialize the secure random number generator (so secure sequences started
    097     * with calls to reseedSecure(long) won't be identical).</li>
    098     * <li>
    099     * This implementation is not synchronized.
    100     * </ul>
    101     * </p>
    102     *
    103     * @version $Revision: 1061496 $ $Date: 2011-01-20 21:32:16 +0100 (jeu. 20 janv. 2011) $
    104     */
    105    public class RandomDataImpl implements RandomData, Serializable {
    106
    107        /** Serializable version identifier */
    108        private static final long serialVersionUID = -626730818244969716L;
    109
    110        /** underlying random number generator */
    111        private RandomGenerator rand = null;
    112
    113        /** underlying secure random number generator */
    114        private SecureRandom secRand = null;
    115
    116        /**
    117         * Construct a RandomDataImpl.
    118         */
    119        public RandomDataImpl() {
    120        }
    121
    122        /**
    123         * Construct a RandomDataImpl using the supplied {@link RandomGenerator} as
    124         * the source of (non-secure) random data.
    125         *
    126         * @param rand
    127         *            the source of (non-secure) random data
    128         * @since 1.1
    129         */
    130        public RandomDataImpl(RandomGenerator rand) {
    131            super();
    132            this.rand = rand;
    133        }
    134
    135        /**
    136         * {@inheritDoc}
    137         * <p>
    138         * <strong>Algorithm Description:</strong> hex strings are generated using a
    139         * 2-step process.
    140         * <ol>
    141         * <li>
    142         * len/2+1 binary bytes are generated using the underlying Random</li>
    143         * <li>
    144         * Each binary byte is translated into 2 hex digits</li>
    145         * </ol>
    146         * </p>
    147         *
    148         * @param len
    149         *            the desired string length.
    150         * @return the random string.
    151         * @throws NotStrictlyPositiveException if {@code len <= 0}.
    152         */
    153        public String nextHexString(int len) {
    154            if (len <= 0) {
    155                throw new NotStrictlyPositiveException(LocalizedFormats.LENGTH, len);
    156            }
    157
    158            // Get a random number generator
    159            RandomGenerator ran = getRan();
    160
    161            // Initialize output buffer
    162            StringBuilder outBuffer = new StringBuilder();
    163
    164            // Get int(len/2)+1 random bytes
    165            byte[] randomBytes = new byte[(len / 2) + 1];
    166            ran.nextBytes(randomBytes);
    167
    168            // Convert each byte to 2 hex digits
    169            for (int i = 0; i < randomBytes.length; i++) {
    170                Integer c = Integer.valueOf(randomBytes[i]);
    171
    172                /*
    173                 * Add 128 to byte value to make interval 0-255 before doing hex
    174                 * conversion. This guarantees <= 2 hex digits from toHexString()
    175                 * toHexString would otherwise add 2^32 to negative arguments.
    176                 */
    177                String hex = Integer.toHexString(c.intValue() + 128);
    178
    179                // Make sure we add 2 hex digits for each byte
    180                if (hex.length() == 1) {
    181                    hex = "0" + hex;
    182                }
    183                outBuffer.append(hex);
    184            }
    185            return outBuffer.toString().substring(0, len);
    186        }
    187
    188        /**
    189         * Generate a random int value uniformly distributed between
    190         * <code>lower</code> and <code>upper</code>, inclusive.
    191         *
    192         * @param lower
    193         *            the lower bound.
    194         * @param upper
    195         *            the upper bound.
    196         * @return the random integer.
    197         * @throws NumberIsTooLargeException if {@code lower >= upper}.
    198         */
    199        public int nextInt(int lower, int upper) {
    200            if (lower >= upper) {
    201                throw new NumberIsTooLargeException(LocalizedFormats.LOWER_BOUND_NOT_BELOW_UPPER_BOUND,
    202                                                    lower, upper, false);
    203            }
    204            double r = getRan().nextDouble();
    205            return (int) ((r * upper) + ((1.0 - r) * lower) + r);
    206        }
    207
    208        /**
    209         * Generate a random long value uniformly distributed between
    210         * <code>lower</code> and <code>upper</code>, inclusive.
    211         *
    212         * @param lower
    213         *            the lower bound.
    214         * @param upper
    215         *            the upper bound.
    216         * @return the random integer.
    217         * @throws NumberIsTooLargeException if {@code lower >= upper}.
    218         */
    219        public long nextLong(long lower, long upper) {
    220            if (lower >= upper) {
    221                throw new NumberIsTooLargeException(LocalizedFormats.LOWER_BOUND_NOT_BELOW_UPPER_BOUND,
    222                                                    lower, upper, false);
    223            }
    224            double r = getRan().nextDouble();
    225            return (long) ((r * upper) + ((1.0 - r) * lower) + r);
    226        }
    227
    228        /**
    229         * {@inheritDoc}
    230         * <p>
    231         * <strong>Algorithm Description:</strong> hex strings are generated in
    232         * 40-byte segments using a 3-step process.
    233         * <ol>
    234         * <li>
    235         * 20 random bytes are generated using the underlying
    236         * <code>SecureRandom</code>.</li>
    237         * <li>
    238         * SHA-1 hash is applied to yield a 20-byte binary digest.</li>
    239         * <li>
    240         * Each byte of the binary digest is converted to 2 hex digits.</li>
    241         * </ol>
    242         * </p>
    243         *
    244         * @param len
    245         *            the length of the generated string
    246         * @return the random string
    247         * @throws NotStrictlyPositiveException if {@code len <= 0}.
    248         */
    249        public String nextSecureHexString(int len) {
    250            if (len <= 0) {
    251                throw new NotStrictlyPositiveException(LocalizedFormats.LENGTH, len);
    252            }
    253
    254            // Get SecureRandom and setup Digest provider
    255            SecureRandom secRan = getSecRan();
    256            MessageDigest alg = null;
    257            try {
    258                alg = MessageDigest.getInstance("SHA-1");
    259            } catch (NoSuchAlgorithmException ex) {
    260                // this should never happen
    261                throw new MathInternalError(ex);
    262            }
    263            alg.reset();
    264
    265            // Compute number of iterations required (40 bytes each)
    266            int numIter = (len / 40) + 1;
    267
    268            StringBuilder outBuffer = new StringBuilder();
    269            for (int iter = 1; iter < numIter + 1; iter++) {
    270                byte[] randomBytes = new byte[40];
    271                secRan.nextBytes(randomBytes);
    272                alg.update(randomBytes);
    273
    274                // Compute hash -- will create 20-byte binary hash
    275                byte hash[] = alg.digest();
    276
    277                // Loop over the hash, converting each byte to 2 hex digits
    278                for (int i = 0; i < hash.length; i++) {
    279                    Integer c = Integer.valueOf(hash[i]);
    280
    281                    /*
    282                     * Add 128 to byte value to make interval 0-255 This guarantees
    283                     * <= 2 hex digits from toHexString() toHexString would
    284                     * otherwise add 2^32 to negative arguments
    285                     */
    286                    String hex = Integer.toHexString(c.intValue() + 128);
    287
    288                    // Keep strings uniform length -- guarantees 40 bytes
    289                    if (hex.length() == 1) {
    290                        hex = "0" + hex;
    291                    }
    292                    outBuffer.append(hex);
    293                }
    294            }
    295            return outBuffer.toString().substring(0, len);
    296        }
    297
    298        /**
    299         * Generate a random int value uniformly distributed between
    300         * <code>lower</code> and <code>upper</code>, inclusive. This algorithm uses
    301         * a secure random number generator.
    302         *
    303         * @param lower
    304         *            the lower bound.
    305         * @param upper
    306         *            the upper bound.
    307         * @return the random integer.
    308         * @throws NumberIsTooLargeException if {@code lower >= upper}.
    309         */
    310        public int nextSecureInt(int lower, int upper) {
    311            if (lower >= upper) {
    312                throw new NumberIsTooLargeException(LocalizedFormats.LOWER_BOUND_NOT_BELOW_UPPER_BOUND,
    313                                                    lower, upper, false);
    314            }
    315            SecureRandom sec = getSecRan();
    316            return lower + (int) (sec.nextDouble() * (upper - lower + 1));
    317        }
    318
    319        /**
    320         * Generate a random long value uniformly distributed between
    321         * <code>lower</code> and <code>upper</code>, inclusive. This algorithm uses
    322         * a secure random number generator.
    323         *
    324         * @param lower
    325         *            the lower bound.
    326         * @param upper
    327         *            the upper bound.
    328         * @return the random integer.
    329         * @throws NumberIsTooLargeException if {@code lower >= upper}.
    330         */
    331        public long nextSecureLong(long lower, long upper) {
    332            if (lower >= upper) {
    333                throw new NumberIsTooLargeException(LocalizedFormats.LOWER_BOUND_NOT_BELOW_UPPER_BOUND,
    334                                                    lower, upper, false);
    335            }
    336            SecureRandom sec = getSecRan();
    337            return lower + (long) (sec.nextDouble() * (upper - lower + 1));
    338        }
    339
    340        /**
    341         * {@inheritDoc}
    342         * <p>
    343         * <strong>Algorithm Description</strong>:
    344         * <ul><li> For small means, uses simulation of a Poisson process
    345         * using Uniform deviates, as described
    346         * <a href="http://irmi.epfl.ch/cmos/Pmmi/interactive/rng7.htm"> here.</a>
    347         * The Poisson process (and hence value returned) is bounded by 1000 * mean.</li>
    348         *
    349         * <li> For large means, uses the rejection algorithm described in <br/>
    350         * Devroye, Luc. (1981).<i>The Computer Generation of Poisson Random Variables</i>
    351         * <strong>Computing</strong> vol. 26 pp. 197-207.</li></ul></p>
    352         *
    353         * @param mean mean of the Poisson distribution.
    354         * @return the random Poisson value.
    355         * @throws NotStrictlyPositiveException if {@code mean <= 0}.
    356         */
    357        public long nextPoisson(double mean) {
    358            if (mean <= 0) {
    359                throw new NotStrictlyPositiveException(LocalizedFormats.MEAN, mean);
    360            }
    361
    362            final RandomGenerator generator = getRan();
    363
    364            final double pivot = 40.0d;
    365            if (mean < pivot) {
    366                double p = FastMath.exp(-mean);
    367                long n = 0;
    368                double r = 1.0d;
    369                double rnd = 1.0d;
    370
    371                while (n < 1000 * mean) {
    372                    rnd = generator.nextDouble();
    373                    r = r * rnd;
    374                    if (r >= p) {
    375                        n++;
    376                    } else {
    377                        return n;
    378                    }
    379                }
    380                return n;
    381            } else {
    382                final double lambda = FastMath.floor(mean);
    383                final double lambdaFractional = mean - lambda;
    384                final double logLambda = FastMath.log(lambda);
    385                final double logLambdaFactorial = MathUtils.factorialLog((int) lambda);
    386                final long y2 = lambdaFractional < Double.MIN_VALUE ? 0 : nextPoisson(lambdaFractional);
    387                final double delta = FastMath.sqrt(lambda * FastMath.log(32 * lambda / FastMath.PI + 1));
    388                final double halfDelta = delta / 2;
    389                final double twolpd = 2 * lambda + delta;
    390                final double a1 = FastMath.sqrt(FastMath.PI * twolpd) * FastMath.exp(1 / 8 * lambda);
    391                final double a2 = (twolpd / delta) * FastMath.exp(-delta * (1 + delta) / twolpd);
    392                final double aSum = a1 + a2 + 1;
    393                final double p1 = a1 / aSum;
    394                final double p2 = a2 / aSum;
    395                final double c1 = 1 / (8 * lambda);
    396
    397                double x = 0;
    398                double y = 0;
    399                double v = 0;
    400                int a = 0;
    401                double t = 0;
    402                double qr = 0;
    403                double qa = 0;
    404                for (;;) {
    405                    final double u = nextUniform(0.0, 1);
    406                    if (u <= p1) {
    407                        final double n = nextGaussian(0d, 1d);
    408                        x = n * FastMath.sqrt(lambda + halfDelta) - 0.5d;
    409                        if (x > delta || x < -lambda) {
    410                            continue;
    411                        }
    412                        y = x < 0 ? FastMath.floor(x) : FastMath.ceil(x);
    413                        final double e = nextExponential(1d);
    414                        v = -e - (n * n / 2) + c1;
    415                    } else {
    416                        if (u > p1 + p2) {
    417                            y = lambda;
    418                            break;
    419                        } else {
    420                            x = delta + (twolpd / delta) * nextExponential(1d);
    421                            y = FastMath.ceil(x);
    422                            v = -nextExponential(1d) - delta * (x + 1) / twolpd;
    423                        }
    424                    }
    425                    a = x < 0 ? 1 : 0;
    426                    t = y * (y + 1) / (2 * lambda);
    427                    if (v < -t && a == 0) {
    428                        y = lambda + y;
    429                        break;
    430                    }
    431                    qr = t * ((2 * y + 1) / (6 * lambda) - 1);
    432                    qa = qr - (t * t) / (3 * (lambda + a * (y + 1)));
    433                    if (v < qa) {
    434                        y = lambda + y;
    435                        break;
    436                    }
    437                    if (v > qr) {
    438                        continue;
    439                    }
    440                    if (v < y * logLambda - MathUtils.factorialLog((int) (y + lambda)) + logLambdaFactorial) {
    441                        y = lambda + y;
    442                        break;
    443                    }
    444                }
    445                return y2 + (long) y;
    446            }
    447        }
    448
    449        /**
    450         * Generate a random value from a Normal (a.k.a. Gaussian) distribution with
    451         * the given mean, <code>mu</code> and the given standard deviation,
    452         * <code>sigma</code>.
    453         *
    454         * @param mu
    455         *            the mean of the distribution
    456         * @param sigma
    457         *            the standard deviation of the distribution
    458         * @return the random Normal value
    459         * @throws NotStrictlyPositiveException if {@code sigma <= 0}.
    460         */
    461        public double nextGaussian(double mu, double sigma) {
    462            if (sigma <= 0) {
    463                throw new NotStrictlyPositiveException(LocalizedFormats.STANDARD_DEVIATION, sigma);
    464            }
    465            return sigma * getRan().nextGaussian() + mu;
    466        }
    467
    468        /**
    469         * Returns a random value from an Exponential distribution with the given
    470         * mean.
    471         * <p>
    472         * <strong>Algorithm Description</strong>: Uses the <a
    473         * href="http://www.jesus.ox.ac.uk/~clifford/a5/chap1/node5.html"> Inversion
    474         * Method</a> to generate exponentially distributed random values from
    475         * uniform deviates.
    476         * </p>
    477         *
    478         * @param mean the mean of the distribution
    479         * @return the random Exponential value
    480         * @throws NotStrictlyPositiveException if {@code mean <= 0}.
    481         */
    482        public double nextExponential(double mean) {
    483            if (mean <= 0.0) {
    484                throw new NotStrictlyPositiveException(LocalizedFormats.MEAN, mean);
    485            }
    486            final RandomGenerator generator = getRan();
    487            double unif = generator.nextDouble();
    488            while (unif == 0.0d) {
    489                unif = generator.nextDouble();
    490            }
    491            return -mean * FastMath.log(unif);
    492        }
    493
    494        /**
    495         * {@inheritDoc}
    496         * <p>
    497         * <strong>Algorithm Description</strong>: scales the output of
    498         * Random.nextDouble(), but rejects 0 values (i.e., will generate another
    499         * random double if Random.nextDouble() returns 0). This is necessary to
    500         * provide a symmetric output interval (both endpoints excluded).
    501         * </p>
    502         *
    503         * @param lower
    504         *            the lower bound.
    505         * @param upper
    506         *            the upper bound.
    507         * @return a uniformly distributed random value from the interval (lower,
    508         *         upper)
    509         * @throws NumberIsTooLargeException if {@code lower >= upper}.
    510         */
    511        public double nextUniform(double lower, double upper) {
    512            if (lower >= upper) {
    513                throw new NumberIsTooLargeException(LocalizedFormats.LOWER_BOUND_NOT_BELOW_UPPER_BOUND,
    514                                                    lower, upper, false);
    515            }
    516            final RandomGenerator generator = getRan();
    517
    518            // ensure nextDouble() isn't 0.0
    519            double u = generator.nextDouble();
    520            while (u <= 0.0) {
    521                u = generator.nextDouble();
    522            }
    523
    524            return lower + u * (upper - lower);
    525        }
    526
    527        /**
    528         * Generates a random value from the {@link BetaDistributionImpl Beta Distribution}.
    529         * This implementation uses {@link #nextInversionDeviate(ContinuousDistribution) inversion}
    530         * to generate random values.
    531         *
    532         * @param alpha first distribution shape parameter
    533         * @param beta second distribution shape parameter
    534         * @return random value sampled from the beta(alpha, beta) distribution
    535         * @throws MathException if an error occurs generating the random value
    536         * @since 2.2
    537         */
    538        public double nextBeta(double alpha, double beta) throws MathException {
    539            return nextInversionDeviate(new BetaDistributionImpl(alpha, beta));
    540        }
    541
    542        /**
    543         * Generates a random value from the {@link BinomialDistributionImpl Binomial Distribution}.
    544         * This implementation uses {@link #nextInversionDeviate(ContinuousDistribution) inversion}
    545         * to generate random values.
    546         *
    547         * @param numberOfTrials number of trials of the Binomial distribution
    548         * @param probabilityOfSuccess probability of success of the Binomial distribution
    549         * @return random value sampled from the Binomial(numberOfTrials, probabilityOfSuccess) distribution
    550         * @throws MathException if an error occurs generating the random value
    551         * @since 2.2
    552         */
    553        public int nextBinomial(int numberOfTrials, double probabilityOfSuccess) throws MathException {
    554            return nextInversionDeviate(new BinomialDistributionImpl(numberOfTrials, probabilityOfSuccess));
    555        }
    556
    557        /**
    558         * Generates a random value from the {@link CauchyDistributionImpl Cauchy Distribution}.
    559         * This implementation uses {@link #nextInversionDeviate(ContinuousDistribution) inversion}
    560         * to generate random values.
    561         *
    562         * @param median the median of the Cauchy distribution
    563         * @param scale the scale parameter of the Cauchy distribution
    564         * @return random value sampled from the Cauchy(median, scale) distribution
    565         * @throws MathException if an error occurs generating the random value
    566         * @since 2.2
    567         */
    568        public double nextCauchy(double median, double scale) throws MathException {
    569            return nextInversionDeviate(new CauchyDistributionImpl(median, scale));
    570        }
    571
    572        /**
    573         * Generates a random value from the {@link ChiSquaredDistributionImpl ChiSquare Distribution}.
    574         * This implementation uses {@link #nextInversionDeviate(ContinuousDistribution) inversion}
    575         * to generate random values.
    576         *
    577         * @param df the degrees of freedom of the ChiSquare distribution
    578         * @return random value sampled from the ChiSquare(df) distribution
    579         * @throws MathException if an error occurs generating the random value
    580         * @since 2.2
    581         */
    582        public double nextChiSquare(double df) throws MathException {
    583            return nextInversionDeviate(new ChiSquaredDistributionImpl(df));
    584        }
    585
    586        /**
    587         * Generates a random value from the {@link FDistributionImpl F Distribution}.
    588         * This implementation uses {@link #nextInversionDeviate(ContinuousDistribution) inversion}
    589         * to generate random values.
    590         *
    591         * @param numeratorDf the numerator degrees of freedom of the F distribution
    592         * @param denominatorDf the denominator degrees of freedom of the F distribution
    593         * @return random value sampled from the F(numeratorDf, denominatorDf) distribution
    594         * @throws MathException if an error occurs generating the random value
    595         * @since 2.2
    596         */
    597        public double nextF(double numeratorDf, double denominatorDf) throws MathException {
    598            return nextInversionDeviate(new FDistributionImpl(numeratorDf, denominatorDf));
    599        }
    600
    601        /**
    602         * Generates a random value from the {@link GammaDistributionImpl Gamma Distribution}.
    603         * This implementation uses {@link #nextInversionDeviate(ContinuousDistribution) inversion}
    604         * to generate random values.
    605         *
    606         * @param shape the median of the Gamma distribution
    607         * @param scale the scale parameter of the Gamma distribution
    608         * @return random value sampled from the Gamma(shape, scale) distribution
    609         * @throws MathException if an error occurs generating the random value
    610         * @since 2.2
    611         */
    612        public double nextGamma(double shape, double scale) throws MathException {
    613            return nextInversionDeviate(new GammaDistributionImpl(shape, scale));
    614        }
    615
    616        /**
    617         * Generates a random value from the {@link HypergeometricDistributionImpl Hypergeometric Distribution}.
    618         * This implementation uses {@link #nextInversionDeviate(IntegerDistribution) inversion}
    619         * to generate random values.
    620         *
    621         * @param populationSize the population size of the Hypergeometric distribution
    622         * @param numberOfSuccesses number of successes in the population of the Hypergeometric distribution
    623         * @param sampleSize the sample size of the Hypergeometric distribution
    624         * @return random value sampled from the Hypergeometric(numberOfSuccesses, sampleSize) distribution
    625         * @throws MathException if an error occurs generating the random value
    626         * @since 2.2
    627         */
    628        public int nextHypergeometric(int populationSize, int numberOfSuccesses, int sampleSize) throws MathException {
    629            return nextInversionDeviate(new HypergeometricDistributionImpl(populationSize, numberOfSuccesses, sampleSize));
    630        }
    631
    632        /**
    633         * Generates a random value from the {@link PascalDistributionImpl Pascal Distribution}.
    634         * This implementation uses {@link #nextInversionDeviate(IntegerDistribution) inversion}
    635         * to generate random values.
    636         *
    637         * @param r the number of successes of the Pascal distribution
    638         * @param p the probability of success of the Pascal distribution
    639         * @return random value sampled from the Pascal(r, p) distribution
    640         * @throws MathException if an error occurs generating the random value
    641         * @since 2.2
    642         */
    643        public int nextPascal(int r, double p) throws MathException {
    644            return nextInversionDeviate(new PascalDistributionImpl(r, p));
    645        }
    646
    647        /**
    648         * Generates a random value from the {@link TDistributionImpl T Distribution}.
    649         * This implementation uses {@link #nextInversionDeviate(ContinuousDistribution) inversion}
    650         * to generate random values.
    651         *
    652         * @param df the degrees of freedom of the T distribution
    653         * @return random value from the T(df) distribution
    654         * @throws MathException if an error occurs generating the random value
    655         * @since 2.2
    656         */
    657        public double nextT(double df) throws MathException {
    658            return nextInversionDeviate(new TDistributionImpl(df));
    659        }
    660
    661        /**
    662         * Generates a random value from the {@link WeibullDistributionImpl Weibull Distribution}.
    663         * This implementation uses {@link #nextInversionDeviate(ContinuousDistribution) inversion}
    664         * to generate random values.
    665         *
    666         * @param shape the shape parameter of the Weibull distribution
    667         * @param scale the scale parameter of the Weibull distribution
    668         * @return random value sampled from the Weibull(shape, size) distribution
    669         * @throws MathException if an error occurs generating the random value
    670         * @since 2.2
    671         */
    672        public double nextWeibull(double shape, double scale) throws MathException {
    673            return nextInversionDeviate(new WeibullDistributionImpl(shape, scale));
    674        }
    675
    676        /**
    677         * Generates a random value from the {@link ZipfDistributionImpl Zipf Distribution}.
    678         * This implementation uses {@link #nextInversionDeviate(IntegerDistribution) inversion}
    679         * to generate random values.
    680         *
    681         * @param numberOfElements the number of elements of the ZipfDistribution
    682         * @param exponent the exponent of the ZipfDistribution
    683         * @return random value sampled from the Zipf(numberOfElements, exponent) distribution
    684         * @throws MathException if an error occurs generating the random value
    685         * @since 2.2
    686         */
    687        public int nextZipf(int numberOfElements, double exponent) throws MathException {
    688            return nextInversionDeviate(new ZipfDistributionImpl(numberOfElements, exponent));
    689        }
    690
    691        /**
    692         * Returns the RandomGenerator used to generate non-secure random data.
    693         * <p>
    694         * Creates and initializes a default generator if null.
    695         * </p>
    696         *
    697         * @return the Random used to generate random data
    698         * @since 1.1
    699         */
    700        private RandomGenerator getRan() {
    701            if (rand == null) {
    702                rand = new JDKRandomGenerator();
    703                rand.setSeed(System.currentTimeMillis());
    704            }
    705            return rand;
    706        }
    707
    708        /**
    709         * Returns the SecureRandom used to generate secure random data.
    710         * <p>
    711         * Creates and initializes if null.
    712         * </p>
    713         *
    714         * @return the SecureRandom used to generate secure random data
    715         */
    716        private SecureRandom getSecRan() {
    717            if (secRand == null) {
    718                secRand = new SecureRandom();
    719                secRand.setSeed(System.currentTimeMillis());
    720            }
    721            return secRand;
    722        }
    723
    724        /**
    725         * Reseeds the random number generator with the supplied seed.
    726         * <p>
    727         * Will create and initialize if null.
    728         * </p>
    729         *
    730         * @param seed
    731         *            the seed value to use
    732         */
    733        public void reSeed(long seed) {
    734            if (rand == null) {
    735                rand = new JDKRandomGenerator();
    736            }
    737            rand.setSeed(seed);
    738        }
    739
    740        /**
    741         * Reseeds the secure random number generator with the current time in
    742         * milliseconds.
    743         * <p>
    744         * Will create and initialize if null.
    745         * </p>
    746         */
    747        public void reSeedSecure() {
    748            if (secRand == null) {
    749                secRand = new SecureRandom();
    750            }
    751            secRand.setSeed(System.currentTimeMillis());
    752        }
    753
    754        /**
    755         * Reseeds the secure random number generator with the supplied seed.
    756         * <p>
    757         * Will create and initialize if null.
    758         * </p>
    759         *
    760         * @param seed
    761         *            the seed value to use
    762         */
    763        public void reSeedSecure(long seed) {
    764            if (secRand == null) {
    765                secRand = new SecureRandom();
    766            }
    767            secRand.setSeed(seed);
    768        }
    769
    770        /**
    771         * Reseeds the random number generator with the current time in
    772         * milliseconds.
    773         */
    774        public void reSeed() {
    775            if (rand == null) {
    776                rand = new JDKRandomGenerator();
    777            }
    778            rand.setSeed(System.currentTimeMillis());
    779        }
    780
    781        /**
    782         * Sets the PRNG algorithm for the underlying SecureRandom instance using
    783         * the Security Provider API. The Security Provider API is defined in <a
    784         * href =
    785         * "http://java.sun.com/j2se/1.3/docs/guide/security/CryptoSpec.html#AppA">
    786         * Java Cryptography Architecture API Specification & Reference.</a>
    787         * <p>
    788         * <strong>USAGE NOTE:</strong> This method carries <i>significant</i>
    789         * overhead and may take several seconds to execute.
    790         * </p>
    791         *
    792         * @param algorithm
    793         *            the name of the PRNG algorithm
    794         * @param provider
    795         *            the name of the provider
    796         * @throws NoSuchAlgorithmException
    797         *             if the specified algorithm is not available
    798         * @throws NoSuchProviderException
    799         *             if the specified provider is not installed
    800         */
    801        public void setSecureAlgorithm(String algorithm, String provider)
    802                throws NoSuchAlgorithmException, NoSuchProviderException {
    803            secRand = SecureRandom.getInstance(algorithm, provider);
    804        }
    805
    806        /**
    807         * Generates an integer array of length <code>k</code> whose entries are
    808         * selected randomly, without repetition, from the integers
    809         * <code>0 through n-1</code> (inclusive).
    810         * <p>
    811         * Generated arrays represent permutations of <code>n</code> taken
    812         * <code>k</code> at a time.
    813         * </p>
    814         * <p>
    815         * <strong>Preconditions:</strong>
    816         * <ul>
    817         * <li> <code>k <= n</code></li>
    818         * <li> <code>n > 0</code></li>
    819         * </ul>
    820         * If the preconditions are not met, an IllegalArgumentException is thrown.
    821         * </p>
    822         * <p>
    823         * Uses a 2-cycle permutation shuffle. The shuffling process is described <a
    824         * href="http://www.maths.abdn.ac.uk/~igc/tch/mx4002/notes/node83.html">
    825         * here</a>.
    826         * </p>
    827         *
    828         * @param n
    829         *            domain of the permutation (must be positive)
    830         * @param k
    831         *            size of the permutation (must satisfy 0 < k <= n).
    832         * @return the random permutation as an int array
    833         * @throws NumberIsTooLargeException if {@code k > n}.
    834         * @throws NotStrictlyPositiveException if {@code k <= 0}.
    835         */
    836        public int[] nextPermutation(int n, int k) {
    837            if (k > n) {
    838                throw new NumberIsTooLargeException(LocalizedFormats.PERMUTATION_EXCEEDS_N,
    839                                                    k, n, true);
    840            }
    841            if (k == 0) {
    842                throw new NotStrictlyPositiveException(LocalizedFormats.PERMUTATION_SIZE,
    843                                                       k);
    844            }
    845
    846            int[] index = getNatural(n);
    847            shuffle(index, n - k);
    848            int[] result = new int[k];
    849            for (int i = 0; i < k; i++) {
    850                result[i] = index[n - i - 1];
    851            }
    852
    853            return result;
    854        }
    855
    856        /**
    857         * Uses a 2-cycle permutation shuffle to generate a random permutation.
    858         * <strong>Algorithm Description</strong>: Uses a 2-cycle permutation
    859         * shuffle to generate a random permutation of <code>c.size()</code> and
    860         * then returns the elements whose indexes correspond to the elements of the
    861         * generated permutation. This technique is described, and proven to
    862         * generate random samples, <a
    863         * href="http://www.maths.abdn.ac.uk/~igc/tch/mx4002/notes/node83.html">
    864         * here</a>
    865         *
    866         * @param c
    867         *            Collection to sample from.
    868         * @param k
    869         *            sample size.
    870         * @return the random sample.
    871         * @throws NumberIsTooLargeException if {@code k > c.size()}.
    872         * @throws NotStrictlyPositiveException if {@code k <= 0}.
    873         */
    874        public Object[] nextSample(Collection<?> c, int k) {
    875            int len = c.size();
    876            if (k > len) {
    877                throw new NumberIsTooLargeException(LocalizedFormats.SAMPLE_SIZE_EXCEEDS_COLLECTION_SIZE,
    878                                                    k, len, true);
    879            }
    880            if (k <= 0) {
    881                throw new NotStrictlyPositiveException(LocalizedFormats.NUMBER_OF_SAMPLES, k);
    882            }
    883
    884            Object[] objects = c.toArray();
    885            int[] index = nextPermutation(len, k);
    886            Object[] result = new Object[k];
    887            for (int i = 0; i < k; i++) {
    888                result[i] = objects[index[i]];
    889            }
    890            return result;
    891        }
    892
    893        /**
    894         * Generate a random deviate from the given distribution using the
    895         * <a href="http://en.wikipedia.org/wiki/Inverse_transform_sampling"> inversion method.</a>
    896         *
    897         * @param distribution Continuous distribution to generate a random value from
    898         * @return a random value sampled from the given distribution
    899         * @throws MathException if an error occurs computing the inverse cumulative distribution function
    900         * @since 2.2
    901         */
    902        public double nextInversionDeviate(ContinuousDistribution distribution) throws MathException {
    903            return distribution.inverseCumulativeProbability(nextUniform(0, 1));
    904
    905        }
    906
    907        /**
    908         * Generate a random deviate from the given distribution using the
    909         * <a href="http://en.wikipedia.org/wiki/Inverse_transform_sampling"> inversion method.</a>
    910         *
    911         * @param distribution Integer distribution to generate a random value from
    912         * @return a random value sampled from the given distribution
    913         * @throws MathException if an error occurs computing the inverse cumulative distribution function
    914         * @since 2.2
    915         */
    916        public int nextInversionDeviate(IntegerDistribution distribution) throws MathException {
    917            final double target = nextUniform(0, 1);
    918            final int glb = distribution.inverseCumulativeProbability(target);
    919            if (distribution.cumulativeProbability(glb) == 1.0d) { // No mass above
    920                return glb;
    921            } else {
    922                return glb + 1;
    923            }
    924        }
    925
    926        // ------------------------Private methods----------------------------------
    927
    928        /**
    929         * Uses a 2-cycle permutation shuffle to randomly re-order the last elements
    930         * of list.
    931         *
    932         * @param list
    933         *            list to be shuffled
    934         * @param end
    935         *            element past which shuffling begins
    936         */
    937        private void shuffle(int[] list, int end) {
    938            int target = 0;
    939            for (int i = list.length - 1; i >= end; i--) {
    940                if (i == 0) {
    941                    target = 0;
    942                } else {
    943                    target = nextInt(0, i);
    944                }
    945                int temp = list[target];
    946                list[target] = list[i];
    947                list[i] = temp;
    948            }
    949        }
    950
    951        /**
    952         * Returns an array representing n.
    953         *
    954         * @param n
    955         *            the natural number to represent
    956         * @return array with entries = elements of n
    957         */
    958        private int[] getNatural(int n) {
    959            int[] natural = new int[n];
    960            for (int i = 0; i < n; i++) {
    961                natural[i] = i;
    962            }
    963            return natural;
    964        }
    965
    966    }