smile.math.random.UniversalGenerator Maven / Gradle / Ivy
/******************************************************************************
* Confidential Proprietary *
* (c) Copyright Haifeng Li 2011, All Rights Reserved *
******************************************************************************/
package smile.math.random;
/**
* The so called "Universal Generator" based on multiplicative congruential
* method, which originally appeared in "Toward a Universal Random Number
* Generator" by Marsaglia, Zaman and Tsang. It was later modified by F. James
* in "A Review of Pseudo-random Number Generators". It passes ALL of the tests
* for random number generators and has a period of 2144. It is
* completely portable (gives bit identical results on all machines with at
* least 24-bit mantissas in the floating point representation).
*
* @author Haifeng Li
*/
public class UniversalGenerator implements RandomNumberGenerator {
/**
* Default seed. DEFAULT_RANDOM_SEED=54217137
*/
private static final int DEFAULT_RANDOM_SEED = 54217137;
/**
* The 46,009,220nd prime number,
* The largest prime less than 9*108. Used as a modulus
* because this version of random() needs a seed between 0
* and 9*108 and BIG_PRIME isn't commensurate
* with any regular period.
* BIG_PRIME = 899999963
*/
private static final int BIG_PRIME = 899999963;
private double c, cd, cm, u[];
private int i97, j97;
/**
* Initialize Random with default seed.
*/
public UniversalGenerator() {
srand(DEFAULT_RANDOM_SEED);
}
/**
* Initialize Random with a specified integer seed
*/
public UniversalGenerator(int seed) {
srand(Math.abs(seed % BIG_PRIME));
}
/**
* Initialize Random with a specified long seed
*/
public UniversalGenerator(long seed) {
srand(Math.abs((int) (seed % BIG_PRIME)));
}
/**
* Initialize the random generator with a seed.
*/
private void srand(int ijkl) {
u = new double[97];
int ij = ijkl / 30082;
int kl = ijkl % 30082;
// Handle the seed range errors
// First random number seed must be between 0 and 31328
// Second seed must have a value between 0 and 30081
if (ij < 0 || ij > 31328 || kl < 0 || kl > 30081) {
ij = ij % 31329;
kl = kl % 30082;
}
int i = ((ij / 177) % 177) + 2;
int j = (ij % 177) + 2;
int k = ((kl / 169) % 178) + 1;
int l = kl % 169;
int m;
double s, t;
for (int ii = 0; ii < 97; ii++) {
s = 0.0;
t = 0.5;
for (int jj = 0; jj < 24; jj++) {
m = (((i * j) % 179) * k) % 179;
i = j;
j = k;
k = m;
l = (53 * l + 1) % 169;
if (((l * m) % 64) >= 32) {
s += t;
}
t *= 0.5;
}
u[ii] = s;
}
c = 362436.0 / 16777216.0;
cd = 7654321.0 / 16777216.0;
cm = 16777213.0 / 16777216.0;
i97 = 96;
j97 = 32;
}
@Override
public double nextDouble() {
double uni;
uni = u[i97] - u[j97];
if (uni < 0.0) {
uni += 1.0;
}
u[i97] = uni;
if (--i97 < 0) {
i97 = 96;
}
if (--j97 < 0) {
j97 = 96;
}
c -= cd;
if (c < 0.0) {
c += cm;
}
uni -= c;
if (uni < 0.0) {
uni++;
}
return uni;
}
@Override
public void nextDoubles(double[] d) {
int n = d.length;
double uni;
for (int i = 0; i < n; i++) {
uni = u[i97] - u[j97];
if (uni < 0.0) {
uni += 1.0;
}
u[i97] = uni;
if (--i97 < 0) {
i97 = 96;
}
if (--j97 < 0) {
j97 = 96;
}
c -= cd;
if (c < 0.0) {
c += cm;
}
uni -= c;
if (uni < 0.0) {
uni += 1.0;
}
d[i] = uni;
}
}
@Override
public int next(int numbits) {
return nextInt() >>> (32 - numbits);
}
@Override
public int nextInt() {
return (int) Math.floor(Integer.MAX_VALUE * (2 * nextDouble() - 1.0));
}
@Override
public int nextInt(int n) {
if (n <= 0) {
throw new IllegalArgumentException("n must be positive");
}
// n is a power of 2
if ((n & -n) == n) {
return (int) ((n * (long) next(31)) >> 31);
}
int bits, val;
do {
bits = next(31);
val = bits % n;
} while (bits - val + (n - 1) < 0);
return val;
}
@Override
public long nextLong() {
return (long) Math.floor(Long.MAX_VALUE * (2 * nextDouble() - 1.0));
}
}