From 555d3301ff2b9b0b43ed2fb0e0e2df968e0e228c Mon Sep 17 00:00:00 2001 From: Norayr Chilingarian Date: Wed, 30 Oct 2013 21:39:40 +0400 Subject: [PATCH] added MersenneTwister.Mod which can be compiled with XDS system. Now will port it to voc --- src/lib/misc/MersenneTwister.Mod | 214 +++++++++++++++++++++++++++++++ 1 file changed, 214 insertions(+) create mode 100644 src/lib/misc/MersenneTwister.Mod diff --git a/src/lib/misc/MersenneTwister.Mod b/src/lib/misc/MersenneTwister.Mod new file mode 100644 index 00000000..53d02a94 --- /dev/null +++ b/src/lib/misc/MersenneTwister.Mod @@ -0,0 +1,214 @@ +<* O2EXTENSIONS + *> +<* IOVERFLOW - *> +MODULE Random; + + IMPORT + Sys:=SYSTEM,Win:=Windows,MathL; + + CONST + (* Period parameter *) + MT19937N*=624; + + (* Period parameters *) + MT19937M=397; + MT19937MATRIX_A =Sys.VAL(SET,9908b0dfH); -- constant vector a + MT19937UPPER_MASK=Sys.VAL(SET,80000000H); -- most significant w-r bits + MT19937LOWER_MASK=Sys.VAL(SET,7fffffffH); -- least significant r bits + + (* Tempering parameters *) + TEMPERING_MASK_B=Sys.VAL(SET,9d2c5680H); + TEMPERING_MASK_C=Sys.VAL(SET,0efc60000H); + + Seed0=4357; + TYPE + tMT19937StateArray=ARRAY MT19937N OF SET; -- the array for the state vector + + VAR + Seed-:LONGINT; + + mt : tMT19937StateArray; + mti: LONGINT; -- mti=MT19937N+1 means mt[] is not initialized + + GaussRandomBuf:LONGREAL; + GaussRandomBufFilled:BOOLEAN; + + (* Initializing the array with a seed *) + PROCEDURE SetSeed*(seed:LONGINT);(* sgenrand_MT19937 *) + CONST + HighBits=Sys.VAL(SET,0ffff0000H); + VAR + i:LONGINT; + BEGIN + Seed:=seed; + FOR i:=0 TO MT19937N-1 DO + mt[i]:=Sys.VAL(SET,seed) * HighBits; + seed:=69069*seed+1; + mt[i]:=mt[i] + (Sys.SHIFT(Sys.VAL(SET,seed) * HighBits,-16)); + seed:=69069*seed+1; + END; + mti := MT19937N; + END SetSeed; + + (* Initialization by array of seeds *) + PROCEDURE SetSeeds*(seed_array-:tMT19937StateArray); (* lsgenrand_MT19937 *) + VAR + i:LONGINT; + BEGIN + FOR i:=0 TO MT19937N-1 DO + mt[i]:=seed_array[i]; + END; + mti:=MT19937N; + END SetSeeds; + + (* random longint (full range) *) + PROCEDURE Int*():LONGINT; (* genrand_MT19937 *) + TYPE + ar=ARRAY 2 OF SET; + VAR + mag01:ARRAY 2 OF SET; + y:SET; + kk:LONGINT; + BEGIN + mag01[0]:={}; + mag01[1]:=MT19937MATRIX_A; + + IF mti>=MT19937N THEN (* generate MT19937N longints at one time *) + IF mti=(MT19937N+1) THEN -- if sgenrand_MT19937() has not been called, + SetSeed(Seed0); -- default initial seed is used + END; + FOR kk:=0 TO MT19937N-MT19937M-1 DO + y:=(mt[kk] * MT19937UPPER_MASK) + (mt[kk+1] * MT19937LOWER_MASK); + mt[kk]:=mt[kk+MT19937M]/Sys.SHIFT(y,-1)/mag01[Sys.VAL(LONGINT,y * {0})]; + END; + FOR kk:=MT19937N-MT19937M TO MT19937N-2 DO + y:=(mt[kk] * MT19937UPPER_MASK) + (mt[kk+1] * MT19937LOWER_MASK); + mt[kk]:=mt[kk+(MT19937M-MT19937N)]/Sys.SHIFT(y,-1)/mag01[Sys.VAL(LONG +INT,y * {0})]; + END; + y:=(mt[MT19937N-1] * MT19937UPPER_MASK) + (mt[0] * MT19937LOWER_MASK); + mt[MT19937N-1]:=mt[MT19937M-1]/Sys.SHIFT(y,-1)/mag01[Sys.VAL(LONGINT,y +* {0})]; + mti:=0; + END; + y:=mt[mti]; INC(mti); + y:=y/Sys.SHIFT(y,-11); + y:=y/(Sys.SHIFT(y,7) * TEMPERING_MASK_B); + y:=y/(Sys.SHIFT(y,15) * TEMPERING_MASK_C); + y:=y/Sys.SHIFT(y,-18); + RETURN Sys.VAL(LONGINT,y); + END Int; + + (*randomization*) + PROCEDURE Randomize*(); (* Randomize_MT19937 *) + VAR + ST:Win.SYSTEMTIME; + BEGIN + Win.GetSystemTime(ST); + SetSeed(((Sys.VAL(LONGINT,ST.wHour)*60+ST.wMinute)*60+ST.wSecond)*1000+S +T.wMilliseconds); + END Randomize; + + (*integer RANDOM with positive range*) + -- bug fixed 21.6.2000. + PROCEDURE IntRange*(Range:LONGINT):LONGINT; (* RandInt_MT19937 *) + TYPE + VAR + X:Sys.CARD64; + BEGIN + X:=Range; + X:=X * Sys.VAL(Sys.CARD64,Int()); + Sys.MOVE(Sys.ADR(X)+4,Sys.ADR(Range),4); + RETURN Range; + END IntRange; + + (*float RANDOM on 0..1 interval*) + PROCEDURE Real*():LONGREAL; (* RandFloat_MT19937 *) + BEGIN + RETURN Sys.VAL(Sys.CARD32,Int())/(1.0*MAX(Sys.CARD32)+1) + END Real; + + PROCEDURE Gauss*(mean,std:LONGREAL):LONGREAL; + VAR + r1,r2,s,result:LONGREAL; + BEGIN + IF GaussRandomBufFilled THEN + result:=GaussRandomBuf*std+mean + ELSE + REPEAT + r1:=2*Real()-1; + r2:=2*Real()-1; + s:=r1*r1+r2*r2; + UNTIL s<1; + s:=MathL.sqrt((-2*MathL.ln(s))/s); + result:=r1*s*std+mean; + GaussRandomBuf:=r2*s; + END; + GaussRandomBufFilled:=NOT GaussRandomBufFilled; + RETURN result + END Gauss; + +BEGIN + Seed:=Seed0; + mti:=MT19937N+1; + GaussRandomBufFilled:=FALSE; +END Random. + +(*---------------------------------------------------------------------- + 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 +----------------------------------------------------------------------*)