Hi all.

I've seen some discussion about random number generators here recently. I came across a quite good one (this being something I use extensively in my work), and thought I'd share it. It is written up in ACM Transactions on Modeling and Computer Simulation 8(1) 3-30, 1990. Supposedly it has a period of -1 + 2^19937. The original was in C; this translation returns the same values as the original. It is called the Mersenne Twister.

enjoy. Toby

(* A Pascal program for MT19937: Real number version([0,1)-interval). Written in C 4/6/1998 by Takuji Nishimura. Translated to GNU Pascal 11/2/1998 by Toby Ewing.

genrand() generates one pseudorandom real number (double) which is uniformly distributed in [0,1)-interval, for each call. sgenrand(seed) sets initial values to the working area of 624 words. Before genrand(), sgenrand(seed) must be called once. (The seed is any 32-bit integer except for 0). Coded by Takuji Nishimura, considering the suggestions by Topher Cooper and Marc Rieffel in July-Aug. 1997.

This library is free software; you can redistribute it and/or modify it under the terms of the GNU Library General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.

This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public License for more details. You should have received a copy of the GNU Library General Public License along with this library; if not, write to the Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura. When you use this, send an email to: matumoto@math.keio.ac.jp with an appropriate reference to your work. *)

Program MersenneTwister;

Const N = 624; {Period parameters} M = 397;

type LI = __longlong__ integer; UL = __unsigned__ LI; {define 64-bit unsigned integer} mtary = array[0..N] of UL; magary = array[0..1] of UL;

var j : integer; mti : integer value N+1; {mti = 1 means mt[N] is not initialized} mt : mtary; {the array for the state vector}

Procedure sgenrand(seed : UL); { Set initial seeds to mt[N] using the generator Line 25 of Table 1 in [KNUTH 1981, The Art of Computer Programming Vol. 2 (2nd Ed.), pp 102]. } begin mt[0] := seed AND 16#ffffffff; for mti := 1 to N do mt[mti] := (69069 * mt[mti-1]) AND 16#ffffffff; end;

Function genrand : real;

var y : UL; kk : __short__ integer; Matrix_A : UL value 16#9908b0df; {constant vector a} Up_Mask : UL value 16#80000000; {most significant w-r bits} Low_Mask : UL value 16#7fffffff; {least significant r bits} T_Mask_B : UL value 16#9d2c5680; {Tempering parameters} T_Mask_C : UL value 16#efc60000; mag01 : magary value (0, Matrix_A);

begin if (mti >= N) then begin {generate N words at one time} if (mti = N+1) then sgenrand(4357);{If sgenrand hasn't been called, use default seed.}

for kk := 0 to N-M do begin y := (mt[kk] AND Up_Mask) OR (mt[kk+1] AND Low_Mask); mt[kk] := mt[kk+M] XOR (y shr 1) XOR mag01[y AND 16#00000001]; end;

for kk := N-M to N-1 do begin y := (mt[kk] AND Up_Mask) OR (mt[kk+1] AND Low_Mask); mt[kk] := mt[kk+(M-N)] XOR (y shr 1) XOR mag01[y AND 16#00000001]; end; y := (mt[N-1] AND Up_Mask) OR (mt[0] AND Low_Mask); mt[N-1] := mt[M-1] XOR (y shr 1) XOR mag01[y AND 16#00000001];

mti := 0; end;

y := mt[mti]; inc(mti); y := y XOR (y shr 11); y := y XOR ((y shl 7) AND T_Mask_B); y := y XOR ((y shl 15) AND T_Mask_C); y := y XOR (y shr 18);

genrand := real(y) * 2.3283064365386963e-10; end;

begin {Output first 1000 generated numbers} sgenrand(4357); {any integer > 0 works as a seed} for j := 1 to 1000 do begin write(genrand:10:8); if ((j mod 8) = 0) then writeln else write(' '); end; end.