From: vladimir@Eng.Sun.COM (Vladimir Ivanovic)
Newsgroups: comp.programming
Subject: Re: Pseudo-Random number generators
Date: 15 Oct 92 23:42:20 GMT
Organization: Sun Microsystems, Inc.
NNTP-Posting-Host: cocteau
In-reply-to: pem@Minsk.DoCS.UU.SE's message of 15 Oct 1992 12:45:02 GMT

>>>>> On 15 Oct 1992 12:45:02 GMT, pem@Minsk.DoCS.UU.SE (Per-Erik Martin) said:
pem> Try these:

pem> Communications of the ACM, October 1988 (v31):
pem> "Random number generators: Good ones are hard to find"
pem> by Stephen K. Park and Keith W. Miller.

pem> Communications of the ACM, October 1991 (v33):
pem> "Random numbers for simulation"
pem> by Pierre L'Ecuyer.

pem> The first one should be read by everyone who deals with random number
pem> generator implementation, and it presents a simple but quite good
pem> generator.
pem> The second article is a second hand reference, I haven't read it myself.

Funny you should mention these... :-)

The June 1988 (v31 #6) issue of the Communications of the ACM has an article
by Pierre L'Ecuyer called, "Efficient and Portable Combined Random Number
Generators".  The following October (v31 #10), Stephen Park and Keith Miller
published "Random Number Generators: Good Ones are Hard to Find."

Both articles are ESSENTIAL reading for anyone interested in random number
generation, and this includes naive users.  The conclusion of both articles is
that horrible random (sic) number generators are used by people who don't know
better, while very good ones take less than 20 lines of Pascal!  Amoung the
horrid generators are some that come with certain systems or are presented in
textbooks! 

Here are the two proposed random number generators.  Even though I have proof
read them twice, I made no guarantees whatsoever about the correctness or the
usability of the code below.  

Here is the Minimal Standard proposed by Park and Miller:

	function Random : real;
	{ Initialize seed with 1..2147483646 }
	{ maxint must be greater than or equal to 2**31 - 1 }
	const
	  a = 16807;
	  m = 2147483647;
	  q = 12773;       { m div a }
	  r = 2836;	   { m mod a }
	var
	  lo, hi, test : integer;
	begin
	  hi := seed div q;
	  lo := seed mod q;
	  test := a * lo - r * hi;
	  if test > 0 then
	     seed := test
	  else
	     seed := test + m;
	 Random := seed / m
       end;

Here is the Portable Combined Generator of L'Ecuyer for 32-bit computers.  It
has a period of roughly 8.12544 x 10**12.

	Function Uniform : real;
	{ Initialize s1 to 1..2147483562; s2 to 1..2147483398 }
	  var
	    Z, k : integer;
	  begin
	    k  := s1 div 53668;
	    s1 := 40014 * (s1 - k * 53668) - k * 12211;
	    if s1 < 0 then s1 := s1 + 2147483563;
	    
	    k  := s2 div 52774;
	    s2 := 40692 * (s2 - k * 52774) - k * 3791;
	    if s2 < 0 then s2 := s2 + 2147483399;

	    Z := s1 - s2;
	    if Z < 1 then Z := Z + 2147483562;

	    Uniform := Z * 4.656613e-10;
	  end;

All are urged to read both articles for a much better presentation.  And you
don't have to be a mathematician to understand them: very little of either
article is not accessible to a competent programmer.

Enjoy!

-- Vladimir

--
Vladimir G. Ivanovic                            Sun Microsystems, Inc
(415) 336-2315                                  MTV12-33
vladimir@Eng.Sun.COM                            2550 Garcia Ave.
{decwrl,hplabs,ucbvax}!sun!Eng!vladimir         Mountain View, CA 94043-1100
                         Disclaimer: I speak for myself.
Newsgroups: comp.programming,comp.lang.c,comp.lang.pascal,comp.lang.objective-c