unit Shared.Math.Random.MT19937;
{$R-} {range checking off}
{$Warnings OFF}
// {$Q-} {overflow checking off}
{----------------------------------------------------------------------
Mersenne Twister: A 623-Dimensionally Equidistributed Uniform
Pseudo-Random Number Generator.
What is Mersenne Twister?
Mersenne Twister(MT) is a pseudorandom number generator developped by
Makoto Matsumoto and Takuji Nishimura (alphabetical order) during
1996-1997. MT has the following merits:
It is designed with consideration on the flaws of various existing
generators.
Far longer period and far higher order of equidistribution than any
other implemented generators. (It is proved that the period is 2^19937-1,
and 623-dimensional equidistribution property is assured.)
Fast generation. (Although it depends on the system, it is reported that
MT is sometimes faster than the standard ANSI-C library in a system
with pipeline and cache memory.)
Efficient use of the memory. (The implemented C-code mt19937.c
consumes only 624 words of working area.)
home page
http://www.math.keio.ac.jp/~matumoto/emt.html
original c source
http://www.math.keio.ac.jp/~nisimura/random/int/mt19937int.c
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, 1999 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.
REFERENCE
M. Matsumoto and T. Nishimura,
"Mersenne Twister: A 623-Dimensionally Equidistributed Uniform
Pseudo-Random Number Generator",
ACM Transactions on Modeling and Computer Simulation,
Vol. 8, No. 1, January 1998, pp 3--30.
Translated to OP and Delphi interface added by Roman Krejci (6.12.1999)
http://www.rksolution.cz/delphi/tips.htm
Revised 21.6.2000: Bug in the function RandInt_MT19937 fixed
----------------------------------------------------------------------}
interface
{Period parameter}
const
MT19937N = 624;
type
TMT19937StateArray =
array [0 .. MT19937N - 1]
of Int64;
// the array for the state vector
procedure sgenrand_MT19937(seed: Int64);
// Initialization by seed
procedure lsgenrand_MT19937(
const seed_array: TMT19937StateArray);
// Initialization by array of seeds
procedure randomize_MT19937;
// randomization
function randInt_MT19937(Range: Int64): Int64;
overload;
// Int64 Random with positive range
function randInt_MT19937(Min, Max: Int64): Int64;
overload;
// Int64 Random with positive range between min and max
function genrand_MT19937: Int64;
// random Int64 (full range);
// function randFloat_MT19937: Double; // float RANDOM on 0..1 interval
implementation
{Period parameters}
const
MT19937M = 397;
MT19937MATRIX_A = $9908B0DF;
// constant vector a
MT19937UPPER_MASK = $80000000;
// most significant w-r bits
MT19937LOWER_MASK = $7FFFFFFF;
// least significant r bits
// Tempering parameters
TEMPERING_MASK_B = $9D2C5680;
TEMPERING_MASK_C = $EFC60000;
type
Int64Rec =
packed record
case Integer
of
0:
(Lo, Hi: Cardinal);
1:
(Cardinals:
array [0 .. 1]
of Cardinal);
2:
(Words:
array [0 .. 3]
of Word);
3:
(Bytes:
array [0 .. 7]
of Byte);
end;
var
mt: TMT19937StateArray;
mti: Int64 = MT19937N + 1;
// mti=MT19937N+1 means mt[] is not initialized
// Initializing the array with a seed
procedure sgenrand_MT19937(seed: Int64);
var
i: Int64;
begin
for i := 0
to MT19937N - 1
do
begin
mt[i] := seed
and $FFFF0000;
seed := 69069 * seed + 1;
mt[i] := mt[i]
or ((seed
and $FFFF0000)
shr 16);
seed := 69069 * seed + 1;
end;
mti := MT19937N;
end;
{
Initialization by "sgenrand_MT19937()" is an example. Theoretically,
there are 2^19937-1 possible states as an intial state.
This function (lsgenrand_MT19937) allows to choose any of 2^19937-1 ones.
Essential bits in "seed_array[]" is following 19937 bits:
(seed_array[0]&MT19937UPPER_MASK), seed_array[1], ..., seed_array[MT19937-1].
(seed_array[0]&MT19937LOWER_MASK) is discarded.
Theoretically,
(seed_array[0]&MT19937UPPER_MASK), seed_array[1], ..., seed_array[MT19937N-1]
can take any values except all zeros.
}
procedure lsgenrand_MT19937(
const seed_array: TMT19937StateArray);
var
i: Int64;
begin
for i := 0
to MT19937N - 1
do
mt[i] := seed_array[i];
mti := MT19937N;
end;
function genrand_MT19937: Int64;
const
mag01:
array [0 .. 1]
of Int64 = (0, MT19937MATRIX_A);
var
y, kk: Int64;
begin
if mti >= MT19937N
then // generate MT19937N Int64 at one time
begin
if mti = (MT19937N + 1)
then // if sgenrand_MT19937() has not been called,
sgenrand_MT19937(4357);
// default initial seed is used
for kk := 0
to MT19937N - MT19937M - 1
do
begin
y := (mt[kk]
and MT19937UPPER_MASK)
or (mt[kk + 1]
and MT19937LOWER_MASK);
mt[kk] := mt[kk + MT19937M]
xor (y
shr 1)
xor mag01[y
and $00000001];
end;
for kk := MT19937N - MT19937M
to MT19937N - 2
do
begin
y := (mt[kk]
and MT19937UPPER_MASK)
or (mt[kk + 1]
and MT19937LOWER_MASK);
mt[kk] := mt[kk + (MT19937M - MT19937N)]
xor (y
shr 1)
xor mag01[y
and $00000001];
end;
y := (mt[MT19937N - 1]
and MT19937UPPER_MASK)
or (mt[0]
and MT19937LOWER_MASK);
mt[MT19937N - 1] := mt[MT19937M - 1]
xor (y
shr 1)
xor mag01[y
and $00000001];
mti := 0;
end;
y := mt[mti];
Inc(mti);
y := y
xor (y
shr 11);
y := y
xor (y
shl 7)
and TEMPERING_MASK_B;
y := y
xor (y
shl 15)
and TEMPERING_MASK_C;
y := y
xor (y
shr 18);
Result := y;
end;
// Delphi interface
procedure randomize_MT19937;
var
OldRandSeed: Integer;
begin
OldRandSeed := System.RandSeed;
// save system RandSeed value
System.Randomize;
// randseed value based on system time is generated
sgenrand_MT19937(System.RandSeed);
// initialize generator state array
System.RandSeed := OldRandSeed;
// restore system RandSeed
end;
// rewritten to delphi code due to 64bit compatibility
// Uwe Raabe - delphipraxis.com
function randInt_MT19937(Range: Int64): Int64;
begin
Result := Int64Rec(Int64(Range) * Cardinal(genrand_MT19937)).Hi;
end;
// bug fixed 21.6.2000.
// asm
// PUSH EAX // Parameter 1 = Range auf Stack sichern
// CALL genrand_MT19937 // Random Int64 aufrufen, Ergebnis steht in EAX (keine Parameter)
// POP EDX // Range vom Stack in EDX holen
// MUL EDX // EAX mit EDX multiplizieren, Ergebnis steht in EDX,EAX
// MOV EAX,EDX // Die h�heren 32-Bits in EDX als Result in EAX zur�ckgeben
// end;
function randInt_MT19937(Min, Max: Int64): Int64;
begin
Result := randInt_MT19937(Max - Min) + Min;
end;
// function randFloat_MT19937: Double;
// const
// Minus32: Double = -32.0;
// asm
// CALL genrand_MT19937
// PUSH 0
// PUSH EAX
// FLD Minus32
// FILD qword ptr [ESP]
// ADD ESP,8
// FSCALE
// FSTP ST(1)
// end;
initialization
// added 25.11.2019
randomize_MT19937;
end.