{RAND4.PAS}
Type
MATRIZ1 = Array[0..2,0..2] of Double;
VETOR3 = Array[1..50,0..5] of Double;
VETOR6 = Array[0..5] of Double;
Const
M1: Double = 4294967087.0;
M2: Double = 4294944443.0;
NORM: Double = 1.0 / (4294967087.0 + 1.0);
A12: Double = 1403580.0;
A13N: Double = 810728.0;
A21: Double  = 527612.0;
A23N: Double = 1370589.0;
TWO17: Double = 131072.0;
TWO53: Double = 9007199254740992.0;
A1P127: MATRIZ1 =  (
       (    2427906178.0, 3580155704.0,  949770784.0),
       (    226153695.0,  1230515664.0,  3580155704.0),
       (    1988835001.0, 986791581.0,   1230515664.0)
       );

A2P127: MATRIZ1 = (
       (    1464411153.0,  277697599.0, 1610723613.0 ),
       (      32183930.0, 1464411153.0, 1022607788.0 ),
       (    2824425944.0,   32183930.0, 2093834863.0 )
       );

Var
CG, BG, IG: VETOR3;
Semente:  VETOR6;
SERIE: Integer;
{-------------------------------------------------------------------}
Function MULTMODM (A, S, C, M: Double): Double;
Var
  V: Double;
  A1: Longint;
Begin
    V := A * S + C;
    If ((V >= TWO53) or (V <= -TWO53)) then 
    Begin
        A1 := trunc(A / TWO17);
        A := A - (A1 * TWO17);
        V  := A1 * S;
        A1 := trunc(V / M);
        V := V - (A1 * M);
        V := V * TWO17 + A * S + C;
    End;
    A1 := trunc(V / M);
    V := V - (A1 * M);
    If (V < 0.0) then  MULTMODM := V + M   else MULTMODM := V;
End;
{-----------------------------------------------------------------}
Procedure MATVECMODM (Var A: MATRIZ1; Var S: array of Double; Var V: Array of Double; M: Double);
Var
   I: Integer;
   X: Array[0..2] of Double;
   ZERO: Double;
Begin
    ZERO := 0.0;
    For I := 0 to 2 do
       Begin
        X[I] := MULTMODM (A[I][0], S[0], ZERO, M);
        X[I] := MULTMODM (A[I][1], S[1], X[I], M);
        X[I] := MULTMODM (A[I][2], S[2], X[I], M);
       End;
     For I := 0 to 2 do
       Begin
        V[I] := X[I];
       End;
End;
{---------------------------------------------------------------------}
Function RAND4 : Double;
Var
   K: Longint;
   P1, P2, U: Double;
Begin
    { Componente 1 }
    P1 := A12 * CG[SERIE][1] - A13N * CG[SERIE][0];
    K := trunc(P1/M1);
    P1 := P1 - (K * M1);
    if (P1 < 0.0) then P1 := P1 + M1;
    CG[SERIE][0] := CG[SERIE][1];
    CG[SERIE][1] := CG[SERIE][2];
    CG[SERIE][2] := P1;
    { Componente 2 }
    P2 := A21 * CG[SERIE][5] - A23N * CG[SERIE][3];
    K := trunc(P2/M2);
    P2 := P2 - (K * M2);
    if (P2 < 0.0) then P2 := P2 + M2;
    CG[SERIE][3] := CG[SERIE][4];
    CG[SERIE][4] := CG[SERIE][5];
    CG[SERIE][5] := P2;
    { Combinacao}
    if (P1 > P2) then U := (P1 - P2) * NORM
        else U := (P1 - P2 + M1) * NORM;
    RAND4 := U;
End;
{-------------------------------------------------------------------------}
Procedure PASSASEMENTE;
Var
   I,J: Integer;
   TEMP: array[0..2] of Double;
Begin
  For I := 0 to 5 do
      Begin
        BG[1][I] :=  Semente[I];
        CG[1][I] :=  Semente[I];
        IG[1][I] :=  Semente[I];
      End;
   For J := 2 to 50 do
     Begin
       For I := 0 to 2 do
          Begin
           TEMP[I] := Semente[I];
          End;
     MATVECMODM (A1P127, TEMP, TEMP, M1);
     For I := 0 to 2 do
       Begin
        Semente[I] := TEMP[I];
       End;
     For I := 0 to 2 do
       Begin
        TEMP[I] := Semente[I + 3];
       End;
     MATVECMODM (A2P127, TEMP, TEMP, M2);
     For I := 0 to 2 do
       Begin
        Semente[I + 3] := TEMP[I];
       End;
     For I := 0 to 5 do
       Begin
         BG[J][I] :=  Semente[I];
         CG[J][I] :=  Semente[I];
         IG[J][I] :=  Semente[I];
       End;
   End;
 End;
{-----------------------------------------------------}