randcard.mod

IMPLEMENTATION MODULE RandCard;

	(********************************************************)
	(*							*)
	(*		Random number generator			*)
	(*							*)
	(*	This version is a more portable (but slightly	*)
	(*	slower) version of the assembly language	*)
	(*	module of the same name.			*)
	(*							*)
	(*  Programmer:		P. Moylan			*)
	(*  Last edited:	18 March 1995			*)
	(*  Status:		Working				*)
	(*							*)
	(********************************************************)

(************************************************************************)
(*									*)
(* The algorithm used is Schrage's method, as described in		*)
(*	Stephen K. Park and Keith W. Miller, "Random Number Generators:	*)
(*	Good ones are hard to find", CACM 31(10), Oct 1988, 1192-1201.	*)
(* A basic property of this particular implementation is that all	*)
(* intermediate results fit into 32 bits (including sign).		*)
(*									*)
(************************************************************************)

FROM LowLevel IMPORT
    (* proc *)	HighWord, MakeLongword, Mul;

(************************************************************************)

CONST
    a = 16807;			(* 7^5 *)
    divisor = modulus DIV a;	(* = 127773 *)
    divisorL = divisor - 65536;	(* low order 16 bits of divisor *)

(************************************************************************)

PROCEDURE Divide (number: LONGCARD;
			VAR (*OUT*) quotient: CARDINAL;
			VAR (*OUT*) remainder: LONGCARD);

    (* Special-purpose division procedure - be careful about using	*)
    (* this in other applications, because it takes advantage of some	*)
    (* foreknowledge of the numbers which arise in this application.	*)
    (* Divides number by divisor=127773.				*)

    VAR step: CARDINAL;

    BEGIN
	remainder := number;
	quotient := 0;
	WHILE remainder >= divisor DO
	    step := ORD(HighWord(remainder)) DIV 2;
	    IF step = 0 THEN step := 1 END(*IF*);
	    INC (quotient, step);

	    (* Here, we take advantage of the fact that	*)
	    (*	divisor := 2^16 + divisorL		*)

	    remainder := remainder - MakeLongword(step, 0)
					- Mul(divisorL,step);
	END (*WHILE*);
    END Divide;

(************************************************************************)

PROCEDURE RandCardinal (): LONGCARD;

    (* Let a = 7^5 = 16807, q = modulus DIV a = 127773, and		*)
    (* r = modulus MOD a = 2836.					*)
    (* We perform the sequence of calculations:				*)
    (*		high := SHORT(seed DIV q);  low := seed MOD q;		*)
    (* 		test := a*low - r*LONG(high);				*)
    (*		IF test > 0 THEN seed := test				*)
    (*		ELSE seed := test + modulus;				*)
    (*		ENDIF;							*)
    (*		RETURN seed;						*)
    (* The division in the first step is complicated by the fact that	*)
    (* our processor does not provide a 32-bit division operation.	*)

    CONST r = modulus MOD a;	(* 2836 *)

    VAR high: CARDINAL;  low: LONGCARD;  test: LONGINT;

    BEGIN
	Divide (seed, high, low);
	test := VAL(LONGINT,a*low) - VAL(LONGINT,Mul(r, high));
	IF test > 0 THEN seed := test
	ELSE seed := test + modulus;
	END (*IF*);
	RETURN seed;
    END RandCardinal;

(****************************************************************)
(*			MODULE INITIALISATION			*)
(****************************************************************)

BEGIN
    seed := 1;
END RandCard.