morgan@ogicse.cse.ogi.edu (Clark O. Morgan) (01/29/91)
Hello. I'm currently using a Sequent Symmetry (3.17) executing BSD4.2 :-) . I need access to a uniformly distributed random number generator over the range (0..1] . The RNG should have a reasonably long period. I found a description of a Linear Congruential method that works fine in theory, but suffers from unsigned long overflow in practice (unless, of course, I implement the RNG is ASM :-( ). Anybody have such a generator that's written in C (can't use Fortran, our machine has no f77 compiler)? Or have pointers to books/archives that I should examine? Thanks in advance.
moss@cs.umass.edu (Eliot Moss) (01/29/91)
There are lots of things to consider in selecting a good random number generator. However, there are a couple of good CACM articles that really helped us, and they can help you, too. Here's the ref to the one we use: "Efficient and Portable Combined Random Number Generators", Pierre L'Ecuyer, Communications of the ACM, Vol. 31, No. 6 (June 1989). -- J. Eliot B. Moss, Assistant Professor Department of Computer and Information Science Lederle Graduate Research Center University of Massachusetts Amherst, MA 01003 (413) 545-4206, 545-1249 (fax); Moss@cs.umass.edu
jesse@altos86.Altos.COM ( Jesse Chisholm) (02/02/91)
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@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 ---------------------------------------------------------moss@cs.umass.edu (Eliot Moss) (02/04/91)
... and here's our version, based on the same source (Note: a previous posting
of mine had some of the citation information wrong, which I have now corrected
below) .... Eliot
------------------------------------------------------------------------
/*
$RCSfile: random.h,v $
$Revision: 1.3 $
$Date: 90/01/21 14:08:34 $
$Author: moss $
$State: Exp $
$Locker: moss $
*/
/* header file for random.c */
# ifndef __RANDOM__
# define __RANDOM__
# ifdef __STDC__
extern double Random(void); /* Random Generator */
extern int RandomInt(int modulus); /* Returns integer 0 <= r < modulus */
extern int TimeRandomize(void); /* Assign seeds using the system time */
extern int SetRandomSeeds(int s1,int s2); /* Assign seeds. */
# else
extern double Random(); /* Random Generator */
extern int RandomInt(); /* Returns integer 0 <= r < modulus */
extern int TimeRandomize(); /* Assign seeds using the system time */
extern int SetRandomSeeds(); /* Assign seeds. */
# endif
# endif
------------------------------------------------------------------------
/*
$RCSfile: random.c,v $
$Revision: 1.3 $
$Date: 90/01/19 06:42:14 $
$Author: nayeri $
$State: Exp $
$Locker: moss $
*/
/*
*
* random.h Random Generator routines for 32-bit machines.
* Author Farshad Nayeri
* Date July 15, 1989
* Description These are a set of library routines for
* random generation. With the current constants
* the period is about 2E18.
* Disclaimer
* The following routines are correct to the best
* of my knowledge, however I will not responsible
* for direct or indirect implications of using these
* routines in case of any error.
* See Also
* "Efficient and Portable Combined Random Number
* Generators" by Pierre L'Ecuyer, Comm. of ACM,
* Vol 31, Number 6, June 1988.
*
*/
#ifdef VMS
# include <time.h>
#else
# include <sys/types.h>
# include <sys/timeb.h>
#endif
#define NORMAL 0
#define ERROR (-1) /* TimeRandomize depeneds on this value */
#define MAX(s1,s2) ((s1) > (s2) ? (s1) : (s2))
#define MIN(s1,s2) ((s1) < (s2) ? (s1) : (s2))
#define MINSEED1 1
#define MAXSEED1 2147483562
#define DIVISOR1 53668 /* constants for random */
#define SEED1KMULT 40014 /* generator */
#define K1MULT 12211
#define MINSEED2 1
#define DIVISOR2 52274
#define MAXSEED2 2147483398
#define SEED2KMULT 40692
#define K2MULT 3791
#define TIMEFACTOR 10000
double normalizer = (1.0 / ((double)MAX(MAXSEED1,MAXSEED2)));
/* normalization factor */
double maxSeed1and2 = MAX(MAXSEED1,MAXSEED2);
static int seed1 = MAXSEED1 / 2 , /* first and second seed*/
seed2 = MAXSEED2 / 2;
/*
*
* DESCRIPTION
* Random will return a random double precision real value
* between 0.0 and 1.0 (See Efficient and Portable Combined
* Random Number Generators by Pierre L'Ecuyer, Comm. of ACM,
* Vol 31, Number 6, June 1988. )
* NOTES
* This is a 32-bit generator.
*/
# ifdef __STDC__
double Random(void)
# else
double Random()
# endif
{
int z,k;
k = seed1 / DIVISOR1;
seed1 = SEED1KMULT * ( seed1 - k * DIVISOR1) - k * K1MULT;
if (seed1 < 0) seed1 += MAXSEED1 + 1;
/* seed1 = seed1 + MAXSEED1 + 1 */
k = seed2 / DIVISOR2;
seed2 = SEED2KMULT * ( seed2 - k * DIVISOR2) - k * K2MULT;
if (seed2 < 0) seed2 += MAXSEED2 + 1;
/* seed2 = seed2 + MAXSEED2 + 1 */
z = seed1 - seed2;
if (z < 1) z += maxSeed1and2;
return ((double) z * normalizer);
}
/*
*
* DESCRIPTION
* RandomInt returns a random integer between 0 and the modulus-1
* Since Random only returns values r > 0, 1-r < 1, so
* this function will always return values between 0 and modulus-1
*
*/
# ifdef __STDC__
int RandomInt(int modulus)
# else
int RandomInt(modulus)
int modulus;
# endif
{
return((int) ((1.0 - Random()) * modulus));
}
/*
*
* DESCRIPTION
* If s1 and s2 are valid, set the seeds, otherwise returns
* error and keep the old seeds.
* RETURN VALUE
* 0 if everything is O.K.
* -1 if there is an error.
*
*/
# ifdef __STDC__
int SetRandomSeeds(int s1,int s2)
# else
int SetRandomSeeds(s1,s2)
int s1;
int s2;
# endif
{
if ((s1>=MINSEED1) && (s1<=MAXSEED1))
seed1 = s1;
else {
return(ERROR);
}
if ((s2>=MINSEED2) && (s2<=MAXSEED2))
seed2 = s2;
else {
return(ERROR);
}
return(NORMAL);
}
/*
*
* DESCRIPTION
* Use the system time to set the seeds.
* RETURN VALUE
* 0 if everything is O.K.
* -1 if there is a problem with time routine.
*
*/
# ifdef __STDC__
int TimeRandomize(void)
# else
int TimeRandomize()
# endif
{
unsigned short theTime = 0;
struct timeb tp;
ftime(&tp);
theTime = tp.millitm / 10; /* millitm is always multiple of 10 */
if ((theTime & 001) == 0)
SetRandomSeeds((int)(MIN(seed1,seed2) + theTime * TIMEFACTOR)
% MAXSEED1 + MINSEED1,
(int)(MAX(seed1,seed2) % MAXSEED2 + MINSEED2));
else
SetRandomSeeds((int)(MAX(seed1,seed2) + theTime * TIMEFACTOR)
% MAXSEED1 + MINSEED1,
(int)(MIN(seed1,seed2) % MAXSEED2 + MINSEED2));
return(NORMAL);
}
------------------------------------------------------------------------
--
J. Eliot B. Moss, Assistant Professor
Department of Computer and Information Science
Lederle Graduate Research Center
University of Massachusetts
Amherst, MA 01003
(413) 545-4206, 545-1249 (fax); Moss@cs.umass.edulan_csse@netrix.nac.dec.com (CSSE LAN Test Account) (02/13/91)
You should also read the article in the Oct 1988 CACM (v.31, no.10) by Stephen K. Park and Keith W. Miller, titled "Random number generators: good ones are hard to find". Among other things, they point out that the common library routines on most systems use modular arithmetic with masking, which is not acceptable, and so they should be replaced. I've implemented their (quite simple) example of a minimally-acceptable uniformly-distributed 32-bit integer on everal systems. It's each, and it works quite well.