Hyperplane Structure in Linear Congruential Generators

A congruential generator is one multiply, one add and one remainder per output. It has a single word of state, an exactly computable period, and a jump-ahead operation that costs a few multiplies regardless of how far ahead the jump goes. Those properties made it the default generator for decades, and they are still the reason it survives inside libraries that use it only for seeding something better. The defect is geometric and easy to see: plot consecutive pairs and the points arrange themselves into a small number of parallel lines.

The recurrence and what its parameters guarantee

The mixed linear congruential generator is

$$x_{n+1} = (a x_n + c) \bmod m, \qquad u_n = x_n / m.$$

The Hull-Dobell conditions give the parameter choices that make the period exactly $m$: $\gcd(c, m) = 1$, every prime factor of $m$ divides $a - 1$, and $a - 1$ is divisible by 4 when $m$ is. With $m = 2^{32}$ this means $c$ odd and $a \equiv 1 \pmod 4$, which yields a period of $2^{32}$ and a state space that a single 32-bit register covers exactly.

Setting $c = 0$ gives the multiplicative variant, which is cheaper because the addition disappears, but the period collapses. For modulus $2^{e}$ with $e \ge 3$, the maximum period of a multiplicative generator is $2^{e-2}$, so a 32-bit modulus caps the period at $2^{30}$ regardless of the multiplier. The low-order bits are worse than the period suggests. In a mixed generator with modulus $2^{e}$, the least significant bit alternates with period 2, the next bit has period 4, and in general bit $j$ (counting from zero) has period at most $2^{j+1}$. Any code that reduces a full-width output modulo a small range, or extracts a few low bits, is reading the part of the word with the shortest cycle.

Reduction modulo a range that does not divide $m$ introduces a second defect. If $n$ does not divide $2^{32}$, then $2^{32} \bmod n$ residues occur $\lfloor 2^{32}/n \rfloor + 1$ times and the rest occur $\lfloor 2^{32}/n \rfloor$ times, so the excess probability on the favored residues is at most $n / 2^{32}$. For a range of one million that is a relative bias near $2 \times 10^{-4}$, invisible in a histogram but comparable to the probability being estimated when the target is a rare outcome with probability around $10^{-4}$. Rejection on the high bits, or scaling by $2^{-32}$ and flooring, avoids the bias at the cost of a loop or a multiply.

Hyperplanes: why pairs fall on lines

Consider the pair $(u_n, u_{n+1})$ for a multiplicative generator. Since $x_{n+1} \equiv a x_n \pmod m$, the pair satisfies $x_{n+1} - a x_n = k m$ for some integer $k$, so every point lies on one of the lines $y = a x + k m$ in the plane. The number of distinct lines visited is bounded, and the bound is not generous. Marsaglia’s result is that $k$-tuples of successive values lie on at most $(k!, m)^{1/k}$ parallel hyperplanes. For $m = 2^{32}$ and $k = 3$ that is about 1600 planes; for $k = 4$ it is about 330. A three-dimensional plot of triples from a 32-bit generator looks like a stack of sheets, and integrals over the unit cube are being evaluated on that stack rather than on the cube.

The structure is a lattice, not an accident of the multiplier. The set of all $k$-tuples reachable from the recurrence is the intersection of a lattice with the cube, and the lattice is generated by the vectors $(1, a, a^2, \dots)$ reduced modulo $m$ against multiples of $m$ in each coordinate. Changing $a$ changes the lattice, so the defect can be minimized but not eliminated: for any multiplier there is a dimension in which the lattice is coarse, and the coarseness grows with the dimension.

The practical consequence is that a generator can pass a one-dimensional test suite and fail a two-dimensional one. Frequency tests on individual values see a nearly uniform histogram; the correlation between successive values is what carries the lattice signature. This is why generator test batteries include serial and gap tests over overlapping tuples, and why a generator validated for one purpose is not automatically validated for another.

Scoring a multiplier with the spectral test

The spectral test measures the lattice directly. For dimension $k$, it takes the shortest nonzero vector of the dual lattice, the lattice of vectors whose dot product with every lattice point is an integer, and reports its length in units of $1/m$. A long shortest vector means the hyperplanes are far apart, which means the lattice is fine and the tuples fill the cube well. The test is scale-free: the score for a given multiplier does not change when the same multiplier is used with a different modulus of the same word size, apart from the normalization.

