# Fast skipping in a linear congruential generator

## Introduction

A linear congruential generator (LCG) is a pseudorandom number generation (PRNG) algorithm. Its behavior is defined by the following recurrence relation:

\begin{align} x_0 &= \text{seed}. \\ x_{k+1} &\equiv (a x_k + b) \mod m \text{, for all k \in \mathbb{N}}. \end{align}

The integer constants $$a$$, $$b$$, $$m$$ are chosen carefully to give good properties. For example: having a full period of $$m$$, having $$m$$ be a power of 2 (which is efficient on computers), or having $$b$$ be 0 (to save an operation). The seed value can be chosen (almost) arbitrarily for each random number generator instance. The sequence of $$x$$ values are the pseudorandom numbers that are generated. They can be used as-is, or they can be filtered and have their distribution altered.

## Fast skipping

Suppose we know $$x_0$$ and we want to determine $$x_n$$ for some $$n \in \mathbb{N}$$. Clearly it is possible to iterate the recurrence formula $$n$$ times to get our desired result, which takes $$\Theta(n)$$ time. But with some clever arithmetic, we can get the result much faster, in only $$\Theta(\log n)$$ time.

Consider the special simple case where $$b = 0$$. Then $$x_n \equiv a^n x_0 \text{ mod } m$$. We can compute $$a^n x_0 \text{ mod } m$$ quickly using the well-known modular exponentiation algorithm, which is exponentiation by squaring with a reduction of each intermediate result modulo $$m$$.

Now for the general case, first consider the recurrence relation $$x_{k+1} = a x_k + b$$, which is the LCG recurrence without the modulo. The closed-form solution[0] to this first-order linear recurrence is $$x_k = a^n x_0 + \frac{a^n - 1}{a - 1} b$$.

Reinserting the modulos, we have $$x_k \equiv ((a^n x_0 \text{ mod } m) + (\frac{a^n - 1}{a - 1} b)) \text{ mod } m$$. The first term, $$a^n x_0 \text{ mod } m$$, is the modular exponentiation discussed above. The second term is trickier. The quotient $$\frac{a^n - 1}{a - 1}$$ is always an integer, because it is equal[1] to $$a^{n-1} + a^{n-2} + \cdots + a^1 + a^0$$. To compute $$\frac{a^n - 1}{a - 1} \text{ mod } m$$, we compute $$(a^n - 1) \text{ mod } (a - 1) m$$, then non-modular-divide[2] by $$a - 1$$.

Altogether, we have: $$x_k \equiv ((a^n x_0 \text{ mod } m) + (\frac{(a^n - 1) \text{ mod } (a -1) m}{a - 1} b)) \text{ mod } m$$

## Backward iteration

By rearranging the recurrence relation, we get $$x_k \: \equiv \: a^{-1} (x_{k+1} - b) \: \equiv \: a^{-1} x_{k+1} + (a^{-1})(-b) \mod m$$. As long as $$a$$ and $$m$$ are coprime, the inverse of $$a$$ exists and can be computed by the extended Euclidean algorithm.

Incidentally, we can do fast backward skipping by using $$(a', b', m)$$ with $$a' \equiv a^{-1} \text{ mod } m$$ and $$b' \equiv -(a^{-1} b) \text{ mod } m$$.

## Source code

To support my claims about fast seeking and backward iteration, here is live Java code that implements the concepts discussed in the article: LcgRandom.java

## Notes

• [0]: Finding this closed-form solution requires some ingenuity, a table of algebraic expansions, or a computer algebra system (CAS). But now that the solution has been stated, you can prove its correctness by induction quite easily.

• [1]: Proof is left as an exercise for the reader.

• [2]: I will not attempt to justify why non-modular division within modular arithmetic works, because I don’t fully understand it myself.