low@melair.UUCP (Rick Low) (12/02/90)
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