third party library: apache-commons-mathandroid-cts-6.0_r9android-cts-6.0_r8android-cts-6.0_r7android-cts-6.0_r6android-cts-6.0_r5android-cts-6.0_r4android-cts-6.0_r32android-cts-6.0_r31android-cts-6.0_r30android-cts-6.0_r3android-cts-6.0_r29android-cts-6.0_r28android-cts-6.0_r27android-cts-6.0_r26android-cts-6.0_r25android-cts-6.0_r24android-cts-6.0_r23android-cts-6.0_r22android-cts-6.0_r21android-cts-6.0_r20android-cts-6.0_r2android-cts-6.0_r19android-cts-6.0_r18android-cts-6.0_r17android-cts-6.0_r16android-cts-6.0_r15android-cts-6.0_r14android-cts-6.0_r13android-cts-6.0_r12android-cts-6.0_r1android-6.0.1_r9android-6.0.1_r81android-6.0.1_r80android-6.0.1_r8android-6.0.1_r79android-6.0.1_r78android-6.0.1_r77android-6.0.1_r74android-6.0.1_r73android-6.0.1_r72android-6.0.1_r70android-6.0.1_r7android-6.0.1_r69android-6.0.1_r68android-6.0.1_r67android-6.0.1_r66android-6.0.1_r65android-6.0.1_r63android-6.0.1_r62android-6.0.1_r61android-6.0.1_r60android-6.0.1_r59android-6.0.1_r58android-6.0.1_r57android-6.0.1_r56android-6.0.1_r55android-6.0.1_r54android-6.0.1_r53android-6.0.1_r52android-6.0.1_r51android-6.0.1_r50android-6.0.1_r5android-6.0.1_r49android-6.0.1_r48android-6.0.1_r47android-6.0.1_r46android-6.0.1_r45android-6.0.1_r43android-6.0.1_r42android-6.0.1_r41android-6.0.1_r40android-6.0.1_r4android-6.0.1_r33android-6.0.1_r32android-6.0.1_r31android-6.0.1_r30android-6.0.1_r3android-6.0.1_r28android-6.0.1_r27android-6.0.1_r26android-6.0.1_r25android-6.0.1_r24android-6.0.1_r22android-6.0.1_r21android-6.0.1_r20android-6.0.1_r18android-6.0.1_r17android-6.0.1_r16android-6.0.1_r13android-6.0.1_r12android-6.0.1_r11android-6.0.1_r10android-6.0.1_r1android-6.0.0_r7android-6.0.0_r6android-6.0.0_r5android-6.0.0_r41android-6.0.0_r4android-6.0.0_r3android-6.0.0_r26android-6.0.0_r25android-6.0.0_r24android-6.0.0_r23android-6.0.0_r2android-6.0.0_r13android-6.0.0_r12android-6.0.0_r11android-6.0.0_r1marshmallow-releasemarshmallow-mr3-releasemarshmallow-mr2-releasemarshmallow-mr1-releasemarshmallow-mr1-devmarshmallow-dr1.6-releasemarshmallow-dr1.5-releasemarshmallow-dr1.5-devmarshmallow-dr-releasemarshmallow-dr-dragon-releasemarshmallow-dr-devmarshmallow-devmarshmallow-cts-release

Change-Id: I52a325624a7f0dd652b362a9840626d6d9f3c42b

26 files changed, 4961 insertions, 0 deletions

