Skip site navigation (1)Skip section navigation (2)
Date:      Sun, 6 Oct 1996 04:56:27 +1000
From:      Bruce Evans <bde@zeta.org.au>
To:        ache@nagual.ru, bde@zeta.org.au
Cc:        current@freebsd.org, freebsd-hackers@freebsd.org, joerg_wunsch@uriah.heep.sax.de
Subject:   Re: I plan to change random() for -current (was Re: rand() and random())
Message-ID:  <199610051856.EAA30329@godzilla.zeta.org.au>

next in thread | raw e-mail | index | archive | help
>> ISO C example:    next = 1103515245 * next         + 12345;
>>                   return (unsigned int)(next / 65536) % 32768;
>> 
>> I.e., it returns bits 16-31 of the current state (right shifted 16).  This
>> is said to be better.  Folklore says that someone broke rand() by not
>> discarding the low bits when ints became 32 bits.
>
>It is not good enough to live due to unrandom nature of original formulae.
>I.e. _all_ bits are unrandom, and lower ones are more unrandom.
>This formulae not good enough even for 16bits.

No, it is probably fine for 16 bits (see Knuth), but probably no good for
random() since random() wants to use all the bits.

>At this moment I worry about random() only, lets consider rand() things
>after it.

OK.

>> This method is also found in the BSD4.4Lite libkern/rand.c.  I guess it
>> can be trusted (as much as the BSD rand.c :-).
>
>So we can safely import it even not from this posting (which can be
>GNU diseased) but from our own libkern/rand.c

Not so safe.  Something might depend on the doemented behaviour that the
sequence after srandom() depends on the seed (not even on the library
version :-).

>Here proposed patch, initial table regenerated to conform default
>seed value (1):

Commit it and see if we get complaints.

!  *      initstate(1, randtbl, 128);

Lost a tab here.

>! 		for (i = 1; i < rand_deg; i++) {
>! 	/*
>! 	 * Compute state[i] = (7^5 * state[i - 1]) mod (2^31 - 1)
>! 	 * wihout overflowing 31 bits:
>! 	 *              (2^31 - 1) = 127773 * (7^5) + 2836
>! 	 * From "Random number generators: good ones are hard to find",
>! 	 * Park and Miller, Communications of the ACM, vol. 31, no. 10,
>! 	 * October 1988, p. 1195.
>! 	 */
>! 			long hi, lo, t;
>! 
>! 			hi = state[i - 1] / 127773;
>! 			lo = state[i - 1] % 127773;
>! 			t = 16807 * lo - 2836 * hi;

This is much slower than the current rand() or random().  On a P133,
measured over 10^7 calls to statically linked libraries:

	rand:                174 nsec
	random:              275 nsec
	random from libkern: 673 nsec

The division in the libkern random() is apparently very expensive,
although gcc optimizes it nicely.  Knuth gives an algorithm for dividing
by (2^n + 1) using only one multiplication, one subtraction, one test
and 2 additions.  The method for (2^n - 1) is an exercise :-).  I think
these methods depend on a double width multiplication, which the i386 has.
The portable C method is slow because it must work without overflowing
32 bits.  We can write the Knuth algorithm in unportable gCc using long
long.  Gcc optimizes multiplications of the form `(long long)x * y'
(where x and y are 32 bit) nicely on the i386.

! 	/*
! 	 * Compute i = (7^5 * state[0]) mod (2^31 - 1)

Perhaps make this a subroutine to avoid duplication and allow easy
changing.  Division takes 10-20 times longer than a subroutine call
on Pentiums.

Bruce



Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?199610051856.EAA30329>