Path: utzoo!utgpu!cunews!mitel!melair!low
From: low@melair.UUCP (Rick Low)
Newsgroups: comp.lang.modula2
Subject: Re: RND Function Needed -- Here's One for FST Modula-2
Message-ID: <407@melair.UUCP>
Date: 1 Dec 90 21:18:42 GMT
References: <1484.2750D9EC@puddle.fidonet.org> <EIVERSON.90Nov27010159@aigyptos.nmsu.edu>
Organization: MEL Defence Systems Ltd., Ottawa, Canada
Lines: 128


Here's a random number generator that I ported to FST Modula-2 a while ago
from the book 'Numerical Recipes in C'.

-----------------------RANDOM.DEF cut here------------------------
(* Random -- generate uniformly distributed random numbers *)
DEFINITION MODULE Random;

(* Ran1 -- return a uniform random deviate in [0.0, 1.0) *)
PROCEDURE Ran1(VAR Idum : INTEGER) : REAL;

END Random.
-----------------------RANDOM.MOD cut here-----------------------
(* Random -- generate uniform random numbers in [0.0, 1.0) *)
IMPLEMENTATION MODULE Random;

(*
 *--------------------------------------------------------------
 * DESCRIPTION:
 * This random number generator is a translation into Modula-2
 * of the routine ran1() in
 *
 * Press, W.H., B.T. Flannery, S.A. Teukolsky, and W.T. Vetterling,
 *	"Numerical Recipes in C -- The Art of Scientific Programming",
 *	New York, Cambridge University Press, 1988.
 *---------------------------------------------------------------
 *)

FROM InOut	IMPORT WriteString;
FROM LMathLib0	IMPORT real, entier;
FROM System	IMPORT ErrorMessage, ErrorLn, Terminate;

CONST
	M1	= 259200;
	IA1	= 7141L;
	IC1	= 54773L;
	RM1	= (1.0/FLOAT(M1));
	M2	= 134456;
	IA2	= 8121L;
	IC2	= 28411L;
	RM2	= (1.0/FLOAT(M2));
	M3	= 243000;
	IA3	= 4561L;
	IC3	= 51349L;

VAR
	iff	: INTEGER;
	ix1, ix2, ix3
		: LONGINT;
	r	: ARRAY [0..97] OF REAL;


(* Ran1 -- return a uniform random deviate in [0.0, 1.0) *)
PROCEDURE Ran1(VAR Idum : INTEGER) : REAL;

VAR
	j	: INTEGER;
	temp	: REAL;

BEGIN
	(* Initialize on first call, even if Idum is not negative. *)
	IF (Idum < 0) OR (iff = 0) THEN
		iff := 1;
		(* Seed first routine... *)
		ix1 := (IC1 - LONG(Idum)) MOD M1;
		ix1 := (IA1 * ix1 + IC1) MOD M1;
		(* and use it to seed the second ...*)
		ix2 := ix1 MOD M2;
		ix1 := (IA1 * ix1 + IC1) MOD M1;
		(* and third routines. *)
		ix3 := ix1 MOD M3;

		(*
		 * Fill the table with sequential uniform deviates
		 * generated by the first two routines.
		 *)
		FOR j := 1 TO 97 DO
			ix1 := (IA1 * ix1 + IC1) MOD M1;
			ix2 := (IA2 * ix2 + IC2) MOD M2;
			(* Combine low and high order pieces. *)
			r[j] := SHORT((real(ix1) + real(ix2) * RM2) * RM1);
		END;
		Idum := 1;
	END; (* initialization *)

	(*
	 * Begin normal (non-initialization) runs here.
	 * Generate the next number for each sequence.
	 *)
	ix1 := (IA1 * ix1 + IC1) MOD M1;
	ix2 := (IA2 * ix2 + IC2) MOD M2;
	ix3 := (IA3 * ix3 + IC3) MOD M3;
	(* Use the third sequence to get an integer between 1 and 97. *)
	j := SHORT(1L + (97L * ix3) DIV M3);

	(*
	 * Handle conditions that should never happen -- limit j to [1,97].
	 * (Press, et. al. abort with an error message here.)
	 *)
	IF ( j < 1) OR (j > 97) THEN 
		ErrorMessage('?? Ran1():  panic');
		ErrorLn;
		Terminate(2);
	END;

	(* Return the table entry and refill it. *)
	temp := r[j];
	r[j] := SHORT((real(ix1) + real(ix2) * RM2) * RM1);
	RETURN(temp);

END Ran1;


(* run-time initialization *)
BEGIN (* Random *)
	iff := 0;
END Random.
--------------------------cut here---------------------

No guarantees, but it works for me!

Rick Low
+1 613 957 9500
mitel!melair!low@uunet.uu.net
-- 
Rick Low
+1 613 957 9500
mitel!melair!low@uunet.uu.net


Brought to you by Super Global Mega Corp .com