diff --git a/src/main/java/net/berack/upo/valpre/rand/Rng.java b/src/main/java/net/berack/upo/valpre/rand/Rng.java
new file mode 100644
index 0000000..37676d0
--- /dev/null
+++ b/src/main/java/net/berack/upo/valpre/rand/Rng.java
@@ -0,0 +1,99 @@
+package net.berack.upo.valpre.rand;
+
+/**
+ * This is an Java library for random number generation.  The use of this
+ * library is recommended as a replacement for the Java class Random,
+ * particularly in simulation applications where the statistical
+ * 'goodness' of the random number generator is important.
+ *
+ * The generator used in this library is a so-called 'Lehmer random number
+ * generator' which returns a pseudo-random number uniformly distributed
+ * between 0.0 and 1.0.  The period is (m - 1) where m = 2,147,483,647 and
+ * the smallest and largest possible values are (1 / m) and 1 - (1 / m)
+ * respectively.  For more details see:
+ *
+ *       "Random Number Generators: Good Ones Are Hard To Find"
+ *                   Steve Park and Keith Miller
+ *              Communications of the ACM, October 1988
+ *
+ * Note that as of 7-11-90 the multiplier used in this library has changed
+ * from the previous "minimal standard" 16807 to a new value of 48271.  To
+ * use this library in its old (16807) form change the constants MULTIPLIER
+ * and CHECK as indicated in the comments.
+ *
+ * Name              : Rng.java  (Random Number Generation - Single Stream)
+ * Authors           : Steve Park & Dave Geyer
+ * Translated by     : Jun Wang & Richard Dutton
+ * Language          : Java
+ * Latest Revision   : 6-10-04
+ *
+ * Program rng       : Section 2.2
+ */
+public class Rng {
+	public final static long CHECK = 399268537L; /* use 1043616065 for the "minimal standard" */
+	public final static long DEFAULT = 123456789L; /* initial seed, use 0 < DEFAULT < MODULUS */
+	public final static long MODULUS = 2147483647; /* DON'T CHANGE THIS VALUE */
+	public final static long MULTIPLIER = 48271; /* use 16807 for the "minimal standard" */
+
+	private long seed = DEFAULT; /* seed is the state of the generator */
+
+	public Rng() {
+	}
+
+	public Rng(long seed) {
+		this.putSeed(seed);
+	}
+
+	/**
+	 * Random is a Lehmer generator that returns a pseudo-random real number
+	 * uniformly distributed between 0.0 and 1.0. The period is (m - 1)
+	 * where m = 2,147,483,647 amd the smallest and largest possible values
+	 * are (1 / m) and 1 - (1 / m) respectively.
+	 */
+	public double random() {
+		this.seed = Rng.newSeed(MODULUS, MULTIPLIER, this.seed);
+		return ((double) this.seed / MODULUS);
+	}
+
+	/**
+	 * Use this (optional) procedure to initialize or reset the state of
+	 * the random number generator according to the following conventions:
+	 * if x > 0 then x is the initial seed (unless too large)
+	 * if x <= 0 then the initial seed is obtained from the system clock
+	 */
+	public void putSeed(long x) {
+		if (x > 0L) {
+			x = x % MODULUS; /* correct if x is too large */
+		} else {
+			x = System.currentTimeMillis();
+			// x = ((unsigned long) time((time_t *) NULL)) % MODULUS;
+		}
+
+		this.seed = x;
+	}
+
+	/**
+	 * Use this (optional) procedure to get the current state of the random
+	 * number generator.
+	 */
+	public long getSeed() {
+		return this.seed;
+	}
+
+	/**
+	 * Use this (optional) procedure to test for a correct implementation.
+	 */
+	public static boolean testRandom() {
+		Rng rng = new Rng(1); /* set initial state to 1 */
+		for (var i = 0; i < 10000; i++)
+			rng.random();
+		return rng.getSeed() == CHECK;
+	}
+
+	public static long newSeed(long modulus, long multiplier, long seed) {
+		var Q = modulus / multiplier;
+		var R = modulus % multiplier;
+		var t = multiplier * (seed % Q) - R * (seed / Q);
+		return t > 0 ? t : (t + modulus);
+	}
+}
diff --git a/src/main/java/net/berack/upo/valpre/rand/Rngs.java b/src/main/java/net/berack/upo/valpre/rand/Rngs.java
new file mode 100644
index 0000000..691248f
--- /dev/null
+++ b/src/main/java/net/berack/upo/valpre/rand/Rngs.java
@@ -0,0 +1,95 @@
+package net.berack.upo.valpre.rand;
+
+/**
+ * This is an Java library for multi-stream random number generation.
+ * The use of this library is recommended as a replacement for the Java
+ * class Random, particularly in simulation applications where the
+ * statistical 'goodness' of the random number generator is important.
+ * The library supplies 256 streams of random numbers; use
+ * selectStream(s) to switch between streams indexed s = 0,1,...,255.
+ *
+ * The streams must be initialized.  The recommended way to do this is by
+ * using the function plantSeeds(x) with the value of x used to initialize
+ * the default stream and all other streams initialized automatically with
+ * values dependent on the value of x.  The following convention is used
+ * to initialize the default stream:
+ *    if x > 0 then x is the state
+ *    if x < 0 then the state is obtained from the system clock
+ *    if x = 0 then the state is to be supplied interactively.
+ *
+ * The generator used in this library is a so-called 'Lehmer random number
+ * generator' which returns a pseudo-random number uniformly distributed
+ * 0.0 and 1.0.  The period is (m - 1) where m = 2,147,483,647 and the
+ * smallest and largest possible values are (1 / m) and 1 - (1 / m)
+ * respectively.  For more details see:
+ *
+ *       "Random Number Generators: Good Ones Are Hard To Find"
+ *                   Steve Park and Keith Miller
+ *              Communications of the ACM, October 1988
+ *
+ * Name            : Rngs.java  (Random Number Generation - Multiple Streams)
+ * Authors         : Steve Park & Dave Geyer
+ * Translated by   : Jun Wang & Richard Dutton
+ * Language        : Java
+ * Latest Revision : 6-10-04
+ */
+class Rngs {
+	private final static int STREAMS = 256; /* # of streams, DON'T CHANGE THIS VALUE */
+	private final static long A256 = 22925; /* jump multiplier, DON'T CHANGE THIS VALUE */
+
+	private Rng[] rngs;
+	private int current = 0;
+
+	public Rngs(long seed) {
+		this.rngs = new Rng[STREAMS];
+		this.plantSeeds(seed);
+	}
+
+	public Rng getRng() {
+		return this.rngs[this.current];
+	}
+
+	public Rng getRng(int stream) {
+		return this.rngs[stream % STREAMS];
+	}
+
+	/**
+	 * Use this function to set the state of all the random number generator
+	 * streams by "planting" a sequence of states (seeds), one per stream,
+	 * with all states dictated by the state of the default stream.
+	 * The sequence of planted states is separated one from the next by
+	 * 8,367,782 calls to Random().
+	 */
+	public void plantSeeds(long seed0) {
+		this.rngs[0] = new Rng(seed0);
+
+		for (int j = 1; j < STREAMS; j++) {
+			seed0 = Rng.newSeed(Rng.MODULUS, A256, seed0);
+			this.rngs[j] = new Rng(seed0);
+		}
+	}
+
+	/**
+	 * Use this function to set the current random number generator
+	 * stream -- that stream from which the next random number will come.
+	 */
+	public void selectStream(int index) {
+		this.current = index % STREAMS;
+	}
+
+	/**
+	 * Use this (optional) function to test for a correct implementation.
+	 */
+	public static boolean testRandom() {
+		var rngs = new Rngs(1);
+		var first = rngs.getRng();
+		for (int i = 0; i < 10000; i++)
+			first.random();
+
+		var x = first.getSeed(); /* get the new state value */
+		var ok = (x == Rng.CHECK); /* and check for correctness */
+
+		x = rngs.getRng(1).getSeed(); /* get the state of stream 1 */
+		return ok && (x == A256); /* x should be the jump multiplier */
+	}
+}
diff --git a/src/main/java/net/berack/upo/valpre/rand/Rvgs.java b/src/main/java/net/berack/upo/valpre/rand/Rvgs.java
new file mode 100644
index 0000000..c315d80
--- /dev/null
+++ b/src/main/java/net/berack/upo/valpre/rand/Rvgs.java
@@ -0,0 +1,209 @@
+package net.berack.upo.valpre.rand;
+
+/* -------------------------------------------------------------------------- 
+ * This is a Java library for generating random variates from six discrete 
+ * distributions
+ *
+ *      Generator         Range (x)     Mean         Variance
+ *
+ *      bernoulli(p)      x = 0,1       p            p*(1-p)
+ *      binomial(n, p)    x = 0,...,n   n*p          n*p*(1-p)
+ *      equilikely(a, b)  x = a,...,b   (a+b)/2      ((b-a+1)*(b-a+1)-1)/12
+ *      Geometric(p)      x = 0,...     p/(1-p)      p/((1-p)*(1-p))
+ *      pascal(n, p)      x = 0,...     n*p/(1-p)    n*p/((1-p)*(1-p))
+ *      poisson(m)        x = 0,...     m            m
+ * 
+ * and seven continuous distributions
+ *
+ *      uniform(a, b)     a < x < b     (a + b)/2    (b - a)*(b - a)/12 
+ *      exponential(m)    x > 0         m            m*m
+ *      erlang(n, b)      x > 0         n*b          n*b*b
+ *      normal(m, s)      all x         m            s*s
+ *      logNormal(a, b)   x > 0            see below
+ *      chiSquare(n)      x > 0         n            2*n 
+ *      student(n)        all x         0  (n > 1)   n/(n - 2)   (n > 2)
+ *
+ * For the a Lognormal(a, b) random variable, the mean and variance are
+ *
+ *                        mean = exp(a + 0.5*b*b)
+ *                    variance = (exp(b*b) - 1) * exp(2*a + b*b)
+ *
+ * Name              : Rvgs.java  (Random Variate GeneratorS)
+ * Authors           : Steve Park & Dave Geyer
+ * Translated by     : Richard Dutton & Jun Wang
+ * Language          : Java
+ * Latest Revision   : 7-1-04
+ * --------------------------------------------------------------------------
+ */
+public class Rvgs {
+
+    private final Rng rng;
+
+    // public Rvgs() {
+    // this.rngs = new Rngs(Rng.DEFAULT);
+    // }
+
+    public Rvgs(Rng rng) {
+        if (rng == null)
+            throw new NullPointerException();
+        this.rng = rng;
+    }
+
+    /**
+     * Returns 1 with probability p or 0 with probability 1 - p.
+     * NOTE: use 0.0 < p < 1.0
+     */
+    public long bernoulli(double p) {
+        return ((this.rng.random() < (1.0 - p)) ? 0 : 1);
+    }
+
+    /**
+     * Returns a binomial distributed integer between 0 and n inclusive.
+     * NOTE: use n > 0 and 0.0 < p < 1.0
+     */
+    public long binomial(long n, double p) {
+        long i, x = 0;
+
+        for (i = 0; i < n; i++)
+            x += bernoulli(p);
+        return (x);
+    }
+
+    /**
+     * Returns an equilikely distributed integer between a and b inclusive.
+     * NOTE: use a < b
+     */
+    public long equilikely(long a, long b) {
+        return (a + (long) ((b - a + 1) * this.rng.random()));
+    }
+
+    /**
+     * Returns a geometric distributed non-negative integer.
+     * NOTE: use 0.0 < p < 1.0
+     */
+    public long geometric(double p) {
+        return ((long) (Math.log(1.0 - this.rng.random()) / Math.log(p)));
+    }
+
+    /**
+     * Returns a Pascal distributed non-negative integer.
+     * NOTE: use n > 0 and 0.0 < p < 1.0
+     */
+    public long pascal(long n, double p) {
+        long i, x = 0;
+
+        for (i = 0; i < n; i++)
+            x += geometric(p);
+        return (x);
+    }
+
+    /**
+     * Returns a Poisson distributed non-negative integer.
+     * NOTE: use m > 0
+     */
+    public long poisson(double m) {
+        double t = 0.0;
+        long x = 0;
+
+        while (t < m) {
+            t += exponential(1.0);
+            x++;
+        }
+        return (x - 1);
+    }
+
+    /**
+     * Returns a uniformly distributed real number between a and b.
+     * NOTE: use a < b
+     */
+    public double uniform(double a, double b) {
+        return (a + (b - a) * this.rng.random());
+    }
+
+    /**
+     * Returns an exponentially distributed positive real number.
+     * NOTE: use m > 0.0
+     */
+    public double exponential(double m) {
+        return (-m * Math.log(1.0 - this.rng.random()));
+    }
+
+    /**
+     * Returns an Erlang distributed positive real number.
+     * NOTE: use n > 0 and b > 0.0
+     */
+    public double erlang(long n, double b) {
+        long i;
+        double x = 0.0;
+
+        for (i = 0; i < n; i++)
+            x += exponential(b);
+        return (x);
+    }
+
+    /**
+     * Returns a normal (Gaussian) distributed real number.
+     * NOTE: use s > 0.0
+     *
+     * Uses a very accurate approximation of the normal idf due to Odeh & Evans,
+     * J. Applied Statistics, 1974, vol 23, pp 96-97.
+     */
+    public double normal(double m, double s) {
+        final double p0 = 0.322232431088;
+        final double q0 = 0.099348462606;
+        final double p1 = 1.0;
+        final double q1 = 0.588581570495;
+        final double p2 = 0.342242088547;
+        final double q2 = 0.531103462366;
+        final double p3 = 0.204231210245e-1;
+        final double q3 = 0.103537752850;
+        final double p4 = 0.453642210148e-4;
+        final double q4 = 0.385607006340e-2;
+        double u, t, p, q, z;
+
+        u = this.rng.random();
+        if (u < 0.5)
+            t = Math.sqrt(-2.0 * Math.log(u));
+        else
+            t = Math.sqrt(-2.0 * Math.log(1.0 - u));
+        p = p0 + t * (p1 + t * (p2 + t * (p3 + t * p4)));
+        q = q0 + t * (q1 + t * (q2 + t * (q3 + t * q4)));
+        if (u < 0.5)
+            z = (p / q) - t;
+        else
+            z = t - (p / q);
+        return (m + s * z);
+    }
+
+    /**
+     * Returns a lognormal distributed positive real number.
+     * NOTE: use b > 0.0
+     */
+    public double logNormal(double a, double b) {
+        return (Math.exp(a + b * normal(0.0, 1.0)));
+    }
+
+    /**
+     * Returns a chi-square distributed positive real number.
+     * NOTE: use n > 0
+     */
+    public double chiSquare(long n) {
+        long i;
+        double z, x = 0.0;
+
+        for (i = 0; i < n; i++) {
+            z = normal(0.0, 1.0);
+            x += z * z;
+        }
+        return (x);
+    }
+
+    /**
+     * Returns a student-t distributed real number.
+     * NOTE: use n > 0
+     */
+    public double student(long n) {
+        return (normal(0.0, 1.0) / Math.sqrt(chiSquare(n) / n));
+    }
+
+}