# Stream-Based Linear Congruential Generator in Java 8

A while ago, I wrote a post on a Stream-Based Linear Congruential Generator in Scala. This post is a similar implementation using the Java 8 Stream API.

Here is a Java 8 stream-based implementation of a pseudo-random number generator using a linear congruential generator (LCG):

package com.michaelcotterell.util;

import java.util.stream.DoubleStream;
import java.util.stream.IntStream;
import java.util.stream.LongStream;

public class Random {

private final long a;    // multiplier
private final long c;    // increment
private final long m;    // modulus
private final long seed; // start value

private final LongStream stream;

public Random(long a, long c, long m, long seed) {
this.a      = a;
this.c      = c;
this.m      = m;
this.seed   = seed % m;
this.stream = LongStream.iterate(this.seed, x -> (a * x + c) % m);
} // Random

public Random(long a, long c, long m) {
this(a, c, m, System.currentTimeMillis());
} // Random

public LongStream longs() {
return stream.map(e -> e);
} // longs

public IntStream ints() {
return stream.mapToInt(e -> (int) (e % (Integer.MAX_VALUE + 1L)));
} // ints

public DoubleStream doubles() {
final double ONE_OVER_M = 1.0 / m;
return stream.mapToDouble(e -> e * ONE_OVER_M);
} // doubles

} // Random

To use this class, you might do something like the following:

Random lcg = new Random(48271, 0, (2L << 31) - 1, System.currentTimeMillis());
Iterator<Double> rand = lcg.doubles().iterator();

for (int i = 0; i < 10; ++i) System.out.println(rand.next()); 

To test different LCGs, you might do the following:

long iters = 100_000;
long seed  = System.currentTimeMillis();

Random[] lcgs = new Random[] {
new Random(48271, 0, (2L << 31) - 1, seed), // MINSTD (updated)
new Random(16807, 0, (2L << 31) - 1, seed), // MINSTD (old)
new Random(65539, 0, (2L << 31),     seed)  // RANDU
};

for (Random lcg : lcgs) {
System.out.println(lcg.doubles().limit(iters).summaryStatistics());
} // for

# Purely Tail Recursive k-th Smallest Element

Finding the k-th smallest (or largest) element of an unordered array is a problem that's been approached by many programmers and computer scientists. There are many methods for solving this problem. To me, one of the more interesting implementations is a purely tail recursive implementation that uses the partition algorithm.

A purely tail recursive function or method is one in which there are no deferred operations. A deferred operation in this case, is one that must wait until a recursive call returns before it can complete. This article goes into this a little detail about the differences between "ordinary" and "pure" tail recursive methods.

Given a function that can swap two elements of an array and the partition function, it is possible to find the k-th smallest element of an array of distinct integers using the following purely tail recursive implementation (the accumulator parameters act in a fashion that is similar to a binary search):

private static int kthMin(int[] a, int left, int right, int k) {
int newPivot = partition(a, left, right, k - 1);
if (newPivot == k - 1) return a[k - 1];
if (newPivot < k - 1) left = newPivot + 1;
if (newPivot > k - 1) right = newPivot - 1;
return kthMin(a, left, right, k);
} // kthMin

If you have any questions, then please post them in the comments.

# Games vs. Simulations

This is not a full blog post. I just liked the differentiation between the terms "gaming" and "simulation" presented in "Applied Mathematical Programming" by Bradley, Hax, and Magnanti (Addison-Wesley, 1977).

[In] gaming [,...] a model is constructed that is an abstract and simplified representation of the real environment. This model provides a responsive mechanism to evaluate the effectiveness of proposed alternatives, which the decision-maker must supply in an organized and sequential fashion. The model is simply a device that allows the decision-maker to test the performance of the various alternatives that seem worthwhile to pursue.

v.s.

Simulation models are similar to gaming models except that all human decision-makers are removed from the modeling process. The model provides the means to evaluate the performance of a number of alternatives, supplied externally to the model by the decision-maker, without allowing for human interactions at intermediate stages of the model computation.

# Multinomial Coefficients in Scala using Fold

Multinomial coefficients are often denoted as $$\left( n;~ k_1, \dots, k_m \right)$$ where $$n = \sum_{i=1}^m k_i$$. Another way to express a multinomial coefficient is using the $${\rm choose}$$ notation commonly employed in combinatorics for binomial coefficients. Given $$n \geq \sum_{i=1}^{m-1} k_i$$, we say that the following are equivalent:

$$\left( n;~ k_1, \dots, k_m \right) ~\mbox{where}~ k_m = n - \sum_{i=1}^{m-1} k_i \equiv {n \choose k_1, \dots, k_{m-1}}$$

For convenience, I'm going to use the first notation.

### Theorem

I'm not sure if this has ever been stated before, but here's an observation that I made a couple months ago:

$$\left( n;~ k_1, \dots, k_m \right) = \left( n-k_m;~ k_1, \dots, k_{m-1} \right) \cdot \frac{1}{k_m !} \prod_{i=1}^{k_m} n-i+1$$

That is, there exists a linear, tail recursive definition for calculating multinomial coefficients.

### Proof of Theorem

Without loss of generality, assume that $$n = \sum_{i=1}^{m} k_i$$. According to the equivalence noted above as well as the Multinomial Theorem, we have:

$$\left( n;~ k_1, \dots, k_m \right) = {n \choose k_1, \dots, k_m} = \frac{ n! }{ k_1 ! \cdot \cdots \cdot k_m ! }$$

(1) If we pull out the term $$k_m$$, we get:

