peters@cubsvax.UUCP (Peter S. Shenkin) (03/01/85)
In attempting to use these functions under 4.1bsd on a VAX/780, I found that
rand() returns strictly alternating even and odd numbers, regardless of
seed!!! So I went and tried the same thing on a 4.2bsd system, and found
the same behavior.
It may be that this has already been hashed out in some
of the above newsgroups, but since I only read net.lang.c among them I will
have missed the discussion. If anyone has a rewritten rand(), or other
advice, I'd appreciate being mailed a fix.
Thanks, Peter S. Shenkin, {philabs,cmcl2!rna}!cubsvax!petersmlm@cmu-cs-cad.ARPA (Michael Mauldin) (03/03/85)
[ generic apology: my mail reply to cubsvax!peters bounced back ]
Rand is a simple linear congruential pseudo-random number generator
(read your Knuth volume II). The low order bits of a linear
congruential sequence are not very random at all. The problem is not
rand, but rather your use of rand. You cannot safely use the result
directly, and you cannot safely use mod (%) to give a smaller range.
You must use div (/) to get the range, or at the very least throw away
the low bits before 'mod'ing.
The right way to use a LCPRNG to get integers frlom 0..n-1 is this:
result = floor (n * rand() / (MAXRAND+1));
Here is a valid call to the system rand():
(int) result = (int) range * ((double) rand () / 2147483648.0);
Here is the salient portion of the system rand:
return ((randx = randx * 1103515245 + 12345) & 0x7fffffff);
Here is my own random number generator which uses Knuth algorithm M and
a few other tricks to make more random numbers. It also contains a
randint(n) which correctly generates pseudo-random integers from 0..n-1.
/*****************************************************************
* rand: A very, very random generator, period approx 6.8064e16.
*
* Uses algorithm M, "Art of Computer Programming", Vol 2. 1969, D.E.Knuth.
*
* Two generators are used to derive the high and low parts of sequence X,
* and another for sequence Y. These three generators were derived by me.
*
* Usage: initialize by calling srand(seed), then rand() returns a random
* integer from 0..2147483647. srand(0) uses the current time as
* the seed. You must call srand() first, because it initializes
* an array which is used by the rand() function.
*
* Author: Michael Mauldin, June 14, 1983.
*****************************************************************/
/* Rand 1, period length 444674 */
# define MUL1 1156
# define OFF1 312342
# define MOD1 1334025
# define RAND1 (seed1=((seed1*MUL1+OFF1)%MOD1))
# define Y RAND1
/* Rand 2, period length 690709 */
# define MUL2 1366
# define OFF2 827291
# define MOD2 1519572
# define RAND2 (seed2=((seed2*MUL2+OFF2)%MOD2))
/* Rand 3, period length 221605 */
# define MUL3 1156
# define OFF3 198273
# define MOD3 1329657
# define RAND3 (seed3=((seed3*MUL3+OFF3)%MOD3))
/*
* RAND2 generates 19 random bits, RAND3 generates 17. The X sequence
* is made up off both, and thus has 31 random bits.
*/
# define X ((RAND2<<13 ^ RAND3>>3) & 017777777777)
# define AUXLEN 97
static int seed1=872978, seed2=518652, seed3=226543, auxtab[AUXLEN];
/****************************************************************
* srand: Initialize the seed (0 => use the current time+pid+uid)
****************************************************************/
srand (seed)
int seed;
{ register int i;
if (seed == 0) seed = (getuid() + getpid() + time(0));
/* Set the three random number seeds */
seed1 = (seed1+seed) % MOD1;
seed2 = (seed2+seed) % MOD2;
seed3 = (seed3+seed) % MOD3;
for (i=AUXLEN; i--; )
auxtab[i] = X;
}
/****************************************************************
* rand: Random integer from 0..2147483647
****************************************************************/
int rand ()
{ register int j, result;
j = AUXLEN * Y / MOD1; /* j random from 0..AUXLEN-1 */
result = auxtab[j];
auxtab[j] = X;
return (result);
}
/****************************************************************
* randint: Return a random integer in a given Range
****************************************************************/
int randint (range)
int range;
{ double rrand;
int result;
rrand = rand();
result = range * (rrand / 2147483648.);
return (result);
}
Share and Enjoy!
Michael L. Mauldin (Fuzzy) Department of Computer Science
Mauldin@CMU-CS-CAD.ARPA Carnegie-Mellon University
(412) 578-3065 Pittsburgh, PA 15213tim@callan.UUCP (Tim Smith) (03/05/85)
In article <320@cubsvax.UUCP> peters@cubsvax.UUCP (Peter S. Shenkin) writes: >In attempting to use these functions under 4.1bsd on a VAX/780, I found that >rand() returns strictly alternating even and odd numbers, regardless of >seed!!! So I went and tried the same thing on a 4.2bsd system, and found >the same behavior. > If BSD rand() is at all like System V rand(), this is not a bug. Rand() is a multiplicitive congruential (sp?) random number generator, which means it works like this: R = A * R + C mod M i+1 i So, depending on the parity of A and C, you have one of the following cases: Ri+1 is always odd Ri+1 is always even Ri+1 is always same parity as Ri or Ri+1 is always opposite parity of Ri In this sort of number generator, the high bits are "more random" than the low bits. Thus, for example, to get a stream or random bits, one would look at the top bit ( assuming M is a power of two... ), not the bottom bit. Look at Knuth, Vol. II for almost all that is known on this topic. -- Duty Now for the Future Tim Smith ihnp4!wlbr!callan!tim or ihnp4!cithep!tim
brooks@lll-crg.ARPA (Eugene D. Brooks III) (03/05/85)
> In attempting to use these functions under 4.1bsd on a VAX/780, I found that > rand() returns strictly alternating even and odd numbers, regardless of > seed!!! So I went and tried the same thing on a 4.2bsd system, and found > the same behavior. > > It may be that this has already been hashed out in some > of the above newsgroups, but since I only read net.lang.c among them I will > have missed the discussion. If anyone has a rewritten rand(), or other > advice, I'd appreciate being mailed a fix. > > Thanks, Peter S. Shenkin, {philabs,cmcl2!rna}!cubsvax!peters rand uses a 32bit multiplicative random number generator. These generators tend to generate a repeating sequence in the lower bits. On the pdp 11 the generator used a long and the returned value was shifted right by 16 bits to get an int. The horribly behaved lower bits never hit the user. On the vax the how 32 bit item is returned. As you noticed those badly behaved bits are there. If you want to use rand shift it down. Shift at least 8 bits out. A better way around the problem is to use the more sophisticated random number generator random() all of its bits are supposed to be random and its really not that much slower than rand(). A simple #define rand random can fix your code.
tim@callan.UUCP (Tim Smith) (03/05/85)
My 'M' in <307@callan.UUCP> is meant to be even, by the way... -- Duty Now for the Future Tim Smith ihnp4!wlbr!callan!tim or ihnp4!cithep!tim
guy@rlgvax.UUCP (Guy Harris) (03/06/85)
> rand uses a 32bit multiplicative random number generator. These generators > tend to generate a repeating sequence in the lower bits. On the pdp 11 > the generator used a long and the returned value was shifted right by 16 > bits to get an int. The horribly behaved lower bits never hit the user. > On the vax the how 32 bit item is returned. Well, to be specific, on the VAX under 4.xBSD the 32 bit item is returned. The System {III,V} systems still return the upper 16 bits. This means that any program that assumes "rand" returns something in the range 0 <= N < MAX16BITINT will work on all UNIX systems on machines with 16-bit "int"s and on all System N systems, but won't work right on 4.xBSD systems. Grumble grumble... -- Guy Harris {seismo,ihnp4,allegra}!rlgvax!guy
alexis@reed.UUCP (Alexis Dimitriadis) (03/09/85)
> A better way around the problem is to use the more > sophisticated random number generator random() all of its bits are supposed > to be random and its really not that much slower than rand(). A simple > #define rand random can fix your code. Our version of 4.2 random() ALWAYS returns an even number the first time it is called after being seeded. That effectively breaks small applications that call random only once. (Or those that seed it for every call). alexis
brooks@lll-crg.ARPA (Eugene D. Brooks III) (03/12/85)
> > A better way around the problem is to use the more > > sophisticated random number generator random() all of its bits are supposed > > to be random and its really not that much slower than rand(). A simple > > #define rand random can fix your code. > > Our version of 4.2 random() ALWAYS returns an even number the first time it > is called after being seeded. That effectively breaks small applications > that call random only once. (Or those that seed it for every call). > > alexis The context of my bug fix was for applications that use many calls to rand and don't want the pattern appearing in the low bits. The question of generating a single number for an application such as fortune is a different can of worms. I don't deal with applications like fortune and such but there is a reasonable way of doing them also. The trick is to get two seeds out of the environment somehow. You can use the time, the number of active jobs, the phase of the moon {don't have /dev/moon on your system? :-)} and as much other crap as you can hash together. Use one seed to seed the random number generator. Use the other to call the generator a few times. Then get your random number from it. Probably the only way to do better is with a shot noise generator. These don't generate lists.
jpk@ur-univax.UUCP (03/15/85)
The same bug in rand() appears on the VMS C compiler (VAX-11 C)! When I reported it to DEC, they had never discovered it. I happened to find it because I was using rand() to generate a random noise graphic, and discovered an awful lot of pattern in a supposedly randomly generated display. What happened was I made two successive calls to rand() to generate an x and y coordinate, so all my x's fell on odd spaces, all my y's on even. Half the points were never chosen on my screen, and a very regular half! The fix was to call MTH$RANDOM on VMS, and random() on unix. The moral: compatible function libraries don't support compatible systems when the original version was bugged! Instead, they support bugged, but compatible, code. Equally wrong on VMS or unix!!! -- Jon Krueger