diff --git a/src/main/java/org/apache/commons/math/random/AbstractWell.java b/src/main/java/org/apache/commons/math/random/AbstractWell.java index 232cd71b2..f85fe6f19 100644 --- a/src/main/java/org/apache/commons/math/random/AbstractWell.java +++ b/src/main/java/org/apache/commons/math/random/AbstractWell.java @@ -38,29 +38,22 @@ public abstract class AbstractWell extends BitsStreamGenerator implements Serial /** Serializable version identifier. */ private static final long serialVersionUID = -8068371019303673353L; - /** Number of bits blocks in the pool. */ - private final int r; - /** Bit mask preserving the first w - p bits in a w bits block. */ - private final int mp; + protected final int mp; /** Bit mask preserving the last p bits in a w bits block. */ - private final int mpTilde; - - /** First parameter of the algorithm. */ - private final int m1; - - /** Second parameter of the algorithm. */ - private final int m2; - - /** Third parameter of the algorithm. */ - private final int m3; + protected final int mpTilde; /** Current index in the bytes pool. */ - private int index; + protected int index; /** Bytes pool. */ - private final int[] v; + protected final int[] v; + protected final int[] iRm1; + protected final int[] iRm2; + protected final int[] i1; + protected final int[] i2; + protected final int[] i3; /** Creates a new random number generator. *

The instance is initialized using the current time as the @@ -99,18 +92,29 @@ public abstract class AbstractWell extends BitsStreamGenerator implements Serial // 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; - this.r = (k + w - 1) / w; + final int r = (k + w - 1) / w; final int p = r * w - k; // set up generator parameters this.mp = (-1) << p; this.mpTilde = ~mp; - this.m1 = m1; - this.m2 = m2; - this.m3 = m3; this.v = new int[r]; this.index = 0; + // set up indirection indices + 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); @@ -171,164 +175,7 @@ public abstract class AbstractWell extends BitsStreamGenerator implements Serial setSeed(new int[] { (int) (seed >>> 32), (int) (seed & 0xffffffffl) }); } - /** Generate next pseudorandom number. - *

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()}.

- *

This implementation is the general WELL algorithm, described in - * a paper by François Panneton, Pierre L'Ecuyer and Makoto Matsumoto - * Improved - * Long-Period Generators Based on Linear Recurrences Modulo 2 ACM - * Transactions on Mathematical Software, 32, 1 (2006).

- * @param bits number of random bits to produce - * @return random bits generated - */ - protected int next(final int bits) { - - final int iRm1 = (index + r - 1) % r; - final int iRm2 = (index + r - 2) % r; - final int i1 = (index + m1) % r; - final int i2 = (index + m2) % r; - final int i3 = (index + m3) % r; - - final int z0 = (mp & v[iRm1]) ^ (mpTilde & v[iRm2]); - final int z1 = t0(v[index]) ^ t1(v[i1]); - final int z2 = t2(v[i2]) ^ t3(v[i3]); - final int z3 = z1 ^ z2; - final int z4 = t4(z0) ^ t5(z1) ^ t6(z2) ^ t7(z3); - - v[index] = z3; - v[iRm1] = z4; - v[iRm2] &= mp; - index = iRm1; - - return z4 >>> (32 - bits); - - } - - /** Apply transform M0 to a bits block. - * @param x bits block to apply transform to - * @return M0(x) - */ - protected int m0(final int x) { - return 0; - } - - /** Apply transform M1 to a bits block. - * @param x bits block to apply transform to - * @return M1(x) - */ - protected int m1(final int x) { - return x; - } - - /** Apply transform M2 to a bits block. - * @param t parameter of the transform - * @param x bits block to apply transform to - * @return M2, t(x) - */ - protected int m2(final int t, final int x) { - return (t >= 0) ? (x >>> t) : (x << -t); - } - - /** Apply transform M3 to a bits block. - * @param t parameter of the transform - * @param x bits block to apply transform to - * @return M3, t(x) - */ - protected int m3(final int t, final int x) { - return x ^ ((t >= 0) ? (x >>> t) : (x << -t)); - } - - /** Apply transform M4 to a bits block. - * @param a parameter of the transform - * @param x bits block to apply transform to - * @return M4, a(x) - */ - protected int m4(final int a, final int x) { - final int shiftedX = x >>> 1; - return ((x & 0x80000000) != 0) ? (shiftedX ^ a) : shiftedX; - } - - /** Apply transform M5 to a bits block. - * @param t first parameter of the transform - * @param b second parameter of the transform - * @param x bits block to apply transform to - * @return M5, t, b(x) - */ - protected int m5(final int t, final int b, final int x) { - // table I of the paper specifies that a left shift for positive t and - // a right shift for negative t, however, reference implementation does - // the opposite (and in fact this transform is used only by Well512a - // with t = -28). Here, we follow the reference implementation with a - // left shift for NEGATIVE t - final int shiftedX = (t >= 0) ? (x >>> t) : (x << -t); - return x ^ (shiftedX & b); - } - - /** Apply transform M6 to a bits block. - * @param q first parameter of the transform - * @param dsMask second parameter of the transform as a bit mask - * @param tMask third parameter of the transform as a bit mask - * @param a fourth parameter of the transform - * @param x bits block to apply transform to - * @return M6, q, s, t, a(x) - */ - protected int m6(final int q, final int dsMask, final int tMask, final int a, final int x) { - final int lShiftedX = x << q; - final int rShiftedX = x >>> (32 - q); - final int z = (lShiftedX ^ rShiftedX) & dsMask; - return ((x & tMask) != 0) ? (z ^ a) : z; - } - - /** Apply transform T0 to a bits block. - * @param vi0 bits block to apply transform to - * @return T0 vi,0 - */ - protected abstract int t0(int vi0); - - /** Apply transform T1 to a bits block. - * @param vim1 bits block to apply transform to - * @return T1 vi,m1 - */ - protected abstract int t1(int vim1); - - /** Apply transform T2 to a bits block. - * @param vim2 bits block to apply transform to - * @return T2 vi,m2 - */ - protected abstract int t2(int vim2); - - /** Apply transform T3 to a bits block. - * @param vim3 bits block to apply transform to - * @return T3 vi,m3 - */ - protected abstract int t3(int vim3); - - /** Apply transform T4 to a bits block. - * @param z0 bits block to apply transform to - * @return T4 z0 - */ - protected abstract int t4(int z0); - - /** Apply transform T5 to a bits block. - * @param z1 bits block to apply transform to - * @return T5 z1 - */ - protected abstract int t5(int z1); - - /** Apply transform T6 to a bits block. - * @param z2 bits block to apply transform to - * @return T6 z2 - */ - protected abstract int t6(int z2); - - /** Apply transform T7 to a bits block. - * @param z3 bits block to apply transform to - * @return T7 z3 - */ - protected abstract int t7(int z3); + /** {@inheritDoc} */ + protected abstract int next(final int bits); } diff --git a/src/main/java/org/apache/commons/math/random/Well1024a.java b/src/main/java/org/apache/commons/math/random/Well1024a.java index d12a2eb3d..a0eb2c8ba 100644 --- a/src/main/java/org/apache/commons/math/random/Well1024a.java +++ b/src/main/java/org/apache/commons/math/random/Well1024a.java @@ -34,7 +34,7 @@ package org.apache.commons.math.random; public class Well1024a extends AbstractWell { /** Serializable version identifier. */ - private static final long serialVersionUID = -5403981908127539981L; + private static final long serialVersionUID = 5680173464174485492L; /** Number of bits in the pool. */ private static final int K = 1024; @@ -79,43 +79,28 @@ public class Well1024a extends AbstractWell { } /** {@inheritDoc} */ - protected int t0(final int vi0) { - return m1(vi0); - } + protected int next(final int bits) { - /** {@inheritDoc} */ - protected int t1(final int vim1) { - return m3(8, vim1); - } + final int indexRm1 = iRm1[index]; + final int indexRm2 = iRm2[index]; - /** {@inheritDoc} */ - protected int t2(final int vim2) { - return m3(-19, vim2); - } + final int v0 = v[index]; + final int vM1 = v[i1[index]]; + final int vM2 = v[i2[index]]; + final int vM3 = v[i3[index]]; - /** {@inheritDoc} */ - protected int t3(final int vim3) { - return m3(-14, vim3); - } + 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)); - /** {@inheritDoc} */ - protected int t4(final int z0) { - return m3(-11, z0); - } + v[index] = z3; + v[indexRm1] = z4; + v[indexRm2] &= mp; + index = indexRm1; - /** {@inheritDoc} */ - protected int t5(final int z1) { - return m3(-7, z1); - } + return z4 >>> (32 - bits); - /** {@inheritDoc} */ - protected int t6(final int z2) { - return m3(-13, z2); } - - /** {@inheritDoc} */ - protected int t7(final int z3) { - return m0(z3); - } - } diff --git a/src/main/java/org/apache/commons/math/random/Well19937a.java b/src/main/java/org/apache/commons/math/random/Well19937a.java index 2b754d356..55a17ea62 100644 --- a/src/main/java/org/apache/commons/math/random/Well19937a.java +++ b/src/main/java/org/apache/commons/math/random/Well19937a.java @@ -34,7 +34,7 @@ package org.apache.commons.math.random; public class Well19937a extends AbstractWell { /** Serializable version identifier. */ - private static final long serialVersionUID = -8052371714518610855L; + private static final long serialVersionUID = -7462102162223815419L; /** Number of bits in the pool. */ private static final int K = 19937; @@ -79,43 +79,28 @@ public class Well19937a extends AbstractWell { } /** {@inheritDoc} */ - protected int t0(final int vi0) { - return m3(-25, vi0); - } + protected int next(final int bits) { - /** {@inheritDoc} */ - protected int t1(final int vim1) { - return m3(27, vim1); - } + final int indexRm1 = iRm1[index]; + final int indexRm2 = iRm2[index]; - /** {@inheritDoc} */ - protected int t2(final int vim2) { - return m2(9, vim2); - } + final int v0 = v[index]; + final int vM1 = v[i1[index]]; + final int vM2 = v[i2[index]]; + final int vM3 = v[i3[index]]; - /** {@inheritDoc} */ - protected int t3(final int vim3) { - return m3(1, vim3); - } + 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)); - /** {@inheritDoc} */ - protected int t4(final int z0) { - return m1(z0); - } + v[index] = z3; + v[indexRm1] = z4; + v[indexRm2] &= mp; + index = indexRm1; - /** {@inheritDoc} */ - protected int t5(final int z1) { - return m3(-9, z1); - } + return z4 >>> (32 - bits); - /** {@inheritDoc} */ - protected int t6(final int z2) { - return m3(-21, z2); } - - /** {@inheritDoc} */ - protected int t7(final int z3) { - return m3(21, z3); - } - } diff --git a/src/main/java/org/apache/commons/math/random/Well19937c.java b/src/main/java/org/apache/commons/math/random/Well19937c.java index aac9109fe..5475549ff 100644 --- a/src/main/java/org/apache/commons/math/random/Well19937c.java +++ b/src/main/java/org/apache/commons/math/random/Well19937c.java @@ -31,23 +31,36 @@ package org.apache.commons.math.random; * @since 2.2 */ -public class Well19937c extends Well19937a { +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. *

The instance is initialized using the current time as the * seed.

*/ 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(seed); + super(K, M1, M2, M3, seed); } /** Creates a new random number generator using an int array seed. @@ -55,29 +68,45 @@ public class Well19937c extends Well19937a { * the seed of the generator will be related to the current time */ public Well19937c(int[] seed) { - super(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(seed); + super(K, M1, M2, M3, seed); } /** {@inheritDoc} */ protected int next(final int bits) { - // compute raw value given by WELL19937a generator - // which is NOT maximally-equidistributed - int z = super.next(32); + 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] &= mp; + index = indexRm1; + // add Matsumoto-Kurita tempering // to get a maximally-equidistributed generator - z = z ^ ((z << 7) & 0xe46e1700); - z = z ^ ((z << 15) & 0x9b868000); + z4 = z4 ^ ((z4 << 7) & 0xe46e1700); + z4 = z4 ^ ((z4 << 15) & 0x9b868000); - return z >>> (32 - bits); + 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 index 6841d008d..4aa55905d 100644 --- a/src/main/java/org/apache/commons/math/random/Well44497a.java +++ b/src/main/java/org/apache/commons/math/random/Well44497a.java @@ -34,7 +34,7 @@ package org.apache.commons.math.random; public class Well44497a extends AbstractWell { /** Serializable version identifier. */ - private static final long serialVersionUID = 5154222742730470272L; + private static final long serialVersionUID = -3859207588353972099L; /** Number of bits in the pool. */ private static final int K = 44497; @@ -79,46 +79,30 @@ public class Well44497a extends AbstractWell { } /** {@inheritDoc} */ - protected int t0(final int vi0) { - return m3(-24, vi0); - } + protected int next(final int bits) { - /** {@inheritDoc} */ - protected int t1(final int vim1) { - return m3(30, vim1); - } + final int indexRm1 = iRm1[index]; + final int indexRm2 = iRm2[index]; - /** {@inheritDoc} */ - protected int t2(final int vim2) { - return m3(-10, vim2); - } + final int v0 = v[index]; + final int vM1 = v[i1[index]]; + final int vM2 = v[i2[index]]; + final int vM3 = v[i3[index]]; - /** {@inheritDoc} */ - protected int t3(final int vim3) { - return m2(-26, vim3); - } + 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; - /** {@inheritDoc} */ - protected int t4(final int z0) { - return m1(z0); - } + v[index] = z3; + v[indexRm1] = z4; + v[indexRm2] &= mp; + index = indexRm1; - /** {@inheritDoc} */ - protected int t5(final int z1) { - return m3(20, z1); - } + return z4 >>> (32 - bits); - /** {@inheritDoc} */ - protected int t6(final int z2) { - // table II of the paper specifies t6 to be m6(9, d14, t5, 0xb729fcec, z2) - // however, the reference implementation uses m6(9, d26, t17, 0xb729fcec, z2). - // Here, we follow the reference implementation - return m6(9, (-1) ^ (0x1 << 26), 0x1 << 17, 0xb729fcec, z2); } - - /** {@inheritDoc} */ - protected int t7(final int z3) { - return m1(z3); - } - } diff --git a/src/main/java/org/apache/commons/math/random/Well44497b.java b/src/main/java/org/apache/commons/math/random/Well44497b.java index 3caa09cc5..2ffea8cc4 100644 --- a/src/main/java/org/apache/commons/math/random/Well44497b.java +++ b/src/main/java/org/apache/commons/math/random/Well44497b.java @@ -31,23 +31,36 @@ package org.apache.commons.math.random; * @since 2.2 */ -public class Well44497b extends Well44497a { +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. *

The instance is initialized using the current time as the * seed.

*/ 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(seed); + super(K, M1, M2, M3, seed); } /** Creates a new random number generator using an int array seed. @@ -55,14 +68,14 @@ public class Well44497b extends Well44497a { * the seed of the generator will be related to the current time */ public Well44497b(int[] seed) { - super(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(seed); + super(K, M1, M2, M3, seed); } /** {@inheritDoc} */ @@ -70,14 +83,33 @@ public class Well44497b extends Well44497a { // compute raw value given by WELL44497a generator // which is NOT maximally-equidistributed - int z = super.next(32); + 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 = (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] &= mp; + index = indexRm1; // add Matsumoto-Kurita tempering // to get a maximally-equidistributed generator - z = z ^ ((z << 7) & 0x93dd1400); - z = z ^ ((z << 15) & 0xfa118000); + z4 = z4 ^ ((z4 << 7) & 0x93dd1400); + z4 = z4 ^ ((z4 << 15) & 0xfa118000); - return z >>> (32 - bits); + 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 index b6ba995bb..99c4a02be 100644 --- a/src/main/java/org/apache/commons/math/random/Well512a.java +++ b/src/main/java/org/apache/commons/math/random/Well512a.java @@ -34,7 +34,7 @@ package org.apache.commons.math.random; public class Well512a extends AbstractWell { /** Serializable version identifier. */ - private static final long serialVersionUID = 8706771840051210473L; + private static final long serialVersionUID = -6104179812103820574L; /** Number of bits in the pool. */ private static final int K = 512; @@ -79,46 +79,27 @@ public class Well512a extends AbstractWell { } /** {@inheritDoc} */ - protected int t0(final int vi0) { - return m3(-16, vi0); - } + protected int next(final int bits) { - /** {@inheritDoc} */ - protected int t1(final int vim1) { - return m3(-15, vim1); - } + final int indexRm1 = iRm1[index]; - /** {@inheritDoc} */ - protected int t2(final int vim2) { - return m3(11, vim2); - } + final int vi = v[index]; + final int vi1 = v[i1[index]]; + final int vi2 = v[i2[index]]; + final int z0 = v[indexRm1]; + // m3: x ^ ((t >= 0) ? (x >>> t) : (x << -t)); - /** {@inheritDoc} */ - protected int t3(final int vim3) { - return m0(vim3); - } + 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)); - /** {@inheritDoc} */ - protected int t4(final int z0) { - return m3(-2, z0); - } + v[index] = z3; + v[indexRm1] = z4; + index = indexRm1; - /** {@inheritDoc} */ - protected int t5(final int z1) { - return m3(-18, z1); - } + return z4 >>> (32 - bits); - /** {@inheritDoc} */ - protected int t6(final int z2) { - // table II of the paper specifies t6 to be m3(-28, z2) - // however, the reference implementation uses m2(-28, z2). - // Here, we follow the reference implementation - return m2(-28, z2); - } - - /** {@inheritDoc} */ - protected int t7(final int z3) { - return m5(-5, 0xda442d24, z3); } }