$$\left( n;~ k_1, \dots, k_m \right) = \frac{1}{k_m !} \cdot \frac{ n! }{ k_1 ! \cdot \cdots \cdot k_{m-1} ! }$$

(2) Now, let's find some $$\zeta$$ such that:

$$\frac{ (n-k_m)! }{ k_1 ! \cdot \cdots \cdot k_{m-1} ! } \cdot \zeta = \frac{ n! }{ k_1 ! \cdot \cdots \cdot k_{m-1} ! }$$

(3) If we isolate $$\zeta$$ on the left-hand side, we observe the following:

$$\zeta = \frac{ n! }{ k_1 ! \cdot \cdots \cdot k_{m-1} ! } \cdot \frac{ k_1 ! \cdot \cdots \cdot k_{m-1} ! }{ (n-k_m)! } = \frac{ n! }{ (n-k_m)! }$$

(4) We know that certain terms cancel in this fraction. It is easy to see that

$$\zeta = \frac{ n! }{ (n-k_m)! } = (n - k_m + 1) \cdot (n - k_m + 2) \cdot \cdots \cdot n = \prod_{i=1}^{k_m} n - k_m + i$$

(5) Since the product in $$\zeta$$ is from $$1$$ to $$k_m$$, we can be further simplify the expression to

$$\zeta = \prod_{i=1}^{k_m} n - k_m + i = \prod_{i=1}^{k_m} n-i+1$$

(6) Replacing $$\zeta$$ in (2) gives us the following:

$$\frac{ (n-k_m)! }{ k_1 ! \cdot \cdots \cdot k_{m-1} ! } \cdot \prod_{i=1}^{k_m} n-i+1 = \frac{ n! }{ k_1 ! \cdot \cdots \cdot k_{m-1} ! }$$

(7) Now we use the result of (6) with (1) to get:

$$\left( n;~ k_1, \dots, k_m \right) = \frac{1}{k_m !} \cdot \frac{ (n-k_m)! }{ k_1 ! \cdot \cdots \cdot k_{m-1} ! } \cdot \prod_{i=1}^{k_m} n-i+1$$

(8) Since the second fraction in (7) can be expressed as a multinomial, we have:

$$\left( n;~ k_1, \dots, k_m \right) = \left( n-k_m;~ k_1, \dots, k_{m-1} \right) \cdot \frac{1}{k_m !} \prod_{i=1}^{k_m} n-i+1$$

Therefore, the theorem is correct.

### Some Notes

This result that we observed in (7) directly relates to the ratios between the coefficients on the face of a multidimensional Pascal triangle.

### Corollary

It is easy to see that the following follows from the theorem:

$$\left( n;~ k_1, \dots, k_m \right) = \left( n-k_m;~ k_1, \dots, k_{m-1} \right) \cdot \prod_{i=1}^{k_m} \frac{n-i+1}{i}$$

### Folding with scala

In Scala, the expression in the corollary can be expressed concisely using the foldLeft operator.This operator applies a binary operator to a start value and all elements of a list, going left to right. If we let our list be the numbers $$1, \dots, k_m$$, then we can foldLeft with the following binary operator:

$$\left(a,b\right) \Rightarrow \frac{a \left(k_m-b+1\right)}{b}$$

For more information on the general fold operator (with which you can pretty much do anything with), see "Introduction to Metamathematics" by Kleene (affiliate link).

### The Code

Here is my recursive implementation of the recurrence relation mentioned in the corollary above:

I've omitted a an iterative implementation of the corollary that I've written for the sake of brevity. If anyone expresses some interest, I may include here or in a future blog post.

object Combinatorics {     /** Computes the multinomial coefficient (n; k_1, .., k_m)      *  where n = the sum of k_1, .., k_m.     *      *  This is a variadic convenience function that allows      *  someone to invoke <code>multi</code> without using an      *  array. Note, however, that the variadic parameter is      *  transformed into an array when this version of      *  <code>multi</code> is invoked.     *      *  @author Michael E. Cotterell <mepcotterell@gmail.com>     *  @param  k items to choose     */    def multi (k: Int*): Long = multi(k.toArray)     /** Computes the multinomial coefficient (n; k_1, .., k_m)      *  where n = the sum of k_1, .., k_m.      *      *  This implementation requires exactly n-many      *  multiplications and m-many recursive calls.     *       *  @see    http://michaelcotterell.com/2013/03/multinomial-coefficients/      *  @author Michael E. Cotterell <mepcotterell@gmail.com>     *  @param  k items to choose     */    def multi (k: Array[Int]): Long =     {        if (k.length == 1) 1L        else {            (1 to k.last).foldLeft(multi(k.slice(0, k.length - 1))){                (prev, i) => prev * (k.sum - i + 1) / i            }        } // if    } // multi     /** Computes the multinomial coefficient (n; k_1, .., k_m)      *  where n = the sum of k_1, .., k_m.     *      *  This implementation requires exactly n-many      *  multiplications and m-many recursive calls. Also, it is      *  experimentally slower than the <code>foldLeft</code>      *  implementation provided by the <code>multi</code>      *  function.     *       *  @see    http://michaelcotterell.com/2013/03/multinomial-coefficients/     *  @author Michael E. Cotterell <mepcotterell@gmail.com>     *  @param  k items to choose     */    def _multi (k: Array[Int]): Long =     {        if (k.length == 1) 1L        else {            var product = _multi(k.slice(0, k.length-1))            for(i <- 1 to k.last) product = product * (k.sum-i+1) / i                        product        } // if    } // _multi } // Combinatorics