Bug in random number generator

[utzoo!decvax!harpo!utah-cs!lepreau]

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 at 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 at 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 at 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 at Berkeley, cithep!ucbvax!cithep!citcsv!kingsley at 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