unit MT19937; {$R-} {range checking 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 longint; // the array for the state vector procedure sgenrand_MT19937(seed: longint); // Initialization by seed procedure lsgenrand_MT19937(const seed_array: tMT19937StateArray); // Initialization by array of seeds procedure randomize_MT19937; // randomization function randInt_MT19937(Range: longint):longint; // integer RANDOM with positive range function genrand_MT19937: longint; // random longint (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; VAR mt : tMT19937StateArray; mti: integer=MT19937N+1; // mti=MT19937N+1 means mt[] is not initialized { Initializing the array with a seed } procedure sgenrand_MT19937(seed: longint); var i: integer; 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: integer; begin for i := 0 to MT19937N-1 do mt[i] := seed_array[i]; mti := MT19937N; end; function genrand_MT19937: longint; const mag01 : array [0..1] of longint =(0, MT19937MATRIX_A); var y: longint; kk: integer; begin if mti >= MT19937N { generate MT19937N longints at one time } then 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: longint; 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; // bug fixed 21.6.2000. Function RandInt_MT19937(Range: longint):longint; // EAX <- Range // Result -> EAX asm PUSH EAX CALL genrand_MT19937 POP EDX MUL EDX MOV EAX,EDX 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; end.