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(a)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.