233 lines
5.8 KiB
C
233 lines
5.8 KiB
C
/* Random number generation using the linear congruential algorithm
|
|
X(n+1) = (a * X(n) + c) mod m
|
|
with a precision of 48 bits.
|
|
|
|
Author: Kriton Kyrimis (kyrimis@cti.gr)
|
|
Code status: Public Domain.
|
|
*/
|
|
|
|
#include <limits.h>
|
|
#include <float.h>
|
|
|
|
/* Parameters for the linear congruential algorithm:
|
|
parm[0..2] is the current value of Xn (internal seed, m.s.word last)
|
|
parm[3..5] is the value of a (m.s.word last)
|
|
parm[6] is the value of c.
|
|
*/
|
|
#define X0 0x1234 /* MSB * Initial value for Xn, obtained using seed48() */
|
|
#define X1 0xABCD /* on SunOS 4.1.3 */
|
|
#define X2 0x330E
|
|
|
|
#define A0 0x0005 /* MSB * Default value for a, taken from the man page */
|
|
#define A1 0xDEEC
|
|
#define A2 0xE66D
|
|
|
|
#define C0 0x000B /* Default value for c, taken from the man page */
|
|
|
|
static unsigned short parm[7] = {
|
|
X2, X1, X0,
|
|
A2, A1, A0,
|
|
C0
|
|
};
|
|
|
|
/* To produce a double random number in [0,1) we get a 32-bit unsigned long
|
|
random number, convert it to double, and divide it by ULONG_MAX + EPSILON.
|
|
(We add the EPSILON to exclude 1.0 from the set of possible results.)
|
|
|
|
We derive EPSILON by noting that for a random value of ULONG_MAX,
|
|
we want to return the smallest double that is less than 1.0.
|
|
Therefore:
|
|
|
|
ULONG_MAX
|
|
--------------------- = (1.0 - DBL_EPSILON)
|
|
(ULONG_MAX + EPSILON)
|
|
|
|
(This is probably overkill.)
|
|
|
|
*/
|
|
|
|
#define EPSILON (double)ULONG_MAX*(1.0/(1.0-DBL_EPSILON)-1.0)
|
|
|
|
|
|
/*--------------------------------------------------------------------------*
|
|
* Parameter initialization functions *
|
|
*--------------------------------------------------------------------------*/
|
|
|
|
/* This function sets the two m.s.words of the internal seed to the value
|
|
supplied by the caller, and the l.s.word of the internal seed to 0x330E.
|
|
*/
|
|
void
|
|
srand48(long seed)
|
|
{
|
|
parm[0] = 0x330E;
|
|
parm[1] = ((unsigned long)seed) & 0xFFFF;
|
|
parm[2] = ((unsigned long)seed >> 16) & 0xFFFF;
|
|
parm[3] = A2;
|
|
parm[4] = A1;
|
|
parm[5] = A0;
|
|
parm[6] = C0;
|
|
}
|
|
|
|
/* This function sets all three words of the internal seed to the value
|
|
supplied by the caller. It returns a pointer to an array containing
|
|
a copy of the previous value of the internal seed.
|
|
*/
|
|
unsigned short *
|
|
seed48(unsigned short *seed)
|
|
{
|
|
static unsigned short oldparm[3];
|
|
unsigned short tmpparm[3];
|
|
|
|
/* Can't assign oldparm[] = parm[] directly, because seed[] may be a pointer
|
|
to oldparm[], obtained from a previous call to seed48 , in which case
|
|
we would destroy the contents of seed[] */
|
|
tmpparm[0] = parm[0];
|
|
tmpparm[1] = parm[1];
|
|
tmpparm[2] = parm[2];
|
|
parm[0] = seed[0];
|
|
parm[1] = seed[1];
|
|
parm[2] = seed[2];
|
|
oldparm[0] = tmpparm[0];
|
|
oldparm[1] = tmpparm[1];
|
|
oldparm[2] = tmpparm[2];
|
|
parm[3] = A2;
|
|
parm[4] = A1;
|
|
parm[5] = A0;
|
|
parm[6] = C0;
|
|
|
|
return oldparm;
|
|
}
|
|
|
|
/* This function sets all seven words of the internal parameters array to the
|
|
values supplied by the caller.
|
|
*/
|
|
void
|
|
lcong48(unsigned short *new_parm)
|
|
{
|
|
parm[0] = new_parm[0];
|
|
parm[1] = new_parm[1];
|
|
parm[2] = new_parm[2];
|
|
parm[3] = new_parm[3];
|
|
parm[4] = new_parm[4];
|
|
parm[5] = new_parm[5];
|
|
parm[6] = new_parm[6];
|
|
}
|
|
|
|
|
|
/*--------------------------------------------------------------------------*
|
|
* Random number generator *
|
|
*--------------------------------------------------------------------------*/
|
|
|
|
/* This function implements the linear congruential algorithm. Thanks to
|
|
gcc's long long ints, implementing 48-bit arithmetic (actually 64-bit,
|
|
truncating the result) is trivial. Limitations of long long int
|
|
implementation in (amiga?) gcc 2.7.0, made me use the kludge with the union
|
|
to convert from short[3] to long long int. (It's probably faster though!)
|
|
|
|
This function takes an array of three shorts (a 48-bit seed, m.s.word
|
|
last) and returns a long between -2**31 and 2**31-1, updating the seed
|
|
(the result is the two m.s.words of the updated seed).
|
|
*/
|
|
|
|
static long
|
|
rng(unsigned short *seed)
|
|
{
|
|
long long int Xn, Xn1, a, c;
|
|
union {
|
|
long long int l;
|
|
unsigned short s[4];
|
|
} i;
|
|
|
|
i.s[0] = 0;
|
|
i.s[1] = seed[2];
|
|
i.s[2] = seed[1];
|
|
i.s[3] = seed[0];
|
|
Xn = i.l;
|
|
|
|
i.s[0] = 0;
|
|
i.s[1] = parm[5];
|
|
i.s[2] = parm[4];
|
|
i.s[3] = parm[3];
|
|
a = i.l;
|
|
|
|
c = (long long int)(parm[6]);
|
|
|
|
Xn1 = a * Xn + c;
|
|
|
|
i.l = Xn1;
|
|
seed[0] = i.s[3];
|
|
seed[1] = i.s[2];
|
|
seed[2] = i.s[1];
|
|
|
|
return (long)((((unsigned long)seed[2]) << 16) + seed[1]);
|
|
}
|
|
|
|
|
|
/*--------------------------------------------------------------------------*
|
|
* Interface functions to the random number generator *
|
|
*--------------------------------------------------------------------------*/
|
|
|
|
/* This function returns a long between 0 and 2**31-1 by calling rng
|
|
with the internal seed, returning the 31 most significant bits of the
|
|
result shifted by one position to the right.
|
|
*/
|
|
long
|
|
lrand48(void)
|
|
{
|
|
return (rng(parm) >> 1) & 0x7FFFFFFF;
|
|
}
|
|
|
|
/* Same as lrand48(), but using an external seed. */
|
|
|
|
long
|
|
nrand48(unsigned short seed[3])
|
|
{
|
|
return (rng(seed) >> 1) & 0x7FFFFFFF;
|
|
}
|
|
|
|
/* This function returns a long between -2**31 and 2**31-1 by calling rng
|
|
with the internal seed.
|
|
*/
|
|
long
|
|
mrand48(void)
|
|
{
|
|
return rng(parm);
|
|
}
|
|
|
|
/* Same as mrand48(), but using an external seed. */
|
|
|
|
long
|
|
jrand48(unsigned short seed[3])
|
|
{
|
|
return rng(seed);
|
|
}
|
|
|
|
/* This function returns a double in the interval [0,1) by calling mrand48()
|
|
and dividing the result by ULONG_MAX + EPSILON. */
|
|
|
|
double
|
|
drand48(void)
|
|
{
|
|
union {
|
|
long l;
|
|
unsigned long u;
|
|
} x;
|
|
|
|
x.l = mrand48();
|
|
return (double)x.u / ((double)(ULONG_MAX) + EPSILON);
|
|
}
|
|
|
|
/* Same as drand48(), but using an external seed. */
|
|
|
|
double
|
|
erand48(unsigned short seed[3])
|
|
{
|
|
union {
|
|
long l;
|
|
unsigned long u;
|
|
} x;
|
|
|
|
x.l = nrand48(seed);
|
|
return (double)x.u / ((double)(ULONG_MAX) + EPSILON);
|
|
}
|