[net.bugs.4bsd] Bug in random number generator

kar (03/15/83)

	The random number generator that came with our 4.1BSD includes the
following line of code (that does the work):

	return((randx = randx * 1103515245 + 12345) & 0x7fffffff);

and the generator on our old V7 system includes this line of code:

	return((randx = randx * 1103515245 + 12345)>>16 & 0x7fff);

The V7 version produced random numbers in the range 0-32767, and the new ver-
sion produces numbers in the range 0-2147483647.

	Unfortunately, the numbers produced by the new version are NOT random.
To verify my claim, try the following program:

	for(;;)
		printf( "%d\n", rand() % 2 );

and see what happens.  If you mod your "random" numbers by any power of 2,
you'll be disappointed with the results.

	The fix?  Go back to the V7 generator.  The algorithm used produces
numbers that are not very random in the low order bits; this is the real
reason for the shift and mask to 16 bits, not the 16 bit architecture of the
PDP-11.

	- Ken Reek, Rochester Institute of Technology
	ucbvax!allegra!rochester!ritcv!kar

lepreau (03/24/83)

No question, the 4.1 rand() is bogus.  If people (Berkeley especially
and Ken Arnold of course) would throw out it out and re-link with a
decent one the nature of a lot of games might change.

Here's a new rand() for the Vax that was writen by the Chris Kingsley,
the same fellow who wrote the super-hot malloc() that Spence posted. I
can't vouch for it beyond that, but he sounds authoritative, and
it couldn't possibly be worse than the current one.  There's even a man
page.  Now if it would only replace the current one in 4.2...
------------------------------------------
>From cithep!citcsv!kingsley@Berkeley Thu Mar  4 20:47:10 1982
Mail-from: ARPANET site UCB-C70 rcvd 4-Mar-82 2228-MST
Date: 4 Mar 1982 20:54:18-PST
From: cithep!citcsv!kingsley at Berkeley
To: cithep!ucbvax!4bsd-bugs@Berkeley
Subject: random number generation.

I have done a little work with rand supplied with the system and I have 
discovered that it is flawed.  The manual page claims that it has a 
period of 2^32 and returns numbers from 0 to 2^31-1.  The code makes it
look like the author thought it was correct, but it is not.  Instead of
masking out the most significant (and also most random) bit, you should
do an unsigned shift to throw out the least significant (and least random)
bit.  I have also found a multiplier that passes Knuth's spectral test
very well.
      I have written a new rand, along with randint(n) which returns 0 to
n-1, and flat() which returns 0.0 to <1.0.  I did it in assembler (Mea
Maxima Culpa!) to use the extended multiply and some bit fiddling.  Would
you like to see it?

Chris Kingsley


>From cithep!citcsv!kingsley@Berkeley Fri Mar  5 20:47:10 1982
Mail-from: ARPANET site UCB-C70 rcvd 5-Mar-82 2037-MST
Date: 5 Mar 1982 19:33:34-PST
From: cithep!citcsv!kingsley at Berkeley
To: cithep!ucbvax!4bsd-bugs@Berkeley, cithep!ucbvax!cithep!citcsv!kingsley@Berkeley

Subject: Re: random number generation.

Yes, I realize that the bottom bits aren't random.  In fact, the bottom
n bits have a period of 2^n.  The rng delivered, though, throws out the
most significant bit to produce a 31 bit number, and claims that it has
a period of 2^32.

For your interest, here's my rand.s (dead languages!),
..data
..align  2
_randx:.long    1
..text
..align      1

..globl  _srand              # set the random seed
_srand: .word   0           # srand(seed) int seed;
	movl    _randx,r0
	movl    4(ap),_randx
	ret
..align  1

..globl  _rand               # give a 31 bit random positive integer
_rand:  .word   0           # int rand()
	jsb     rand
	bicl2   $1,r0
	rotl    $-1,r0,r0
	ret

