[M3devel] Mersenne Twister in Modula-3, Error?

[Ken Durocher]

I implemented the Mersenne Twister in Modula-3 some time ago, however when I
test it out (the same test data that the authors of the Mersenne Twister
used), the output is different.  I tested this code on both 32 and 64 bit
machines, and they both give me the same results, but the results are
different from that of the original Mersenne Twister code.

Here's the code:

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(): REAL;
PROCEDURE genrandReal2(): REAL;
PROCEDURE genrandReal3(): REAL;

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] 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(): REAL =
  BEGIN
    RETURN FLOAT(genrandInt32(), REAL) * (1.0 / 4294967295.0);
  END genrandReal1;

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

PROCEDURE genrandReal3(): REAL =
  BEGIN
    RETURN (FLOAT(genrandInt32(), REAL) + 0.5) * (1.0 / 4294967296.0);
  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.

And here is the first few numbers it outputs:

3499211612 581869302 3890346734 3586334585 ...

Compare that to the first few numbers of output from the original Mersenne
Twister code, in C:

1067595299 955945823 477289528 4107218783 ...

Anyone know why they are different? Is the code wrong, or is it an
implementation difference, or what?
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://m3lists.elegosoft.com/pipermail/m3devel/attachments/20111018/900f2289/attachment-0001.html>