edu.cornell.lassp.houle.RngPack.Ranlux Maven / Gradle / Ivy
package edu.cornell.lassp.houle.RngPack;
import java.util.*;
import java.io.Serializable;
//
// RngPack 1.1a by Paul Houle
// http://www.honeylocust.com/RngPack/
//
/**
*
* RANLUX is an advanced pseudo-random number generator based on the
* RCARRY algorithm proposed in 1991 by Marsaglia and Zaman.
* RCARRY
* used a subtract-and-borrow algorithm with a period on the order of
* 10171 but still had detectable correlations between
* numbers. Martin Luescher proposed the RANLUX
* algorithm in 1993; RANLUX generates pseudo-random numbers using
* RCARRY but throws away numbers to destroy correlations.
* Thus RANLUX trades execution speed for quality; by choosing a
* larger luxury setting one gets better random numbers slower.
* By the tests availible at the time it was proposed,
* RANLUX at the default luxury setting appears to be a
* significant advance quality over previous
* generators.
*
*
*
*
*
*
*
* LUXURY LEVELS
*
*
*
* level p
* 0 24
* equivalent to the original RCARRY of Marsaglia
* and Zaman, very long period, but fails many tests.
* 1 48 considerable improvement in quality over level 0,
* now passes the gap test, but still fails spectral test.
*
* 2 97 passes all known tests, but theoretically still
* defective.
* 3 223
* DEFAULT VALUE. Any theoretically possible
* correlations have very small chance of being observed.
* 4 389 highest possible luxury, all 24 bits chaotic.
*
*
*
*
*
* VALIDATION
*
*
* The Java version of RANLUX has been verified against published
* values of numbers 1-5 and 101-105 produced by the reference implementation
* of RANLUX for the following initial conditions:
*
*
* - Default initialization:
Ranlux()
* - Initialization with:
Ranlux(0,0)
* - Initialization with:
Ranlux(389,1)
* - Initialization with:
Ranlux(75,0)
*
* References:
*
* -
* M. Luscher, Computer Physics Communications 79 (1994) 100
*
-
* F. James, Computer Physics Communications 79 (1994) 111
*
- About RANLUX random number generator: Excerpts from discussion in the Usenet news groups
*
- Miller's FORTRAN 90 implementation of RANLUX with test code
*
*
*
*
*
* Source code is available.
*
* @author Paul Houle (E-mail: [email protected])
* @version 1.1a
*/
public class Ranlux extends RandomSeedable implements Serializable {
/**
* Maximum luxury level: maxlev=4
*/
public static final int maxlev=4;
/**
* Default luxury level: lxdflt=3
*/
public static final int lxdflt=3;
static final int igiga = 1000000000, jsdflt = 314159265;
static final int twop12 = 4096,itwo24 = 1 << 24,icons = 2147483563;
static final int ndskip[]={0,24,73,199,365};
int iseeds[],isdext[],next[];
int luxlev=lxdflt,nskip,inseed,jseed;
int in24 = 0, kount = 0, mkount = 0, i24 = 24, j24 = 10;
float seeds[],carry= (float) 0.0,twom24,twom12;
boolean diagOn=false;
/**
* Default initialization of RANLUX. Uses default seed
* jsdflt=314159265 and luxury level 3.
*/
public Ranlux() {
init_arrays();
rluxdef();
};
/**
* Initialize RANLUX with specified luxury level
* and seed.
*
* @param lux luxury level from 0-4.
* @param ins seed, a positive integer.
*
*/
public Ranlux(int lux,int ins) {
init_arrays();
rluxgo(lux,Math.abs(ins));
};
/**
* Initialize RANLUX with specified luxury level
* and seed.
*
* @param lux luxury level from 0-4.
* @param ins seed, a positive long.
*
*/
public Ranlux(int lux,long ins) {
init_arrays();
rluxgo(lux,Math.abs((int) (ins%Integer.MAX_VALUE)));
};
/**
*
* Initialize RANLUX with default luxury level
* and a specified seed.
*
* @param ins seed, a positive integer
*/
public Ranlux(int ins) {
init_arrays();
rluxgo(lxdflt,Math.abs(ins));
};
/**
*
* Initialize RANLUX with default luxury level
* and a specified seed.
*
* @param ins seed, a positive integer
*/
public Ranlux(long ins) {
init_arrays();
rluxgo(lxdflt,Math.abs((int) (ins%Integer.MAX_VALUE)));
};
/**
*
* Initialize RANLUX with specified luxury level
* and a Date object. Can be used to conveniently initialize RANLUX
* from the clock,
*
*
* RandomElement e = Ranlux(4,new Date());
*
*
* @param lux luxury level from 0-4.
* @param d date used to generate seed
*
*/
public Ranlux(int lux,Date d) {
init_arrays();
rluxgo(lux,(int) (ClockSeed(d)%Integer.MAX_VALUE) );
};
/**
*
* Initialize RANLUX with default luxury level
* and a Date object. Can be used to conveniently initialize RANLUX
* from the clock,
*
*
* RandomElement e = Ranlux(new Date());
*
*
* @param d date used to generate seed
*
*/
public Ranlux(Date d) {
init_arrays();
rluxgo(lxdflt,(int) (ClockSeed(d)%Integer.MAX_VALUE) );
};
/**
* Turns diagnostic messages on and off. If setDiag(true) is
* called, RANLUX will print diagnostic information to
* System.err
*
* @param b diagnostic message status
*/
public void setDiag(boolean b) {
diagOn=b;
};
/**
*
* The random number generator.
*
* @return a pseudo-random double in the range (0,1)
*/
public final double raw() {
int i,k,lp;
float uni,out;
uni=seeds[j24]-seeds[i24]-carry;
if (uni < (float) 0.0) {
uni=uni+ (float) 1.0;
carry = twom24;
} else carry = (float) 0.0;
seeds[i24]=uni;
i24=next[i24];
j24=next[j24];
out=uni;
if (uni=igiga) {
mkount++;
kount -= igiga;
};
return out;
};
private void init_arrays() {
/*
*
* converted from fortran: fortran arrays start at 1, java arrays start at
* 0. Here we take the low road to compatibility... We declare arrays that
* are one bigger than the fortran code. This wastes three ints and a double;
* If you're porting this to a Commodore 64 you might need the space.
*
*/
iseeds = new int[24+1];
isdext = new int[25+1];
next = new int[24+1];
seeds = new float[24+1];
};
private void rluxdef()
{
int lp,i,k;
jseed = jsdflt;
inseed = jseed;
diag("RANLUX DEFAULT INITIALIZATION: "+jseed);
luxlev = lxdflt;
nskip = ndskip[luxlev];
lp = nskip + 24;
in24 = 0;
kount = 0;
mkount = 0;
diag("RANLUX DEFAULT LUXURY LEVEL = "+luxlev+" p = "+lp);
twom24 = (float) 1.0;
for(i = 1;i <= 24;i++)
{
twom24 = twom24 * (float) 0.5;
k = jseed / 53668;
jseed = 40014 * (jseed-k*53668) - k * 12211;
if (jseed<0) jseed = jseed + icons;
iseeds[i] = jseed % itwo24;
};
twom12 = twom24 * (float) 4096.0;
for(i = 1 ; i<=24;i++)
{
seeds[i] = iseeds[i] * twom24;
next[i] = i - 1;
};
next[1] = 24;
i24 = 24;
j24 = 10;
carry = (float) 0.0;
if (seeds[24] == 0.0) carry = twom24;
};
private final void rluxgo(int lux,int ins)
{
int ilx, i, iouter, isk, k, inner, izip, izip2;
float uni;
if (lux<00)
{
luxlev = lxdflt;
}
else if (lux<=maxlev)
{
luxlev = lux;
}
else if (lux<24 || lux>2000)
{
luxlev = maxlev;
diag("RANLUX ILLEGAL LUXURY RLUXGO: "+lux);
} else
{
luxlev = lux;
for(ilx=0;ilx<=maxlev;ilx++)
if (lux == ndskip[ilx]+24)
luxlev=ilx;
};
if (luxlev <= maxlev) {
nskip = ndskip[luxlev];
diag("RANLUX LUXURY LEVEL SET BY RLUXGO : "+luxlev+" P= "+(nskip+24));
} else {
nskip = luxlev-24;
diag("RANLUX P-VALUE SET BY RLUXGO TO: "+luxlev);
};
in24=0;
if (ins<0)
diag("Illegal initialization by RLUXGO, negative input seed");
if (ins>0) {
jseed = ins;
diag("RANLUX INITIALIZED BY RLUXGO FROM SEED "+jseed);
} else
{
jseed=jsdflt;
diag("RANLUX INITIALIZED BY RLUXGO FROM DEFAULT SEED");
};
inseed = jseed;
twom24 = (float) 1.0;
for(i=1;i<=24;i++) {
twom24 = twom24 * (float) 0.5;
k = jseed / 53668;
jseed = 40014 * (jseed-k*53668) - k * 12211;
if (jseed < 0) jseed = jseed + icons;
iseeds[i] = jseed % itwo24;
};
twom12=twom24*4096;
for(i=1;i<=24;i++) {
seeds[i] = iseeds[i] * twom24;
next[i] = i - 1;
};
next[1]=24;
i24=24;
j24=10;
carry = (float) 0.0;
if (seeds[24] == 0.0) carry = twom24;
kount = 0;
mkount = 0;
};
private void diag(String s)
{
if (diagOn)
System.err.println(s);
};
};