..globl  _randint            # give a random positive integer from 0 to n-1
_randint:.word  0xc         # int randint(n) int n;
	jsb     rand
	emul    4(ap),r0,$0,r2
	tstl    r0
	jgeq    L1
	addl3   4(ap),r3,r0
	jbr     retins
L1:     movl    r3,r0
retins: ret
..align  1

# compute the next 32 bit random number
rand:   mull3   $505360173,_randx,r0
	addl2   $907633385,r0
	movl    r0,_randx
	rsb
..align 1

..globl  _flat               # give a random double fro 0. to <1.
_flat:  .word   0xc         # double flat()
	jsb     rand
	movl    r0,r2
	movf    $0f1.0,r0
	extzv   $25,$7,r2,r3
	insv    r3,$0,$7,r0
	extzv   $9,$16,r2,r3
	insv    r3,$16,$16,r0
	extzv   $0,$9,r2,r1
	subd2   $0d1.0,r0
	ret


(the actual generator is the routine rand, the global routines just 
do range conversion)
and the man page for it,

..TH RAND 3 "Caltech Homegrown"
..SH NAME
srand, rand, randint, flat \- random number generator
..SH SYNOPSIS
..nf
..B srand(seed)
..B int seed;
..PP
..B rand()
..PP
..B randint(n)
..B int n;
..PP
..B double flat()
..fi
..SH DESCRIPTION
These subroutines use a
multiplicative congruential
random number generator
with period
..if t 2\u\s732\s0\d
..if n 2^32
to return successive pseudo-random
numbers.  The numbers range over
..if t 2\u\s731\s10\d\-1.
..if n 2^31-1
for
..I rand ,
0 to n-1 for randint, and 0 to <1.0 for flat.
..PP
..I Srand
sets the random seed and returns its previous value.
..PP
The generator has been tested with the spectral test in Knuth.  That test 
checks for independance of n-tuples produced by a multiplicative congruential 
generator.
The results are:
pairs of numbers are independant to one part in 65742; 3 numbers,
one part in 1580; 4 numbers, one part in 247; 5 numbers, one part in 71.9;
and 6 numbers, one part in 37.8.  (The code for the test is in
/usr/local/src/spectr.bc)


Enjoy!

Chris Kingsley

smb (03/24/83)

As Knuth points out (Art of Computer Programming, Vol. 2 -- Seminumerical
Algorithms) the low-order bits of *any* linear congruential random number
generator are not particularly random.  The proper way to view the output
is as a fraction between 0 and 1.  Thus, to get a random bit, one should
never take rand()%2; rather, one should use rand() > (MAXINT/2) (or some
such).

Of course, none of this means that rand() is a good random number generator --
it's not.


		--Steve

hall (03/25/83)

#R:ritcv:-19100:uiucdcs:8200012:000:529
uiucdcs!hall    Mar 24 21:00:00 1983

I suggest replacing the random number generator with one based
on the linear congruential method described in Knuth, "The Art
Of Computer Programming, Vol 2: Seminumerical Algorithms".
There is an EXCELLENT article in the Bell System Technical
Journal, Vol. 61, No. 8, October, 1982 by C.S. Roberts that
describes just what you want. The appropriate constants are
listed for a generator with period 2**48; test data is also
tabulated to check your implementation. I recommend this
article highly.
--John
(...pur-ee!uiucdcs!hall)

mark (03/25/83)

It is curious that V7's rand produces 15 bit random integers,
but 32V (and hence 4.1BSD) produce 31 bit numbers.  This makes
it interesting trying to write portable programs.  USG UNIX
rand's produce 15 bit numbers, so the proper direction does
indeed seem to be toward 15 bits.

Let me add my two cents worth to the comment about the low
order bits.  It's certainly convenient to use the % operator
to get a random integer in a given range.  But it's simply
not very random.  Here is a macro that will do the job
much more randomly on 4.1BSD (and, presumably, on other UNIX
systems)

#define randint(n) ( ((long) rand() & 0177777) * n / 0177777)
n is the range you want, e.g. randint(10) will produce a random
integer in the range 0 through 9.