diff --git a/src/main/java/org/apache/commons/math/random/AbstractRandomGenerator.java b/src/main/java/org/apache/commons/math/random/AbstractRandomGenerator.java
new file mode 100644
index 0000000..afec63e
--- /dev/null
+++ b/src/main/java/org/apache/commons/math/random/AbstractRandomGenerator.java
@@ -0,0 +1,272 @@
+/*
+ * 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.math.random;
+
+import org.apache.commons.math.exception.NotStrictlyPositiveException;
+import org.apache.commons.math.util.FastMath;
+
+/**
+ * Abstract class implementing the {@link  RandomGenerator} interface.
+ * Default implementations for all methods other than {@link #nextDouble()} and
+ * {@link #setSeed(long)} are provided.
+ * <p>
+ * All data generation methods are based on {@code code nextDouble()}.
+ * Concrete implementations <strong>must</strong> override
+ * this method and <strong>should</strong> provide better / more
+ * performant implementations of the other methods if the underlying PRNG
+ * supplies them.</p>
+ *
+ * @since 1.1
+ * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $
+ */
+public abstract class AbstractRandomGenerator implements RandomGenerator {
+
+    /**
+     * Cached random normal value.  The default implementation for
+     * {@link #nextGaussian} generates pairs of values and this field caches the
+     * second value so that the full algorithm is not executed for every
+     * activation.  The value {@code Double.NaN} signals that there is
+     * no cached value.  Use {@link #clear} to clear the cached value.
+     */
+    private double cachedNormalDeviate = Double.NaN;
+
+    /**
+     * Construct a RandomGenerator.
+     */
+    public AbstractRandomGenerator() {
+        super();
+
+    }
+
+    /**
+     * Clears the cache used by the default implementation of
+     * {@link #nextGaussian}. Implemementations that do not override the
+     * default implementation of {@code nextGaussian} should call this
+     * method in the implementation of {@link #setSeed(long)}
+     */
+    public void clear() {
+        cachedNormalDeviate = Double.NaN;
+    }
+
+    /** {@inheritDoc} */
+    public void setSeed(int seed) {
+        setSeed((long) seed);
+    }
+
+    /** {@inheritDoc} */
+    public void setSeed(int[] seed) {
+        // the following number is the largest prime that fits in 32 bits (it is 2^32 - 5)
+        final long prime = 4294967291l;
+
+        long combined = 0l;
+        for (int s : seed) {
+            combined = combined * prime + s;
+        }
+        setSeed(combined);
+    }
+
+    /**
+     * Sets the seed of the underyling random number generator using a
+     * {@code long} seed.  Sequences of values generated starting with the
+     * same seeds should be identical.
+     * <p>
+     * Implementations that do not override the default implementation of
+     * {@code nextGaussian} should include a call to {@link #clear} in the
+     * implementation of this method.</p>
+     *
+     * @param seed the seed value
+     */
+    public abstract void setSeed(long seed);
+
+    /**
+     * Generates random bytes and places them into a user-supplied
+     * byte array.  The number of random bytes produced is equal to
+     * the length of the byte array.
+     * <p>
+     * The default implementation fills the array with bytes extracted from
+     * random integers generated using {@link #nextInt}.</p>
+     *
+     * @param bytes the non-null byte array in which to put the
+     * random bytes
+     */
+    public void nextBytes(byte[] bytes) {
+        int bytesOut = 0;
+        while (bytesOut < bytes.length) {
+          int randInt = nextInt();
+          for (int i = 0; i < 3; i++) {
+              if ( i > 0) {
+                  randInt = randInt >> 8;
+              }
+              bytes[bytesOut++] = (byte) randInt;
+              if (bytesOut == bytes.length) {
+                  return;
+              }
+          }
+        }
+    }
+
+     /**
+     * Returns the next pseudorandom, uniformly distributed {@code int}
+     * value from this random number generator's sequence.
+     * All 2<font size="-1"><sup>32</sup></font> possible {@code int} values
+     * should be produced with  (approximately) equal probability.
+     * <p>
+     * The default implementation provided here returns
+     * <pre>
+     * <code>(int) (nextDouble() * Integer.MAX_VALUE)</code>
+     * </pre></p>
+     *
+     * @return the next pseudorandom, uniformly distributed {@code int}
+     *  value from this random number generator's sequence
+     */
+    public int nextInt() {
+        return (int) (nextDouble() * Integer.MAX_VALUE);
+    }
+
+    /**
+     * Returns a pseudorandom, uniformly distributed {@code int} value
+     * between 0 (inclusive) and the specified value (exclusive), drawn from
+     * this random number generator's sequence.
+     * <p>
+     * The default implementation returns
+     * <pre>
+     * <code>(int) (nextDouble() * n</code>
+     * </pre></p>
+     *
+     * @param n the bound on the random number to be returned.  Must be
+     * positive.
+     * @return  a pseudorandom, uniformly distributed {@code int}
+     * value between 0 (inclusive) and n (exclusive).
+     * @throws NotStrictlyPositiveException if {@code n <= 0}.
+     */
+    public int nextInt(int n) {
+        if (n <= 0 ) {
+            throw new NotStrictlyPositiveException(n);
+        }
+        int result = (int) (nextDouble() * n);
+        return result < n ? result : n - 1;
+    }
+
+     /**
+     * Returns the next pseudorandom, uniformly distributed {@code long}
+     * value from this random number generator's sequence.  All
+     * 2<font size="-1"><sup>64</sup></font> possible {@code long} values
+     * should be produced with (approximately) equal probability.
+     * <p>
+     * The default implementation returns
+     * <pre>
+     * <code>(long) (nextDouble() * Long.MAX_VALUE)</code>
+     * </pre></p>
+     *
+     * @return  the next pseudorandom, uniformly distributed {@code long}
+     *value from this random number generator's sequence
+     */
+    public long nextLong() {
+        return (long) (nextDouble() * Long.MAX_VALUE);
+    }
+
+    /**
+     * Returns the next pseudorandom, uniformly distributed
+     * {@code boolean} value from this random number generator's
+     * sequence.
+     * <p>
+     * The default implementation returns
+     * <pre>
+     * <code>nextDouble() <= 0.5</code>
+     * </pre></p>
+     *
+     * @return  the next pseudorandom, uniformly distributed
+     * {@code boolean} value from this random number generator's
+     * sequence
+     */
+    public boolean nextBoolean() {
+        return nextDouble() <= 0.5;
+    }
+
+     /**
+     * Returns the next pseudorandom, uniformly distributed {@code float}
+     * value between {@code 0.0} and {@code 1.0} from this random
+     * number generator's sequence.
+     * <p>
+     * The default implementation returns
+     * <pre>
+     * <code>(float) nextDouble() </code>
+     * </pre></p>
+     *
+     * @return  the next pseudorandom, uniformly distributed {@code float}
+     * value between {@code 0.0} and {@code 1.0} from this
+     * random number generator's sequence
+     */
+    public float nextFloat() {
+        return (float) nextDouble();
+    }
+
+    /**
+     * Returns the next pseudorandom, uniformly distributed
+     * {@code double} value between {@code 0.0} and
+     * {@code 1.0} from this random number generator's sequence.
+     * <p>
+     * This method provides the underlying source of random data used by the
+     * other methods.</p>
+     *
+     * @return  the next pseudorandom, uniformly distributed
+     *  {@code double} value between {@code 0.0} and
+     *  {@code 1.0} from this random number generator's sequence
+     */
+    public abstract double nextDouble();
+
+    /**
+     * Returns the next pseudorandom, Gaussian ("normally") distributed
+     * {@code double} value with mean {@code 0.0} and standard
+     * deviation {@code 1.0} from this random number generator's sequence.
+     * <p>
+     * The default implementation uses the <em>Polar Method</em>
+     * due to G.E.P. Box, M.E. Muller and G. Marsaglia, as described in
+     * D. Knuth, <u>The Art of Computer Programming</u>, 3.4.1C.</p>
+     * <p>
+     * The algorithm generates a pair of independent random values.  One of
+     * these is cached for reuse, so the full algorithm is not executed on each
+     * activation.  Implementations that do not override this method should
+     * make sure to call {@link #clear} to clear the cached value in the
+     * implementation of {@link #setSeed(long)}.</p>
+     *
+     * @return  the next pseudorandom, Gaussian ("normally") distributed
+     * {@code double} value with mean {@code 0.0} and
+     * standard deviation {@code 1.0} from this random number
+     *  generator's sequence
+     */
+    public double nextGaussian() {
+        if (!Double.isNaN(cachedNormalDeviate)) {
+            double dev = cachedNormalDeviate;
+            cachedNormalDeviate = Double.NaN;
+            return dev;
+        }
+        double v1 = 0;
+        double v2 = 0;
+        double s = 1;
+        while (s >=1 ) {
+            v1 = 2 * nextDouble() - 1;
+            v2 = 2 * nextDouble() - 1;
+            s = v1 * v1 + v2 * v2;
+        }
+        if (s != 0) {
+            s = FastMath.sqrt(-2 * FastMath.log(s) / s);
+        }
+        cachedNormalDeviate = v2 * s;
+        return v1 * s;
+    }
+}
diff --git a/src/main/java/org/apache/commons/math/random/AbstractWell.java b/src/main/java/org/apache/commons/math/random/AbstractWell.java
new file mode 100644
index 0000000..96e18a2
--- /dev/null
+++ b/src/main/java/org/apache/commons/math/random/AbstractWell.java
@@ -0,0 +1,186 @@
+/*
+ * 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.math.random;
+
+import java.io.Serializable;
+
+
+/** This abstract class implements the WELL class of pseudo-random number generator
+ * from Fran&ccedil;ois Panneton, Pierre L'Ecuyer and Makoto Matsumoto.
+
+ * <p>This generator is described in a paper by Fran&ccedil;ois Panneton,
+ * Pierre L'Ecuyer and Makoto Matsumoto <a
+ * href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng.pdf">Improved
+ * Long-Period Generators Based on Linear Recurrences Modulo 2</a> ACM
+ * Transactions on Mathematical Software, 32, 1 (2006). The errata for the paper
+ * are in <a href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng-errata.txt">wellrng-errata.txt</a>.</p>
+
+ * @see <a href="http://www.iro.umontreal.ca/~panneton/WELLRNG.html">WELL Random number generator</a>
+ * @version $Revision: 1003892 $ $Date: 2010-10-02 23:28:56 +0200 (sam. 02 oct. 2010) $
+ * @since 2.2
+
+ */
+public abstract class AbstractWell extends BitsStreamGenerator implements Serializable {
+
+    /** Serializable version identifier. */
+    private static final long serialVersionUID = -817701723016583596L;
+
+    /** Current index in the bytes pool. */
+    protected int index;
+
+    /** Bytes pool. */
+    protected final int[] v;
+
+    /** Index indirection table giving for each index its predecessor taking table size into account. */
+    protected final int[] iRm1;
+
+    /** Index indirection table giving for each index its second predecessor taking table size into account. */
+    protected final int[] iRm2;
+
+    /** Index indirection table giving for each index the value index + m1 taking table size into account. */
+    protected final int[] i1;
+
+    /** Index indirection table giving for each index the value index + m2 taking table size into account. */
+    protected final int[] i2;
+
+    /** Index indirection table giving for each index the value index + m3 taking table size into account. */
+    protected final int[] i3;
+
+    /** Creates a new random number generator.
+     * <p>The instance is initialized using the current time as the
+     * seed.</p>
+     * @param k number of bits in the pool (not necessarily a multiple of 32)
+     * @param m1 first parameter of the algorithm
+     * @param m2 second parameter of the algorithm
+     * @param m3 third parameter of the algorithm
+     */
+    protected AbstractWell(final int k, final int m1, final int m2, final int m3) {
+        this(k, m1, m2, m3, System.currentTimeMillis());
+    }
+
+    /** Creates a new random number generator using a single int seed.
+     * @param k number of bits in the pool (not necessarily a multiple of 32)
+     * @param m1 first parameter of the algorithm
+     * @param m2 second parameter of the algorithm
+     * @param m3 third parameter of the algorithm
+     * @param seed the initial seed (32 bits integer)
+     */
+    protected AbstractWell(final int k, final int m1, final int m2, final int m3, final int seed) {
+        this(k, m1, m2, m3, new int[] { seed });
+    }
+
+    /** Creates a new random number generator using an int array seed.
+     * @param k number of bits in the pool (not necessarily a multiple of 32)
+     * @param m1 first parameter of the algorithm
+     * @param m2 second parameter of the algorithm
+     * @param m3 third parameter of the algorithm
+     * @param seed the initial seed (32 bits integers array), if null
+     * the seed of the generator will be related to the current time
+     */
+    protected AbstractWell(final int k, final int m1, final int m2, final int m3, final int[] seed) {
+
+        // the bits pool contains k bits, k = r w - p where r is the number
+        // of w bits blocks, w is the block size (always 32 in the original paper)
+        // and p is the number of unused bits in the last block
+        final int w = 32;
+        final int r = (k + w - 1) / w;
+        this.v      = new int[r];
+        this.index  = 0;
+
+        // precompute indirection index tables. These tables are used for optimizing access
+        // they allow saving computations like "(j + r - 2) % r" with costly modulo operations
+        iRm1 = new int[r];
+        iRm2 = new int[r];
+        i1   = new int[r];
+        i2   = new int[r];
+        i3   = new int[r];
+        for (int j = 0; j < r; ++j) {
+            iRm1[j] = (j + r - 1) % r;
+            iRm2[j] = (j + r - 2) % r;
+            i1[j]   = (j + m1)    % r;
+            i2[j]   = (j + m2)    % r;
+            i3[j]   = (j + m3)    % r;
+        }
+
+        // initialize the pool content
+        setSeed(seed);
+
+    }
+
+    /** Creates a new random number generator using a single long seed.
+     * @param k number of bits in the pool (not necessarily a multiple of 32)
+     * @param m1 first parameter of the algorithm
+     * @param m2 second parameter of the algorithm
+     * @param m3 third parameter of the algorithm
+     * @param seed the initial seed (64 bits integer)
+     */
+    protected AbstractWell(final int k, final int m1, final int m2, final int m3, final long seed) {
+        this(k, m1, m2, m3, new int[] { (int) (seed >>> 32), (int) (seed & 0xffffffffl) });
+    }
+
+    /** Reinitialize the generator as if just built with the given int seed.
+     * <p>The state of the generator is exactly the same as a new
+     * generator built with the same seed.</p>
+     * @param seed the initial seed (32 bits integer)
+     */
+    @Override
+    public void setSeed(final int seed) {
+        setSeed(new int[] { seed });
+    }
+
+    /** Reinitialize the generator as if just built with the given int array seed.
+     * <p>The state of the generator is exactly the same as a new
+     * generator built with the same seed.</p>
+     * @param seed the initial seed (32 bits integers array), if null
+     * the seed of the generator will be related to the current time
+     */
+    @Override
+    public void setSeed(final int[] seed) {
+
+        if (seed == null) {
+            setSeed(System.currentTimeMillis());
+            return;
+        }
+
+        System.arraycopy(seed, 0, v, 0, Math.min(seed.length, v.length));
+
+        if (seed.length < v.length) {
+            for (int i = seed.length; i < v.length; ++i) {
+                final long l = v[i - seed.length];
+                v[i] = (int) ((1812433253l * (l ^ (l >> 30)) + i) & 0xffffffffL);
+            }
+        }
+
+        index = 0;
+
+    }
+
+    /** Reinitialize the generator as if just built with the given long seed.
+     * <p>The state of the generator is exactly the same as a new
+     * generator built with the same seed.</p>
+     * @param seed the initial seed (64 bits integer)
+     */
+    @Override
+    public void setSeed(final long seed) {
+        setSeed(new int[] { (int) (seed >>> 32), (int) (seed & 0xffffffffl) });
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    protected abstract int next(final int bits);
+
+}
diff --git a/src/main/java/org/apache/commons/math/random/BitsStreamGenerator.java b/src/main/java/org/apache/commons/math/random/BitsStreamGenerator.java
new file mode 100644
index 0000000..8979473
--- /dev/null
+++ b/src/main/java/org/apache/commons/math/random/BitsStreamGenerator.java
@@ -0,0 +1,153 @@
+/*
+ * 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.math.random;
+
+import org.apache.commons.math.exception.NotStrictlyPositiveException;
+import org.apache.commons.math.util.FastMath;
+
+/** Base class for random number generators that generates bits streams.
+
+ * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $
+ * @since 2.0
+
+ */
+public abstract class BitsStreamGenerator implements RandomGenerator {
+
+    /** Next gaussian. */
+    private double nextGaussian;
+
+    /** Creates a new random number generator.
+     */
+    public BitsStreamGenerator() {
+        nextGaussian = Double.NaN;
+    }
+
+    /** {@inheritDoc} */
+    public abstract void setSeed(int seed);
+
+    /** {@inheritDoc} */
+    public abstract void setSeed(int[] seed);
+
+    /** {@inheritDoc} */
+    public abstract void setSeed(long seed);
+
+    /** Generate next pseudorandom number.
+     * <p>This method is the core generation algorithm. It is used by all the
+     * public generation methods for the various primitive types {@link
+     * #nextBoolean()}, {@link #nextBytes(byte[])}, {@link #nextDouble()},
+     * {@link #nextFloat()}, {@link #nextGaussian()}, {@link #nextInt()},
+     * {@link #next(int)} and {@link #nextLong()}.</p>
+     * @param bits number of random bits to produce
+     * @return random bits generated
+     */
+    protected abstract int next(int bits);
+
+    /** {@inheritDoc} */
+    public boolean nextBoolean() {
+        return next(1) != 0;
+    }
+
+    /** {@inheritDoc} */
+    public void nextBytes(byte[] bytes) {
+        int i = 0;
+        final int iEnd = bytes.length - 3;
+        while (i < iEnd) {
+            final int random = next(32);
+            bytes[i]     = (byte) (random & 0xff);
+            bytes[i + 1] = (byte) ((random >>  8) & 0xff);
+            bytes[i + 2] = (byte) ((random >> 16) & 0xff);
+            bytes[i + 3] = (byte) ((random >> 24) & 0xff);
+            i += 4;
+        }
+        int random = next(32);
+        while (i < bytes.length) {
+            bytes[i++] = (byte) (random & 0xff);
+            random     = random >> 8;
+        }
+    }
+
+    /** {@inheritDoc} */
+    public double nextDouble() {
+        final long high = ((long) next(26)) << 26;
+        final int  low  = next(26);
+        return (high | low) * 0x1.0p-52d;
+    }
+
+    /** {@inheritDoc} */
+    public float nextFloat() {
+        return next(23) * 0x1.0p-23f;
+    }
+
+    /** {@inheritDoc} */
+    public double nextGaussian() {
+
+        final double random;
+        if (Double.isNaN(nextGaussian)) {
+            // generate a new pair of gaussian numbers
+            final double x = nextDouble();
+            final double y = nextDouble();
+            final double alpha = 2 * FastMath.PI * x;
+            final double r      = FastMath.sqrt(-2 * FastMath.log(y));
+            random       = r * FastMath.cos(alpha);
+            nextGaussian = r * FastMath.sin(alpha);
+        } else {
+            // use the second element of the pair already generated
+            random = nextGaussian;
+            nextGaussian = Double.NaN;
+        }
+
+        return random;
+
+    }
+
+    /** {@inheritDoc} */
+    public int nextInt() {
+        return next(32);
+    }
+
+    /** {@inheritDoc} */
+    public int nextInt(int n) throws IllegalArgumentException {
+
+        if (n < 1) {
+            throw new NotStrictlyPositiveException(n);
+        }
+
+        // find bit mask for n
+        int mask = n;
+        mask |= mask >> 1;
+        mask |= mask >> 2;
+        mask |= mask >> 4;
+        mask |= mask >> 8;
+        mask |= mask >> 16;
+
+        while (true) {
+            final int random = next(32) & mask;
+            if (random < n) {
+                return random;
+            }
+        }
+
+    }
+
+    /** {@inheritDoc} */
+    public long nextLong() {
+        final long high  = ((long) next(32)) << 32;
+        final long  low  = ((long) next(32)) & 0xffffffffL;
+        return high | low;
+    }
+
+}
diff --git a/src/main/java/org/apache/commons/math/random/CorrelatedRandomVectorGenerator.java b/src/main/java/org/apache/commons/math/random/CorrelatedRandomVectorGenerator.java
new file mode 100644
index 0000000..f1df51d
--- /dev/null
+++ b/src/main/java/org/apache/commons/math/random/CorrelatedRandomVectorGenerator.java
@@ -0,0 +1,304 @@
+/*
+ * 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.math.random;
+
+import org.apache.commons.math.DimensionMismatchException;
+import org.apache.commons.math.linear.MatrixUtils;
+import org.apache.commons.math.linear.NotPositiveDefiniteMatrixException;
+import org.apache.commons.math.linear.RealMatrix;
+import org.apache.commons.math.util.FastMath;
+
+/**
+ * A {@link RandomVectorGenerator} that generates vectors with with
+ * correlated components.
+ * <p>Random vectors with correlated components are built by combining
+ * the uncorrelated components of another random vector in such a way that
+ * the resulting correlations are the ones specified by a positive
+ * definite covariance matrix.</p>
+ * <p>The main use for correlated random vector generation is for Monte-Carlo
+ * simulation of physical problems with several variables, for example to
+ * generate error vectors to be added to a nominal vector. A particularly
+ * interesting case is when the generated vector should be drawn from a <a
+ * href="http://en.wikipedia.org/wiki/Multivariate_normal_distribution">
+ * Multivariate Normal Distribution</a>. The approach using a Cholesky
+ * decomposition is quite usual in this case. However, it can be extended
+ * to other cases as long as the underlying random generator provides
+ * {@link NormalizedRandomGenerator normalized values} like {@link
+ * GaussianRandomGenerator} or {@link UniformRandomGenerator}.</p>
+ * <p>Sometimes, the covariance matrix for a given simulation is not
+ * strictly positive definite. This means that the correlations are
+ * not all independent from each other. In this case, however, the non
+ * strictly positive elements found during the Cholesky decomposition
+ * of the covariance matrix should not be negative either, they
+ * should be null. Another non-conventional extension handling this case
+ * is used here. Rather than computing <code>C = U<sup>T</sup>.U</code>
+ * where <code>C</code> is the covariance matrix and <code>U</code>
+ * is an upper-triangular matrix, we compute <code>C = B.B<sup>T</sup></code>
+ * where <code>B</code> is a rectangular matrix having
+ * more rows than columns. The number of columns of <code>B</code> is
+ * the rank of the covariance matrix, and it is the dimension of the
+ * uncorrelated random vector that is needed to compute the component
+ * of the correlated vector. This class handles this situation
+ * automatically.</p>
+ *
+ * @version $Revision: 1043908 $ $Date: 2010-12-09 12:53:14 +0100 (jeu. 09 déc. 2010) $
+ * @since 1.2
+ */
+
+public class CorrelatedRandomVectorGenerator
+    implements RandomVectorGenerator {
+
+    /** Mean vector. */
+    private final double[] mean;
+
+    /** Underlying generator. */
+    private final NormalizedRandomGenerator generator;
+
+    /** Storage for the normalized vector. */
+    private final double[] normalized;
+
+    /** Permutated Cholesky root of the covariance matrix. */
+    private RealMatrix root;
+
+    /** Rank of the covariance matrix. */
+    private int rank;
+
+    /** Simple constructor.
+     * <p>Build a correlated random vector generator from its mean
+     * vector and covariance matrix.</p>
+     * @param mean expected mean values for all components
+     * @param covariance covariance matrix
+     * @param small diagonal elements threshold under which  column are
+     * considered to be dependent on previous ones and are discarded
+     * @param generator underlying generator for uncorrelated normalized
+     * components
+     * @exception IllegalArgumentException if there is a dimension
+     * mismatch between the mean vector and the covariance matrix
+     * @exception NotPositiveDefiniteMatrixException if the
+     * covariance matrix is not strictly positive definite
+     * @exception DimensionMismatchException if the mean and covariance
+     * arrays dimensions don't match
+     */
+    public CorrelatedRandomVectorGenerator(double[] mean,
+                                           RealMatrix covariance, double small,
+                                           NormalizedRandomGenerator generator)
+    throws NotPositiveDefiniteMatrixException, DimensionMismatchException {
+
+        int order = covariance.getRowDimension();
+        if (mean.length != order) {
+            throw new DimensionMismatchException(mean.length, order);
+        }
+        this.mean = mean.clone();
+
+        decompose(covariance, small);
+
+        this.generator = generator;
+        normalized = new double[rank];
+
+    }
+
+    /** Simple constructor.
+     * <p>Build a null mean random correlated vector generator from its
+     * covariance matrix.</p>
+     * @param covariance covariance matrix
+     * @param small diagonal elements threshold under which  column are
+     * considered to be dependent on previous ones and are discarded
+     * @param generator underlying generator for uncorrelated normalized
+     * components
+     * @exception NotPositiveDefiniteMatrixException if the
+     * covariance matrix is not strictly positive definite
+     */
+    public CorrelatedRandomVectorGenerator(RealMatrix covariance, double small,
+                                           NormalizedRandomGenerator generator)
+    throws NotPositiveDefiniteMatrixException {
+
+        int order = covariance.getRowDimension();
+        mean = new double[order];
+        for (int i = 0; i < order; ++i) {
+            mean[i] = 0;
+        }
+
+        decompose(covariance, small);
+
+        this.generator = generator;
+        normalized = new double[rank];
+
+    }
+
+    /** Get the underlying normalized components generator.
+     * @return underlying uncorrelated components generator
+     */
+    public NormalizedRandomGenerator getGenerator() {
+        return generator;
+    }
+
+    /** Get the root of the covariance matrix.
+     * The root is the rectangular matrix <code>B</code> such that
+     * the covariance matrix is equal to <code>B.B<sup>T</sup></code>
+     * @return root of the square matrix
+     * @see #getRank()
+     */
+    public RealMatrix getRootMatrix() {
+        return root;
+    }
+
+    /** Get the rank of the covariance matrix.
+     * The rank is the number of independent rows in the covariance
+     * matrix, it is also the number of columns of the rectangular
+     * matrix of the decomposition.
+     * @return rank of the square matrix.
+     * @see #getRootMatrix()
+     */
+    public int getRank() {
+        return rank;
+    }
+
+    /** Decompose the original square matrix.
+     * <p>The decomposition is based on a Choleski decomposition
+     * where additional transforms are performed:
+     * <ul>
+     *   <li>the rows of the decomposed matrix are permuted</li>
+     *   <li>columns with the too small diagonal element are discarded</li>
+     *   <li>the matrix is permuted</li>
+     * </ul>
+     * This means that rather than computing M = U<sup>T</sup>.U where U
+     * is an upper triangular matrix, this method computed M=B.B<sup>T</sup>
+     * where B is a rectangular matrix.
+     * @param covariance covariance matrix
+     * @param small diagonal elements threshold under which  column are
+     * considered to be dependent on previous ones and are discarded
+     * @exception NotPositiveDefiniteMatrixException if the
+     * covariance matrix is not strictly positive definite
+     */
+    private void decompose(RealMatrix covariance, double small)
+    throws NotPositiveDefiniteMatrixException {
+
+        int order = covariance.getRowDimension();
+        double[][] c = covariance.getData();
+        double[][] b = new double[order][order];
+
+        int[] swap  = new int[order];
+        int[] index = new int[order];
+        for (int i = 0; i < order; ++i) {
+            index[i] = i;
+        }
+
+        rank = 0;
+        for (boolean loop = true; loop;) {
+
+            // find maximal diagonal element
+            swap[rank] = rank;
+            for (int i = rank + 1; i < order; ++i) {
+                int ii  = index[i];
+                int isi = index[swap[i]];
+                if (c[ii][ii] > c[isi][isi]) {
+                    swap[rank] = i;
+                }
+            }
+
+
+            // swap elements
+            if (swap[rank] != rank) {
+                int tmp = index[rank];
+                index[rank] = index[swap[rank]];
+                index[swap[rank]] = tmp;
+            }
+
+            // check diagonal element
+            int ir = index[rank];
+            if (c[ir][ir] < small) {
+
+                if (rank == 0) {
+                    throw new NotPositiveDefiniteMatrixException();
+                }
+
+                // check remaining diagonal elements
+                for (int i = rank; i < order; ++i) {
+                    if (c[index[i]][index[i]] < -small) {
+                        // there is at least one sufficiently negative diagonal element,
+                        // the covariance matrix is wrong
+                        throw new NotPositiveDefiniteMatrixException();
+                    }
+                }
+
+                // all remaining diagonal elements are close to zero,
+                // we consider we have found the rank of the covariance matrix
+                ++rank;
+                loop = false;
+
+            } else {
+
+                // transform the matrix
+                double sqrt = FastMath.sqrt(c[ir][ir]);
+                b[rank][rank] = sqrt;
+                double inverse = 1 / sqrt;
+                for (int i = rank + 1; i < order; ++i) {
+                    int ii = index[i];
+                    double e = inverse * c[ii][ir];
+                    b[i][rank] = e;
+                    c[ii][ii] -= e * e;
+                    for (int j = rank + 1; j < i; ++j) {
+                        int ij = index[j];
+                        double f = c[ii][ij] - e * b[j][rank];
+                        c[ii][ij] = f;
+                        c[ij][ii] = f;
+                    }
+                }
+
+                // prepare next iteration
+                loop = ++rank < order;
+
+            }
+
+        }
+
+        // build the root matrix
+        root = MatrixUtils.createRealMatrix(order, rank);
+        for (int i = 0; i < order; ++i) {
+            for (int j = 0; j < rank; ++j) {
+                root.setEntry(index[i], j, b[i][j]);
+            }
+        }
+
+    }
+
+    /** Generate a correlated random vector.
+     * @return a random vector as an array of double. The returned array
+     * is created at each call, the caller can do what it wants with it.
+     */
+    public double[] nextVector() {
+
+        // generate uncorrelated vector
+        for (int i = 0; i < rank; ++i) {
+            normalized[i] = generator.nextNormalizedDouble();
+        }
+
+        // compute correlated vector
+        double[] correlated = new double[mean.length];
+        for (int i = 0; i < correlated.length; ++i) {
+            correlated[i] = mean[i];
+            for (int j = 0; j < rank; ++j) {
+                correlated[i] += root.getEntry(i, j) * normalized[j];
+            }
+        }
+
+        return correlated;
+
+    }
+
+}
diff --git a/src/main/java/org/apache/commons/math/random/EmpiricalDistribution.java b/src/main/java/org/apache/commons/math/random/EmpiricalDistribution.java
new file mode 100644
index 0000000..7f08a06
--- /dev/null
+++ b/src/main/java/org/apache/commons/math/random/EmpiricalDistribution.java
@@ -0,0 +1,133 @@
+/*
+ * 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.math.random;
+
+import java.io.IOException;
+import java.io.File;
+import java.net.URL;
+import java.util.List;
+
+import org.apache.commons.math.stat.descriptive.StatisticalSummary;
+import org.apache.commons.math.stat.descriptive.SummaryStatistics;
+
+/**
+ * Represents an <a href="http://random.mat.sbg.ac.at/~ste/dipl/node11.html">
+ * empirical probability distribution</a> -- a probability distribution derived
+ * from observed data without making any assumptions about the functional form
+ * of the population distribution that the data come from.<p>
+ * Implementations of this interface maintain data structures, called
+ * <i>distribution digests</i>, that describe empirical distributions and
+ * support the following operations: <ul>
+ * <li>loading the distribution from a file of observed data values</li>
+ * <li>dividing the input data into "bin ranges" and reporting bin frequency
+ *     counts (data for histogram)</li>
+ * <li>reporting univariate statistics describing the full set of data values
+ *     as well as the observations within each bin</li>
+ * <li>generating random values from the distribution</li>
+ * </ul>
+ * Applications can use <code>EmpiricalDistribution</code> implementations to
+ * build grouped frequency histograms representing the input data or to
+ * generate random values "like" those in the input file -- i.e., the values
+ * generated will follow the distribution of the values in the file.</p>
+ *
+ * @version $Revision: 817128 $ $Date: 2009-09-21 03:30:53 +0200 (lun. 21 sept. 2009) $
+ */
+public interface EmpiricalDistribution {
+
+    /**
+     * Computes the empirical distribution from the provided
+     * array of numbers.
+     *
+     * @param dataArray the data array
+     */
+    void load(double[] dataArray);
+
+    /**
+     * Computes the empirical distribution from the input file.
+     *
+     * @param file the input file
+     * @throws IOException if an IO error occurs
+     */
+    void load(File file) throws IOException;
+
+    /**
+     * Computes the empirical distribution using data read from a URL.
+     *
+     * @param url url of the input file
+     * @throws IOException if an IO error occurs
+     */
+    void load(URL url) throws IOException;
+
+    /**
+     * Generates a random value from this distribution.
+     * <strong>Preconditions:</strong><ul>
+     * <li>the distribution must be loaded before invoking this method</li></ul>
+     * @return the random value.
+     *
+     * @throws IllegalStateException if the distribution has not been loaded
+     */
+    double getNextValue() throws IllegalStateException;
+
+
+    /**
+     * Returns a
+     * {@link org.apache.commons.math.stat.descriptive.StatisticalSummary}
+     * describing this distribution.
+     * <strong>Preconditions:</strong><ul>
+     * <li>the distribution must be loaded before invoking this method</li>
+     * </ul>
+     *
+     * @return the sample statistics
+     * @throws IllegalStateException if the distribution has not been loaded
+     */
+    StatisticalSummary getSampleStats() throws IllegalStateException;
+
+    /**
+     * Property indicating whether or not the distribution has been loaded.
+     *
+     * @return true if the distribution has been loaded
+     */
+    boolean isLoaded();
+
+     /**
+     * Returns the number of bins.
+     *
+     * @return the number of bins
+     */
+    int getBinCount();
+
+    /**
+     * Returns a list of
+     * {@link org.apache.commons.math.stat.descriptive.SummaryStatistics}
+     * containing statistics describing the values in each of the bins.  The
+     * List is indexed on the bin number.
+     *
+     * @return List of bin statistics
+     */
+    List<SummaryStatistics> getBinStats();
+
+    /**
+     * Returns the array of upper bounds for the bins.  Bins are: <br/>
+     * [min,upperBounds[0]],(upperBounds[0],upperBounds[1]],...,
+     *  (upperBounds[binCount-2], upperBounds[binCount-1] = max].
+     *
+     * @return array of bin upper bounds
+     */
+    double[] getUpperBounds();
+
+}
diff --git a/src/main/java/org/apache/commons/math/random/EmpiricalDistributionImpl.java b/src/main/java/org/apache/commons/math/random/EmpiricalDistributionImpl.java
new file mode 100644
index 0000000..a3ae8a7
--- /dev/null
+++ b/src/main/java/org/apache/commons/math/random/EmpiricalDistributionImpl.java
@@ -0,0 +1,479 @@
+/*
+ * 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.math.random;
+
+import java.io.BufferedReader;
+import java.io.File;
+import java.io.FileReader;
+import java.io.IOException;
+import java.io.InputStreamReader;
+import java.io.Serializable;
+import java.net.URL;
+import java.util.ArrayList;
+import java.util.List;
+
+import org.apache.commons.math.MathRuntimeException;
+import org.apache.commons.math.exception.util.LocalizedFormats;
+import org.apache.commons.math.stat.descriptive.StatisticalSummary;
+import org.apache.commons.math.stat.descriptive.SummaryStatistics;
+import org.apache.commons.math.util.FastMath;
+
+/**
+ * Implements <code>EmpiricalDistribution</code> interface.  This implementation
+ * uses what amounts to the
+ * <a href="http://nedwww.ipac.caltech.edu/level5/March02/Silverman/Silver2_6.html">
+ * Variable Kernel Method</a> with Gaussian smoothing:<p>
+ * <strong>Digesting the input file</strong>
+ * <ol><li>Pass the file once to compute min and max.</li>
+ * <li>Divide the range from min-max into <code>binCount</code> "bins."</li>
+ * <li>Pass the data file again, computing bin counts and univariate
+ *     statistics (mean, std dev.) for each of the bins </li>
+ * <li>Divide the interval (0,1) into subintervals associated with the bins,
+ *     with the length of a bin's subinterval proportional to its count.</li></ol>
+ * <strong>Generating random values from the distribution</strong><ol>
+ * <li>Generate a uniformly distributed value in (0,1) </li>
+ * <li>Select the subinterval to which the value belongs.
+ * <li>Generate a random Gaussian value with mean = mean of the associated
+ *     bin and std dev = std dev of associated bin.</li></ol></p><p>
+ *<strong>USAGE NOTES:</strong><ul>
+ *<li>The <code>binCount</code> is set by default to 1000.  A good rule of thumb
+ *    is to set the bin count to approximately the length of the input file divided
+ *    by 10. </li>
+ *<li>The input file <i>must</i> be a plain text file containing one valid numeric
+ *    entry per line.</li>
+ * </ul></p>
+ *
+ * @version $Revision: 1003886 $ $Date: 2010-10-02 23:04:44 +0200 (sam. 02 oct. 2010) $
+ */
+public class EmpiricalDistributionImpl implements Serializable, EmpiricalDistribution {
+
+    /** Serializable version identifier */
+    private static final long serialVersionUID = 5729073523949762654L;
+
+    /** List of SummaryStatistics objects characterizing the bins */
+    private final List<SummaryStatistics> binStats;
+
+    /** Sample statistics */
+    private SummaryStatistics sampleStats = null;
+
+    /** Max loaded value */
+    private double max = Double.NEGATIVE_INFINITY;
+
+    /** Min loaded value */
+    private double min = Double.POSITIVE_INFINITY;
+
+    /** Grid size */
+    private double delta = 0d;
+
+    /** number of bins */
+    private final int binCount;
+
+    /** is the distribution loaded? */
+    private boolean loaded = false;
+
+    /** upper bounds of subintervals in (0,1) "belonging" to the bins */
+    private double[] upperBounds = null;
+
+    /** RandomData instance to use in repeated calls to getNext() */
+    private final RandomData randomData = new RandomDataImpl();
+
+    /**
+     * Creates a new EmpiricalDistribution with the default bin count.
+     */
+    public EmpiricalDistributionImpl() {
+        binCount = 1000;
+        binStats = new ArrayList<SummaryStatistics>();
+    }
+
+    /**
+     * Creates a new EmpiricalDistribution  with the specified bin count.
+     *
+     * @param binCount number of bins
+     */
+    public EmpiricalDistributionImpl(int binCount) {
+        this.binCount = binCount;
+        binStats = new ArrayList<SummaryStatistics>();
+    }
+
+     /**
+     * Computes the empirical distribution from the provided
+     * array of numbers.
+     *
+     * @param in the input data array
+     */
+    public void load(double[] in) {
+        DataAdapter da = new ArrayDataAdapter(in);
+        try {
+            da.computeStats();
+            fillBinStats(in);
+        } catch (IOException e) {
+            throw new MathRuntimeException(e);
+        }
+        loaded = true;
+
+    }
+
+    /**
+     * Computes the empirical distribution using data read from a URL.
+     * @param url  url of the input file
+     *
+     * @throws IOException if an IO error occurs
+     */
+    public void load(URL url) throws IOException {
+        BufferedReader in =
+            new BufferedReader(new InputStreamReader(url.openStream()));
+        try {
+            DataAdapter da = new StreamDataAdapter(in);
+            da.computeStats();
+            if (sampleStats.getN() == 0) {
+                throw MathRuntimeException.createEOFException(LocalizedFormats.URL_CONTAINS_NO_DATA,
+                                                              url);
+            }
+            in = new BufferedReader(new InputStreamReader(url.openStream()));
+            fillBinStats(in);
+            loaded = true;
+        } finally {
+           try {
+               in.close();
+           } catch (IOException ex) {
+               // ignore
+           }
+        }
+    }
+
+    /**
+     * Computes the empirical distribution from the input file.
+     *
+     * @param file the input file
+     * @throws IOException if an IO error occurs
+     */
+    public void load(File file) throws IOException {
+        BufferedReader in = new BufferedReader(new FileReader(file));
+        try {
+            DataAdapter da = new StreamDataAdapter(in);
+            da.computeStats();
+            in = new BufferedReader(new FileReader(file));
+            fillBinStats(in);
+            loaded = true;
+        } finally {
+            try {
+                in.close();
+            } catch (IOException ex) {
+                // ignore
+            }
+        }
+    }
+
+    /**
+     * Provides methods for computing <code>sampleStats</code> and
+     * <code>beanStats</code> abstracting the source of data.
+     */
+    private abstract class DataAdapter{
+
+        /**
+         * Compute bin stats.
+         *
+         * @throws IOException  if an error occurs computing bin stats
+         */
+        public abstract void computeBinStats() throws IOException;
+
+        /**
+         * Compute sample statistics.
+         *
+         * @throws IOException if an error occurs computing sample stats
+         */
+        public abstract void computeStats() throws IOException;
+
+    }
+
+    /**
+     * Factory of <code>DataAdapter</code> objects. For every supported source
+     * of data (array of doubles, file, etc.) an instance of the proper object
+     * is returned.
+     */
+    private class DataAdapterFactory{
+        /**
+         * Creates a DataAdapter from a data object
+         *
+         * @param in object providing access to the data
+         * @return DataAdapter instance
+         */
+        public DataAdapter getAdapter(Object in) {
+            if (in instanceof BufferedReader) {
+                BufferedReader inputStream = (BufferedReader) in;
+                return new StreamDataAdapter(inputStream);
+            } else if (in instanceof double[]) {
+                double[] inputArray = (double[]) in;
+                return new ArrayDataAdapter(inputArray);
+            } else {
+                throw MathRuntimeException.createIllegalArgumentException(
+                      LocalizedFormats.INPUT_DATA_FROM_UNSUPPORTED_DATASOURCE,
+                      in.getClass().getName(),
+                      BufferedReader.class.getName(), double[].class.getName());
+            }
+        }
+    }
+    /**
+     * <code>DataAdapter</code> for data provided through some input stream
+     */
+    private class StreamDataAdapter extends DataAdapter{
+
+        /** Input stream providing access to the data */
+        private BufferedReader inputStream;
+
+        /**
+         * Create a StreamDataAdapter from a BufferedReader
+         *
+         * @param in BufferedReader input stream
+         */
+        public StreamDataAdapter(BufferedReader in){
+            super();
+            inputStream = in;
+        }
+
+        /** {@inheritDoc} */
+        @Override
+        public void computeBinStats() throws IOException {
+            String str = null;
+            double val = 0.0d;
+            while ((str = inputStream.readLine()) != null) {
+                val = Double.parseDouble(str);
+                SummaryStatistics stats = binStats.get(findBin(val));
+                stats.addValue(val);
+            }
+
+            inputStream.close();
+            inputStream = null;
+        }
+
+        /** {@inheritDoc} */
+        @Override
+        public void computeStats() throws IOException {
+            String str = null;
+            double val = 0.0;
+            sampleStats = new SummaryStatistics();
+            while ((str = inputStream.readLine()) != null) {
+                val = Double.valueOf(str).doubleValue();
+                sampleStats.addValue(val);
+            }
+            inputStream.close();
+            inputStream = null;
+        }
+    }
+
+    /**
+     * <code>DataAdapter</code> for data provided as array of doubles.
+     */
+    private class ArrayDataAdapter extends DataAdapter {
+
+        /** Array of input  data values */
+        private double[] inputArray;
+
+        /**
+         * Construct an ArrayDataAdapter from a double[] array
+         *
+         * @param in double[] array holding the data
+         */
+        public ArrayDataAdapter(double[] in){
+            super();
+            inputArray = in;
+        }
+
+        /** {@inheritDoc} */
+        @Override
+        public void computeStats() throws IOException {
+            sampleStats = new SummaryStatistics();
+            for (int i = 0; i < inputArray.length; i++) {
+                sampleStats.addValue(inputArray[i]);
+            }
+        }
+
+        /** {@inheritDoc} */
+        @Override
+        public void computeBinStats() throws IOException {
+            for (int i = 0; i < inputArray.length; i++) {
+                SummaryStatistics stats =
+                    binStats.get(findBin(inputArray[i]));
+                stats.addValue(inputArray[i]);
+            }
+        }
+    }
+
+    /**
+     * Fills binStats array (second pass through data file).
+     *
+     * @param in object providing access to the data
+     * @throws IOException  if an IO error occurs
+     */
+    private void fillBinStats(Object in) throws IOException {
+        // Set up grid
+        min = sampleStats.getMin();
+        max = sampleStats.getMax();
+        delta = (max - min)/(Double.valueOf(binCount)).doubleValue();
+
+        // Initialize binStats ArrayList
+        if (!binStats.isEmpty()) {
+            binStats.clear();
+        }
+        for (int i = 0; i < binCount; i++) {
+            SummaryStatistics stats = new SummaryStatistics();
+            binStats.add(i,stats);
+        }
+
+        // Filling data in binStats Array
+        DataAdapterFactory aFactory = new DataAdapterFactory();
+        DataAdapter da = aFactory.getAdapter(in);
+        da.computeBinStats();
+
+        // Assign upperBounds based on bin counts
+        upperBounds = new double[binCount];
+        upperBounds[0] =
+        ((double) binStats.get(0).getN()) / (double) sampleStats.getN();
+        for (int i = 1; i < binCount-1; i++) {
+            upperBounds[i] = upperBounds[i-1] +
+            ((double) binStats.get(i).getN()) / (double) sampleStats.getN();
+        }
+        upperBounds[binCount-1] = 1.0d;
+    }
+
+    /**
+     * Returns the index of the bin to which the given value belongs
+     *
+     * @param value  the value whose bin we are trying to find
+     * @return the index of the bin containing the value
+     */
+    private int findBin(double value) {
+        return FastMath.min(
+                FastMath.max((int) FastMath.ceil((value- min) / delta) - 1, 0),
+                binCount - 1);
+        }
+
+    /**
+     * Generates a random value from this distribution.
+     *
+     * @return the random value.
+     * @throws IllegalStateException if the distribution has not been loaded
+     */
+    public double getNextValue() throws IllegalStateException {
+
+        if (!loaded) {
+            throw MathRuntimeException.createIllegalStateException(LocalizedFormats.DISTRIBUTION_NOT_LOADED);
+        }
+
+        // Start with a uniformly distributed random number in (0,1)
+        double x = FastMath.random();
+
+        // Use this to select the bin and generate a Gaussian within the bin
+        for (int i = 0; i < binCount; i++) {
+           if (x <= upperBounds[i]) {
+               SummaryStatistics stats = binStats.get(i);
+               if (stats.getN() > 0) {
+                   if (stats.getStandardDeviation() > 0) {  // more than one obs
+                        return randomData.nextGaussian
+                            (stats.getMean(),stats.getStandardDeviation());
+                   } else {
+                       return stats.getMean(); // only one obs in bin
+                   }
+               }
+           }
+        }
+        throw new MathRuntimeException(LocalizedFormats.NO_BIN_SELECTED);
+    }
+
+    /**
+     * Returns a {@link StatisticalSummary} describing this distribution.
+     * <strong>Preconditions:</strong><ul>
+     * <li>the distribution must be loaded before invoking this method</li></ul>
+     *
+     * @return the sample statistics
+     * @throws IllegalStateException if the distribution has not been loaded
+     */
+    public StatisticalSummary getSampleStats() {
+        return sampleStats;
+    }
+
+    /**
+     * Returns the number of bins.
+     *
+     * @return the number of bins.
+     */
+    public int getBinCount() {
+        return binCount;
+    }
+
+    /**
+     * Returns a List of {@link SummaryStatistics} instances containing
+     * statistics describing the values in each of the bins.  The list is
+     * indexed on the bin number.
+     *
+     * @return List of bin statistics.
+     */
+    public List<SummaryStatistics> getBinStats() {
+        return binStats;
+    }
+
+    /**
+     * <p>Returns a fresh copy of the array of upper bounds for the bins.
+     * Bins are: <br/>
+     * [min,upperBounds[0]],(upperBounds[0],upperBounds[1]],...,
+     *  (upperBounds[binCount-2], upperBounds[binCount-1] = max].</p>
+     *
+     * <p>Note: In versions 1.0-2.0 of commons-math, this method
+     * incorrectly returned the array of probability generator upper
+     * bounds now returned by {@link #getGeneratorUpperBounds()}.</p>
+     *
+     * @return array of bin upper bounds
+     * @since 2.1
+     */
+    public double[] getUpperBounds() {
+        double[] binUpperBounds = new double[binCount];
+        binUpperBounds[0] = min + delta;
+        for (int i = 1; i < binCount - 1; i++) {
+            binUpperBounds[i] = binUpperBounds[i-1] + delta;
+        }
+        binUpperBounds[binCount - 1] = max;
+        return binUpperBounds;
+    }
+
+    /**
+     * <p>Returns a fresh copy of the array of upper bounds of the subintervals
+     * of [0,1] used in generating data from the empirical distribution.
+     * Subintervals correspond to bins with lengths proportional to bin counts.</p>
+     *
+     * <p>In versions 1.0-2.0 of commons-math, this array was (incorrectly) returned
+     * by {@link #getUpperBounds()}.</p>
+     *
+     * @since 2.1
+     * @return array of upper bounds of subintervals used in data generation
+     */
+    public double[] getGeneratorUpperBounds() {
+        int len = upperBounds.length;
+        double[] out = new double[len];
+        System.arraycopy(upperBounds, 0, out, 0, len);
+        return out;
+    }
+
+    /**
+     * Property indicating whether or not the distribution has been loaded.
+     *
+     * @return true if the distribution has been loaded
+     */
+    public boolean isLoaded() {
+        return loaded;
+    }
+}
diff --git a/src/main/java/org/apache/commons/math/random/GaussianRandomGenerator.java b/src/main/java/org/apache/commons/math/random/GaussianRandomGenerator.java
new file mode 100644
index 0000000..6bd9775
--- /dev/null
+++ b/src/main/java/org/apache/commons/math/random/GaussianRandomGenerator.java
@@ -0,0 +1,47 @@
+/*
+ * 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.math.random;
+
+/**
+ * This class is a gaussian normalized random generator for scalars.
+ * <p>This class is a simple wrapper around the {@link
+ * RandomGenerator#nextGaussian} method.</p>
+ * @version $Revision: 811827 $ $Date: 2009-09-06 17:32:50 +0200 (dim. 06 sept. 2009) $
+ * @since 1.2
+ */
+
+public class GaussianRandomGenerator implements NormalizedRandomGenerator {
+
+    /** Underlying generator. */
+    private final RandomGenerator generator;
+
+    /** Create a new generator.
+     * @param generator underlying random generator to use
+     */
+    public GaussianRandomGenerator(final RandomGenerator generator) {
+        this.generator = generator;
+    }
+
+    /** Generate a random scalar with null mean and unit standard deviation.
+     * @return a random scalar with null mean and unit standard deviation
+     */
+    public double nextNormalizedDouble() {
+        return generator.nextGaussian();
+    }
+
+}
diff --git a/src/main/java/org/apache/commons/math/random/JDKRandomGenerator.java b/src/main/java/org/apache/commons/math/random/JDKRandomGenerator.java
new file mode 100644
index 0000000..573ef61
--- /dev/null
+++ b/src/main/java/org/apache/commons/math/random/JDKRandomGenerator.java
@@ -0,0 +1,50 @@
+/*
+ * 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.math.random;
+
+import java.util.Random;
+
+/**
+ * Extension of <code>java.util.Random</code> to implement
+ * {@link RandomGenerator}.
+ *
+ * @since 1.1
+ * @version $Revision: 811685 $ $Date: 2009-09-05 19:36:48 +0200 (sam. 05 sept. 2009) $
+ */
+public class JDKRandomGenerator extends Random implements RandomGenerator {
+
+    /** Serializable version identifier. */
+    private static final long serialVersionUID = -7745277476784028798L;
+
+    /** {@inheritDoc} */
+    public void setSeed(int seed) {
+        setSeed((long) seed);
+    }
+
+    /** {@inheritDoc} */
+    public void setSeed(int[] seed) {
+        // the following number is the largest prime that fits in 32 bits (it is 2^32 - 5)
+        final long prime = 4294967291l;
+
+        long combined = 0l;
+        for (int s : seed) {
+            combined = combined * prime + s;
+        }
+        setSeed(combined);
+    }
+
+}
diff --git a/src/main/java/org/apache/commons/math/random/MersenneTwister.java b/src/main/java/org/apache/commons/math/random/MersenneTwister.java
new file mode 100644
index 0000000..6a6fadd
--- /dev/null
+++ b/src/main/java/org/apache/commons/math/random/MersenneTwister.java
@@ -0,0 +1,259 @@
+/*
+ * 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.math.random;
+
+import java.io.Serializable;
+
+import org.apache.commons.math.util.FastMath;
+
+
+/** This class implements a powerful pseudo-random number generator
+ * developed by Makoto Matsumoto and Takuji Nishimura during
+ * 1996-1997.
+
+ * <p>This generator features an extremely long period
+ * (2<sup>19937</sup>-1) and 623-dimensional equidistribution up to 32
+ * bits accuracy. The home page for this generator is located at <a
+ * href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html">
+ * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html</a>.</p>
+
+ * <p>This generator is described in a paper by Makoto Matsumoto and
+ * Takuji Nishimura in 1998: <a
+ * href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/mt.pdf">Mersenne
+ * Twister: A 623-Dimensionally Equidistributed Uniform Pseudo-Random
+ * Number Generator</a>, ACM Transactions on Modeling and Computer
+ * Simulation, Vol. 8, No. 1, January 1998, pp 3--30</p>
+
+ * <p>This class is mainly a Java port of the 2002-01-26 version of
+ * the generator written in C by Makoto Matsumoto and Takuji
+ * Nishimura. Here is their original copyright:</p>
+
+ * <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0">
+ * <tr><td>Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
+ *     All rights reserved.</td></tr>
+
+ * <tr><td>Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * <ol>
+ *   <li>Redistributions of source code must retain the above copyright
+ *       notice, this list of conditions and the following disclaimer.</li>
+ *   <li>Redistributions in binary form must reproduce the above copyright
+ *       notice, this list of conditions and the following disclaimer in the
+ *       documentation and/or other materials provided with the distribution.</li>
+ *   <li>The names of its contributors may not be used to endorse or promote
+ *       products derived from this software without specific prior written
+ *       permission.</li>
+ * </ol></td></tr>
+
+ * <tr><td><strong>THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
+ * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
+ * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+ * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+ * DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS
+ * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
+ * OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
+ * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
+ * USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
+ * DAMAGE.</strong></td></tr>
+ * </table>
+
+ * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $
+ * @since 2.0
+
+ */
+public class MersenneTwister extends BitsStreamGenerator implements Serializable {
+
+    /** Serializable version identifier. */
+    private static final long serialVersionUID = 8661194735290153518L;
+
+    /** Size of the bytes pool. */
+    private static final int   N     = 624;
+
+    /** Period second parameter. */
+    private static final int   M     = 397;
+
+    /** X * MATRIX_A for X = {0, 1}. */
+    private static final int[] MAG01 = { 0x0, 0x9908b0df };
+
+    /** Bytes pool. */
+    private int[] mt;
+
+    /** Current index in the bytes pool. */
+    private int   mti;
+
+    /** Creates a new random number generator.
+     * <p>The instance is initialized using the current time as the
+     * seed.</p>
+     */
+    public MersenneTwister() {
+        mt = new int[N];
+        setSeed(System.currentTimeMillis());
+    }
+
+    /** Creates a new random number generator using a single int seed.
+     * @param seed the initial seed (32 bits integer)
+     */
+    public MersenneTwister(int seed) {
+        mt = new int[N];
+        setSeed(seed);
+    }
+
+    /** Creates a new random number generator using an int array seed.
+     * @param seed the initial seed (32 bits integers array), if null
+     * the seed of the generator will be related to the current time
+     */
+    public MersenneTwister(int[] seed) {
+        mt = new int[N];
+        setSeed(seed);
+    }
+
+    /** Creates a new random number generator using a single long seed.
+     * @param seed the initial seed (64 bits integer)
+     */
+    public MersenneTwister(long seed) {
+        mt = new int[N];
+        setSeed(seed);
+    }
+
+    /** Reinitialize the generator as if just built with the given int seed.
+     * <p>The state of the generator is exactly the same as a new
+     * generator built with the same seed.</p>
+     * @param seed the initial seed (32 bits integer)
+     */
+    @Override
+    public void setSeed(int seed) {
+        // we use a long masked by 0xffffffffL as a poor man unsigned int
+        long longMT = seed;
+        mt[0]= (int) longMT;
+        for (mti = 1; mti < N; ++mti) {
+            // See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier.
+            // initializer from the 2002-01-09 C version by Makoto Matsumoto
+            longMT = (1812433253l * (longMT ^ (longMT >> 30)) + mti) & 0xffffffffL;
+            mt[mti]= (int) longMT;
+        }
+    }
+
+    /** Reinitialize the generator as if just built with the given int array seed.
+     * <p>The state of the generator is exactly the same as a new
+     * generator built with the same seed.</p>
+     * @param seed the initial seed (32 bits integers array), if null
+     * the seed of the generator will be related to the current time
+     */
+    @Override
+    public void setSeed(int[] seed) {
+
+        if (seed == null) {
+            setSeed(System.currentTimeMillis());
+            return;
+        }
+
+        setSeed(19650218);
+        int i = 1;
+        int j = 0;
+
+        for (int k = FastMath.max(N, seed.length); k != 0; k--) {
+            long l0 = (mt[i] & 0x7fffffffl)   | ((mt[i]   < 0) ? 0x80000000l : 0x0l);
+            long l1 = (mt[i-1] & 0x7fffffffl) | ((mt[i-1] < 0) ? 0x80000000l : 0x0l);
+            long l  = (l0 ^ ((l1 ^ (l1 >> 30)) * 1664525l)) + seed[j] + j; // non linear
+            mt[i]   = (int) (l & 0xffffffffl);
+            i++; j++;
+            if (i >= N) {
+                mt[0] = mt[N - 1];
+                i = 1;
+            }
+            if (j >= seed.length) {
+                j = 0;
+            }
+        }
+
+        for (int k = N - 1; k != 0; k--) {
+            long l0 = (mt[i] & 0x7fffffffl)   | ((mt[i]   < 0) ? 0x80000000l : 0x0l);
+            long l1 = (mt[i-1] & 0x7fffffffl) | ((mt[i-1] < 0) ? 0x80000000l : 0x0l);
+            long l  = (l0 ^ ((l1 ^ (l1 >> 30)) * 1566083941l)) - i; // non linear
+            mt[i]   = (int) (l & 0xffffffffL);
+            i++;
+            if (i >= N) {
+                mt[0] = mt[N - 1];
+                i = 1;
+            }
+        }
+
+        mt[0] = 0x80000000; // MSB is 1; assuring non-zero initial array
+
+    }
+
+    /** Reinitialize the generator as if just built with the given long seed.
+     * <p>The state of the generator is exactly the same as a new
+     * generator built with the same seed.</p>
+     * @param seed the initial seed (64 bits integer)
+     */
+    @Override
+    public void setSeed(long seed) {
+        setSeed(new int[] { (int) (seed >>> 32), (int) (seed & 0xffffffffl) });
+    }
+
+    /** Generate next pseudorandom number.
+     * <p>This method is the core generation algorithm. It is used by all the
+     * public generation methods for the various primitive types {@link
+     * #nextBoolean()}, {@link #nextBytes(byte[])}, {@link #nextDouble()},
+     * {@link #nextFloat()}, {@link #nextGaussian()}, {@link #nextInt()},
+     * {@link #next(int)} and {@link #nextLong()}.</p>
+     * @param bits number of random bits to produce
+     * @return random bits generated
+     */
+    @Override
+    protected int next(int bits) {
+
+        int y;
+
+        if (mti >= N) { // generate N words at one time
+            int mtNext = mt[0];
+            for (int k = 0; k < N - M; ++k) {
+                int mtCurr = mtNext;
+                mtNext = mt[k + 1];
+                y = (mtCurr & 0x80000000) | (mtNext & 0x7fffffff);
+                mt[k] = mt[k + M] ^ (y >>> 1) ^ MAG01[y & 0x1];
+            }
+            for (int k = N - M; k < N - 1; ++k) {
+                int mtCurr = mtNext;
+                mtNext = mt[k + 1];
+                y = (mtCurr & 0x80000000) | (mtNext & 0x7fffffff);
+                mt[k] = mt[k + (M - N)] ^ (y >>> 1) ^ MAG01[y & 0x1];
+            }
+            y = (mtNext & 0x80000000) | (mt[0] & 0x7fffffff);
+            mt[N - 1] = mt[M - 1] ^ (y >>> 1) ^ MAG01[y & 0x1];
+
+            mti = 0;
+        }
+
+        y = mt[mti++];
+
+        // tempering
+        y ^=  y >>> 11;
+        y ^= (y <<   7) & 0x9d2c5680;
+        y ^= (y <<  15) & 0xefc60000;
+        y ^=  y >>> 18;
+
+        return y >>> (32 - bits);
+
+    }
+
+}
diff --git a/src/main/java/org/apache/commons/math/random/NormalizedRandomGenerator.java b/src/main/java/org/apache/commons/math/random/NormalizedRandomGenerator.java
new file mode 100644
index 0000000..0f13b81
--- /dev/null
+++ b/src/main/java/org/apache/commons/math/random/NormalizedRandomGenerator.java
@@ -0,0 +1,38 @@
+/*
+ * 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.math.random;
+
+/**
+ * This interface represent a normalized random generator for
+ * scalars.
+ * Normalized generator provide null mean and unit standard deviation scalars.
+ * @version $Revision: 811786 $ $Date: 2009-09-06 11:36:08 +0200 (dim. 06 sept. 2009) $
+ * @since 1.2
+ */
+public interface NormalizedRandomGenerator {
+
+  /** Generate a random scalar with null mean and unit standard deviation.
+   * <p>This method does <strong>not</strong> specify the shape of the
+   * distribution, it is the implementing class that provides it. The
+   * only contract here is to generate numbers with null mean and unit
+   * standard deviation.</p>
+   * @return a random scalar with null mean and unit standard deviation
+   */
+  double nextNormalizedDouble();
+
+}
diff --git a/src/main/java/org/apache/commons/math/random/RandomAdaptor.java b/src/main/java/org/apache/commons/math/random/RandomAdaptor.java
new file mode 100644
index 0000000..8091a48
--- /dev/null
+++ b/src/main/java/org/apache/commons/math/random/RandomAdaptor.java
@@ -0,0 +1,198 @@
+/*
+ * 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.math.random;
+
+import java.util.Random;
+
+/**
+ * Extension of <code>java.util.Random</code> wrapping a
+ * {@link RandomGenerator}.
+ *
+ * @since 1.1
+ * @version $Revision: 1003886 $ $Date: 2010-10-02 23:04:44 +0200 (sam. 02 oct. 2010) $
+ */
+public class RandomAdaptor extends Random implements RandomGenerator {
+
+    /** Serializable version identifier. */
+    private static final long serialVersionUID = 2306581345647615033L;
+
+    /** Wrapped randomGenerator instance */
+    private final RandomGenerator randomGenerator;
+
+    /**
+     * Prevent instantiation without a generator argument
+     */
+    @SuppressWarnings("unused")
+    private RandomAdaptor() { randomGenerator = null; }
+
+    /**
+     * Construct a RandomAdaptor wrapping the supplied RandomGenerator.
+     *
+     * @param randomGenerator  the wrapped generator
+     */
+    public RandomAdaptor(RandomGenerator randomGenerator) {
+        this.randomGenerator = randomGenerator;
+    }
+
+    /**
+     * Factory method to create a <code>Random</code> using the supplied
+     * <code>RandomGenerator</code>.
+     *
+     * @param randomGenerator  wrapped RandomGenerator instance
+     * @return a Random instance wrapping the RandomGenerator
+     */
+    public static Random createAdaptor(RandomGenerator randomGenerator) {
+        return new RandomAdaptor(randomGenerator);
+    }
+
+    /**
+     * Returns the next pseudorandom, uniformly distributed
+     * <code>boolean</code> value from this random number generator's
+     * sequence.
+     *
+     * @return  the next pseudorandom, uniformly distributed
+     * <code>boolean</code> value from this random number generator's
+     * sequence
+     */
+    @Override
+    public boolean nextBoolean() {
+        return randomGenerator.nextBoolean();
+    }
+
+     /**
+     * Generates random bytes and places them into a user-supplied
+     * byte array.  The number of random bytes produced is equal to
+     * the length of the byte array.
+     *
+     * @param bytes the non-null byte array in which to put the
+     * random bytes
+     */
+    @Override
+    public void nextBytes(byte[] bytes) {
+        randomGenerator.nextBytes(bytes);
+    }
+
+     /**
+     * Returns the next pseudorandom, uniformly distributed
+     * <code>double</code> value between <code>0.0</code> and
+     * <code>1.0</code> from this random number generator's sequence.
+     *
+     * @return  the next pseudorandom, uniformly distributed
+     *  <code>double</code> value between <code>0.0</code> and
+     *  <code>1.0</code> from this random number generator's sequence
+     */
+    @Override
+    public double nextDouble() {
+        return randomGenerator.nextDouble();
+    }
+
+    /**
+     * Returns the next pseudorandom, uniformly distributed <code>float</code>
+     * value between <code>0.0</code> and <code>1.0</code> from this random
+     * number generator's sequence.
+     *
+     * @return  the next pseudorandom, uniformly distributed <code>float</code>
+     * value between <code>0.0</code> and <code>1.0</code> from this
+     * random number generator's sequence
+     */
+    @Override
+    public float nextFloat() {
+        return randomGenerator.nextFloat();
+    }
+
+    /**
+     * Returns the next pseudorandom, Gaussian ("normally") distributed
+     * <code>double</code> value with mean <code>0.0</code> and standard
+     * deviation <code>1.0</code> from this random number generator's sequence.
+     *
+     * @return  the next pseudorandom, Gaussian ("normally") distributed
+     * <code>double</code> value with mean <code>0.0</code> and
+     * standard deviation <code>1.0</code> from this random number
+     *  generator's sequence
+     */
+    @Override
+    public double nextGaussian() {
+        return randomGenerator.nextGaussian();
+    }
+
+     /**
+     * Returns the next pseudorandom, uniformly distributed <code>int</code>
+     * value from this random number generator's sequence.
+     * All 2<font size="-1"><sup>32</sup></font> possible <tt>int</tt> values
+     * should be produced with  (approximately) equal probability.
+     *
+     * @return the next pseudorandom, uniformly distributed <code>int</code>
+     *  value from this random number generator's sequence
+     */
+    @Override
+    public int nextInt() {
+        return randomGenerator.nextInt();
+    }
+
+    /**
+     * Returns a pseudorandom, uniformly distributed <tt>int</tt> value
+     * between 0 (inclusive) and the specified value (exclusive), drawn from
+     * this random number generator's sequence.
+     *
+     * @param n the bound on the random number to be returned.  Must be
+     * positive.
+     * @return  a pseudorandom, uniformly distributed <tt>int</tt>
+     * value between 0 (inclusive) and n (exclusive).
+     * @throws IllegalArgumentException  if n is not positive.
+     */
+    @Override
+    public int nextInt(int n) {
+        return randomGenerator.nextInt(n);
+    }
+
+    /**
+     * Returns the next pseudorandom, uniformly distributed <code>long</code>
+     * value from this random number generator's sequence.  All
+     * 2<font size="-1"><sup>64</sup></font> possible <tt>long</tt> values
+     * should be produced with (approximately) equal probability.
+     *
+     * @return  the next pseudorandom, uniformly distributed <code>long</code>
+     *value from this random number generator's sequence
+     */
+    @Override
+    public long nextLong() {
+        return randomGenerator.nextLong();
+    }
+
+    /** {@inheritDoc} */
+    public void setSeed(int seed) {
+        if (randomGenerator != null) {  // required to avoid NPE in constructor
+            randomGenerator.setSeed(seed);
+        }
+    }
+
+    /** {@inheritDoc} */
+    public void setSeed(int[] seed) {
+        if (randomGenerator != null) {  // required to avoid NPE in constructor
+            randomGenerator.setSeed(seed);
+        }
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    public void setSeed(long seed) {
+        if (randomGenerator != null) {  // required to avoid NPE in constructor
+            randomGenerator.setSeed(seed);
+        }
+    }
+
+}
diff --git a/src/main/java/org/apache/commons/math/random/RandomData.java b/src/main/java/org/apache/commons/math/random/RandomData.java
new file mode 100644
index 0000000..0fc5136
--- /dev/null
+++ b/src/main/java/org/apache/commons/math/random/RandomData.java
@@ -0,0 +1,272 @@
+/*
+ * 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.math.random;
+import java.util.Collection;
+
+/**
+ * Random data generation utilities.
+ * @version $Revision: 780975 $ $Date: 2009-06-02 11:05:37 +0200 (mar. 02 juin 2009) $
+ */
+public interface RandomData {
+    /**
+     * Generates a random string of hex characters of length
+     * <code>len</code>.
+     * <p>
+     * The generated string will be random, but not cryptographically
+     * secure. To generate cryptographically secure strings, use
+     * <code>nextSecureHexString</code></p>
+     * <p>
+     * <strong>Preconditions</strong>:<ul>
+     * <li><code>len > 0</code> (otherwise an IllegalArgumentException
+     *     is thrown.)</li>
+     * </ul></p>
+     *
+     * @param len the length of the string to be generated
+     * @return random string of hex characters of length <code>len</code>
+     */
+    String nextHexString(int len);
+
+    /**
+     * Generates a uniformly distributed random integer between
+     * <code>lower</code> and <code>upper</code> (endpoints included).
+     * <p>
+     * The generated integer will be random, but not cryptographically secure.
+     * To generate cryptographically secure integer sequences, use
+     * <code>nextSecureInt</code>.</p>
+     * <p>
+     * <strong>Preconditions</strong>:<ul>
+     * <li><code>lower < upper</code> (otherwise an IllegalArgumentException
+     *     is thrown.)</li>
+     * </ul></p>
+     *
+     * @param lower lower bound for generated integer
+     * @param upper upper bound for generated integer
+     * @return a random integer greater than or equal to <code>lower</code>
+     * and less than or equal to <code>upper</code>.
+     */
+    int nextInt(int lower, int upper);
+
+    /**
+     * Generates a uniformly distributed random long integer between
+     * <code>lower</code> and <code>upper</code> (endpoints included).
+     * <p>
+     * The generated long integer values will be random, but not
+     * cryptographically secure.
+     * To generate cryptographically secure sequences of longs, use
+     * <code>nextSecureLong</code></p>
+     * <p>
+     * <strong>Preconditions</strong>:<ul>
+     * <li><code>lower < upper</code> (otherwise an IllegalArgumentException
+     *     is thrown.)</li>
+     * </ul></p>
+     *
+     * @param lower lower bound for generated integer
+     * @param upper upper bound for generated integer
+     * @return a random integer greater than or equal to <code>lower</code>
+     * and less than or equal to <code>upper</code>.
+     */
+    long nextLong(long lower, long upper);
+
+    /**
+     * Generates a random string of hex characters from a secure random
+     * sequence.
+     * <p>
+     * If cryptographic security is not required,
+     * use <code>nextHexString()</code>.</p>
+     * <p>
+     * <strong>Preconditions</strong>:<ul>
+     * <li><code>len > 0</code> (otherwise an IllegalArgumentException
+     *     is thrown.)</li>
+     * </ul></p>
+     * @param len length of return string
+     * @return the random hex string
+     */
+    String nextSecureHexString(int len);
+
+    /**
+     * Generates a uniformly distributed random integer between
+     * <code>lower</code> and <code>upper</code> (endpoints included)
+     * from a secure random sequence.
+     * <p>
+     * Sequences of integers generated using this method will be
+     * cryptographically secure. If cryptographic security is not required,
+     * <code>nextInt</code> should be used instead of this method.</p>
+     * <p>
+     * <strong>Definition</strong>:
+     * <a href="http://en.wikipedia.org/wiki/Cryptographically_secure_pseudo-random_number_generator">
+     * Secure Random Sequence</a></p>
+     * <p>
+     * <strong>Preconditions</strong>:<ul>
+     * <li><code>lower < upper</code> (otherwise an IllegalArgumentException
+     *     is thrown.)</li>
+     * </ul></p>
+     *
+     * @param lower lower bound for generated integer
+     * @param upper upper bound for generated integer
+     * @return a random integer greater than or equal to <code>lower</code>
+     * and less than or equal to <code>upper</code>.
+     */
+    int nextSecureInt(int lower, int upper);
+
+    /**
+     * Generates a random long integer between <code>lower</code>
+     * and <code>upper</code> (endpoints included).
+     * <p>
+     * Sequences of long values generated using this method will be
+     * cryptographically secure. If cryptographic security is not required,
+     * <code>nextLong</code> should be used instead of this method.</p>
+     * <p>
+     * <strong>Definition</strong>:
+     * <a href="http://en.wikipedia.org/wiki/Cryptographically_secure_pseudo-random_number_generator">
+     * Secure Random Sequence</a></p>
+     * <p>
+     * <strong>Preconditions</strong>:<ul>
+     * <li><code>lower < upper</code> (otherwise an IllegalArgumentException
+     *     is thrown.)</li>
+     * </ul></p>
+     *
+     * @param lower lower bound for generated integer
+     * @param upper upper bound for generated integer
+     * @return a long integer greater than or equal to <code>lower</code>
+     * and less than or equal to <code>upper</code>.
+     */
+    long nextSecureLong(long lower, long upper);
+
+    /**
+     * Generates a random value from the Poisson distribution with
+     * the given mean.
+     * <p>
+     * <strong>Definition</strong>:
+     * <a href="http://www.itl.nist.gov/div898/handbook/eda/section3/eda366j.htm">
+     * Poisson Distribution</a></p>
+     * <p>
+     * <strong>Preconditions</strong>: <ul>
+     * <li>The specified mean <i>must</i> be positive (otherwise an
+     *     IllegalArgumentException is thrown.)</li>
+     * </ul></p>
+     * @param mean Mean of the distribution
+     * @return poisson deviate with the specified mean
+     */
+    long nextPoisson(double mean);
+
+    /**
+     * Generates a random value from the
+     * Normal (or Gaussian) distribution with the given mean
+     * and standard deviation.
+     * <p>
+     * <strong>Definition</strong>:
+     * <a href="http://www.itl.nist.gov/div898/handbook/eda/section3/eda3661.htm">
+     * Normal Distribution</a></p>
+     * <p>
+     * <strong>Preconditions</strong>: <ul>
+     * <li><code>sigma > 0</code> (otherwise an IllegalArgumentException
+     *     is thrown.)</li>
+     * </ul></p>
+     * @param mu Mean of the distribution
+     * @param sigma Standard deviation of the distribution
+     * @return random value from Gaussian distribution with mean = mu,
+     * standard deviation = sigma
+     */
+    double nextGaussian(double mu, double sigma);
+
+    /**
+     * Generates a random value from the exponential distribution
+     * with expected value = <code>mean</code>.
+     * <p>
+     * <strong>Definition</strong>:
+     * <a href="http://www.itl.nist.gov/div898/handbook/eda/section3/eda3667.htm">
+     * Exponential Distribution</a></p>
+     * <p>
+     * <strong>Preconditions</strong>: <ul>
+     * <li><code>mu >= 0</code> (otherwise an IllegalArgumentException
+     *     is thrown.)</li>
+     * </ul></p>
+     * @param mean Mean of the distribution
+     * @return random value from exponential distribution
+     */
+    double nextExponential(double mean);
+
+    /**
+     * Generates a uniformly distributed random value from the open interval
+     * (<code>lower</code>,<code>upper</code>) (i.e., endpoints excluded).
+     * <p>
+     * <strong>Definition</strong>:
+     * <a href="http://www.itl.nist.gov/div898/handbook/eda/section3/eda3662.htm">
+     * Uniform Distribution</a> <code>lower</code> and
+     * <code>upper - lower</code> are the
+     * <a href = "http://www.itl.nist.gov/div898/handbook/eda/section3/eda364.htm">
+     * location and scale parameters</a>, respectively.</p>
+     * <p>
+     * <strong>Preconditions</strong>:<ul>
+     * <li><code>lower < upper</code> (otherwise an IllegalArgumentException
+     *     is thrown.)</li>
+     * </ul></p>
+     *
+     * @param lower lower endpoint of the interval of support
+     * @param upper upper endpoint of the interval of support
+     * @return uniformly distributed random value between lower
+     * and upper (exclusive)
+     */
+    double nextUniform(double lower, double upper);
+
+    /**
+     * Generates an integer array of length <code>k</code> whose entries
+     * are selected randomly, without repetition, from the integers <code>
+     * 0 through n-1</code> (inclusive).
+     * <p>
+     * Generated arrays represent permutations
+     * of <code>n</code> taken <code>k</code> at a time.</p>
+     * <p>
+     * <strong>Preconditions:</strong><ul>
+     * <li> <code>k <= n</code></li>
+     * <li> <code>n > 0</code> </li>
+     * </ul>
+     * If the preconditions are not met, an IllegalArgumentException is
+     * thrown.</p>
+     *
+     * @param n domain of the permutation
+     * @param k size of the permutation
+     * @return random k-permutation of n
+     */
+    int[] nextPermutation(int n, int k);
+
+    /**
+     * Returns an array of <code>k</code> objects selected randomly
+     * from the Collection <code>c</code>.
+     * <p>
+     * Sampling from <code>c</code>
+     * is without replacement; but if <code>c</code> contains identical
+     * objects, the sample may include repeats.  If all elements of <code>
+     * c</code> are distinct, the resulting object array represents a
+     * <a href="http://rkb.home.cern.ch/rkb/AN16pp/node250.html#SECTION0002500000000000000000">
+     * Simple Random Sample</a> of size
+     * <code>k</code> from the elements of <code>c</code>.</p>
+     * <p>
+     * <strong>Preconditions:</strong><ul>
+     * <li> k must be less than or equal to the size of c </li>
+     * <li> c must not be empty </li>
+     * </ul>
+     * If the preconditions are not met, an IllegalArgumentException is
+     * thrown.</p>
+     *
+     * @param c collection to be sampled
+     * @param k size of the sample
+     * @return random sample of k elements from c
+     */
+    Object[] nextSample(Collection<?> c, int k);
+}
diff --git a/src/main/java/org/apache/commons/math/random/RandomDataImpl.java b/src/main/java/org/apache/commons/math/random/RandomDataImpl.java
new file mode 100644
index 0000000..e9ccab7
--- /dev/null
+++ b/src/main/java/org/apache/commons/math/random/RandomDataImpl.java
@@ -0,0 +1,966 @@
+/*
+ * 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.math.random;
+
+import java.io.Serializable;
+import java.security.MessageDigest;
+import java.security.NoSuchAlgorithmException;
+import java.security.NoSuchProviderException;
+import java.security.SecureRandom;
+import java.util.Collection;
+
+import org.apache.commons.math.MathException;
+import org.apache.commons.math.distribution.BetaDistributionImpl;
+import org.apache.commons.math.distribution.BinomialDistributionImpl;
+import org.apache.commons.math.distribution.CauchyDistributionImpl;
+import org.apache.commons.math.distribution.ChiSquaredDistributionImpl;
+import org.apache.commons.math.distribution.ContinuousDistribution;
+import org.apache.commons.math.distribution.FDistributionImpl;
+import org.apache.commons.math.distribution.GammaDistributionImpl;
+import org.apache.commons.math.distribution.HypergeometricDistributionImpl;
+import org.apache.commons.math.distribution.IntegerDistribution;
+import org.apache.commons.math.distribution.PascalDistributionImpl;
+import org.apache.commons.math.distribution.TDistributionImpl;
+import org.apache.commons.math.distribution.WeibullDistributionImpl;
+import org.apache.commons.math.distribution.ZipfDistributionImpl;
+import org.apache.commons.math.exception.MathInternalError;
+import org.apache.commons.math.exception.NotStrictlyPositiveException;
+import org.apache.commons.math.exception.NumberIsTooLargeException;
+import org.apache.commons.math.exception.util.LocalizedFormats;
+import org.apache.commons.math.util.FastMath;
+import org.apache.commons.math.util.MathUtils;
+
+/**
+ * Implements the {@link RandomData} interface using a {@link RandomGenerator}
+ * instance to generate non-secure data and a {@link java.security.SecureRandom}
+ * instance to provide data for the <code>nextSecureXxx</code> methods. If no
+ * <code>RandomGenerator</code> is provided in the constructor, the default is
+ * to use a generator based on {@link java.util.Random}. To plug in a different
+ * implementation, either implement <code>RandomGenerator</code> directly or
+ * extend {@link AbstractRandomGenerator}.
+ * <p>
+ * Supports reseeding the underlying pseudo-random number generator (PRNG). The
+ * <code>SecurityProvider</code> and <code>Algorithm</code> used by the
+ * <code>SecureRandom</code> instance can also be reset.
+ * </p>
+ * <p>
+ * For details on the default PRNGs, see {@link java.util.Random} and
+ * {@link java.security.SecureRandom}.
+ * </p>
+ * <p>
+ * <strong>Usage Notes</strong>:
+ * <ul>
+ * <li>
+ * Instance variables are used to maintain <code>RandomGenerator</code> and
+ * <code>SecureRandom</code> instances used in data generation. Therefore, to
+ * generate a random sequence of values or strings, you should use just
+ * <strong>one</strong> <code>RandomDataImpl</code> instance repeatedly.</li>
+ * <li>
+ * The "secure" methods are *much* slower. These should be used only when a
+ * cryptographically secure random sequence is required. A secure random
+ * sequence is a sequence of pseudo-random values which, in addition to being
+ * well-dispersed (so no subsequence of values is an any more likely than other
+ * subsequence of the the same length), also has the additional property that
+ * knowledge of values generated up to any point in the sequence does not make
+ * it any easier to predict subsequent values.</li>
+ * <li>
+ * When a new <code>RandomDataImpl</code> is created, the underlying random
+ * number generators are <strong>not</strong> initialized. If you do not
+ * explicitly seed the default non-secure generator, it is seeded with the
+ * current time in milliseconds on first use. The same holds for the secure
+ * generator. If you provide a <code>RandomGenerator</code> to the constructor,
+ * however, this generator is not reseeded by the constructor nor is it reseeded
+ * on first use.</li>
+ * <li>
+ * The <code>reSeed</code> and <code>reSeedSecure</code> methods delegate to the
+ * corresponding methods on the underlying <code>RandomGenerator</code> and
+ * <code>SecureRandom</code> instances. Therefore, <code>reSeed(long)</code>
+ * fully resets the initial state of the non-secure random number generator (so
+ * that reseeding with a specific value always results in the same subsequent
+ * random sequence); whereas reSeedSecure(long) does <strong>not</strong>
+ * reinitialize the secure random number generator (so secure sequences started
+ * with calls to reseedSecure(long) won't be identical).</li>
+ * <li>
+ * This implementation is not synchronized.
+ * </ul>
+ * </p>
+ *
+ * @version $Revision: 1061496 $ $Date: 2011-01-20 21:32:16 +0100 (jeu. 20 janv. 2011) $
+ */
+public class RandomDataImpl implements RandomData, Serializable {
+
+    /** Serializable version identifier */
+    private static final long serialVersionUID = -626730818244969716L;
+
+    /** underlying random number generator */
+    private RandomGenerator rand = null;
+
+    /** underlying secure random number generator */
+    private SecureRandom secRand = null;
+
+    /**
+     * Construct a RandomDataImpl.
+     */
+    public RandomDataImpl() {
+    }
+
+    /**
+     * Construct a RandomDataImpl using the supplied {@link RandomGenerator} as
+     * the source of (non-secure) random data.
+     *
+     * @param rand
+     *            the source of (non-secure) random data
+     * @since 1.1
+     */
+    public RandomDataImpl(RandomGenerator rand) {
+        super();
+        this.rand = rand;
+    }
+
+    /**
+     * {@inheritDoc}
+     * <p>
+     * <strong>Algorithm Description:</strong> hex strings are generated using a
+     * 2-step process.
+     * <ol>
+     * <li>
+     * len/2+1 binary bytes are generated using the underlying Random</li>
+     * <li>
+     * Each binary byte is translated into 2 hex digits</li>
+     * </ol>
+     * </p>
+     *
+     * @param len
+     *            the desired string length.
+     * @return the random string.
+     * @throws NotStrictlyPositiveException if {@code len <= 0}.
+     */
+    public String nextHexString(int len) {
+        if (len <= 0) {
+            throw new NotStrictlyPositiveException(LocalizedFormats.LENGTH, len);
+        }
+
+        // Get a random number generator
+        RandomGenerator ran = getRan();
+
+        // Initialize output buffer
+        StringBuilder outBuffer = new StringBuilder();
+
+        // Get int(len/2)+1 random bytes
+        byte[] randomBytes = new byte[(len / 2) + 1];
+        ran.nextBytes(randomBytes);
+
+        // Convert each byte to 2 hex digits
+        for (int i = 0; i < randomBytes.length; i++) {
+            Integer c = Integer.valueOf(randomBytes[i]);
+
+            /*
+             * Add 128 to byte value to make interval 0-255 before doing hex
+             * conversion. This guarantees <= 2 hex digits from toHexString()
+             * toHexString would otherwise add 2^32 to negative arguments.
+             */
+            String hex = Integer.toHexString(c.intValue() + 128);
+
+            // Make sure we add 2 hex digits for each byte
+            if (hex.length() == 1) {
+                hex = "0" + hex;
+            }
+            outBuffer.append(hex);
+        }
+        return outBuffer.toString().substring(0, len);
+    }
+
+    /**
+     * Generate a random int value uniformly distributed between
+     * <code>lower</code> and <code>upper</code>, inclusive.
+     *
+     * @param lower
+     *            the lower bound.
+     * @param upper
+     *            the upper bound.
+     * @return the random integer.
+     * @throws NumberIsTooLargeException if {@code lower >= upper}.
+     */
+    public int nextInt(int lower, int upper) {
+        if (lower >= upper) {
+            throw new NumberIsTooLargeException(LocalizedFormats.LOWER_BOUND_NOT_BELOW_UPPER_BOUND,
+                                                lower, upper, false);
+        }
+        double r = getRan().nextDouble();
+        return (int) ((r * upper) + ((1.0 - r) * lower) + r);
+    }
+
+    /**
+     * Generate a random long value uniformly distributed between
+     * <code>lower</code> and <code>upper</code>, inclusive.
+     *
+     * @param lower
+     *            the lower bound.
+     * @param upper
+     *            the upper bound.
+     * @return the random integer.
+     * @throws NumberIsTooLargeException if {@code lower >= upper}.
+     */
+    public long nextLong(long lower, long upper) {
+        if (lower >= upper) {
+            throw new NumberIsTooLargeException(LocalizedFormats.LOWER_BOUND_NOT_BELOW_UPPER_BOUND,
+                                                lower, upper, false);
+        }
+        double r = getRan().nextDouble();
+        return (long) ((r * upper) + ((1.0 - r) * lower) + r);
+    }
+
+    /**
+     * {@inheritDoc}
+     * <p>
+     * <strong>Algorithm Description:</strong> hex strings are generated in
+     * 40-byte segments using a 3-step process.
+     * <ol>
+     * <li>
+     * 20 random bytes are generated using the underlying
+     * <code>SecureRandom</code>.</li>
+     * <li>
+     * SHA-1 hash is applied to yield a 20-byte binary digest.</li>
+     * <li>
+     * Each byte of the binary digest is converted to 2 hex digits.</li>
+     * </ol>
+     * </p>
+     *
+     * @param len
+     *            the length of the generated string
+     * @return the random string
+     * @throws NotStrictlyPositiveException if {@code len <= 0}.
+     */
+    public String nextSecureHexString(int len) {
+        if (len <= 0) {
+            throw new NotStrictlyPositiveException(LocalizedFormats.LENGTH, len);
+        }
+
+        // Get SecureRandom and setup Digest provider
+        SecureRandom secRan = getSecRan();
+        MessageDigest alg = null;
+        try {
+            alg = MessageDigest.getInstance("SHA-1");
+        } catch (NoSuchAlgorithmException ex) {
+            // this should never happen
+            throw new MathInternalError(ex);
+        }
+        alg.reset();
+
+        // Compute number of iterations required (40 bytes each)
+        int numIter = (len / 40) + 1;
+
+        StringBuilder outBuffer = new StringBuilder();
+        for (int iter = 1; iter < numIter + 1; iter++) {
+            byte[] randomBytes = new byte[40];
+            secRan.nextBytes(randomBytes);
+            alg.update(randomBytes);
+
+            // Compute hash -- will create 20-byte binary hash
+            byte hash[] = alg.digest();
+
+            // Loop over the hash, converting each byte to 2 hex digits
+            for (int i = 0; i < hash.length; i++) {
+                Integer c = Integer.valueOf(hash[i]);
+
+                /*
+                 * Add 128 to byte value to make interval 0-255 This guarantees
+                 * <= 2 hex digits from toHexString() toHexString would
+                 * otherwise add 2^32 to negative arguments
+                 */
+                String hex = Integer.toHexString(c.intValue() + 128);
+
+                // Keep strings uniform length -- guarantees 40 bytes
+                if (hex.length() == 1) {
+                    hex = "0" + hex;
+                }
+                outBuffer.append(hex);
+            }
+        }
+        return outBuffer.toString().substring(0, len);
+    }
+
+    /**
+     * Generate a random int value uniformly distributed between
+     * <code>lower</code> and <code>upper</code>, inclusive. This algorithm uses
+     * a secure random number generator.
+     *
+     * @param lower
+     *            the lower bound.
+     * @param upper
+     *            the upper bound.
+     * @return the random integer.
+     * @throws NumberIsTooLargeException if {@code lower >= upper}.
+     */
+    public int nextSecureInt(int lower, int upper) {
+        if (lower >= upper) {
+            throw new NumberIsTooLargeException(LocalizedFormats.LOWER_BOUND_NOT_BELOW_UPPER_BOUND,
+                                                lower, upper, false);
+        }
+        SecureRandom sec = getSecRan();
+        return lower + (int) (sec.nextDouble() * (upper - lower + 1));
+    }
+
+    /**
+     * Generate a random long value uniformly distributed between
+     * <code>lower</code> and <code>upper</code>, inclusive. This algorithm uses
+     * a secure random number generator.
+     *
+     * @param lower
+     *            the lower bound.
+     * @param upper
+     *            the upper bound.
+     * @return the random integer.
+     * @throws NumberIsTooLargeException if {@code lower >= upper}.
+     */
+    public long nextSecureLong(long lower, long upper) {
+        if (lower >= upper) {
+            throw new NumberIsTooLargeException(LocalizedFormats.LOWER_BOUND_NOT_BELOW_UPPER_BOUND,
+                                                lower, upper, false);
+        }
+        SecureRandom sec = getSecRan();
+        return lower + (long) (sec.nextDouble() * (upper - lower + 1));
+    }
+
+    /**
+     * {@inheritDoc}
+     * <p>
+     * <strong>Algorithm Description</strong>:
+     * <ul><li> For small means, uses simulation of a Poisson process
+     * using Uniform deviates, as described
+     * <a href="http://irmi.epfl.ch/cmos/Pmmi/interactive/rng7.htm"> here.</a>
+     * The Poisson process (and hence value returned) is bounded by 1000 * mean.</li>
+     *
+     * <li> For large means, uses the rejection algorithm described in <br/>
+     * Devroye, Luc. (1981).<i>The Computer Generation of Poisson Random Variables</i>
+     * <strong>Computing</strong> vol. 26 pp. 197-207.</li></ul></p>
+     *
+     * @param mean mean of the Poisson distribution.
+     * @return the random Poisson value.
+     * @throws NotStrictlyPositiveException if {@code mean <= 0}.
+     */
+    public long nextPoisson(double mean) {
+        if (mean <= 0) {
+            throw new NotStrictlyPositiveException(LocalizedFormats.MEAN, mean);
+        }
+
+        final RandomGenerator generator = getRan();
+
+        final double pivot = 40.0d;
+        if (mean < pivot) {
+            double p = FastMath.exp(-mean);
+            long n = 0;
+            double r = 1.0d;
+            double rnd = 1.0d;
+
+            while (n < 1000 * mean) {
+                rnd = generator.nextDouble();
+                r = r * rnd;
+                if (r >= p) {
+                    n++;
+                } else {
+                    return n;
+                }
+            }
+            return n;
+        } else {
+            final double lambda = FastMath.floor(mean);
+            final double lambdaFractional = mean - lambda;
+            final double logLambda = FastMath.log(lambda);
+            final double logLambdaFactorial = MathUtils.factorialLog((int) lambda);
+            final long y2 = lambdaFractional < Double.MIN_VALUE ? 0 : nextPoisson(lambdaFractional);
+            final double delta = FastMath.sqrt(lambda * FastMath.log(32 * lambda / FastMath.PI + 1));
+            final double halfDelta = delta / 2;
+            final double twolpd = 2 * lambda + delta;
+            final double a1 = FastMath.sqrt(FastMath.PI * twolpd) * FastMath.exp(1 / 8 * lambda);
+            final double a2 = (twolpd / delta) * FastMath.exp(-delta * (1 + delta) / twolpd);
+            final double aSum = a1 + a2 + 1;
+            final double p1 = a1 / aSum;
+            final double p2 = a2 / aSum;
+            final double c1 = 1 / (8 * lambda);
+
+            double x = 0;
+            double y = 0;
+            double v = 0;
+            int a = 0;
+            double t = 0;
+            double qr = 0;
+            double qa = 0;
+            for (;;) {
+                final double u = nextUniform(0.0, 1);
+                if (u <= p1) {
+                    final double n = nextGaussian(0d, 1d);
+                    x = n * FastMath.sqrt(lambda + halfDelta) - 0.5d;
+                    if (x > delta || x < -lambda) {
+                        continue;
+                    }
+                    y = x < 0 ? FastMath.floor(x) : FastMath.ceil(x);
+                    final double e = nextExponential(1d);
+                    v = -e - (n * n / 2) + c1;
+                } else {
+                    if (u > p1 + p2) {
+                        y = lambda;
+                        break;
+                    } else {
+                        x = delta + (twolpd / delta) * nextExponential(1d);
+                        y = FastMath.ceil(x);
+                        v = -nextExponential(1d) - delta * (x + 1) / twolpd;
+                    }
+                }
+                a = x < 0 ? 1 : 0;
+                t = y * (y + 1) / (2 * lambda);
+                if (v < -t && a == 0) {
+                    y = lambda + y;
+                    break;
+                }
+                qr = t * ((2 * y + 1) / (6 * lambda) - 1);
+                qa = qr - (t * t) / (3 * (lambda + a * (y + 1)));
+                if (v < qa) {
+                    y = lambda + y;
+                    break;
+                }
+                if (v > qr) {
+                    continue;
+                }
+                if (v < y * logLambda - MathUtils.factorialLog((int) (y + lambda)) + logLambdaFactorial) {
+                    y = lambda + y;
+                    break;
+                }
+            }
+            return y2 + (long) y;
+        }
+    }
+
+    /**
+     * Generate a random value from a Normal (a.k.a. Gaussian) distribution with
+     * the given mean, <code>mu</code> and the given standard deviation,
+     * <code>sigma</code>.
+     *
+     * @param mu
+     *            the mean of the distribution
+     * @param sigma
+     *            the standard deviation of the distribution
+     * @return the random Normal value
+     * @throws NotStrictlyPositiveException if {@code sigma <= 0}.
+     */
+    public double nextGaussian(double mu, double sigma) {
+        if (sigma <= 0) {
+            throw new NotStrictlyPositiveException(LocalizedFormats.STANDARD_DEVIATION, sigma);
+        }
+        return sigma * getRan().nextGaussian() + mu;
+    }
+
+    /**
+     * Returns a random value from an Exponential distribution with the given
+     * mean.
+     * <p>
+     * <strong>Algorithm Description</strong>: Uses the <a
+     * href="http://www.jesus.ox.ac.uk/~clifford/a5/chap1/node5.html"> Inversion
+     * Method</a> to generate exponentially distributed random values from
+     * uniform deviates.
+     * </p>
+     *
+     * @param mean the mean of the distribution
+     * @return the random Exponential value
+     * @throws NotStrictlyPositiveException if {@code mean <= 0}.
+     */
+    public double nextExponential(double mean) {
+        if (mean <= 0.0) {
+            throw new NotStrictlyPositiveException(LocalizedFormats.MEAN, mean);
+        }
+        final RandomGenerator generator = getRan();
+        double unif = generator.nextDouble();
+        while (unif == 0.0d) {
+            unif = generator.nextDouble();
+        }
+        return -mean * FastMath.log(unif);
+    }
+
+    /**
+     * {@inheritDoc}
+     * <p>
+     * <strong>Algorithm Description</strong>: scales the output of
+     * Random.nextDouble(), but rejects 0 values (i.e., will generate another
+     * random double if Random.nextDouble() returns 0). This is necessary to
+     * provide a symmetric output interval (both endpoints excluded).
+     * </p>
+     *
+     * @param lower
+     *            the lower bound.
+     * @param upper
+     *            the upper bound.
+     * @return a uniformly distributed random value from the interval (lower,
+     *         upper)
+     * @throws NumberIsTooLargeException if {@code lower >= upper}.
+     */
+    public double nextUniform(double lower, double upper) {
+        if (lower >= upper) {
+            throw new NumberIsTooLargeException(LocalizedFormats.LOWER_BOUND_NOT_BELOW_UPPER_BOUND,
+                                                lower, upper, false);
+        }
+        final RandomGenerator generator = getRan();
+
+        // ensure nextDouble() isn't 0.0
+        double u = generator.nextDouble();
+        while (u <= 0.0) {
+            u = generator.nextDouble();
+        }
+
+        return lower + u * (upper - lower);
+    }
+
+    /**
+     * Generates a random value from the {@link BetaDistributionImpl Beta Distribution}.
+     * This implementation uses {@link #nextInversionDeviate(ContinuousDistribution) inversion}
+     * to generate random values.
+     *
+     * @param alpha first distribution shape parameter
+     * @param beta second distribution shape parameter
+     * @return random value sampled from the beta(alpha, beta) distribution
+     * @throws MathException if an error occurs generating the random value
+     * @since 2.2
+     */
+    public double nextBeta(double alpha, double beta) throws MathException {
+        return nextInversionDeviate(new BetaDistributionImpl(alpha, beta));
+    }
+
+    /**
+     * Generates a random value from the {@link BinomialDistributionImpl Binomial Distribution}.
+     * This implementation uses {@link #nextInversionDeviate(ContinuousDistribution) inversion}
+     * to generate random values.
+     *
+     * @param numberOfTrials number of trials of the Binomial distribution
+     * @param probabilityOfSuccess probability of success of the Binomial distribution
+     * @return random value sampled from the Binomial(numberOfTrials, probabilityOfSuccess) distribution
+     * @throws MathException if an error occurs generating the random value
+     * @since 2.2
+     */
+    public int nextBinomial(int numberOfTrials, double probabilityOfSuccess) throws MathException {
+        return nextInversionDeviate(new BinomialDistributionImpl(numberOfTrials, probabilityOfSuccess));
+    }
+
+    /**
+     * Generates a random value from the {@link CauchyDistributionImpl Cauchy Distribution}.
+     * This implementation uses {@link #nextInversionDeviate(ContinuousDistribution) inversion}
+     * to generate random values.
+     *
+     * @param median the median of the Cauchy distribution
+     * @param scale the scale parameter of the Cauchy distribution
+     * @return random value sampled from the Cauchy(median, scale) distribution
+     * @throws MathException if an error occurs generating the random value
+     * @since 2.2
+     */
+    public double nextCauchy(double median, double scale) throws MathException {
+        return nextInversionDeviate(new CauchyDistributionImpl(median, scale));
+    }
+
+    /**
+     * Generates a random value from the {@link ChiSquaredDistributionImpl ChiSquare Distribution}.
+     * This implementation uses {@link #nextInversionDeviate(ContinuousDistribution) inversion}
+     * to generate random values.
+     *
+     * @param df the degrees of freedom of the ChiSquare distribution
+     * @return random value sampled from the ChiSquare(df) distribution
+     * @throws MathException if an error occurs generating the random value
+     * @since 2.2
+     */
+    public double nextChiSquare(double df) throws MathException {
+        return nextInversionDeviate(new ChiSquaredDistributionImpl(df));
+    }
+
+    /**
+     * Generates a random value from the {@link FDistributionImpl F Distribution}.
+     * This implementation uses {@link #nextInversionDeviate(ContinuousDistribution) inversion}
+     * to generate random values.
+     *
+     * @param numeratorDf the numerator degrees of freedom of the F distribution
+     * @param denominatorDf the denominator degrees of freedom of the F distribution
+     * @return random value sampled from the F(numeratorDf, denominatorDf) distribution
+     * @throws MathException if an error occurs generating the random value
+     * @since 2.2
+     */
+    public double nextF(double numeratorDf, double denominatorDf) throws MathException {
+        return nextInversionDeviate(new FDistributionImpl(numeratorDf, denominatorDf));
+    }
+
+    /**
+     * Generates a random value from the {@link GammaDistributionImpl Gamma Distribution}.
+     * This implementation uses {@link #nextInversionDeviate(ContinuousDistribution) inversion}
+     * to generate random values.
+     *
+     * @param shape the median of the Gamma distribution
+     * @param scale the scale parameter of the Gamma distribution
+     * @return random value sampled from the Gamma(shape, scale) distribution
+     * @throws MathException if an error occurs generating the random value
+     * @since 2.2
+     */
+    public double nextGamma(double shape, double scale) throws MathException {
+        return nextInversionDeviate(new GammaDistributionImpl(shape, scale));
+    }
+
+    /**
+     * Generates a random value from the {@link HypergeometricDistributionImpl Hypergeometric Distribution}.
+     * This implementation uses {@link #nextInversionDeviate(IntegerDistribution) inversion}
+     * to generate random values.
+     *
+     * @param populationSize the population size of the Hypergeometric distribution
+     * @param numberOfSuccesses number of successes in the population of the Hypergeometric distribution
+     * @param sampleSize the sample size of the Hypergeometric distribution
+     * @return random value sampled from the Hypergeometric(numberOfSuccesses, sampleSize) distribution
+     * @throws MathException if an error occurs generating the random value
+     * @since 2.2
+     */
+    public int nextHypergeometric(int populationSize, int numberOfSuccesses, int sampleSize) throws MathException {
+        return nextInversionDeviate(new HypergeometricDistributionImpl(populationSize, numberOfSuccesses, sampleSize));
+    }
+
+    /**
+     * Generates a random value from the {@link PascalDistributionImpl Pascal Distribution}.
+     * This implementation uses {@link #nextInversionDeviate(IntegerDistribution) inversion}
+     * to generate random values.
+     *
+     * @param r the number of successes of the Pascal distribution
+     * @param p the probability of success of the Pascal distribution
+     * @return random value sampled from the Pascal(r, p) distribution
+     * @throws MathException if an error occurs generating the random value
+     * @since 2.2
+     */
+    public int nextPascal(int r, double p) throws MathException {
+        return nextInversionDeviate(new PascalDistributionImpl(r, p));
+    }
+
+    /**
+     * Generates a random value from the {@link TDistributionImpl T Distribution}.
+     * This implementation uses {@link #nextInversionDeviate(ContinuousDistribution) inversion}
+     * to generate random values.
+     *
+     * @param df the degrees of freedom of the T distribution
+     * @return random value from the T(df) distribution
+     * @throws MathException if an error occurs generating the random value
+     * @since 2.2
+     */
+    public double nextT(double df) throws MathException {
+        return nextInversionDeviate(new TDistributionImpl(df));
+    }
+
+    /**
+     * Generates a random value from the {@link WeibullDistributionImpl Weibull Distribution}.
+     * This implementation uses {@link #nextInversionDeviate(ContinuousDistribution) inversion}
+     * to generate random values.
+     *
+     * @param shape the shape parameter of the Weibull distribution
+     * @param scale the scale parameter of the Weibull distribution
+     * @return random value sampled from the Weibull(shape, size) distribution
+     * @throws MathException if an error occurs generating the random value
+     * @since 2.2
+     */
+    public double nextWeibull(double shape, double scale) throws MathException {
+        return nextInversionDeviate(new WeibullDistributionImpl(shape, scale));
+    }
+
+    /**
+     * Generates a random value from the {@link ZipfDistributionImpl Zipf Distribution}.
+     * This implementation uses {@link #nextInversionDeviate(IntegerDistribution) inversion}
+     * to generate random values.
+     *
+     * @param numberOfElements the number of elements of the ZipfDistribution
+     * @param exponent the exponent of the ZipfDistribution
+     * @return random value sampled from the Zipf(numberOfElements, exponent) distribution
+     * @throws MathException if an error occurs generating the random value
+     * @since 2.2
+     */
+    public int nextZipf(int numberOfElements, double exponent) throws MathException {
+        return nextInversionDeviate(new ZipfDistributionImpl(numberOfElements, exponent));
+    }
+
+    /**
+     * Returns the RandomGenerator used to generate non-secure random data.
+     * <p>
+     * Creates and initializes a default generator if null.
+     * </p>
+     *
+     * @return the Random used to generate random data
+     * @since 1.1
+     */
+    private RandomGenerator getRan() {
+        if (rand == null) {
+            rand = new JDKRandomGenerator();
+            rand.setSeed(System.currentTimeMillis());
+        }
+        return rand;
+    }
+
+    /**
+     * Returns the SecureRandom used to generate secure random data.
+     * <p>
+     * Creates and initializes if null.
+     * </p>
+     *
+     * @return the SecureRandom used to generate secure random data
+     */
+    private SecureRandom getSecRan() {
+        if (secRand == null) {
+            secRand = new SecureRandom();
+            secRand.setSeed(System.currentTimeMillis());
+        }
+        return secRand;
+    }
+
+    /**
+     * Reseeds the random number generator with the supplied seed.
+     * <p>
+     * Will create and initialize if null.
+     * </p>
+     *
+     * @param seed
+     *            the seed value to use
+     */
+    public void reSeed(long seed) {
+        if (rand == null) {
+            rand = new JDKRandomGenerator();
+        }
+        rand.setSeed(seed);
+    }
+
+    /**
+     * Reseeds the secure random number generator with the current time in
+     * milliseconds.
+     * <p>
+     * Will create and initialize if null.
+     * </p>
+     */
+    public void reSeedSecure() {
+        if (secRand == null) {
+            secRand = new SecureRandom();
+        }
+        secRand.setSeed(System.currentTimeMillis());
+    }
+
+    /**
+     * Reseeds the secure random number generator with the supplied seed.
+     * <p>
+     * Will create and initialize if null.
+     * </p>
+     *
+     * @param seed
+     *            the seed value to use
+     */
+    public void reSeedSecure(long seed) {
+        if (secRand == null) {
+            secRand = new SecureRandom();
+        }
+        secRand.setSeed(seed);
+    }
+
+    /**
+     * Reseeds the random number generator with the current time in
+     * milliseconds.
+     */
+    public void reSeed() {
+        if (rand == null) {
+            rand = new JDKRandomGenerator();
+        }
+        rand.setSeed(System.currentTimeMillis());
+    }
+
+    /**
+     * Sets the PRNG algorithm for the underlying SecureRandom instance using
+     * the Security Provider API. The Security Provider API is defined in <a
+     * href =
+     * "http://java.sun.com/j2se/1.3/docs/guide/security/CryptoSpec.html#AppA">
+     * Java Cryptography Architecture API Specification & Reference.</a>
+     * <p>
+     * <strong>USAGE NOTE:</strong> This method carries <i>significant</i>
+     * overhead and may take several seconds to execute.
+     * </p>
+     *
+     * @param algorithm
+     *            the name of the PRNG algorithm
+     * @param provider
+     *            the name of the provider
+     * @throws NoSuchAlgorithmException
+     *             if the specified algorithm is not available
+     * @throws NoSuchProviderException
+     *             if the specified provider is not installed
+     */
+    public void setSecureAlgorithm(String algorithm, String provider)
+            throws NoSuchAlgorithmException, NoSuchProviderException {
+        secRand = SecureRandom.getInstance(algorithm, provider);
+    }
+
+    /**
+     * Generates an integer array of length <code>k</code> whose entries are
+     * selected randomly, without repetition, from the integers
+     * <code>0 through n-1</code> (inclusive).
+     * <p>
+     * Generated arrays represent permutations of <code>n</code> taken
+     * <code>k</code> at a time.
+     * </p>
+     * <p>
+     * <strong>Preconditions:</strong>
+     * <ul>
+     * <li> <code>k <= n</code></li>
+     * <li> <code>n > 0</code></li>
+     * </ul>
+     * If the preconditions are not met, an IllegalArgumentException is thrown.
+     * </p>
+     * <p>
+     * Uses a 2-cycle permutation shuffle. The shuffling process is described <a
+     * href="http://www.maths.abdn.ac.uk/~igc/tch/mx4002/notes/node83.html">
+     * here</a>.
+     * </p>
+     *
+     * @param n
+     *            domain of the permutation (must be positive)
+     * @param k
+     *            size of the permutation (must satisfy 0 < k <= n).
+     * @return the random permutation as an int array
+     * @throws NumberIsTooLargeException if {@code k > n}.
+     * @throws NotStrictlyPositiveException if {@code k <= 0}.
+     */
+    public int[] nextPermutation(int n, int k) {
+        if (k > n) {
+            throw new NumberIsTooLargeException(LocalizedFormats.PERMUTATION_EXCEEDS_N,
+                                                k, n, true);
+        }
+        if (k == 0) {
+            throw new NotStrictlyPositiveException(LocalizedFormats.PERMUTATION_SIZE,
+                                                   k);
+        }
+
+        int[] index = getNatural(n);
+        shuffle(index, n - k);
+        int[] result = new int[k];
+        for (int i = 0; i < k; i++) {
+            result[i] = index[n - i - 1];
+        }
+
+        return result;
+    }
+
+    /**
+     * Uses a 2-cycle permutation shuffle to generate a random permutation.
+     * <strong>Algorithm Description</strong>: Uses a 2-cycle permutation
+     * shuffle to generate a random permutation of <code>c.size()</code> and
+     * then returns the elements whose indexes correspond to the elements of the
+     * generated permutation. This technique is described, and proven to
+     * generate random samples, <a
+     * href="http://www.maths.abdn.ac.uk/~igc/tch/mx4002/notes/node83.html">
+     * here</a>
+     *
+     * @param c
+     *            Collection to sample from.
+     * @param k
+     *            sample size.
+     * @return the random sample.
+     * @throws NumberIsTooLargeException if {@code k > c.size()}.
+     * @throws NotStrictlyPositiveException if {@code k <= 0}.
+     */
+    public Object[] nextSample(Collection<?> c, int k) {
+        int len = c.size();
+        if (k > len) {
+            throw new NumberIsTooLargeException(LocalizedFormats.SAMPLE_SIZE_EXCEEDS_COLLECTION_SIZE,
+                                                k, len, true);
+        }
+        if (k <= 0) {
+            throw new NotStrictlyPositiveException(LocalizedFormats.NUMBER_OF_SAMPLES, k);
+        }
+
+        Object[] objects = c.toArray();
+        int[] index = nextPermutation(len, k);
+        Object[] result = new Object[k];
+        for (int i = 0; i < k; i++) {
+            result[i] = objects[index[i]];
+        }
+        return result;
+    }
+
+    /**
+     * Generate a random deviate from the given distribution using the
+     * <a href="http://en.wikipedia.org/wiki/Inverse_transform_sampling"> inversion method.</a>
+     *
+     * @param distribution Continuous distribution to generate a random value from
+     * @return a random value sampled from the given distribution
+     * @throws MathException if an error occurs computing the inverse cumulative distribution function
+     * @since 2.2
+     */
+    public double nextInversionDeviate(ContinuousDistribution distribution) throws MathException {
+        return distribution.inverseCumulativeProbability(nextUniform(0, 1));
+
+    }
+
+    /**
+     * Generate a random deviate from the given distribution using the
+     * <a href="http://en.wikipedia.org/wiki/Inverse_transform_sampling"> inversion method.</a>
+     *
+     * @param distribution Integer distribution to generate a random value from
+     * @return a random value sampled from the given distribution
+     * @throws MathException if an error occurs computing the inverse cumulative distribution function
+     * @since 2.2
+     */
+    public int nextInversionDeviate(IntegerDistribution distribution) throws MathException {
+        final double target = nextUniform(0, 1);
+        final int glb = distribution.inverseCumulativeProbability(target);
+        if (distribution.cumulativeProbability(glb) == 1.0d) { // No mass above
+            return glb;
+        } else {
+            return glb + 1;
+        }
+    }
+
+    // ------------------------Private methods----------------------------------
+
+    /**
+     * Uses a 2-cycle permutation shuffle to randomly re-order the last elements
+     * of list.
+     *
+     * @param list
+     *            list to be shuffled
+     * @param end
+     *            element past which shuffling begins
+     */
+    private void shuffle(int[] list, int end) {
+        int target = 0;
+        for (int i = list.length - 1; i >= end; i--) {
+            if (i == 0) {
+                target = 0;
+            } else {
+                target = nextInt(0, i);
+            }
+            int temp = list[target];
+            list[target] = list[i];
+            list[i] = temp;
+        }
+    }
+
+    /**
+     * Returns an array representing n.
+     *
+     * @param n
+     *            the natural number to represent
+     * @return array with entries = elements of n
+     */
+    private int[] getNatural(int n) {
+        int[] natural = new int[n];
+        for (int i = 0; i < n; i++) {
+            natural[i] = i;
+        }
+        return natural;
+    }
+
+}
diff --git a/src/main/java/org/apache/commons/math/random/RandomGenerator.java b/src/main/java/org/apache/commons/math/random/RandomGenerator.java
new file mode 100644
index 0000000..0b86c2a
--- /dev/null
+++ b/src/main/java/org/apache/commons/math/random/RandomGenerator.java
@@ -0,0 +1,148 @@
+/*
+ * 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.math.random;
+
+
+/**
+ * Interface extracted from <code>java.util.Random</code>.  This interface is
+ * implemented by {@link AbstractRandomGenerator}.
+ *
+ * @since 1.1
+ * @version $Revision: 949750 $ $Date: 2010-05-31 16:06:04 +0200 (lun. 31 mai 2010) $
+ */
+public interface RandomGenerator {
+
+    /**
+     * Sets the seed of the underlying random number generator using an
+     * <code>int</code> seed.
+     * <p>Sequences of values generated starting with the same seeds
+     * should be identical.
+     * </p>
+     * @param seed the seed value
+     */
+    void setSeed(int seed);
+
+    /**
+     * Sets the seed of the underlying random number generator using an
+     * <code>int</code> array seed.
+     * <p>Sequences of values generated starting with the same seeds
+     * should be identical.
+     * </p>
+     * @param seed the seed value
+     */
+    void setSeed(int[] seed);
+
+    /**
+     * Sets the seed of the underlying random number generator using a
+     * <code>long</code> seed.
+     * <p>Sequences of values generated starting with the same seeds
+     * should be identical.
+     * </p>
+     * @param seed the seed value
+     */
+    void setSeed(long seed);
+
+    /**
+     * Generates random bytes and places them into a user-supplied
+     * byte array.  The number of random bytes produced is equal to
+     * the length of the byte array.
+     *
+     * @param bytes the non-null byte array in which to put the
+     * random bytes
+     */
+    void nextBytes(byte[] bytes);
+
+    /**
+     * Returns the next pseudorandom, uniformly distributed <code>int</code>
+     * value from this random number generator's sequence.
+     * All 2<font size="-1"><sup>32</sup></font> possible <tt>int</tt> values
+     * should be produced with  (approximately) equal probability.
+     *
+     * @return the next pseudorandom, uniformly distributed <code>int</code>
+     *  value from this random number generator's sequence
+     */
+    int nextInt();
+
+    /**
+     * Returns a pseudorandom, uniformly distributed <tt>int</tt> value
+     * between 0 (inclusive) and the specified value (exclusive), drawn from
+     * this random number generator's sequence.
+     *
+     * @param n the bound on the random number to be returned.  Must be
+     * positive.
+     * @return  a pseudorandom, uniformly distributed <tt>int</tt>
+     * value between 0 (inclusive) and n (exclusive).
+     * @throws IllegalArgumentException  if n is not positive.
+     */
+    int nextInt(int n);
+
+    /**
+     * Returns the next pseudorandom, uniformly distributed <code>long</code>
+     * value from this random number generator's sequence.  All
+     * 2<font size="-1"><sup>64</sup></font> possible <tt>long</tt> values
+     * should be produced with (approximately) equal probability.
+     *
+     * @return  the next pseudorandom, uniformly distributed <code>long</code>
+     *value from this random number generator's sequence
+     */
+    long nextLong();
+
+    /**
+     * Returns the next pseudorandom, uniformly distributed
+     * <code>boolean</code> value from this random number generator's
+     * sequence.
+     *
+     * @return  the next pseudorandom, uniformly distributed
+     * <code>boolean</code> value from this random number generator's
+     * sequence
+     */
+    boolean nextBoolean();
+
+    /**
+     * Returns the next pseudorandom, uniformly distributed <code>float</code>
+     * value between <code>0.0</code> and <code>1.0</code> from this random
+     * number generator's sequence.
+     *
+     * @return  the next pseudorandom, uniformly distributed <code>float</code>
+     * value between <code>0.0</code> and <code>1.0</code> from this
+     * random number generator's sequence
+     */
+    float nextFloat();
+
+    /**
+     * Returns the next pseudorandom, uniformly distributed
+     * <code>double</code> value between <code>0.0</code> and
+     * <code>1.0</code> from this random number generator's sequence.
+     *
+     * @return  the next pseudorandom, uniformly distributed
+     *  <code>double</code> value between <code>0.0</code> and
+     *  <code>1.0</code> from this random number generator's sequence
+     */
+    double nextDouble();
+
+    /**
+     * Returns the next pseudorandom, Gaussian ("normally") distributed
+     * <code>double</code> value with mean <code>0.0</code> and standard
+     * deviation <code>1.0</code> from this random number generator's sequence.
+     *
+     * @return  the next pseudorandom, Gaussian ("normally") distributed
+     * <code>double</code> value with mean <code>0.0</code> and
+     * standard deviation <code>1.0</code> from this random number
+     *  generator's sequence
+     */
+    double nextGaussian();
+}
diff --git a/src/main/java/org/apache/commons/math/random/RandomVectorGenerator.java b/src/main/java/org/apache/commons/math/random/RandomVectorGenerator.java
new file mode 100644
index 0000000..15abbd7
--- /dev/null
+++ b/src/main/java/org/apache/commons/math/random/RandomVectorGenerator.java
@@ -0,0 +1,35 @@
+/*
+ * 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.math.random;
+
+
+/** This interface represents a random generator for whole vectors.
+ *
+ * @since 1.2
+ * @version $Revision: 811786 $ $Date: 2009-09-06 11:36:08 +0200 (dim. 06 sept. 2009) $
+ *
+ */
+
+public interface RandomVectorGenerator {
+
+    /** Generate a random vector.
+     * @return a random vector as an array of double.
+     */
+    double[] nextVector();
+
+}
diff --git a/src/main/java/org/apache/commons/math/random/UncorrelatedRandomVectorGenerator.java b/src/main/java/org/apache/commons/math/random/UncorrelatedRandomVectorGenerator.java
new file mode 100644
index 0000000..d365f9b
--- /dev/null
+++ b/src/main/java/org/apache/commons/math/random/UncorrelatedRandomVectorGenerator.java
@@ -0,0 +1,93 @@
+/*
+ * 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.math.random;
+
+import java.util.Arrays;
+
+import org.apache.commons.math.exception.DimensionMismatchException;
+
+/**
+ * A {@link RandomVectorGenerator} that generates vectors with uncorrelated
+ * components. Components of generated vectors follow (independent) Gaussian
+ * distributions, with parameters supplied in the constructor.
+ *
+ * @version $Revision: 962515 $ $Date: 2010-07-09 15:15:28 +0200 (ven. 09 juil. 2010) $
+ * @since 1.2
+ */
+
+public class UncorrelatedRandomVectorGenerator
+  implements RandomVectorGenerator {
+
+    /** Underlying scalar generator. */
+    private final NormalizedRandomGenerator generator;
+
+    /** Mean vector. */
+    private final double[] mean;
+
+    /** Standard deviation vector. */
+    private final double[] standardDeviation;
+
+  /** Simple constructor.
+   * <p>Build an uncorrelated random vector generator from
+   * its mean and standard deviation vectors.</p>
+   * @param mean expected mean values for each component
+   * @param standardDeviation standard deviation for each component
+   * @param generator underlying generator for uncorrelated normalized
+   * components
+   */
+  public UncorrelatedRandomVectorGenerator(double[] mean,
+                                           double[] standardDeviation,
+                                           NormalizedRandomGenerator generator) {
+    if (mean.length != standardDeviation.length) {
+        throw new DimensionMismatchException(mean.length, standardDeviation.length);
+    }
+    this.mean              = mean.clone();
+    this.standardDeviation = standardDeviation.clone();
+    this.generator = generator;
+  }
+
+  /** Simple constructor.
+   * <p>Build a null mean random and unit standard deviation
+   * uncorrelated vector generator</p>
+   * @param dimension dimension of the vectors to generate
+   * @param generator underlying generator for uncorrelated normalized
+   * components
+   */
+  public UncorrelatedRandomVectorGenerator(int dimension,
+                                           NormalizedRandomGenerator generator) {
+    mean              = new double[dimension];
+    standardDeviation = new double[dimension];
+    Arrays.fill(standardDeviation, 1.0);
+    this.generator = generator;
+  }
+
+  /** Generate an uncorrelated random vector.
+   * @return a random vector as a newly built array of double
+   */
+  public double[] nextVector() {
+
+    double[] random = new double[mean.length];
+    for (int i = 0; i < random.length; ++i) {
+      random[i] = mean[i] + standardDeviation[i] * generator.nextNormalizedDouble();
+    }
+
+    return random;
+
+  }
+
+}
diff --git a/src/main/java/org/apache/commons/math/random/UniformRandomGenerator.java b/src/main/java/org/apache/commons/math/random/UniformRandomGenerator.java
new file mode 100644
index 0000000..54492d0
--- /dev/null
+++ b/src/main/java/org/apache/commons/math/random/UniformRandomGenerator.java
@@ -0,0 +1,61 @@
+/*
+ * 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.math.random;
+
+import org.apache.commons.math.util.FastMath;
+
+/**
+ * This class implements a normalized uniform random generator.
+ * <p>Since it is a normalized random generator, it generates values
+ * from a uniform distribution with mean equal to 0 and standard
+ * deviation equal to 1. Generated values fall in the range
+ * [-&#x0221A;3, +&#x0221A;3].</p>
+ *
+ * @since 1.2
+ *
+ * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $
+ */
+
+public class UniformRandomGenerator implements NormalizedRandomGenerator {
+
+    /** Serializable version identifier. */
+    private static final long serialVersionUID = 1569292426375546027L;
+
+    /** Square root of three. */
+    private static final double SQRT3 = FastMath.sqrt(3.0);
+
+    /** Underlying generator. */
+    private final RandomGenerator generator;
+
+    /** Create a new generator.
+     * @param generator underlying random generator to use
+     */
+    public UniformRandomGenerator(RandomGenerator generator) {
+        this.generator = generator;
+    }
+
+    /** Generate a random scalar with null mean and unit standard deviation.
+     * <p>The number generated is uniformly distributed between -&sqrt;(3)
+     * and +&sqrt;(3).</p>
+     * @return a random scalar with null mean and unit standard deviation
+     */
+    public double nextNormalizedDouble() {
+        return SQRT3 * (2 * generator.nextDouble() - 1.0);
+    }
+
+}
diff --git a/src/main/java/org/apache/commons/math/random/UnitSphereRandomVectorGenerator.java b/src/main/java/org/apache/commons/math/random/UnitSphereRandomVectorGenerator.java
new file mode 100644
index 0000000..cac8f18
--- /dev/null
+++ b/src/main/java/org/apache/commons/math/random/UnitSphereRandomVectorGenerator.java
@@ -0,0 +1,84 @@
+/*
+ * 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.math.random;
+
+import org.apache.commons.math.util.FastMath;
+
+
+/**
+ * Generate random vectors isotropically located on the surface of a sphere.
+ *
+ * @since 2.1
+ * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $
+ */
+
+public class UnitSphereRandomVectorGenerator
+    implements RandomVectorGenerator {
+    /**
+     * RNG used for generating the individual components of the vectors.
+     */
+    private final RandomGenerator rand;
+    /**
+     * Space dimension.
+     */
+    private final int dimension;
+
+    /**
+     * @param dimension Space dimension.
+     * @param rand RNG for the individual components of the vectors.
+     */
+    public UnitSphereRandomVectorGenerator(final int dimension,
+                                           final RandomGenerator rand) {
+        this.dimension = dimension;
+        this.rand = rand;
+    }
+    /**
+     * Create an object that will use a default RNG ({@link MersenneTwister}),
+     * in order to generate the individual components.
+     *
+     * @param dimension Space dimension.
+     */
+    public UnitSphereRandomVectorGenerator(final int dimension) {
+        this(dimension, new MersenneTwister());
+    }
+
+    /** {@inheritDoc} */
+    public double[] nextVector() {
+
+        final double[] v = new double[dimension];
+
+        double normSq;
+        do {
+            normSq = 0;
+            for (int i = 0; i < dimension; i++) {
+                final double comp = 2 * rand.nextDouble() - 1;
+                v[i] = comp;
+                normSq += comp * comp;
+            }
+        } while (normSq > 1);
+
+        final double f = 1 / FastMath.sqrt(normSq);
+        for (int i = 0; i < dimension; i++) {
+            v[i] *= f;
+        }
+
+        return v;
+
+    }
+
+}
diff --git a/src/main/java/org/apache/commons/math/random/ValueServer.java b/src/main/java/org/apache/commons/math/random/ValueServer.java
new file mode 100644
index 0000000..9146e69
--- /dev/null
+++ b/src/main/java/org/apache/commons/math/random/ValueServer.java
@@ -0,0 +1,385 @@
+/*
+ * 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.math.random;
+import java.io.BufferedReader;
+import java.io.IOException;
+import java.io.InputStreamReader;
+import java.net.MalformedURLException;
+import java.net.URL;
+
+import org.apache.commons.math.MathRuntimeException;
+import org.apache.commons.math.exception.util.LocalizedFormats;
+
+/**
+ * Generates values for use in simulation applications.
+ * <p>
+ * How values are generated is determined by the <code>mode</code>
+ * property.</p>
+ * <p>
+ * Supported <code>mode</code> values are: <ul>
+ * <li> DIGEST_MODE -- uses an empirical distribution </li>
+ * <li> REPLAY_MODE -- replays data from <code>valuesFileURL</code></li>
+ * <li> UNIFORM_MODE -- generates uniformly distributed random values with
+ *                      mean = <code>mu</code> </li>
+ * <li> EXPONENTIAL_MODE -- generates exponentially distributed random values
+ *                         with mean = <code>mu</code></li>
+ * <li> GAUSSIAN_MODE -- generates Gaussian distributed random values with
+ *                       mean = <code>mu</code> and
+ *                       standard deviation = <code>sigma</code></li>
+ * <li> CONSTANT_MODE -- returns <code>mu</code> every time.</li></ul></p>
+ *
+ * @version $Revision: 1003886 $ $Date: 2010-10-02 23:04:44 +0200 (sam. 02 oct. 2010) $
+ *
+ */
+public class ValueServer {
+
+    /** Use empirical distribution.  */
+    public static final int DIGEST_MODE = 0;
+
+    /** Replay data from valuesFilePath. */
+    public static final int REPLAY_MODE = 1;
+
+    /** Uniform random deviates with mean = &mu;. */
+    public static final int UNIFORM_MODE = 2;
+
+    /** Exponential random deviates with mean = &mu;. */
+    public static final int EXPONENTIAL_MODE = 3;
+
+    /** Gaussian random deviates with mean = &mu;, std dev = &sigma;. */
+    public static final int GAUSSIAN_MODE = 4;
+
+    /** Always return mu */
+    public static final int CONSTANT_MODE = 5;
+
+    /** mode determines how values are generated. */
+    private int mode = 5;
+
+    /** URI to raw data values. */
+    private URL valuesFileURL = null;
+
+    /** Mean for use with non-data-driven modes. */
+    private double mu = 0.0;
+
+    /** Standard deviation for use with GAUSSIAN_MODE. */
+    private double sigma = 0.0;
+
+    /** Empirical probability distribution for use with DIGEST_MODE. */
+    private EmpiricalDistribution empiricalDistribution = null;
+
+    /** File pointer for REPLAY_MODE. */
+    private BufferedReader filePointer = null;
+
+    /** RandomDataImpl to use for random data generation. */
+    private final RandomData randomData;
+
+    // Data generation modes ======================================
+
+    /** Creates new ValueServer */
+    public ValueServer() {
+        randomData = new RandomDataImpl();
+    }
+
+    /**
+     * Construct a ValueServer instance using a RandomData as its source
+     * of random data.
+     *
+     * @param randomData the RandomData instance used to source random data
+     * @since 1.1
+     */
+    public ValueServer(RandomData randomData) {
+        this.randomData = randomData;
+    }
+
+    /**
+     * Returns the next generated value, generated according
+     * to the mode value (see MODE constants).
+     *
+     * @return generated value
+     * @throws IOException in REPLAY_MODE if a file I/O error occurs
+     */
+    public double getNext() throws IOException {
+        switch (mode) {
+            case DIGEST_MODE: return getNextDigest();
+            case REPLAY_MODE: return getNextReplay();
+            case UNIFORM_MODE: return getNextUniform();
+            case EXPONENTIAL_MODE: return getNextExponential();
+            case GAUSSIAN_MODE: return getNextGaussian();
+            case CONSTANT_MODE: return mu;
+            default: throw MathRuntimeException.createIllegalStateException(
+                    LocalizedFormats.UNKNOWN_MODE,
+                    mode,
+                    "DIGEST_MODE",   DIGEST_MODE,   "REPLAY_MODE",      REPLAY_MODE,
+                    "UNIFORM_MODE",  UNIFORM_MODE,  "EXPONENTIAL_MODE", EXPONENTIAL_MODE,
+                    "GAUSSIAN_MODE", GAUSSIAN_MODE, "CONSTANT_MODE",    CONSTANT_MODE);
+        }
+    }
+
+    /**
+     * Fills the input array with values generated using getNext() repeatedly.
+     *
+     * @param values array to be filled
+     * @throws IOException in REPLAY_MODE if a file I/O error occurs
+     */
+    public void fill(double[] values) throws IOException {
+        for (int i = 0; i < values.length; i++) {
+            values[i] = getNext();
+        }
+    }
+
+    /**
+     * Returns an array of length <code>length</code> with values generated
+     * using getNext() repeatedly.
+     *
+     * @param length length of output array
+     * @return array of generated values
+     * @throws IOException in REPLAY_MODE if a file I/O error occurs
+     */
+    public double[] fill(int length) throws IOException {
+        double[] out = new double[length];
+        for (int i = 0; i < length; i++) {
+            out[i] = getNext();
+        }
+        return out;
+    }
+
+    /**
+     * Computes the empirical distribution using values from the file
+     * in <code>valuesFileURL</code>, using the default number of bins.
+     * <p>
+     * <code>valuesFileURL</code> must exist and be
+     * readable by *this at runtime.</p>
+     * <p>
+     * This method must be called before using <code>getNext()</code>
+     * with <code>mode = DIGEST_MODE</code></p>
+     *
+     * @throws IOException if an I/O error occurs reading the input file
+     */
+    public void computeDistribution() throws IOException {
+        empiricalDistribution = new EmpiricalDistributionImpl();
+        empiricalDistribution.load(valuesFileURL);
+    }
+
+    /**
+     * Computes the empirical distribution using values from the file
+     * in <code>valuesFileURL</code> and <code>binCount</code> bins.
+     * <p>
+     * <code>valuesFileURL</code> must exist and be readable by this process
+     * at runtime.</p>
+     * <p>
+     * This method must be called before using <code>getNext()</code>
+     * with <code>mode = DIGEST_MODE</code></p>
+     *
+     * @param binCount the number of bins used in computing the empirical
+     * distribution
+     * @throws IOException if an error occurs reading the input file
+     */
+    public void computeDistribution(int binCount)
+            throws IOException {
+        empiricalDistribution = new EmpiricalDistributionImpl(binCount);
+        empiricalDistribution.load(valuesFileURL);
+        mu = empiricalDistribution.getSampleStats().getMean();
+        sigma = empiricalDistribution.getSampleStats().getStandardDeviation();
+    }
+
+    /** Getter for property mode.
+     * @return Value of property mode.
+     */
+    public int getMode() {
+        return mode;
+    }
+
+    /** Setter for property mode.
+     * @param mode New value of property mode.
+     */
+    public void setMode(int mode) {
+        this.mode = mode;
+    }
+
+    /**
+     * Getter for <code>valuesFileURL<code>
+     * @return Value of property valuesFileURL.
+     */
+    public URL getValuesFileURL() {
+        return valuesFileURL;
+    }
+
+    /**
+     * Sets the <code>valuesFileURL</code> using a string URL representation
+     * @param url String representation for new valuesFileURL.
+     * @throws MalformedURLException if url is not well formed
+     */
+    public void setValuesFileURL(String url) throws MalformedURLException {
+        this.valuesFileURL = new URL(url);
+    }
+
+    /**
+     * Sets the <code>valuesFileURL</code>
+     * @param url New value of property valuesFileURL.
+     */
+    public void setValuesFileURL(URL url) {
+        this.valuesFileURL = url;
+    }
+
+    /** Getter for property empiricalDistribution.
+     * @return Value of property empiricalDistribution.
+     */
+    public EmpiricalDistribution getEmpiricalDistribution() {
+        return empiricalDistribution;
+    }
+
+    /**
+     * Resets REPLAY_MODE file pointer to the beginning of the <code>valuesFileURL</code>.
+     *
+     * @throws IOException if an error occurs opening the file
+     */
+    public void resetReplayFile() throws IOException {
+        if (filePointer != null) {
+            try {
+                filePointer.close();
+                filePointer = null;
+            } catch (IOException ex) {
+                // ignore
+            }
+        }
+        filePointer = new BufferedReader(new InputStreamReader(valuesFileURL.openStream()));
+    }
+
+    /**
+     * Closes <code>valuesFileURL</code> after use in REPLAY_MODE.
+     *
+     * @throws IOException if an error occurs closing the file
+     */
+    public void closeReplayFile() throws IOException {
+        if (filePointer != null) {
+            filePointer.close();
+            filePointer = null;
+        }
+    }
+
+    /** Getter for property mu.
+     * @return Value of property mu.
+     */
+    public double getMu() {
+        return mu;
+    }
+
+    /** Setter for property mu.
+     * @param mu New value of property mu.
+     */
+    public void setMu(double mu) {
+        this.mu = mu;
+    }
+
+    /** Getter for property sigma.
+     * @return Value of property sigma.
+     */
+    public double getSigma() {
+        return sigma;
+    }
+
+    /** Setter for property sigma.
+     * @param sigma New value of property sigma.
+     */
+    public void setSigma(double sigma) {
+        this.sigma = sigma;
+    }
+
+    //------------- private methods ---------------------------------
+
+    /**
+     * Gets a random value in DIGEST_MODE.
+     * <p>
+     * <strong>Preconditions</strong>: <ul>
+     * <li>Before this method is called, <code>computeDistribution()</code>
+     * must have completed successfully; otherwise an
+     * <code>IllegalStateException</code> will be thrown</li></ul></p>
+     *
+     * @return next random value from the empirical distribution digest
+     */
+    private double getNextDigest() {
+        if ((empiricalDistribution == null) ||
+            (empiricalDistribution.getBinStats().size() == 0)) {
+            throw MathRuntimeException.createIllegalStateException(LocalizedFormats.DIGEST_NOT_INITIALIZED);
+        }
+        return empiricalDistribution.getNextValue();
+    }
+
+    /**
+     * Gets next sequential value from the <code>valuesFileURL</code>.
+     * <p>
+     * Throws an IOException if the read fails.</p>
+     * <p>
+     * This method will open the <code>valuesFileURL</code> if there is no
+     * replay file open.</p>
+     * <p>
+     * The <code>valuesFileURL</code> will be closed and reopened to wrap around
+     * from EOF to BOF if EOF is encountered. EOFException (which is a kind of
+     * IOException) may still be thrown if the <code>valuesFileURL</code> is
+     * empty.</p>
+     *
+     * @return next value from the replay file
+     * @throws IOException if there is a problem reading from the file
+     * @throws NumberFormatException if an invalid numeric string is
+     *   encountered in the file
+     */
+    private double getNextReplay() throws IOException {
+        String str = null;
+        if (filePointer == null) {
+            resetReplayFile();
+        }
+        if ((str = filePointer.readLine()) == null) {
+            // we have probably reached end of file, wrap around from EOF to BOF
+            closeReplayFile();
+            resetReplayFile();
+            if ((str = filePointer.readLine()) == null) {
+                throw MathRuntimeException.createEOFException(LocalizedFormats.URL_CONTAINS_NO_DATA,
+                                                              valuesFileURL);
+            }
+        }
+        return Double.valueOf(str).doubleValue();
+    }
+
+    /**
+     * Gets a uniformly distributed random value with mean = mu.
+     *
+     * @return random uniform value
+     */
+    private double getNextUniform() {
+        return randomData.nextUniform(0, 2 * mu);
+    }
+
+    /**
+     * Gets an exponentially distributed random value with mean = mu.
+     *
+     * @return random exponential value
+     */
+    private double getNextExponential() {
+        return randomData.nextExponential(mu);
+    }
+
+    /**
+     * Gets a Gaussian distributed random value with mean = mu
+     * and standard deviation = sigma.
+     *
+     * @return random Gaussian value
+     */
+    private double getNextGaussian() {
+        return randomData.nextGaussian(mu, sigma);
+    }
+
+}
diff --git a/src/main/java/org/apache/commons/math/random/Well1024a.java b/src/main/java/org/apache/commons/math/random/Well1024a.java
new file mode 100644
index 0000000..8406ed5
--- /dev/null
+++ b/src/main/java/org/apache/commons/math/random/Well1024a.java
@@ -0,0 +1,106 @@
+/*
+ * 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.math.random;
+
+
+/** This class implements the WELL1024a pseudo-random number generator
+ * from Fran&ccedil;ois Panneton, Pierre L'Ecuyer and Makoto Matsumoto.
+
+ * <p>This generator is described in a paper by Fran&ccedil;ois Panneton,
+ * Pierre L'Ecuyer and Makoto Matsumoto <a
+ * href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng.pdf">Improved
+ * Long-Period Generators Based on Linear Recurrences Modulo 2</a> ACM
+ * Transactions on Mathematical Software, 32, 1 (2006). The errata for the paper
+ * are in <a href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng-errata.txt">wellrng-errata.txt</a>.</p>
+
+ * @see <a href="http://www.iro.umontreal.ca/~panneton/WELLRNG.html">WELL Random number generator</a>
+ * @version $Revision: 1003892 $ $Date: 2010-10-02 23:28:56 +0200 (sam. 02 oct. 2010) $
+ * @since 2.2
+
+ */
+public class Well1024a extends AbstractWell {
+
+    /** Serializable version identifier. */
+    private static final long serialVersionUID = 5680173464174485492L;
+
+    /** Number of bits in the pool. */
+    private static final int K = 1024;
+
+    /** First parameter of the algorithm. */
+    private static final int M1 = 3;
+
+    /** Second parameter of the algorithm. */
+    private static final int M2 = 24;
+
+    /** Third parameter of the algorithm. */
+    private static final int M3 = 10;
+
+    /** Creates a new random number generator.
+     * <p>The instance is initialized using the current time as the
+     * seed.</p>
+     */
+    public Well1024a() {
+        super(K, M1, M2, M3);
+    }
+
+    /** Creates a new random number generator using a single int seed.
+     * @param seed the initial seed (32 bits integer)
+     */
+    public Well1024a(int seed) {
+        super(K, M1, M2, M3, seed);
+    }
+
+    /** Creates a new random number generator using an int array seed.
+     * @param seed the initial seed (32 bits integers array), if null
+     * the seed of the generator will be related to the current time
+     */
+    public Well1024a(int[] seed) {
+        super(K, M1, M2, M3, seed);
+    }
+
+    /** Creates a new random number generator using a single long seed.
+     * @param seed the initial seed (64 bits integer)
+     */
+    public Well1024a(long seed) {
+        super(K, M1, M2, M3, seed);
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    protected int next(final int bits) {
+
+        final int indexRm1 = iRm1[index];
+
+        final int v0       = v[index];
+        final int vM1      = v[i1[index]];
+        final int vM2      = v[i2[index]];
+        final int vM3      = v[i3[index]];
+
+        final int z0 = v[indexRm1];
+        final int z1 = v0  ^ (vM1 ^ (vM1 >>> 8));
+        final int z2 = (vM2 ^ (vM2 << 19)) ^ (vM3 ^ (vM3 << 14));
+        final int z3 = z1      ^ z2;
+        final int z4 = (z0 ^ (z0 << 11)) ^ (z1 ^ (z1 << 7)) ^ (z2 ^ (z2 << 13));
+
+        v[index]     = z3;
+        v[indexRm1]  = z4;
+        index        = indexRm1;
+
+        return z4 >>> (32 - bits);
+
+    }
+}
diff --git a/src/main/java/org/apache/commons/math/random/Well19937a.java b/src/main/java/org/apache/commons/math/random/Well19937a.java
new file mode 100644
index 0000000..97f9bfd
--- /dev/null
+++ b/src/main/java/org/apache/commons/math/random/Well19937a.java
@@ -0,0 +1,108 @@
+/*
+ * 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.math.random;
+
+
+/** This class implements the WELL19937a pseudo-random number generator
+ * from Fran&ccedil;ois Panneton, Pierre L'Ecuyer and Makoto Matsumoto.
+
+ * <p>This generator is described in a paper by Fran&ccedil;ois Panneton,
+ * Pierre L'Ecuyer and Makoto Matsumoto <a
+ * href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng.pdf">Improved
+ * Long-Period Generators Based on Linear Recurrences Modulo 2</a> ACM
+ * Transactions on Mathematical Software, 32, 1 (2006). The errata for the paper
+ * are in <a href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng-errata.txt">wellrng-errata.txt</a>.</p>
+
+ * @see <a href="http://www.iro.umontreal.ca/~panneton/WELLRNG.html">WELL Random number generator</a>
+ * @version $Revision: 1003892 $ $Date: 2010-10-02 23:28:56 +0200 (sam. 02 oct. 2010) $
+ * @since 2.2
+
+ */
+public class Well19937a extends AbstractWell {
+
+    /** Serializable version identifier. */
+    private static final long serialVersionUID = -7462102162223815419L;
+
+    /** Number of bits in the pool. */
+    private static final int K = 19937;
+
+    /** First parameter of the algorithm. */
+    private static final int M1 = 70;
+
+    /** Second parameter of the algorithm. */
+    private static final int M2 = 179;
+
+    /** Third parameter of the algorithm. */
+    private static final int M3 = 449;
+
+    /** Creates a new random number generator.
+     * <p>The instance is initialized using the current time as the
+     * seed.</p>
+     */
+    public Well19937a() {
+        super(K, M1, M2, M3);
+    }
+
+    /** Creates a new random number generator using a single int seed.
+     * @param seed the initial seed (32 bits integer)
+     */
+    public Well19937a(int seed) {
+        super(K, M1, M2, M3, seed);
+    }
+
+    /** Creates a new random number generator using an int array seed.
+     * @param seed the initial seed (32 bits integers array), if null
+     * the seed of the generator will be related to the current time
+     */
+    public Well19937a(int[] seed) {
+        super(K, M1, M2, M3, seed);
+    }
+
+    /** Creates a new random number generator using a single long seed.
+     * @param seed the initial seed (64 bits integer)
+     */
+    public Well19937a(long seed) {
+        super(K, M1, M2, M3, seed);
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    protected int next(final int bits) {
+
+        final int indexRm1 = iRm1[index];
+        final int indexRm2 = iRm2[index];
+
+        final int v0       = v[index];
+        final int vM1      = v[i1[index]];
+        final int vM2      = v[i2[index]];
+        final int vM3      = v[i3[index]];
+
+        final int z0 = (0x80000000 & v[indexRm1]) ^ (0x7FFFFFFF & v[indexRm2]);
+        final int z1 = (v0 ^ (v0 << 25))  ^ (vM1 ^ (vM1 >>> 27));
+        final int z2 = (vM2 >>> 9) ^ (vM3 ^ (vM3 >>> 1));
+        final int z3 = z1      ^ z2;
+        final int z4 = z0 ^ (z1 ^ (z1 << 9)) ^ (z2 ^ (z2 << 21)) ^ (z3 ^ (z3 >>> 21));
+
+        v[index]     = z3;
+        v[indexRm1]  = z4;
+        v[indexRm2] &= 0x80000000;
+        index        = indexRm1;
+
+        return z4 >>> (32 - bits);
+
+    }
+}
diff --git a/src/main/java/org/apache/commons/math/random/Well19937c.java b/src/main/java/org/apache/commons/math/random/Well19937c.java
new file mode 100644
index 0000000..53029c1
--- /dev/null
+++ b/src/main/java/org/apache/commons/math/random/Well19937c.java
@@ -0,0 +1,115 @@
+/*
+ * 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.math.random;
+
+
+/** This class implements the WELL19937c pseudo-random number generator
+ * from Fran&ccedil;ois Panneton, Pierre L'Ecuyer and Makoto Matsumoto.
+
+ * <p>This generator is described in a paper by Fran&ccedil;ois Panneton,
+ * Pierre L'Ecuyer and Makoto Matsumoto <a
+ * href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng.pdf">Improved
+ * Long-Period Generators Based on Linear Recurrences Modulo 2</a> ACM
+ * Transactions on Mathematical Software, 32, 1 (2006). The errata for the paper
+ * are in <a href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng-errata.txt">wellrng-errata.txt</a>.</p>
+
+ * @see <a href="http://www.iro.umontreal.ca/~panneton/WELLRNG.html">WELL Random number generator</a>
+ * @version $Revision: 1003892 $ $Date: 2010-10-02 23:28:56 +0200 (sam. 02 oct. 2010) $
+ * @since 2.2
+
+ */
+public class Well19937c extends AbstractWell {
+
+    /** Serializable version identifier. */
+    private static final long serialVersionUID = -7203498180754925124L;
+
+    /** Number of bits in the pool. */
+    private static final int K = 19937;
+
+    /** First parameter of the algorithm. */
+    private static final int M1 = 70;
+
+    /** Second parameter of the algorithm. */
+    private static final int M2 = 179;
+
+    /** Third parameter of the algorithm. */
+    private static final int M3 = 449;
+
+    /** Creates a new random number generator.
+     * <p>The instance is initialized using the current time as the
+     * seed.</p>
+     */
+    public Well19937c() {
+        super(K, M1, M2, M3);
+    }
+
+    /** Creates a new random number generator using a single int seed.
+     * @param seed the initial seed (32 bits integer)
+     */
+    public Well19937c(int seed) {
+        super(K, M1, M2, M3, seed);
+    }
+
+    /** Creates a new random number generator using an int array seed.
+     * @param seed the initial seed (32 bits integers array), if null
+     * the seed of the generator will be related to the current time
+     */
+    public Well19937c(int[] seed) {
+        super(K, M1, M2, M3, seed);
+    }
+
+    /** Creates a new random number generator using a single long seed.
+     * @param seed the initial seed (64 bits integer)
+     */
+    public Well19937c(long seed) {
+        super(K, M1, M2, M3, seed);
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    protected int next(final int bits) {
+
+        final int indexRm1 = iRm1[index];
+        final int indexRm2 = iRm2[index];
+
+        final int v0       = v[index];
+        final int vM1      = v[i1[index]];
+        final int vM2      = v[i2[index]];
+        final int vM3      = v[i3[index]];
+
+        final int z0 = (0x80000000 & v[indexRm1]) ^ (0x7FFFFFFF & v[indexRm2]);
+        final int z1 = (v0 ^ (v0 << 25))  ^ (vM1 ^ (vM1 >>> 27));
+        final int z2 = (vM2 >>> 9) ^ (vM3 ^ (vM3 >>> 1));
+        final int z3 = z1      ^ z2;
+        int z4 = z0 ^ (z1 ^ (z1 << 9)) ^ (z2 ^ (z2 << 21)) ^ (z3 ^ (z3 >>> 21));
+
+        v[index]     = z3;
+        v[indexRm1]  = z4;
+        v[indexRm2] &= 0x80000000;
+        index        = indexRm1;
+
+
+        // add Matsumoto-Kurita tempering
+        // to get a maximally-equidistributed generator
+        z4 = z4 ^ ((z4 <<  7) & 0xe46e1700);
+        z4 = z4 ^ ((z4 << 15) & 0x9b868000);
+
+        return z4 >>> (32 - bits);
+
+    }
+
+}
diff --git a/src/main/java/org/apache/commons/math/random/Well44497a.java b/src/main/java/org/apache/commons/math/random/Well44497a.java
new file mode 100644
index 0000000..70d672c
--- /dev/null
+++ b/src/main/java/org/apache/commons/math/random/Well44497a.java
@@ -0,0 +1,111 @@
+/*
+ * 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.math.random;
+
+
+/** This class implements the WELL44497a pseudo-random number generator
+ * from Fran&ccedil;ois Panneton, Pierre L'Ecuyer and Makoto Matsumoto.
+
+ * <p>This generator is described in a paper by Fran&ccedil;ois Panneton,
+ * Pierre L'Ecuyer and Makoto Matsumoto <a
+ * href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng.pdf">Improved
+ * Long-Period Generators Based on Linear Recurrences Modulo 2</a> ACM
+ * Transactions on Mathematical Software, 32, 1 (2006). The errata for the paper
+ * are in <a href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng-errata.txt">wellrng-errata.txt</a>.</p>
+
+ * @see <a href="http://www.iro.umontreal.ca/~panneton/WELLRNG.html">WELL Random number generator</a>
+ * @version $Revision: 1003892 $ $Date: 2010-10-02 23:28:56 +0200 (sam. 02 oct. 2010) $
+ * @since 2.2
+
+ */
+public class Well44497a extends AbstractWell {
+
+    /** Serializable version identifier. */
+    private static final long serialVersionUID = -3859207588353972099L;
+
+    /** Number of bits in the pool. */
+    private static final int K = 44497;
+
+    /** First parameter of the algorithm. */
+    private static final int M1 = 23;
+
+    /** Second parameter of the algorithm. */
+    private static final int M2 = 481;
+
+    /** Third parameter of the algorithm. */
+    private static final int M3 = 229;
+
+    /** Creates a new random number generator.
+     * <p>The instance is initialized using the current time as the
+     * seed.</p>
+     */
+    public Well44497a() {
+        super(K, M1, M2, M3);
+    }
+
+    /** Creates a new random number generator using a single int seed.
+     * @param seed the initial seed (32 bits integer)
+     */
+    public Well44497a(int seed) {
+        super(K, M1, M2, M3, seed);
+    }
+
+    /** Creates a new random number generator using an int array seed.
+     * @param seed the initial seed (32 bits integers array), if null
+     * the seed of the generator will be related to the current time
+     */
+    public Well44497a(int[] seed) {
+        super(K, M1, M2, M3, seed);
+    }
+
+    /** Creates a new random number generator using a single long seed.
+     * @param seed the initial seed (64 bits integer)
+     */
+    public Well44497a(long seed) {
+        super(K, M1, M2, M3, seed);
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    protected int next(final int bits) {
+
+        final int indexRm1 = iRm1[index];
+        final int indexRm2 = iRm2[index];
+
+        final int v0       = v[index];
+        final int vM1      = v[i1[index]];
+        final int vM2      = v[i2[index]];
+        final int vM3      = v[i3[index]];
+
+        // the values below include the errata of the original article
+        final int z0       = (0xFFFF8000 & v[indexRm1]) ^ (0x00007FFF & v[indexRm2]);
+        final int z1       = (v0 ^ (v0 << 24))  ^ (vM1 ^ (vM1 >>> 30));
+        final int z2       = (vM2 ^ (vM2 << 10)) ^ (vM3 << 26);
+        final int z3       = z1      ^ z2;
+        final int z2Prime  = ((z2 << 9) ^ (z2 >>> 23)) & 0xfbffffff;
+        final int z2Second = ((z2 & 0x00020000) != 0) ? (z2Prime ^ 0xb729fcec) : z2Prime;
+        final int z4       = z0 ^ (z1 ^ (z1 >>> 20)) ^ z2Second ^ z3;
+
+        v[index]     = z3;
+        v[indexRm1]  = z4;
+        v[indexRm2] &= 0xFFFF8000;
+        index        = indexRm1;
+
+        return z4 >>> (32 - bits);
+
+    }
+}
diff --git a/src/main/java/org/apache/commons/math/random/Well44497b.java b/src/main/java/org/apache/commons/math/random/Well44497b.java
new file mode 100644
index 0000000..28b5dbb
--- /dev/null
+++ b/src/main/java/org/apache/commons/math/random/Well44497b.java
@@ -0,0 +1,119 @@
+/*
+ * 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.math.random;
+
+
+/** This class implements the WELL44497b pseudo-random number generator
+ * from Fran&ccedil;ois Panneton, Pierre L'Ecuyer and Makoto Matsumoto.
+
+ * <p>This generator is described in a paper by Fran&ccedil;ois Panneton,
+ * Pierre L'Ecuyer and Makoto Matsumoto <a
+ * href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng.pdf">Improved
+ * Long-Period Generators Based on Linear Recurrences Modulo 2</a> ACM
+ * Transactions on Mathematical Software, 32, 1 (2006). The errata for the paper
+ * are in <a href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng-errata.txt">wellrng-errata.txt</a>.</p>
+
+ * @see <a href="http://www.iro.umontreal.ca/~panneton/WELLRNG.html">WELL Random number generator</a>
+ * @version $Revision: 1003892 $ $Date: 2010-10-02 23:28:56 +0200 (sam. 02 oct. 2010) $
+ * @since 2.2
+
+ */
+public class Well44497b extends AbstractWell {
+
+    /** Serializable version identifier. */
+    private static final long serialVersionUID = 4032007538246675492L;
+
+    /** Number of bits in the pool. */
+    private static final int K = 44497;
+
+    /** First parameter of the algorithm. */
+    private static final int M1 = 23;
+
+    /** Second parameter of the algorithm. */
+    private static final int M2 = 481;
+
+    /** Third parameter of the algorithm. */
+    private static final int M3 = 229;
+
+    /** Creates a new random number generator.
+     * <p>The instance is initialized using the current time as the
+     * seed.</p>
+     */
+    public Well44497b() {
+        super(K, M1, M2, M3);
+    }
+
+    /** Creates a new random number generator using a single int seed.
+     * @param seed the initial seed (32 bits integer)
+     */
+    public Well44497b(int seed) {
+        super(K, M1, M2, M3, seed);
+    }
+
+    /** Creates a new random number generator using an int array seed.
+     * @param seed the initial seed (32 bits integers array), if null
+     * the seed of the generator will be related to the current time
+     */
+    public Well44497b(int[] seed) {
+        super(K, M1, M2, M3, seed);
+    }
+
+    /** Creates a new random number generator using a single long seed.
+     * @param seed the initial seed (64 bits integer)
+     */
+    public Well44497b(long seed) {
+        super(K, M1, M2, M3, seed);
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    protected int next(final int bits) {
+
+        // compute raw value given by WELL44497a generator
+        // which is NOT maximally-equidistributed
+        final int indexRm1 = iRm1[index];
+        final int indexRm2 = iRm2[index];
+
+        final int v0       = v[index];
+        final int vM1      = v[i1[index]];
+        final int vM2      = v[i2[index]];
+        final int vM3      = v[i3[index]];
+
+        // the values below include the errata of the original article
+        final int z0       = (0xFFFF8000 & v[indexRm1]) ^ (0x00007FFF & v[indexRm2]);
+        final int z1       = (v0 ^ (v0 << 24))  ^ (vM1 ^ (vM1 >>> 30));
+        final int z2       = (vM2 ^ (vM2 << 10)) ^ (vM3 << 26);
+        final int z3       = z1      ^ z2;
+        final int z2Prime  = ((z2 << 9) ^ (z2 >>> 23)) & 0xfbffffff;
+        final int z2Second = ((z2 & 0x00020000) != 0) ? (z2Prime ^ 0xb729fcec) : z2Prime;
+        int z4             = z0 ^ (z1 ^ (z1 >>> 20)) ^ z2Second ^ z3;
+
+        v[index]     = z3;
+        v[indexRm1]  = z4;
+        v[indexRm2] &= 0xFFFF8000;
+        index        = indexRm1;
+
+        // add Matsumoto-Kurita tempering
+        // to get a maximally-equidistributed generator
+        z4 = z4 ^ ((z4 <<  7) & 0x93dd1400);
+        z4 = z4 ^ ((z4 << 15) & 0xfa118000);
+
+        return z4 >>> (32 - bits);
+
+    }
+
+}
diff --git a/src/main/java/org/apache/commons/math/random/Well512a.java b/src/main/java/org/apache/commons/math/random/Well512a.java
new file mode 100644
index 0000000..c822c9d
--- /dev/null
+++ b/src/main/java/org/apache/commons/math/random/Well512a.java
@@ -0,0 +1,107 @@
+/*
+ * 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.math.random;
+
+
+/** This class implements the WELL512a pseudo-random number generator
+ * from Fran&ccedil;ois Panneton, Pierre L'Ecuyer and Makoto Matsumoto.
+
+ * <p>This generator is described in a paper by Fran&ccedil;ois Panneton,
+ * Pierre L'Ecuyer and Makoto Matsumoto <a
+ * href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng.pdf">Improved
+ * Long-Period Generators Based on Linear Recurrences Modulo 2</a> ACM
+ * Transactions on Mathematical Software, 32, 1 (2006). The errata for the paper
+ * are in <a href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng-errata.txt">wellrng-errata.txt</a>.</p>
+
+ * @see <a href="http://www.iro.umontreal.ca/~panneton/WELLRNG.html">WELL Random number generator</a>
+ * @version $Revision: 1003892 $ $Date: 2010-10-02 23:28:56 +0200 (sam. 02 oct. 2010) $
+ * @since 2.2
+
+ */
+public class Well512a extends AbstractWell {
+
+    /** Serializable version identifier. */
+    private static final long serialVersionUID = -6104179812103820574L;
+
+    /** Number of bits in the pool. */
+    private static final int K = 512;
+
+    /** First parameter of the algorithm. */
+    private static final int M1 = 13;
+
+    /** Second parameter of the algorithm. */
+    private static final int M2 = 9;
+
+    /** Third parameter of the algorithm. */
+    private static final int M3 = 5;
+
+    /** Creates a new random number generator.
+     * <p>The instance is initialized using the current time as the
+     * seed.</p>
+     */
+    public Well512a() {
+        super(K, M1, M2, M3);
+    }
+
+    /** Creates a new random number generator using a single int seed.
+     * @param seed the initial seed (32 bits integer)
+     */
+    public Well512a(int seed) {
+        super(K, M1, M2, M3, seed);
+    }
+
+    /** Creates a new random number generator using an int array seed.
+     * @param seed the initial seed (32 bits integers array), if null
+     * the seed of the generator will be related to the current time
+     */
+    public Well512a(int[] seed) {
+        super(K, M1, M2, M3, seed);
+    }
+
+    /** Creates a new random number generator using a single long seed.
+     * @param seed the initial seed (64 bits integer)
+     */
+    public Well512a(long seed) {
+        super(K, M1, M2, M3, seed);
+    }
+
+    /** {@inheritDoc} */
+    @Override
+    protected int next(final int bits) {
+
+        final int indexRm1 = iRm1[index];
+
+        final int vi = v[index];
+        final int vi1 = v[i1[index]];
+        final int vi2 = v[i2[index]];
+        final int z0 = v[indexRm1];
+
+        // the values below include the errata of the original article
+        final int z1 = (vi ^ (vi << 16))   ^ (vi1 ^ (vi1 << 15));
+        final int z2 = vi2 ^ (vi2 >>> 11);
+        final int z3 = z1 ^ z2;
+        final int z4 = (z0 ^ (z0 << 2)) ^ (z1 ^ (z1 << 18)) ^ (z2 << 28) ^ (z3 ^ ((z3 << 5) & 0xda442d24));
+
+        v[index] = z3;
+        v[indexRm1]  = z4;
+        index    = indexRm1;
+
+        return z4 >>> (32 - bits);
+
+    }
+
+}
diff --git a/src/main/java/org/apache/commons/math/random/package.html b/src/main/java/org/apache/commons/math/random/package.html
new file mode 100644
index 0000000..9978ca0
--- /dev/null
+++ b/src/main/java/org/apache/commons/math/random/package.html
@@ -0,0 +1,132 @@
+<html>
+<!--
+   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.
+  -->
+    <!-- $Revision: 1054186 $ $Date: 2011-01-01 03:28:46 +0100 (sam. 01 janv. 2011) $ -->
+    <body>
+      <p>Random number and random data generators.</p>
+      <p>Commons-math provides a few pseudo random number generators. The top level interface is RandomGenerator.
+      It is implemented by three classes:
+      <ul>
+        <li>{@link org.apache.commons.math.random.JDKRandomGenerator JDKRandomGenerator}
+            that extends the JDK provided generator</li>
+        <li>AbstractRandomGenerator as a helper for users generators</li>
+        <li>BitStreamGenerator which is an abstract class for several generators and
+            which in turn is extended by:
+            <ul>
+              <li>{@link org.apache.commons.math.random.MersenneTwister MersenneTwister}</li>
+              <li>{@link org.apache.commons.math.random.Well512a Well512a}</li>
+              <li>{@link org.apache.commons.math.random.Well1024a Well1024a}</li>
+              <li>{@link org.apache.commons.math.random.Well19937a Well19937a}</li>
+              <li>{@link org.apache.commons.math.random.Well19937c Well19937c}</li>
+              <li>{@link org.apache.commons.math.random.Well44497a Well44497a}</li>
+              <li>{@link org.apache.commons.math.random.Well44497b Well44497b}</li>
+            </ul>
+          </li>
+        </ul>
+      </p>
+
+      <p>
+      The JDK provided generator is a simple one that can be used only for very simple needs.
+      The Mersenne Twister is a fast generator with very good properties well suited for
+      Monte-Carlo simulation. It is equidistributed for generating vectors up to dimension 623
+      and has a huge period: 2<sup>19937</sup> - 1 (which is a Mersenne prime). This generator
+      is described in a paper by Makoto Matsumoto and Takuji Nishimura in 1998: <a
+      href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/mt.pdf">Mersenne Twister:
+      A 623-Dimensionally Equidistributed Uniform Pseudo-Random Number Generator</a>, ACM
+      Transactions on Modeling and Computer Simulation, Vol. 8, No. 1, January 1998, pp 3--30.
+      The WELL generators are a family of generators with period ranging from 2<sup>512</sup> - 1
+      to 2<sup>44497</sup> - 1 (this last one is also a Mersenne prime) with even better properties
+      than Mersenne Twister. These generators are described in a paper by Fran&ccedil;ois Panneton,
+      Pierre L'Ecuyer and Makoto Matsumoto <a
+      href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng.pdf">Improved Long-Period
+      Generators Based on Linear Recurrences Modulo 2</a> ACM Transactions on Mathematical Software,
+      32, 1 (2006). The errata for the paper are in <a
+      href="http://www.iro.umontreal.ca/~lecuyer/myftp/papers/wellrng-errata.txt">wellrng-errata.txt</a>.
+      </p>
+
+      <p>
+      For simple sampling, any of these generators is sufficient. For Monte-Carlo simulations the
+      JDK generator does not have any of the good mathematical properties of the other generators,
+      so it should be avoided. The Mersenne twister and WELL generators have equidistribution properties
+      proven according to their bits pool size which is directly linked to their period (all of them
+      have maximal period, i.e. a generator with size n pool has a period 2<sup>n</sup>-1). They also
+      have equidistribution properties for 32 bits blocks up to s/32 dimension where s is their pool size.
+      So WELL19937c for exemple is equidistributed up to dimension 623 (19937/32). This means a Monte-Carlo
+      simulation generating a vector of n variables at each iteration has some guarantees on the properties
+      of the vector as long as its dimension does not exceed the limit. However, since we use bits from two
+      successive 32 bits generated integers to create one double, this limit is smaller when the variables are
+      of type double. so for Monte-Carlo simulation where less the 16 doubles are generated at each round,
+      WELL1024 may be sufficient. If a larger number of doubles are needed a generator with a larger pool
+      would be useful.
+      </p>
+
+      <p>
+      The WELL generators are more modern then MersenneTwister (the paper describing than has been published
+      in 2006 instead of 1998) and fix some of its (few) drawbacks. If initialization array contains many
+      zero bits, MersenneTwister may take a very long time (several hundreds of thousands of iterations to
+      reach a steady state with a balanced number of zero and one in its bits pool). So the WELL generators
+      are better to <i>escape zeroland</i> as explained by the WELL generators creators. The Well19937a and
+      Well44497a generator are not maximally equidistributed (i.e. there are some dimensions or bits blocks
+      size for which they are not equidistributed). The Well512a, Well1024a, Well19937c and Well44497b are
+      maximally equidistributed for blocks size up to 32 bits (they should behave correctly also for double
+      based on more than 32 bits blocks, but equidistribution is not proven at these blocks sizes).
+      </p>
+
+      <p>
+      The MersenneTwister generator uses a 624 elements integer array, so it consumes less than 2.5 kilobytes.
+      The WELL generators use 6 integer arrays with a size equal to the pool size, so for example the
+      WELL44497b generator uses about 33 kilobytes. This may be important if a very large number of
+      generator instances were used at the same time.
+      </p>
+
+      <p>
+      All generators are quite fast. As an example, here are some comparisons, obtained on a 64 bits JVM on a
+      linux computer with a 2008 processor (AMD phenom Quad 9550 at 2.2 GHz). The generation rate for
+      MersenneTwister was about 27 millions doubles per second (remember we generate two 32 bits integers for
+      each double). Generation rates for other PRNG, relative to MersenneTwister:
+      </p>
+
+      <p>
+        <table border="1" align="center">
+          <tr BGCOLOR="#CCCCFF"><td colspan="2"><font size="+2">Example of performances</font></td></tr>
+          <tr BGCOLOR="#EEEEFF"><font size="+1"><td>Name</td><td>generation rate (relative to MersenneTwister)</td></font></tr>
+          <tr><td>{@link org.apache.commons.math.random.MersenneTwister MersenneTwister}</td><td>1</td></tr>
+          <tr><td>{@link org.apache.commons.math.random.JDKRandomGenerator JDKRandomGenerator}</td><td>between 0.96 and 1.16</td></tr>
+          <tr><td>{@link org.apache.commons.math.random.Well512a Well512a}</td><td>between 0.85 and 0.88</td></tr>
+          <tr><td>{@link org.apache.commons.math.random.Well1024a Well1024a}</td><td>between 0.63 and 0.73</td></tr>
+          <tr><td>{@link org.apache.commons.math.random.Well19937a Well19937a}</td><td>between 0.70 and 0.71</td></tr>
+          <tr><td>{@link org.apache.commons.math.random.Well19937c Well19937c}</td><td>between 0.57 and 0.71</td></tr>
+          <tr><td>{@link org.apache.commons.math.random.Well44497a Well44497a}</td><td>between 0.69 and 0.71</td></tr>
+          <tr><td>{@link org.apache.commons.math.random.Well44497b Well44497b}</td><td>between 0.65 and 0.71</td></tr>
+        </table>
+      </p>
+
+      <p>
+      So for most simulation problems, the better generators like {@link
+      org.apache.commons.math.random.Well19937c Well19937c} and {@link
+      org.apache.commons.math.random.Well44497b Well44497b} are probably very good choices.
+      </p>
+
+      <p>
+      Note that <em>none</em> of these generators are suitable for cryptography. They are devoted
+      to simulation, and to generate very long series with strong properties on the series as a whole
+      (equidistribution, no correlation ...). They do not attempt to create small series but with
+      very strong properties of unpredictability as needed in cryptography.
+      </p>
+
+    </body>
+</html>