INTERFACE MT19937ar;

IMPORT Word;

PROCEDURE initGenrand(s: Word.T);
PROCEDURE initByArray(key: ARRAY OF Word.T; length: INTEGER);

PROCEDURE genrandInt32(): Word.T;
PROCEDURE genrandInt31(): INTEGER;
PROCEDURE genrandReal1(): LONGREAL;
PROCEDURE genrandReal2(): LONGREAL;
PROCEDURE genrandReal3(): LONGREAL;

END MT19937ar.

(**********************************)

MODULE MT19937ar;

IMPORT Word;
FROM Word IMPORT And, Xor, Or, RightShift, LeftShift;

CONST
  n = 624;
  m = 397;
  matrix_a = 16_9908b0df;
  upper_mask = 16_80000000;
  lower_mask = 16_7fffffff;

VAR
  mt: ARRAY [0..n-1] OF Word.T;
  mti := n + 1;

PROCEDURE initGenrand(s: Word.T) =
  BEGIN
    mt[0] := And(s, 16_ffffffff);
    FOR mti := 1 TO n - 1 DO
      mt[mti] := 1812433253 * (Xor(mt[mti-1], RightShift(mt[mti-1], 30))) + mti;
      mt[mti] := And(mt[mti], 16_ffffffff);
    END;
  END initGenrand;

PROCEDURE initByArray(key: ARRAY OF Word.T; length: INTEGER) =
  VAR i, j, k: INTEGER;
  BEGIN
    i := 1; j := 0;
    initGenrand(19650218);
    IF n > length THEN k := n ELSE k := length; END;

    WHILE k # 0 DO
      mt[i] := Xor(mt[i], Xor(mt[i-1], RightShift(mt[i-1], 30)) * 1664525) +
                   key[j] + j;
      mt[i] := And(mt[i], 16_ffffffff);
      INC(i); INC(j);

      IF i >= n THEN
        mt[0] := mt[n-1];
        i := 1;
      END;

      IF j >= length THEN
        j := 0;
      END;
      DEC(k);
    END;

    k := n - 1;
    WHILE k # 0 DO
      mt[i] := Xor(mt[i], Xor(mt[i-1], RightShift(mt[i-1], 30)) * 1566083941) - i;
      mt[i] := And(mt[i], 16_ffffffff);
      INC(i);

      IF i >= n THEN
        mt[0] := mt[n-1];
        i := 1;
      END;
      DEC(k);
    END;

    mt[0] := 16_80000000;
  END initByArray;

PROCEDURE genrandInt32(): Word.T =
  VAR y: Word.T;
      k: INTEGER := 0;
      mag01 := ARRAY [0..1] OF Word.T {0, matrix_a};

  BEGIN
    IF mti >= n THEN
      IF mti = (n + 1) THEN
        initGenrand(5489);
      END;

      WHILE k < (n - m) DO
        y := Or(And(mt[k], upper_mask), And(mt[k+1], lower_mask));
        mt[k] := Xor(Xor(mt[k+m], RightShift(y, 1)), mag01[And(y, 1)]);
        INC(k);
      END;

      WHILE k < (n - 1) DO
        y := Or(And(mt[k], upper_mask), And(mt[k+1], lower_mask));
        mt[k] := Xor(Xor(mt[k+(m-n)], RightShift(y, 1)), mag01[And(y, 1)]);
        INC(k);
      END;

      y := Or(And(mt[n-1], upper_mask), And(mt[0], lower_mask));
      mt[n-1] := Xor(Xor(mt[m-1], RightShift(y, 1)), mag01[And(y, 1)]);
      mti := 0;
    END;

    y := mt[mti];
    INC(mti);

    y := Xor(y, RightShift(y, 11));
    y := Xor(y, And(LeftShift(y, 7), 16_9d2c5680));
    y := Xor(y, And(LeftShift(y, 15), 16_efc60000));
    y := Xor(y, RightShift(y, 18));

    RETURN y;
  END genrandInt32;

PROCEDURE genrandInt31(): INTEGER =
  BEGIN
    RETURN RightShift(genrandInt32(), 1);
  END genrandInt31;

PROCEDURE genrandReal1(): LONGREAL =
  BEGIN
    RETURN FLOAT(genrandInt32(), LONGREAL) * (1.0D0 / 4294967295.0D0);
  END genrandReal1;

PROCEDURE genrandReal2(): LONGREAL =
  BEGIN
    RETURN FLOAT(genrandInt32(), LONGREAL) * (1.0D0 / 4294967296.0D0);
  END genrandReal2;

PROCEDURE genrandReal3(): LONGREAL =
  BEGIN
    RETURN (FLOAT(genrandInt32(), LONGREAL) + 0.5D0) * (1.0D0 / 4294967296.0D0);
  END genrandReal3;

BEGIN
END MT19937ar.

(**********************************)

MODULE Main;

IMPORT IO, Fmt, MT19937ar;

VAR init := ARRAY [0..3] OF INTEGER {16_123, 16_234, 16_345, 16_456};
    len := 4;

BEGIN
  MT19937ar.initByArray(init, len);
  IO.Put("Random integers: ");
  FOR i := 1 TO 999 DO
    IO.Put(Fmt.Unsigned(MT19937ar.genrandInt32(), base := 10));
    IO.Put(" ");
    IF (i MOD 5) = 4 THEN IO.Put("\n"); END;
  END;
  IO.Put("\n");
END Main.