Project Nayuki

Fast skipping in a linear congruential generator


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 \left[\left(a^n x_0 \text{ mod } m\right) + \left(\frac{(a^n - 1) \text{ mod } (a -1) m}{a - 1} b\right)\right] \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 runnable code that implements the concepts discussed in the article:


More info