IMPLEMENTATION MODULE clcg4; FROM MYINOUT IMPORT WriteInt, WriteLn, WriteString; FROM NUM IMPORT WriteLR; (****************************************************************************** Environnement : MCS Modula-2 for Solaris. Programmeur : Pierre L'Ecuyer *******************************************************************************) (* Ref. : Modula-2 Implementation of a variant of AECL random number package, *) (* with new parameter values. *) VAR aw, avw : ARRAY [1..4] OF LONGINT; (* a_j^W and a_j^{VW}. *) a, m : ARRAY [1..4] OF LONGINT; Ig, Lg, Cg : ARRAY [1..4] OF ARRAY Gen OF LONGINT; (* Init, last, and current seeds. *) i, j : [1..4]; (* Index of a component. *) (******************************************************************************) PROCEDURE MultMod (a, s, m : LONGINT) : LONGINT; (* Suppose que 0 < a < m et 0 < s < m. Retourne (a*s) MOD m. *) CONST H = 32768; (* Doit etre egal a la valeur de 2**d.*) VAR a0, a1, q, qh, rh, k, p : LONGINT; BEGIN IF a < H THEN a0 := a; p := 0 ELSE a1 := a DIV H; a0 := a - H * a1; qh := m DIV H; rh := m - H * qh; IF a1 >= H THEN (* A2 = 1. *) a1 := a1 - H; k := s DIV qh; p := H * (s - k * qh) - k * rh; WHILE p < 0 DO p := p + m END ELSE (* p = (A2 * s * h) MOD m. *) p := 0 END; IF a1 <> 0 THEN q := m DIV a1; k := s DIV q; p := p - k * (m - a1 * q); IF p > 0 THEN p := p - m END; p := p + a1 * (s - k * q); WHILE p < 0 DO p := p + m END END; (* p = ((A2 * h + a1) * s) MOD m. *) k := p DIV qh; p := H * (p - k * qh) - k * rh; WHILE p < 0 DO p := p + m END END; (* p = ((A2 * h + a1) * h * s) MOD m *) IF a0 <> 0 THEN q := m DIV a0; k := s DIV q; p := p - k * (m - a0 * q); IF p > 0 THEN p := p - m END; p := p + a0 * (s - k * q); WHILE p < 0 DO p := p + m END END; RETURN p END MultMod; (******************************************************************************) PROCEDURE SetSeed ( g : Gen ; S1, S2, S3, S4 : LONGINT ); BEGIN Ig [1, g] := S1; Ig [2, g] := S2; Ig [3, g] := S3; Ig [4, g] := S4; InitGenerator (g, InitialSeed) END SetSeed; (******************************************************************************) PROCEDURE WriteState ( g : Gen); VAR i : INTEGER; BEGIN WriteLn; WriteString("State of generator g="); WriteInt(g, 2); WriteLn; FOR i:=1 TO 4 DO WriteString(" Cg["); WriteInt(i, 2); WriteString("] = "); WriteInt(Cg[i, g], 12); WriteLn; END; END WriteState; (******************************************************************************) PROCEDURE GetState ( g : Gen ; VAR S1, S2, S3, S4 : LONGINT ); BEGIN S1 := Cg [1, g]; S2 := Cg [2, g]; S3 := Cg [3, g]; S4 := Cg [4, g] END GetState; (******************************************************************************) PROCEDURE InitGenerator ( g : Gen ; Where : SeedType ); BEGIN FOR j := 1 TO 4 DO CASE Where OF InitialSeed : Lg [j, g] := Ig [j, g] | NewSeed : Lg [j, g] := MultMod (aw[j], Lg [j, g], m[j]) | LastSeed : END; Cg [j, g] := Lg [j, g]; END END InitGenerator; (******************************************************************************) PROCEDURE SetInitialSeed (S1, S2, S3, S4 : LONGINT); VAR g : Gen; BEGIN Ig [1, 0] := S1; Ig [2, 0] := S2; Ig [3, 0] := S3; Ig [4, 0] := S4; InitGenerator (0, InitialSeed); FOR g := 1 TO Maxgen DO FOR j := 1 TO 4 DO Ig [j, g] := MultMod (avw[j], Ig [j, g-1], m[j]) END; InitGenerator (g, InitialSeed) END END SetInitialSeed; (******************************************************************************) PROCEDURE Init (v, w : LONGINT); BEGIN FOR j := 1 TO 4 DO aw [j] := a [j]; FOR i := 1 TO w DO aw [j] := MultMod (aw[j], aw[j], m[j]) END; avw [j] := aw [j]; FOR i := 1 TO v DO avw [j] := MultMod (avw[j], avw[j], m[j]) END; END; SetInitialSeed (11111111, 22222222, 33333333, 44444444); (* Default seed for subsequence (generator) number 0. *) END Init; (******************************************************************************) PROCEDURE GenVal (g : Gen) : LONGREAL; VAR k, s1, s2, s3, s4 : LONGINT; U : LONGREAL; BEGIN U := 0.0; s1 := Cg [1, g]; k := s1 DIV 46693; s1 := 45991 * (s1 - k * 46693) - k * 25884; IF s1 < 0 THEN s1 := s1 + 2147483647 END; Cg [1, g] := s1; U := U + 4.65661287524579692E-10 * VAL (LONGREAL, s1); s2 := Cg [2, g]; k := s2 DIV 10339; s2 := 207707 * (s2 - k * 10339) - k * 870; IF s2 < 0 THEN s2 := s2 + 2147483543 END; Cg [2, g] := s2; U := U - 4.65661310075985993E-10 * VAL (LONGREAL, s2); IF U < 0.0 THEN U := U + 1.0 END; s3 := Cg [3, g]; k := s3 DIV 15499; s3 := 138556 * (s3 - k * 15499) - k * 3979; IF s3 < 0 THEN s3 := s3 + 2147483423 END; Cg [3, g] := s3; U := U + 4.65661336096842131E-10 * VAL (LONGREAL, s3); IF U > 1.0 THEN U := U - 1.0 END; s4 := Cg [4, g]; k := s4 DIV 43218; s4 := 49689 * (s4 - k * 43218) - k * 24121; IF s4 < 0 THEN s4 := s4 + 2147483323 END; Cg [4, g] := s4; U := U - 4.65661357780891134E-10 * VAL (LONGREAL, s4); IF U < 0.0 THEN U := U + 1.0 END; RETURN U; END GenVal; (******************************************************************************) PROCEDURE InitDefault; BEGIN Init (31, 41) END InitDefault; (******************************************************************************) BEGIN m[1] := 2147483647; a[1] := 45991; m[2] := 2147483543; a[2] := 207707; m[3] := 2147483423; a[3] := 138556; m[4] := 2147483323; a[4] := 49689; InitDefault; END clcg4.