Computing the score for $k = 2$ is a continued-fraction expansion of $a/m$. For higher dimensions it is a lattice reduction problem, and the classical algorithms do it in time polynomial in the dimension for fixed $k$. Knuth’s treatment in Volume 2 of The Art of Computer Programming gives the reduction and the tables of good multipliers for moduli near $2^{32}$ and $2^{35}$, which is why those specific constants appear in so many old codes.

A multiplier that scores well in low dimensions can still score badly in higher ones, and the score has to be read against the dimension that the application actually uses. A simulation that consumes uniforms in blocks of 100 and forms a statistic from each block is sensitive to structure in dimension 100, where every 32-bit congruential generator is hopeless. The number of hyperplanes at that dimension is small enough that the block statistic is a deterministic function of a few lattice coordinates.

State size and what a long period costs

The Mersenne Twister replaced the congruential family as the default because it fixes the lattice problem in the dimensions that matter while keeping a long period. Its state is 624 32-bit words, 2496 bytes, and its period is $2^{19937} - 1$, a Mersenne prime. The recurrence is linear over GF(2), and the generator is 623-dimensionally equidistributed: over the full period, every 32-bit pattern appears equally often in each of the first 623 coordinates of the output tuple. Equidistribution at that dimension is the property the congruential generators cannot approach.

The costs are concrete. The state does not fit in a few registers, so a generator per thread multiplies the footprint by 2.5 KB. The output is 32 bits, so a double-precision uniform needs two draws and a shift, which is why variants exist that emit 53-bit doubles directly. And the linearity over GF(2) is exploitable: the generator passes the standard statistical batteries but fails tests designed to detect linear dependence, and it is trivially predictable from a small number of outputs.

For parallel work the state size becomes a correctness issue rather than a performance one. Two threads that seed the same generator with the same value produce identical streams; two threads that seed with adjacent values produce streams whose initial outputs are correlated through the seeding recurrence. The reliable constructions either jump ahead by a fixed offset in the period, which requires a generator with an efficient jump operation, or abandon the sequence model entirely.

Splitting a stream without correlating it

Counter-based generators change the interface. Instead of carrying a state forward, they compute an output as a keyed bijection of a counter: output $= f(\text{key}, \text{counter})$, where $f$ is a small number of multiply-xor rounds. The counter is the position, the key is the stream identifier, and there is no state to advance. Skipping ahead is assigning a counter. Two streams are independent if their keys differ, and the guarantee is only as good as the mixing function, so keys are usually derived by hashing a master seed together with a stream index rather than by adding one.

The trade is that a counter-based generator cannot be used where a single sequential stream with a documented period is required, because the period is the width of the counter and the sequence is not a single cycle in the classical sense. It also costs more per output than a one-word congruential step, since the mixing function is several multiplies. What it buys is reproducibility under parallel decomposition: the value at position $n$ of stream $k$ does not depend on how the work was partitioned, so a run split across a different number of threads produces bit-identical results.

def jump(a, c, m, k, x):
    """Advance the LCG x -> (a*x + c) mod m by k steps in O(log k)."""
    am, cm = 1, 0                      # accumulated map: x -> am*x + cm
    while k:
        if k & 1:
            am, cm = am * a % m, (am * c + cm) % m
        a, c = a * a % m, (a * c + c) % m   # square the base map
        k >>= 1
    return (am * x + cm) % m

Jump-ahead makes subsequence splitting practical: give thread $i$ the starting state $T^{i \cdot L}(x_0)$ for a block length $L$, and the streams are disjoint segments of one cycle. The catch is that the segments are not independent in the statistical sense, only non-overlapping, and the lattice structure of the underlying generator applies within each segment. For a generator with a 32-bit state, disjoint segments of a bad lattice are still bad; for a generator with a long period and good equidistribution, the segments behave acceptably.

Whatever the generator, the uniforms it produces enter the sampling machinery with their defects intact. The inverse-transform construction maps a correlated tuple to a correlated tuple of the target distribution, and the variance of a Monte Carlo estimator is computed as if the draws were independent, so any residual correlation shows up as an understated error bar rather than as an obvious bug.