RandomLib.mod.m2pp
662 lines (586 with data), 21.3 kB
IMPLEMENTATION MODULE RandomLib;
(*========================================================================*)
(* HINWEIS : Bitte nur die Datei RandomLib.mod.m2pp editieren *)
(*========================================================================*)
(* Es sind 2 Versionen enthalten die mit *)
(* *)
(* m2pp -D __{Parameter}__ < RandomLib.mod.m2pp > RandomLib.mod *)
(* *)
(* mit Parameter = {XDS|GM2} erzeugt werden koennen. *)
(* *)
(* XDS : Der Export von Objekten aus lokalen Modulen wird mit *)
(* "IMPORT" angezeigt - das scheint nicht ganz "regelkonform" *)
(* GM2 : Normales "EXPORT" wird genutzt (ELSE-Fall). *)
(* *)
(* ansonsten gibt es keine Aenderungen am Quellcode *)
(*------------------------------------------------------------------------*)
(* Dieses Modul stellt verschiedene uniforme Pseudozufallszahlen- *)
(* generatoren zur Verfuegung. *)
(* This module providing several uniform pseudorandom number generators. *)
(*------------------------------------------------------------------------*)
(* Letzte Aenderung: *)
(* *)
(* 03.06.93, MRi: Erstellen der ersten Version von RanMar / RanMarInit *)
(* 13.10.15, MRi: Ersetzen von LFLOAT(LTRUNC( durch LowLong.intpart( *)
(* 21.01.16, MRi: Einfuehren von InitZufallSysT, *)
(* Verbesserung der Doku im Definitionsmodul *)
(* 06.09.15, MRi: Uebersetzen von ACM 133 von Algol60 nach Modula-2 *)
(* 02.12.15, MRi: Einfuegen von InitRandom fuer ACM 133 *)
(* 22.01.16, MRi: Einfuegen von ACM133 in Modul Random *)
(* 22.01.16, MRi: Erstellen der ersten Version des CERN Library RanLux *)
(* Random number generators - uebersetzt von Fortran 77 *)
(* nach Modula-2. *)
(* 23.01.16, MRi: Kleinere Anpassungen an RanLux, EA hinzugefuegt *)
(* RanLux,RLuxGo gegen F77 Version mit 8*10**6 generierten *)
(* Zufallszahlen erfolgreich getestet. *)
(* 23.01.16, MRi: Integrieren von RanLux in Random, umbenennen des Moduls *)
(* in RandomLib, testen der Restart-Optionen von RanLux *)
(* 04.04.16, MRi: Initialisierungsroutinen von Random und Zufall ueber- *)
(* arbeitet um sicherzustellen dass nicht mit 0 *)
(* initialisiert wird. *)
(* 12.04.16, MRi: Initialisierungsroutinen von Random und Zufall ueber- *)
(* arbeitet damit nach kuerzerer Zeit eine neue Zufalls- *)
(* zahlenreihe abgerufen werden kann. *)
(*------------------------------------------------------------------------*)
(* Offene Punkte *)
(* *)
(* - Wieso funktioniert ein Aufruf von RLuxGo mit K1,K2 # 0 nicht bei *)
(* luxlev 3 & 4 - dann schlaegt auch der Restart mit RLuxAt/RLuxGo *)
(* fehl (merkwuergerweise sind die ersten Werte noch richtig), nicht *)
(* aber der mit RLuxUt/RLuxGo&RLuxIn. Wenn K1+K2=0 ist alles gut. *)
(* Dir Fortran-Routine verhaelt sich aber auch so ... *)
(* - Herausnahme der debug-optionen in RanLux (ersetzen durch *)
(* Praeprozessor-direktiven). *)
(*------------------------------------------------------------------------*)
(* For procedure Random see the ACM licence terms at *)
(* https://www.acm.org/publications/policies/software-copyright-notice, *)
(* all other routines are under the LGPL. *)
(* *)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: RandomLib.mod.m2pp,v 1.7 2018/01/16 09:19:51 mriedl Exp mriedl $ *)
FROM SysClock IMPORT maxSecondParts,DateTime,CanGetClock,GetClock;
IMPORT LowLong;
IMPORT LongMath;
IMPORT Errors;
IMPORT TIO;
MODULE ACM133;
IMPORT maxSecondParts,DateTime,CanGetClock,GetClock;
<* IF (__XDS__) THEN *>
IMPORT Random,InitRandom; (* EXPORT *)
<* ELSE *>
EXPORT Random,InitRandom;
<* END *>
VAR X : ARRAY [1..2] OF INTEGER;
PROCEDURE Random(a,b : LONGREAL) : LONGREAL;
(* ACM Algorithm 133 *)
CONST mult = 5;
mscale = 1048576;
topscale = 32768;
mrecip = 3.435973E-10; (* 2.0**(-35) *)
toprecip = 3.051757E-05; (* 2.0**(-15) *)
BEGIN
X[1] := mult * X[1];
X[2] := mult * X[2];
(* Add overflow from X[2] and remove it from X[2] *)
IF (X[2] >= mscale) THEN
X[1] := X[1] + X[2] DIV mscale;
X[2] := X[2] MOD mscale;
END;
(* Possibly need to tidy up X[1] *)
IF (X[1] >= topscale) THEN
X[1] := X[1] MOD topscale;
END;
(* Generate a real value in the required range *)
RETURN (VAL(LONGREAL,X[1])*toprecip + VAL(LONGREAL,X[2])*mrecip) *
(b-a) + a;
END Random;
PROCEDURE InitRandom();
VAR dt : DateTime;
T,frac : CARDINAL;
x : REAL;
BEGIN
IF CanGetClock() THEN
GetClock(dt);
frac := dt.fractions + 1;
IF (frac = 1) THEN frac:=maxSecondParts; END;
(*
T:=(((dt.hour + 1)*60 + (dt.minute + 1))*60 +
(dt.second + 1))*maxSecondParts + frac;
*)
T:=((dt.second + 1))*maxSecondParts + frac;
x:=VAL(REAL,T); (* / VAL(REAL,24*360*maxSecondParts); *)
WHILE (x >= 10.0) DO x:=x / 3.14152927; END;
WHILE (x < 10.0) DO x:=x * 10.00; END;
X[1]:=TRUNC(x);
X[2]:=TRUNC(4096.0*(x-VAL(REAL,X[1])));
ELSE
X[1]:=VAL(INTEGER,Random(1.0,111.0)) DIV VAL(INTEGER,Random(1.0,111.0));
X[2]:=VAL(INTEGER,Random(5.0,555.0)) MOD VAL(INTEGER,Random(5.0,555.0));
END;
END InitRandom;
BEGIN
X[1]:=13;
X[2]:=19;
END ACM133;
MODULE RandomMar;
IMPORT Errors;
<* IF (__XDS__) THEN *>
IMPORT RanMarInit,RanMar; (* EXPORT *)
<* ELSE *>
EXPORT RanMarInit,RanMar;
<* END *>
VAR U : ARRAY [1..97] OF LONGREAL;
c,cd,cm : LONGREAL;
i97,j97 : CARDINAL;
PROCEDURE RanMarInit(ij,kl : CARDINAL);
VAR i,j,k,l,m : CARDINAL;
ii,jj : CARDINAL;
s,t : LONGREAL;
BEGIN
IF (ij > 31328) OR (kl > 30081) THEN
Errors.ErrOut("(ij IN {0..31328}) oder ");
Errors.FatalError("(kl IN {0..30081}) nicht erfuellt (Random.RanMarInit)");
END;
i := ((ij DIV 177) MOD 177) + 2;
j := ( ij MOD 177) + 2;
k := ((kl DIV 169) MOD 178) + 1;
l := ( kl MOD 169);
FOR ii:= 1 TO 97 DO
s := 0.0;
t := 0.5;
FOR jj:=1 TO 24 DO
m := ((i*j MOD 179)*k) MOD 179;
i := j;
j := k;
k := m;
l := (53*l+1) MOD 169;
IF (((l*m) MOD 64) >= 32) THEN s:=s + t; END;
t := 0.5 * t;
END;
U[ii] := s;
END;
c := 362436.0 / 16777216.0;
cd := 7654321.0 / 16777216.0;
cm := 16777213.0 / 16777216.0;
i97 := 97;
j97 := 33;
END RanMarInit;
PROCEDURE RanMar() : LONGREAL;
VAR uni : LONGREAL;
BEGIN
uni := U[i97] - U[j97];
IF ( uni < 0.0 ) THEN uni:=uni + 1.0; END;
U[i97] := uni;
DEC(i97); IF (i97 = 0) THEN i97 := 97; END;
DEC(j97); IF (j97 = 0) THEN j97 := 97; END;
c := c - cd;
IF (c < 0.0) THEN c:=c + cm; END;
uni := uni - c;
IF (uni < 0.0) THEN uni:=uni + 1.0; END;
RETURN uni;
END RanMar;
BEGIN
RanMarInit(1802,9373);
END RandomMar;
MODULE ZUFALL;
IMPORT LowLong;
IMPORT maxSecondParts,DateTime,CanGetClock,GetClock;
IMPORT TIO;
<* IF (__XDS__) THEN *>
IMPORT InitZufallSysT,InitZufall,Zufall; (* EXPORT *)
<* ELSE *>
EXPORT InitZufallSysT,InitZufall,Zufall;
<* END *>
CONST Order = 13;
VAR XV : ARRAY [1..Order] OF LONGREAL;
PROCEDURE InitZufall;
VAR i : CARDINAL;
BEGIN
FOR i:=3 TO Order DO XV[i]:=0.0; END;
XV[1]:=0.14159265E+00;
XV[2]:=0.70710678E+00;
END InitZufall;
PROCEDURE InitZufallSysT;
VAR i : CARDINAL;
dt : DateTime;
T,frac : CARDINAL;
x : REAL;
BEGIN
IF CanGetClock() THEN
GetClock(dt);
frac := dt.fractions + 1;
IF (frac = 1) THEN frac:=maxSecondParts; END;
T:= (dt.second + 1)*maxSecondParts + frac;
x:=VAL(REAL,T);
WHILE (x >= 10.0) DO x:=x / 3.1415927; END;
WHILE (x < 10.0) DO x:=x*10.0; END;
XV[1]:=VAL(LONGREAL,TRUNC(x));
XV[2]:=VAL(LONGREAL,TRUNC(1000.0*(x-VAL(REAL,XV[1]))));
XV[1]:=0.14159265*XV[1];
XV[2]:=0.70710678*XV[2];
FOR i:=3 TO Order DO XV[i]:=0.0; END;
ELSE
InitZufall;
END;
END InitZufallSysT;
PROCEDURE Zufall() : LONGREAL;
VAR i : CARDINAL;
BEGIN
FOR i:=1 TO Order-1 DO
XV[i+1]:=(XV[i+1]+XV[i]) - LowLong.intpart(XV[i+1]+XV[i]);
END;
RETURN XV[Order];
END Zufall;
BEGIN
InitZufall();
END ZUFALL;
PROCEDURE IRandom(a,b : CARDINAL) : CARDINAL; (* MR, 12.02.2015 *)
VAR min : CARDINAL;
BEGIN
IF (b >= a) THEN min:=a; a:=b-a; ELSE min:=b; a:=a-b; END;
RETURN VAL(CARDINAL,(LongMath.round(VAL(LONGREAL,a)*Zufall()))) + min;
END IRandom;
MODULE LuxaryRandom;
IMPORT Errors;
IMPORT TIO;
IMPORT SeedVektor;
<* IF (__XDS__) THEN *>
IMPORT RLuxGo,RLuxAt,RLuxUt,RLuxIn,RanLux; (* EXPORT *)
<* ELSE *>
EXPORT RLuxGo,RLuxAt,RLuxUt,RLuxIn,RanLux;
<* END *>
(*------------------------------------------------------------------------*)
(* All code derived from CERN Program Library SUBROUTINE RANLUX *)
(* Translated from FORTRAN77 to Modula-2 by Michael Riedl, Jan. 2016 *)
(*------------------------------------------------------------------------*)
CONST debug = TRUE;
CONST MaxLev = 4;
LXdflt = 3;
TwoP12 = 4096.0;
JSdflt = 314159265;
IGiga = 1000000000;
Itwo24 = 16777216; (* 2**24 *)
ICons = 2147483563;
VAR ndskip : ARRAY[0..MaxLev] OF INTEGER; (* save *)
seeds : ARRAY[1..24] OF LONGREAL; (* save *)
iseeds : ARRAY[1..24] OF INTEGER; (* save *)
next : ARRAY[1..24] OF INTEGER; (* save *)
notyet : BOOLEAN; (* save *)
luxlev : INTEGER; (* save *)
i24,j24,in24 : INTEGER; (* save *)
kount,mkount : INTEGER; (* save *)
nskip : INTEGER; (* save *)
twom24,twom12 : LONGREAL; (* save *)
carry : LONGREAL; (* save *)
uni : LONGREAL;
inseed,jseed : INTEGER;
PROCEDURE RLuxGo(Lux,InS,K1,K2 : INTEGER);
VAR izip,izip2 : INTEGER;
inner,iouter : INTEGER;
i,k,ilx,isk : INTEGER;
BEGIN
IF (Lux < 0) THEN
luxlev := LXdflt
ELSIF (Lux <= MaxLev) THEN
luxlev := Lux
ELSIF (Lux < 24) OR (Lux > 2000) THEN
luxlev := MaxLev;
Errors.ErrOut('RanLux ILLEGAL Luxury RLuxGo (not in 24 .. 2000) ');
ELSE
luxlev := Lux;
FOR ilx:=0 TO MaxLev DO
IF (Lux = (ndskip[ilx] + 24)) THEN
luxlev := ilx;
END;
END;
END;
IF (luxlev <= MaxLev) THEN
nskip := ndskip[luxlev];
IF debug THEN
TIO.WrStr(' RanLux luxury level set by RLuxGo : ');
TIO.WrInt(luxlev,2); TIO.WrStr(', P = '); TIO.WrInt(nskip+24,1);
TIO.WrLn;
END;
ELSE
nskip := luxlev - 24;
IF debug THEN
TIO.WrStr(' RanLux P-value set by RLuxGo to :');
TIO.WrInt(luxlev,1); TIO.WrLn;
END;
END;
in24 := 0;
IF (InS < 0) THEN
Errors.ErrOut('Illegal initialization by RLuxGo, negative input seed');
END;
IF (InS > 0) THEN
jseed := InS;
IF debug THEN
TIO.WrStr(' RanLux initialized by RLuxGo from seeds ');
TIO.WrInt(jseed,1); TIO.WrChar(",");
TIO.WrInt(K1,1); TIO.WrChar(",");
TIO.WrInt(K2,1); TIO.WrLn;
END;
ELSE
jseed := JSdflt;
IF debug THEN
TIO.WrStr(' RanLux initialized by RLuxGo form default seed');
TIO.WrLn;
END;
END;
inseed := jseed;
notyet:=FALSE;
twom24 := 1.0;
FOR i:=1 TO 24 DO
twom24 := twom24*0.5;
k := jseed DIV 53668;
jseed:=40014*(jseed - k*53668) - k*12211;
IF (jseed < 0) THEN
jseed:=jseed + ICons;
END;
iseeds[i]:= (jseed MOD Itwo24);
END;
twom12 := twom24*4096.0;
FOR i:=1 TO 24 DO;
seeds[i]:=VAL(LONGREAL,iseeds[i])*twom24;
next[i]:=i - 1;
END;
next[1]:=24;
i24 := 24;
j24 := 10;
carry := 0.0;
IF (seeds[24] = 0.0) THEN
carry := twom24;
END;
(* If restarting at a break point, skip K1 + IGiga*K2 *)
(* Note that this is the number of numbers delivered to *)
(* the user PLUS the number skipped (if luxury .GT. 0). *)
kount := K1;
mkount := K2;
IF (K1+K2 # 0) THEN
FOR iouter:=1 TO K2+1 DO
inner := IGiga;
IF (iouter = (K2+1)) THEN
inner := K1;
END;
FOR isk:=1 TO inner DO
uni:=seeds[j24] - seeds[i24] - carry;
IF (uni < 0.0) THEN
uni:=uni + 1.0;
carry := twom24;
ELSE
carry := 0.0;
END;
seeds[i24]:=uni;
i24 := next[i24];
j24 := next[j24];
END;
END;
(* Get the right value of IN24 by direct calculation *)
in24 := kount MOD (nskip+24);
IF (mkount > 0) THEN
izip := IGiga MOD (nskip+24);
izip2 := mkount*izip + in24;
in24 := izip2 MOD (nskip+24);
END;
(* Now IN24 had better be between zero and 23 inclusive *)
IF (in24 > 23) THEN
Errors.ErrOut('Error in restarting with RLuxGo:');
Errors.WriteString(' The values ');
Errors.WriteInt(InS); Errors.WriteChar(",");
Errors.WriteInt(K1); Errors.WriteChar(",");
Errors.WriteInt(K2);
Errors.WriteString(' cannot occur at luxury level ');
Errors.WriteInt(luxlev);
Errors.WriteLn; Errors.WriteLn;
in24 := 0;
END;
END; (* IF (K1+K2 # 0) *)
END RLuxGo;
PROCEDURE RLuxAt(VAR Lout,InOut,K1,K2 : INTEGER);
BEGIN
Lout:=luxlev;
InOut:=inseed;
K1:=kount;
K2:=mkount;
END RLuxAt;
PROCEDURE RLuxUt(VAR ISDext : SeedVektor);
VAR i : INTEGER;
BEGIN
FOR i:=1 TO 24 DO
ISDext[i] := TRUNC(seeds[i]*TwoP12*TwoP12);
END;
ISDext[25] := i24 + 100*j24 + 10000*in24 + 1000000*luxlev;
IF (carry > 0.0) THEN
ISDext[25] := -ISDext[25];
END;
END RLuxUt;
PROCEDURE RLuxIn(VAR ISDext : SeedVektor);
VAR i,isd : INTEGER;
BEGIN
notyet:=FALSE;
twom24 := 1.0;
FOR i:=1 TO 24 DO
next[i] := i-1;
twom24 := twom24*0.5;
END;
next[1] := 24;
twom12 := twom24*4096.0;
IF debug THEN
TIO.WrStr(' Full initialization of RanLux with 25 integers:');
TIO.WrLn;
FOR i:=1 TO 25 DO
TIO.WrInt(ISDext[i],1);
IF ((i MOD 5) = 0) THEN TIO.WrLn; ELSE TIO.WrChar(','); END;
END;
TIO.WrLn;
END;
FOR i:=1 TO 24 DO
seeds[i] := VAL(LONGREAL,ISDext[i])*twom24;
END;
carry := 0.0;
IF (ISDext[25] < 0) THEN
carry := twom24;
END;
isd := ABS(ISDext[25]);
i24 := (isd MOD 100);
isd := isd DIV 100;
j24 := (isd MOD 100);
isd := isd DIV 100;
in24 := (isd MOD 100);
isd := isd DIV 100;
luxlev := isd;
IF (luxlev <= MaxLev) THEN
nskip := ndskip[luxlev];
IF debug THEN
TIO.WrStr(' RanLux Luxury level set by RLuxIn to : ');
TIO.WrInt(luxlev,1); TIO.WrLn;
END;
ELSIF (luxlev >= 24) THEN
nskip := luxlev-24;
IF debug THEN
TIO.WrStr(' RanLux P-value set by RLuxIn to : ');
TIO.WrInt(luxlev,1); TIO.WrLn;
END;
ELSE
nskip := ndskip[MaxLev];
IF debug THEN
TIO.WrStr(' RanLux illegal Luxury RLuxIn : ');
TIO.WrInt(luxlev,1); TIO.WrLn;
END;
luxlev := MaxLev;
END;
inseed := -1;
END RLuxIn;
PROCEDURE RanLux(VAR RVec : ARRAY OF LONGREAL;
LenV : INTEGER);
VAR i,k,ivec,isk : INTEGER;
BEGIN
(* NOTYET is .TRUE. if no initialization has been performed yet. *)
(* Default Initialization by Multiplicative Congruential *)
IF (notyet) THEN
notyet := FALSE;
jseed := JSdflt;
inseed := JSdflt;
IF debug THEN
TIO.WrStr(' RanLux default initialization : ');
TIO.WrInt(jseed,1); TIO.WrLn;
END;
luxlev := LXdflt;
nskip := ndskip[luxlev];
in24 := 0;
kount := 0;
mkount := 0;
IF debug THEN
TIO.WrStr(' RanLux default Luxury level = ');
TIO.WrInt(luxlev,1);
TIO.WrStr(', p ='); TIO.WrInt(nskip+24,1); TIO.WrLn;
END;
twom24 := 1.0;
FOR i:=1 TO 24 DO
twom24 := twom24*0.5;
k := jseed DIV 53668;
jseed := 40014*(jseed-k*53668) - k*12211;
IF (jseed < 0) THEN
jseed:=jseed + ICons;
END;
iseeds[i] := (jseed MOD Itwo24);
END;
twom12 := twom24*4096.0;
FOR i:=1 TO 24 DO
seeds[i] := VAL(LONGREAL,iseeds[i])*twom24;
next[i] := i-1;
END;
next[1] := 24;
i24 := 24;
j24 := 10;
carry := 0.0;
IF (seeds[24] = 0.0) THEN
carry := twom24;
END;
END; (* IF notyet *)
(* The Generator proper: "Subtract-with-borrow", as proposed by *)
(* Marsaglia and Zaman, Florida State University, March, 1989 *)
FOR ivec:=0 TO LenV-1 DO
uni := seeds[j24] - seeds[i24] - carry;
IF (uni < 0.0) THEN
uni:=uni + 1.0;
carry := twom24
ELSE
carry := 0.0;
END;
seeds[i24] := uni;
i24 := next[i24];
j24 := next[j24];
RVec[ivec] := uni;
(* small numbers (with less than 12 "significant" bits) are "padded". *)
IF (uni < twom12) THEN
RVec[ivec]:=RVec[ivec] + twom24*seeds[j24];
(* and zero is forbidden in case someone takes a logarithm *)
IF (RVec[ivec] = 0.0) THEN
RVec[ivec] := twom24*twom24;
END;
END;
(* Skipping to luxury. As proposed by Martin Luscher. *)
in24:=in24 + 1;
IF (in24 = 24) THEN
in24 := 0;
kount:=kount + nskip;
FOR isk:=1 TO nskip DO
uni := seeds[j24] - seeds[i24] - carry;
IF (uni < 0.0) THEN
uni:=uni + 1.0;
carry := twom24;
ELSE
carry := 0.0;
END;
seeds[i24] := uni;
i24 := next[i24];
j24 := next[j24];
END;
END;
END;
kount:=kount + LenV;
IF (kount >= IGiga) THEN
mkount:=mkount + 1;
kount:=kount - IGiga;
END;
END RanLux;
BEGIN
notyet := TRUE; (* save *)
luxlev := LXdflt; (* save xxx *)
in24 := 0; (* save *)
kount := 0; (* save xxx *)
mkount := 0; (* save xxx *)
inseed := JSdflt; (* xxx *)
i24 := 24; (* save *)
j24 := 10; (* save *)
carry := 0.0; (* save *)
ndskip[0] := 0; (* save *)
ndskip[1] := 24;
ndskip[2] := 73;
ndskip[3] := 199;
ndskip[4] := 365;
END LuxaryRandom;
END RandomLib.