looking for uniformly distributed (0..1] random number generator
Jesse Chisholm
jesse at altos86.Altos.COM
Sat Feb 2 10:56:06 AEST 1991
Yes. Attached is a module originally written for DOS but
tweaked for UNIX that generates nicely linear uniform numbers
with a period of 2147483646. The math is done in 32-bit ints
in such a way as to avoid the overflows. It only work for
modulo, not division (sigh).
It is a work in progress, but is functional as is.
I'm sure you can convert this long to a double in the range you
want.
========== "A witty saying proves nothing." -- Voltaire ==========
jesse at Altos86.Altos.COM Tel:1-408-432-6200x4810 Fax:1-408-434-0273
--- cut here ---------------------------------------------------------
/*
** This program originally written by Jesse Chisholm 90.08.08
** for use in the DOS world, MicroSoft C 5.1, i386.
**
** Tweaked to compile under unix cc because someone asked
** on the usenet for a linear random number generator with
** a decent period.
**
** This module is in the PUBLIC DOMAIN as of 91.01.01!
** It may be transmitted via any known means, and used
** for any purpose that is not illegal, immoral, or fattening.
*/
#ifdef TEST
#include <stdio.h>
#endif /* TEST */
/*
** added because unix didn't seem to have it
*/
typedef struct {
long quot;
long rem;
} ldiv_t;
/*
** Prime Multiplicative Congruential Generator
**
** Prime == 2147483647
** Standard root == 16807
** Equation new = old * root mod prime;
**
** For alternate sequences:
**
** 1st Choice Primitive root == 48271
** 2nd Choice Primitive root == 69621
**
*/
#define PRIME 2147483647L
#define PROOT 16807L /* This gives the minimum standard */
#define PRQUO 127773L /* PRIME / PROOT */
#define PRREM 2836L /* PRIME % PROOT */
#define PRCHK 1043618065L /* 10000th seed after 1L */
#ifdef NOTSTANDARD
#define PROOT1 48271L /* one of 16807, 48271, 397204094 */
#define PRQUO1 44488L /* PRIME / PROOT1 */
#define PRREM1 3399L /* PRIME % PROOT1 */
#define PRCHK1 399268537L /* 10000th seed after 1L */
#define PROOT2 69621L /* one of 39373, 69621, 630360016 */
#define PRQUO2 30845L /* PRIME / PROOT2 */
#define PRREM2 23902L /* PRIME % PROOT2 */
#define PRCHK2 190055451L /* 10000th seed after 1L */
#endif /* NOTSTANDARD */
#define NORMAL 4 /* this many uniforms at a time */
/*
**
** Using only the primitive root (16807) and the prime (2147483647) gives
** a sequence of psuedo random numbers in the range 1..2147483646 that
** is exactly 2147483646 numbers long before sequence repeat.
**
** Entry points:
**
** void initran(long prime, long root); use alternate numbers
** long _rand(); the basics
** long setran(long seed); start sequence
** long urand(long lo, long hi); uniform in range
** long nrand(long lo, long hi); normal in range
**
*/
/*
** static storage for running random number.
*/
static long _l_prim_ = PRIME;
static long _l_root_ = PROOT;
static long _i_seed_ = 1L; /* 1..prime-1 */
static long _i_quot_ = PRQUO; /* prime div root */
static long _i_rmnd_ = PRREM; /* prime mod root */
/*
** provided because unix didn't seem to have this routine
** this not as efficient as assembler would have been
*/
ldiv_t ldiv(n, d)
long n, d;
{
ldiv_t t;
t.quot = n / d;
t.rem = n % d;
return(t);
}
/*
** initran(prime, root)
**
** define which prime and root are to be used for the random number
** generation from this seed on.
*/
void
initran(prime, root)
long prime, root;
{
ldiv_t temp;
if ((_l_prim_ != prime) || (_l_root_ != root)) {
temp = ldiv(_l_prim_ = prime, _l_root_ = root);
_i_quot_ = temp.quot;
_i_rmnd_ = temp.rem;
}
}
/*
** setran(seed)
**
** use the argument as the start of the sequence
**
** if seed == 0, then use time of day as seed
*/
long
setran(seed)
long seed;
{
long time();
for(_i_seed_ = seed;
_i_seed_ == 0L;
_i_seed_ = time((long*)0) % _l_prim_)
continue;
_i_seed_ &= 0x7FFFFFFF; /* mask off sign bit, if any */
return(_i_seed_);
}
/*
** _rand()
**
** returns UNIFORM random numbers in the range 1..2147483646
**
** This is the basic routine.
**
** new = old * root % prime;
**
** The table below shows some convenient values
**
** prime primitive root alternate root range
**
** 17 5 6 1..16
** 257 19 17 1..256
** 32749 182 128 1..32748
** 65537 32749 32719 1..65536
** 2147483647 16807 39373 1..2147483646
** 2147483647 48271 69621 1..2147483646
** 2147483647 397204094 630360016 1..2147483646
**
** This is the best single primitive root for 2^31-1 according to
** Pierre L'Ecuyer; "Random Numbers for Simulation"; CACM 10/90.
**
** 2147483647 41358 1..2147483646
**
**
** He also lists the prime 4611685301167870637
** and the primitive root 1968402271571654650
** as being pretty good.
**
** There are many other primitive roots available for these
** or other primes. These were chosen as esthetically pleasing.
** Also chosen so that `x*y % p != 1' for roots x and y.
**
** A number (x) is a primitive root of a prime (p)
** iff in the equation:
**
** t = x^i % p
**
** when (i == p-1) and (t == 1)
**
** for all (1 <= i < p-1) then (t >= 2)
**
** The number of primitive roots for a given prime (p) is the number
** of numbers less that (p-1) that are relatively prime to (p-1).
**
** Note: that is how many, not which ones.
**
** if x is a primitive root of p, and y is relatively prime to p-1
** then (x^y % p) is also a primitive root of p.
** also (x^(p-1-y) % p) is a primitive root of p.
**
** The basic algorythm:
** new = old * root % prime;
** has the drawback that intermediate values are larger than the
** 4 byte register size. So an alternate way of implementing
** this expression is used so that no intermediate value exceeds
** the register size.
*/
#ifdef LMULDIV
/***** _lmuldiv isn't working. sigh. needs more thought *****/
/*
** primitive. return.quot == a*b/c return.rem == a*b%c
*/
ldiv_t
_lmuldiv(a, b, c)
long a, b, c;
{
ldiv_t t, answer;
long alpha, beta, gamma, delta;
long cq, cr;
alpha = a * b;
answer = ldiv(alpha, c);
return(answer);
if ((a == 0L) || (b == 0L)) {
answer.quot = answer.rem = 0L;
printf("\n%ld * %ld / %ld == %ld .+. %ld",
a, b, c, answer.quot, answer.rem);
return(answer);
}
if (c == 0L) {
answer.quot = 0L;
answer.rem = -1L;
printf("\n%ld * %ld / %ld == %ld .+. %ld",
a, b, c, answer.quot, answer.rem);
return(answer);
}
t = ldiv(c, a);
cq = t.quot;
cr = t.rem;
t = ldiv(b, cq);
alpha = a * t.rem;
beta = cr * t.quot;
gamma = alpha - beta;
if (gamma > 0L) {
answer.rem = gamma;
answer.quot = 0L;
} else {
if (cr < cq) {
answer.rem = gamma + c;
answer.quot = 1L;
} else {
delta = ((-gamma / c) + 1L);
answer.rem = gamma + c*delta;
answer.quot = t.quot - delta;
}
}
printf("\n%ld * %ld / %ld == %ld .+. %ld",
a, b, c, answer.quot, answer.rem);
return(answer);
}
#endif /* LMULDIV */
long
_rand()
{
long test;
ldiv_t temp;
long alpha, beta, delta;
if (_i_seed_ == 0L)
setran(0L);
/*
** implement multiplication and modulo in such a way that
** intermediate values all fit in a long int.
*/
/* _i_seed_ = _i_seed_ * PROOT % PRIME */
temp = ldiv(_i_seed_, _i_quot_);
alpha = _l_root_ * temp.rem;
beta = _i_rmnd_ * temp.quot;
/*
** normalize the intermediate values
*/
if (alpha > beta) {
_i_seed_ = alpha - beta;
} else {
if (_i_rmnd_ < _i_quot_) {
_i_seed_ = alpha - beta + _l_prim_;
} else {
/*
** TODO: implement (a*z div m) in such a way that intermediate
** values fit in long int.
*/
/*
** delta = temp.quot - ((_l_root_ * _i_seed_) / _l_prim_);
*/
/*
** delta = temp.quot - _lmuldiv(_l_root_, _i_seed_, _l_prim_);
*/
delta = (((beta - alpha) / _l_prim_) + 1L);
_i_seed_ = alpha - beta + _l_prim_*delta;
}
}
return( _i_seed_ );
}
/*
** returns UNIFORM random numbers in the range
*/
int urand(lo, hi)
int lo, hi;
{
if (lo > hi)
return(urand(hi,lo));
if (lo == hi)
return(lo);
return( (int)((_rand() % ((long)(hi - lo + 1))) + (long)lo) );
}
/*
** returns NORMAL random numbers in the range
*/
int nrand(lo, hi)
int lo, hi;
{
long t;
int i;
for(t=i=0; i<NORMAL; i++)
t += (long)urand(lo,hi);
i = (int)((t + (NORMAL/2)) / NORMAL);
return(i);
}
#ifdef TEST
main()
{
int num;
long dd, ee;
ldiv_t t, u, v;
printf("Testing validity of arithmatic in this implementation\n");
printf("of the MINIMAL STANDARD RANDOM NUMBER GENERATOR\n\n");
initran(PRIME, PROOT);
setran(1L);
ee = 1L;
for(num=0; num<10000; num++) {
if ((num % 100) == 0) {
printf("\rIteration %d /10000", num);
}
#ifdef LMULDIV
t = _lmuldiv(PROOT, ee, PRIME);
if (t.quot == 0L) {
if (t.rem != (long)((unsigned long)PROOT*(unsigned long)ee)) {
printf("\n(%ld*%ld)/%ld = %ld#%ld ??? %ld\n",
PROOT, ee, PRIME,
t.quot, t.rem, PROOT*ee);
}
} else if (t.rem == 0L) {
u = _lmuldiv(PRIME, t.quot, PROOT);
if ((u.quot != ee) || (u.rem != 0L)) {
printf("\n(%ld*%ld)/%ld = %ld#%ld ???\n",
PROOT, ee, PRIME, t.quot, t.rem);
printf(" (%ld*%ld)/%ld = %ld#%ld ???\n",
PRIME, t.quot, PROOT, u.quot, u.rem);
}
} else {
u = _lmuldiv(PRIME, t.quot, PROOT);
v = _lmuldiv(t.rem, 1L, PROOT);
if (((u.rem+v.rem) != PROOT)
|| ((u.quot+v.quot) != (ee-1L))) {
printf("\n(%ld*%ld)/%ld = %ld#%ld ???\n",
PROOT, ee, PRIME, t.quot, t.rem);
printf(" (%ld*%ld)/%ld = %ld#%ld ???\n",
PRIME, t.quot, PROOT, u.quot, u.rem);
printf(" (%ld*%ld)/%ld = %ld#%ld ???\n",
t.rem, 1L, PROOT, v.quot, v.rem);
}
}
ee = t.rem;
#endif /* LMULDIV */
dd = _rand();
#ifdef LMULDIV
if (ee != dd) {
printf("\rIteration %d ee=%ld, dd=%ld\n", num, ee, dd);
}
#endif /* LMULDIV */
}
printf("\rSeed should be %ld and is %ld.\n", PRCHK, dd);
printf("\n%s\n\n\n", PRCHK == dd ? "All is correct!" : "** ERROR **");
#ifdef LMULDIV
if (PRCHK != ee)
printf("**ERROR** in _lmuldiv, %ld != %ld\n", PRCHK, ee);
#endif /* LMULDIV */
}
#endif /* TEST */
--- cut here ---------------------------------------------------------
More information about the Comp.unix.programmer
mailing list