IMPLEMENTATION MODULE RAND; (****************************************************************************** Environnement : SOLARIS 2.3 ------------- Programmeurs : ------------ Pierre L'Ecuyer, Pierre Poulin, Gaetan Perron, Jean Belanger, Denis Alain Transportabilite : ---------------- On suppose que les entiers sont representes sur au moins 32 bits. *******************************************************************************) (* Ref. : L'Ecuyer, P. et Cote, S. (1987). "A Random Number Package with *) (* Splitting Facilities." Rapport no. DIUL-RR-8705, *) (* Departement d'informatique , Universite Laval. *) CONST Maxg = 32; (* Maxg peut etre initialisee a une autre valeur positive, *) (* pourvu que Maxg * (2 ** (v + w)) < 2.3058E18 et les constantes *) (* a1W, a2W, a1VW and a2VW soient recalculees . *) (* Notez que 2.3058E18 > (2**60). *) M1 = 2147483563; M2 = 2147483399; A1 = 40014; A2 = 40692; A1W = 1033780774; (* = A1**(2**w) MOD M1 *) A2W = 1494757890; (* = A2**(2**w) MOD M1 *) A1VW = 2082007225; (* = A1**(2**(v+w)) MOD M1 *) A2VW = 784306273; (* = A2**(2**(v+w)) MOD M2 *) UnSurM1 = 4.656613057E-10; (* 1 / M1 *) VAR Ig1, Ig2, Lg1, Lg2, Cg1, Cg2 : ARRAY Gen OF LONGINT; Antithetic : ARRAY Gen OF BOOLEAN; GenNumber : Gen; (******************************************************************************) 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 GetState ( g : Gen ; VAR S1, S2 : LONGINT ); BEGIN S1 := Cg1 [g]; S2 := Cg2 [g] END GetState; (******************************************************************************) PROCEDURE SetSeed ( g : Gen ; S1, S2 : LONGINT ); BEGIN Ig1 [g] := S1; Ig2 [g] := S2; InitGenerator (g, InitialSeed) END SetSeed ; (******************************************************************************) PROCEDURE AdvanceState ( g : Gen ; k : LONGCARD ); VAR B1, B2 : LONGINT; I : LONGCARD; BEGIN B1 := A1; B2 := A2; FOR I := 1 TO k DO B1 := MultMod (B1, B1, M1); B2 := MultMod (B2, B2, M2); END ; (* B1 = A1 ** k and B2 = A2 ** k. *) SetSeed (g, MultMod (B1, Cg1[g], M1), MultMod (B2, Cg2[g], M2)) END AdvanceState; (******************************************************************************) PROCEDURE InitGenerator ( g : Gen ; Where : SeedType ); BEGIN CASE Where OF InitialSeed : Lg1 [g] := Ig1 [g]; Lg2 [g] := Ig2 [g] | NewSeed : Lg1 [g] := MultMod (A1W, Lg1 [g], M1); Lg2 [g] := MultMod (A2W, Lg2 [g], M2) | LastSeed : END; Cg1 [g] := Lg1 [g]; Cg2 [g] := Lg2 [g] END InitGenerator; (******************************************************************************) PROCEDURE SetInitialSeed (S1, S2 : LONGINT ); VAR g : Gen; BEGIN Ig1 [1] := S1; Ig2 [1] := S2; InitGenerator (1, InitialSeed); FOR g := 2 TO Maxg DO Ig1 [g] := MultMod (A1VW, Ig1 [g-1], M1); Ig2 [g] := MultMod (A2VW, Ig2 [g-1], M2); InitGenerator (g, InitialSeed) END END SetInitialSeed; (******************************************************************************) PROCEDURE SetAntithetic ( g : Gen; B : BOOLEAN ); BEGIN Antithetic [g] := B END SetAntithetic; (******************************************************************************) PROCEDURE Random ( g : Gen ) : LONGCARD; VAR Z, k, s1, s2 : LONGINT; BEGIN s1 := Cg1 [g]; s2 := Cg2 [g]; k := s1 DIV 53668; s1 := 40014 * (s1 - k * 53668) - k * 12211; IF s1 < 0 THEN s1 := s1 + 2147483563 END; k := s2 DIV 52774; s2 := 40692 * (s2 - k * 52774) - k * 3791; IF s2 < 0 THEN s2 := s2 + 2147483399 END; Cg1 [g] := s1; Cg2 [g] := s2; Z := s1 - s2; IF Z < 1 THEN Z := Z + 2147483562 END; IF Antithetic [g] THEN Z := 2147483563 - Z END; RETURN LONGCARD (Z) END Random; (******************************************************************************) PROCEDURE Uniform01 ( g : Gen ) : LONGREAL; BEGIN RETURN LFLOAT (Random (g)) * UnSurM1; END Uniform01; (******************************************************************************) BEGIN (* Initialisation *) FOR GenNumber := 1 TO Maxg DO Antithetic [GenNumber] := FALSE END; SetInitialSeed (1234567890, 123456789) END RAND.