--- a
+++ b/RandomLib.mod.m2pp
@@ -0,0 +1,657 @@
+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).                                           *)
+  (*------------------------------------------------------------------------*)
+  (* 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.