From:	CSBVAX::CSBVAX::MRGATE::"SMTP::CRVAX.SRI.COM::RELAY-INFO-VAX"  7-MAR-1989 00:02
To:	MRGATE::"ARISIA::EVERHART"
Subj:	Random Number Generators


Received: From KL.SRI.COM by CRVAX.SRI.COM with TCP; Mon,  6 MAR 89 14:04:18 PDT
Received: from circus.llnl.gov by KL.SRI.COM with TCP; Mon, 6 Mar 89 13:48:09 PST
Date: Mon, 6 Mar 89 13:41 PST
From: "Vladimir Ivanovic 415.423.7786" <IVANOVIC%VAXR@circus.llnl.gov>
Subject: Random Number Generators
To: info-vax@kl.sri.COM
X-VMS-To: INFO-VAX,IVANOVIC

    There has been some recent traffic about random number generators.
    ( If this posting is a repetition, it's due to the fact I'm 200+
    messages behind! Sorry. )
    
    After reading the two articles mentioned below, I would NEVER use
    anything except one of the random number generators recommended
    in those articles.  The benefit of using the ones presented here
    is that their period is greater than the "standard" generator presented
    in the Park and Miller article.
    
--------------------------------------------------------------------------------
unit RandomNumberGenerators;

{ Two random number generators are used here.  Uniform16 generates 16-bit
random numbers, while Uniform32 generates 32-bit numbers.  Init16 and Init32
initialize the appropriate seeds for the two generators. } 

{ The algorithms for Uniform16 and Uniform32 are implementations taken from the
excellent article "Efficient and Portable Combined Random Number Generators" by
Pierre L'Ecuyer in the June 1988 ( v 31 #6 ) issue of the Communications of the
ACM . Uniform16 has a period of ~6.95E12, while Uniform32 has a period of 
~2.30584E18.  Both pass numerous, stringent tests for randomness with flying
colors.  Another recent article by Stephen Park and Keith Miller, "Random
Number Generators: Good Ones are Hard to Find", in the October 1988 issue of
the Communications of the ACM is also highly recommended. } 

interface

	type
	{ Not strictly needed, but included to make things obvious. }
		int16 = integer;  {        -32,768 .. 32,767        }
		int32 = longint;  { -2,147,483,648 .. 2,147,483,647 }

	var
		s1, s2: int32;			{ seeds for the 32-bit version }
		s3, s4, s5: int16;		{ seeds for the 16-bit version }

	{ These procedures ask the user to initialize the seeds. }
	procedure Init16;
	procedure Init32;

	{ These two functions return a new random number each time they're called. }
	function Uniform16: real;
	function Uniform32: real;


implementation

	procedure Init16;
	{ Initializes the three seeds needed for the 16-bit version of Uniform. }
	{ The bounds for initializing the seeds are taken from L'Ecuyer's article. }
		var
			ok: boolean;
	begin
		ok := false;
		repeat
			writeln('Enter an integer in the range [ 1, 32362 ] for seed #1.');
			write('>  ');
			readln(s3);
			writeln;
			if (s3 >= 1) and (s3 <= 32362) then ok := true;
		until ok;

		ok := false;
		repeat
			writeln('Enter an integer in the range [ 1, 31726 ] for seed #2.');
			write('>  ');
			readln(s4);
			writeln;
			if (s4 >= 1) and (s4 <= 31726) then ok := true;
		until ok;

		ok := false;
		repeat
			writeln('Enter an integer in the range [ 1, 31656 ] for seed #3.');
			write('>  ');
			readln(s5);
			writeln;
			if (s5 >= 1) and (s5 <= 31726) then ok := true;
		until ok;

	end;


	procedure Init32;
	{ Initializes the two seeds needed for the 32-bit version of Uniform. }
	{ The bounds for initializing the seeds are taken from L'Ecuyer's article. }
		var
			ok: boolean;
	begin
		ok := false;
		repeat
			writeln('Enter an integer in the range [ 1 .. 2,147,483,562 ] for seed #1.');
			write('>  ');
			readln(s1);
			writeln;
			if (s1 >= 1) and (s1 <= 2147483562) then ok := true;
		until ok;

		ok := false;
		repeat
			writeln('Enter an integer in the range [ 1 .. 2,147,483,398 ] for seed #2.');
			write('>  ');
			readln(s2);
			writeln;
			if (s2 >= 1) and (s2 <= 2147483398) then ok := true;
		until ok;

	end;


	function Uniform16: real;
	{ See the L'Ecuyer article for an informative discussion of why things 
	are done this way. } 
		var
			z, k: int16;
	begin
		k := s3 div 206;
		s3 := 157 * (s3 - k * 206) - k * 21;
		if s3 < 0 then s3 := s3 + 32363;

		k := s4 div 217;
		s4 := 146 * (s4 - k * 217) - k * 45;
		if s4 < 0 then s4 := s4 + 31727;

		k := s5 div 222;
		s5 := 142 * (s5 - k * 222) - k * 133;
		if s5 < 0 then s5 := s5 + 31657;

		z := s3 - s4;
		if z > 706 then z := z - 32362;
		z := z + s5;
		if z < 1 then z := z + 32362;

		Uniform16 := z * 3.0899E-5;
	end;


	function Uniform32: real;
	{ See the L'Ecuyer article for an informative discussion of why things are
    done this way. }
		var
			z, k: int32;
	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;

		Uniform32 := z * 4.656613e-10;
	end;

end.
   

