The Entropy Problem: A Study of Pseudorandom Number Generator Theory
A technical survey of pseudorandom number generator theory and statistical validation
Abstract
The mathematical integrity of any application relying on pseudorandomness depends critically on properties of the underlying pseudorandom number generator (PRNG) that are frequently underspecified in practice. This article surveys the structural properties of F₂-linear generators and their scramblers, including algebraic normal form analysis, equidistribution theory, and the word polynomial framework for computing characteristic polynomials of non-commutative linear engines. We examine the BigCrush battery of TestU01 and the Hamming-weight dependency (HWD) test, drawing directly on Blackman and Vigna [1,2]. We document the statistical weaknesses of the Mersenne Twister relative to the current state of the art, and address the seeding, parallel simulation, and floating-point conversion problems as independent sources of error. All numerical results are drawn from the primary literature with explicit source citations.
Keywords: pseudorandom number generation, F₂-linear generators, xoshiro, xoroshiro, scrambler, algebraic normal form, BigCrush, Hamming-weight dependencies
1. Introduction
Randomness is more naturally defined by absence — the lack of any detectable pattern or predictability — than by any positive property it possesses [13]. A pseudorandom number generator is a deterministic program designed to produce output sequences that are statistically indistinguishable from truly random sequences under a defined set of tests. In any application relying on pseudorandomness, the PRNG is the sole source of stochasticity in the model, so the statistical properties of the PRNG directly constrain the validity of any conclusion derived from it.
The research literature on PRNG quality has advanced substantially since the introduction of the Mersenne Twister in 1998 [5]. The Mersenne Twister, despite its documented statistical weaknesses, remains the default PRNG in Python, R, MATLAB, and a large fraction of general-purpose software. This article presents the current state of the research, drawing principally on Blackman and Vigna [1,2], L'Ecuyer and Simard [6], and associated primary literature. All claims are cited to their source; no results beyond what the cited literature establishes are asserted.
2. The Mathematical Structure of Pseudorandom Number Generators
2.1 Formal Definition
Following Blackman and Vigna [1], a w-bit PRNG is defined by a finite state space S, a transition function τ : S → S, and an output function φ : S → {0,1}ʷ. Given an initial state σ ∈ S (the seed), the generator produces:
φ(σ), φ(τ(σ)), φ(τ²(σ)), φ(τ³(σ)), …
The sequence is entirely determined by σ. Critical properties for practical use are: full period (the generator visits every non-zero state before repeating), statistical uniformity of the output distribution, low correlation between successive outputs, and reproducibility across platforms.
2.2 F₂-Linear Generators
In F₂-linear generators, the state space S = (F₂)ⁿ is an n-dimensional vector space over F₂ = {0,1}, and the transition function τ is an F₂-linear transformation representable as a matrix M ∈ F₂ⁿˣⁿ:
τ(s) = Ms (mod 2)
The generator has full period 2ⁿ − 1 if and only if the characteristic polynomial p(x) = det(M − xI) is primitive over F₂ [1]. Primitivity requires that p(x) be irreducible and that x have maximum order in F₂[x]/p(x). The weight of p(x) is a quality indicator: Compagner [4] establishes that weight close to half the degree is optimal.
Every bit of the state satisfies the linear recurrence defined by p(x). This linear structure is detectable by the binary rank test (MatrixRank) and the linear complexity test (LinearComp) in BigCrush [6]. All F₂-linear generators fail these two tests without exception. Blackman and Vigna note that the community considers the advantage of provable mathematical properties to outweigh these linearity test failures [1].
2.3 Scramblers: Theory and Algebraic Normal Form Analysis
Blackman and Vigna [1] introduce scramblers: nonlinear functions applied to the state of the linear engine before output, designed to mask linear artifacts while preserving the period and equidistribution properties of the underlying engine. They analyze four scramblers using the algebraic normal form (ANF) representation of Boolean functions and the splitting field of the characteristic polynomial.
2.3.1 Representation by Generating Functions and Splitting Field
Every bit of a linear engine satisfies a linear recurrence with characteristic polynomial p(x). Blackman and Vigna [1] represent each bit's output sequence as a rational generating function f(z). In the splitting field E = (Z/2Z)[α]/p(α), where α, α², α⁴, …, α²ⁿ⁻¹ are all roots of p, the j-th output has the explicit form:
bⱼ = Σᶜ βᶜ(αᶜ)ⱼ, summed over the cyclotomic coset C = {1, 2, 4, …, 2ⁿ⁻¹}, with βᶜ ∈ E
This representation makes it possible to compute the linear complexity of any Boolean function of the state bits [1]. Blackman and Vigna use this to establish that the linear complexity of scrambler output bits grows rapidly with bit index b, via the standard bound on linear complexity of a filtered LFSR (Proposition 9.2 of [1]):
U(n, d) = Σᵈⱼ₌₁ C(n, j)
where d is the polynomial degree of the Boolean function in ANF and n is the state size.
2.3.2 The + Scrambler
The + scrambler adds two words of the state array in Z/2ʷZ. The closed-form ANF of bit b is in ANF with exactly 2ᵇ + 1 monomials of maximum degree b + 1 [1]. Blackman and Vigna note that combining bits does not lose polynomial degree: any linear combination of + scrambler output bits has linear complexity equal to that of the bit of highest degree.
Bit | xoroshiro128+ (n=128) | xoshiro256+ (n=256) | xoshiro512+ (n=512) |
0 (lowest) | 8,256 | 32,896 | 131,328 |
1 | 349,632 | 2,796,416 | 22,370,048 |
2 | 11,017,632 | 177,589,056 | 2,852,247,168 |
3 | 275,584,032 | 8,987,138,112 | 290,367,762,560 |
4 | ~6.9 × 10⁹ | ~2.2 × 10¹¹ | ~7.3 × 10¹² |
Source: Table 11 of Blackman and Vigna [1].
Blackman and Vigna classify + as a weak scrambler: the lowest output bit retains detectable linearity, though higher bits become undetectable in any feasible test [1].
2.3.3 The * Scrambler
The * scrambler multiplies a chosen state word by an odd constant, recommended close to φ·2ʷ (φ = the golden ratio), following Knuth [1]. Blackman and Vigna prove (Theorem 9.4 of [1]) that for multiplier 3x, the number of monomials of bit b equals the count of binary palindromes smaller than 2ᵇ. Blackman and Vigna classify * as a weak scrambler: it eliminates Hamming-weight dependencies better than + due to richer monomial structure, but retains linear artifacts in MatrixRank when output is bit-reversed [1].
2.3.4 The ++ Scrambler
The ++ scrambler takes two words x, y from the state and computes z = ρᵣ(x + y) + x, where ρᵣ is left rotation by r bits. Blackman and Vigna [1] recommend r ∈ {17, 23, 41, 47} for w = 64. By Proposition 9.3 of [1], the ++ scrambler achieves polynomial degree w − r + 1 in the lowest output bit.
Bit | xoroshiro128++ (n=128) | xoshiro256++ (n=256) | xoshiro512++ (n=512) | xoroshiro1024++ (n=1024) |
0 | ~10³⁶ | ~3×10⁴⁸ | ~10⁶⁸ | ~9×10⁷⁴ |
1 | ~2×10³⁶ | ~2×10⁴⁹ | ~10⁶⁹ | ~2×10⁷⁶ |
2 | ~3×10³⁶ | ~10⁵⁰ | ~9×10⁶⁹ | ~4×10⁷⁷ |
3 | ~4×10³⁶ | ~4×10⁵⁰ | ~8×10⁷⁰ | ~10⁷⁸ |
Source: Table 12 of Blackman and Vigna [1].
Blackman and Vigna classify ++ as a strong scrambler; xoshiro256++ and xoroshiro128++ pass all BigCrush and HWD tests [1,2].
2.3.5 The ** Scrambler
The ** scrambler applies multiply-rotate-multiply: z = ρᵣ(x·S)·T. Blackman and Vigna [1] identify optimal parameters for w = 64 as S = 5, T = 9, r = 7, giving polynomial degree 56 for the lowest bit.
Bit | xoroshiro128** (n=128) | xoshiro256** (n=256) | xoshiro512** (n=512) | xoroshiro1024** (n=1024) |
0 | ~3×10³⁷ | ~2×10⁵⁷ | ~4×10⁷⁵ | ~10⁹³ |
1 | ~4×10³⁷ | ~7×10⁵⁷ | ~3×10⁷⁶ | ~2×10⁹⁴ |
2 | ~6×10³⁷ | ~3×10⁵⁸ | ~2×10⁷⁷ | ~3×10⁹⁵ |
3 | ~7×10³⁷ | ~9×10⁵⁸ | ~2×10⁷⁸ | ~6×10⁹⁶ |
Source: Table 13 of Blackman and Vigna [1].
Blackman and Vigna classify ** as a strong scrambler; xoshiro256** and xoroshiro128** pass all BigCrush and HWD tests [1,2].
2.4 The xoshiro and xoroshiro Linear Engines
2.4.1 xoroshiro
Updates two words of the state cyclically. Blackman and Vigna [1] apply Sothanaphan's theorem on non-commuting block matrices [3] to compute the column-determinant word polynomial, since the transition matrix blocks do not commute under standard determinant rules.
2.4.2 xoshiro
Updates all state words simultaneously. Parameters a = 17, b = 45 yield a primitive characteristic polynomial of weight 115 for xoshiro256, close to the Compagner optimum of n/2 = 128 [1,4]. The reference C implementation (Figure 4 of [1]):
const uint64_t result_plus = s[0] + s[3];
const uint64_t result_plusplus = rotl(s[0] + s[3], 23) + s[0];
const uint64_t result_starstar = rotl(s[1] * 5, 7) * 9;
const uint64_t t = s[1] << 17;
s[2] ^= s[0]; s[3] ^= s[1]; s[1] ^= s[2]; s[0] ^= s[3];
s[2] ^= t; s[3] = rotl(s[3], 45);
2.5 Equidistribution
A generator is d-dimensionally equidistributed if every d-tuple of consecutive outputs appears with the theoretically expected frequency over the full period [1]. Blackman and Vigna prove several equidistribution results for + and ++ scrambled generators (Propositions 7.1–7.6 of [1]); generators using * or ** are k-dimensionally equidistributed (the maximum possible), since these scramblers remap sequences bijectively. xoshiro256** is therefore 4-dimensionally equidistributed; xoshiro256++ is 3-dimensionally equidistributed [1].
2.6 Escaping Zeroland
When a linear engine is initialized with low Hamming weight, it requires many iterations before the distribution of ones approaches 0.5. Panneton, L'Ecuyer, and Matsumoto [10] document this "zeroland" problem; for MT19937 the warm-up period is approximately one million iterations. Blackman and Vigna [1] show xoshiro256++ and xoroshiro128++ reach equilibrium in fewer than 100 iterations from any non-zero state seeded via SplitMix64. This warm-up behavior can be observed directly by deliberately initializing a generator's state with a single set bit and watching the bias gradually diffuse across the state over many iterations before settling into balanced output [13].
2.7 Jump Functions
Because xoshiro and xoroshiro are F₂-linear, their state evolution is representable as matrix exponentiation over F₂. A jump function advances the state by J steps in O(n) time without iterating individually [1,11]. Blackman and Vigna provide jump polynomials for xoshiro256 for skips of 2¹²⁸ and 2¹⁹² steps [1]. The Mersenne Twister provides no jump function.
3. The Mersenne Twister: Documented Properties and Limitations
3.1 Overview
MT19937 [5] has period 2¹⁹⁹³⁷ − 1 and is 623-dimensionally equidistributed. It passed DIEHARD at introduction and remains the default PRNG in Python, R, and MATLAB. Vigna [12] has published a detailed critique of the family.
3.2 Statistical Test Failures
MT19937-64 has two systematic BigCrush failures [6]: MatrixRank and LinearComp, in both forward and reversed output — a direct consequence of its F₂-linear structure and sparse transition matrix [1].
3.3 Zeroland
Panneton et al. [10] document that MT19937 requires ~1 million iterations to escape zeroland; the first ~10⁶ outputs from a freshly seeded instance are non-uniform.
3.4 Speed
Generator | ns / 64-bit output | cycles / byte |
MT19937-64 | 2.19 | 0.94 |
xoshiro256++ | 0.86 | 0.34 |
xoshiro256** | 0.84 | 0.33 |
xoroshiro128++ | 0.90 | 0.38 |
xoroshiro128+ | 0.80 | 0.36 |
SplitMix64 | 1.14 | 0.46 |
Source: Table 1 of Blackman and Vigna [1], Intel Core i7-8700B @ 3.20GHz.
4. Statistical Validation: BigCrush and the HWD Test
4.1 TestU01 and BigCrush
TestU01 [6] is a C library comprising 106 statistical tests. Blackman and Vigna [1] define systematic failure as a p-value outside [0.001, 0.999] across all tested seeds. Passing BigCrush is necessary but not sufficient for generator quality, as the HWD test demonstrates [2].
4.2 The Hamming-Weight Dependency Test
Blackman and Vigna [2] introduce HWD to detect bias in Hamming-weight distributions of consecutive output words, a bias class BigCrush misses in several generators, since F₂-linear generators with sparse matrices tend to map low-Hamming-weight states to other low-Hamming-weight states.
4.2.1 Construction
For window size k and word width w, each word is classified by Hamming weight relative to a threshold, producing a base-3 signature over k outputs. The test tracks the average Hamming weight of the following word per signature, normalizes to standard normal variables, and applies a Kronecker-power transform of a small unitary matrix [2]. The full test (w=64, k=8) processes one terabyte of output in ~3 minutes on reference hardware [2].
4.3 HWD Test Results
Generator | BigCrush | HWD failure at | Faulty signature |
xorshift128 | Fails MR, LC | 8 × 10⁸ bytes | 00000021 |
xorshift128+ | Passes all | 6 × 10⁹ bytes | 00000012 (transitional) |
xoroshiro128+ | Passes all | 5 × 10¹² bytes | 00000012 |
dSFMT (521 bits) | Passes all | 6 × 10¹² bytes | 1001000100100010 |
Tiny MT (127 bits) | Passes all | 8 × 10¹³ bytes | 00000202 (w=32) |
MT (521 bits) | Fails MR, LC | 4 × 10¹⁰ bytes | 1000000100000000 |
MT (607 bits) | Fails MR, LC | 4×10⁸–4×10¹⁰ bytes | 1000000001000000000 |
SFMT (607 bits) | Fails MR, LC | 4 × 10⁸ bytes | 001000001000 |
WELL512 | Passes all | 3 × 10¹⁵ bytes | 2001002200000000 |
xoshiro256++ | Passes all | None to 10¹⁵ bytes | — |
xoshiro256** | Passes all | None to 10¹⁵ bytes | — |
Source: Table 3 of Blackman and Vigna [2].
Blackman and Vigna [2] note this is, to their knowledge, the first demonstration of bias in Tiny Mersenne Twister, dSFMT (521 bits), and WELL512 beyond linearity tests — establishing that BigCrush passage alone is insufficient evidence of quality. xorshift128+, the PRNG embedded in the JavaScript engines of Chrome, Firefox, Safari, and Edge [1], exhibits an HWD failure at approximately 6 gigabytes of output.
5. The Seeding Problem
5.1 Correlation on Similar Seeds
F₂-linear generators must not be initialized by passing the seed directly into the state array: similar seeds produce correlated states, since the linear transition function maps nearby states to nearby successors until differences propagate through the full state [1].
5.2 SplitMix64
The recommended initializer, from the LXM family [11], uses multiplication in Z/2⁶⁴Z rather than F₂-linear operations, so single-bit seed differences produce uncorrelated outputs [1]:
x += 0x9E3779B97F4A7C15;
x ^= x >> 30; x *= 0xBF58476D1CE4E5B9;
x ^= x >> 27; x *= 0x94D049BB133111EB;
output = x ^ (x >> 31);
To seed a 256-bit xoshiro256 state, call SplitMix64 four times, filling one word each [1].
5.3 The Reseeding Error
A common mistake is reseeding a new generator instance per batch rather than seeding once and using the transition function thereafter; system-time-based reseeding at high frequency produces correlated sequences. This exact failure mode is easy to demonstrate directly: repeatedly resetting a generator to the same seed reproduces the same output sequence every time, since the seed is the only true "input" the generator has once running [13].
5.4 Reproducibility
For audit purposes, results must be bit-identical on re-run given the same seed and implementation version, requiring stored seeds and version-controlled implementations.
6. Parallel Simulation and Non-Overlapping Sequences
6.1 The Overlap Probability
For period P and n threads each consuming L outputs, Pr[overlap] ≤ n²L/P [1]. For xoshiro256 (P = 2²⁵⁶ − 1) this is negligible at any practical scale; the binding constraint is seed construction, since sequentially-seeded parallel instances remain subject to the correlation-on-similar-seeds problem.
6.2 Jump-Based Parallelization
Seed a master instance via SplitMix64; for thread i, apply i jump operations of size J = 2¹²⁸ from the master state. This guarantees non-overlapping subsequences by construction [1], with negligible per-thread overhead.
7. Floating-Point Conversion
7.1 The Standard Conversion
Converting a 64-bit integer x to a double u ∈ [0,1) requires the top 53 bits [1]:
u = (x >> 11) * 0x1.0p-53;
This produces a distribution uniform over all 2⁵³ dyadic rationals in [0,1), and prefers the high bits of x — appropriate for + and ++ scramblers, whose highest bits carry the greatest linear complexity.
7.2 The Half-Values Problem
An alternative bit-manipulation ("union trick") conversion generates only 52 significant bits, producing half the representable doubles in [0,1); Blackman and Vigna note dSFMT shares this limitation. The multiplicative conversion above avoids it.
8. Comparison of Current Generators
8.1 Generator Properties
Generator | Period | State | Speed (ns/64b) | BigCrush | HWD | Jump fn |
MT19937-64 | 2¹⁹⁹³⁷−1 | 19937 bits | 2.19 | Fails MR, LC | Not tested | No |
xorshift128+ | 2¹²⁸−1 | 128 bits | 0.72 | Passes | 6×10⁹ B | No |
xoroshiro128++ | 2¹²⁸−1 | 128 bits | 0.90 | Passes | Clean | Yes |
xoshiro256++ | 2²⁵⁶−1 | 256 bits | 0.86 | Passes | Clean | Yes |
xoshiro256** | 2²⁵⁶−1 | 256 bits | 0.84 | Passes | Clean | Yes |
PCG-DXSM | 2¹²⁸ | 128 bits | ~0.90 | Passes | Not formally tested | Partial |
8.2 Adoption
xoshiro256++: Rust SmallRng, Julia (default), GNU Fortran, Lua, Facebook Folly. xoshiro256**: .NET 6+ default. xorshift128+: Chrome, Firefox, Safari, Edge JS engines [1]. PCG-DXSM: NumPy default (v1.17+), Go math/rand/v2 [7]. MT19937: Python, R, MATLAB defaults.
8.3 The xoshiro vs PCG Debate
Vigna has argued the LCG foundation of PCG introduces directional bias; O'Neill has argued xoshiro's F₂-linear structure produces state-space clustering near all-zero states. Both families pass all currently published empirical tests at feasible scale; xoshiro256++ has the advantage of formal HWD testing [2], while PCG-DXSM has extensive deployment as NumPy's default [7].
9. Implementation Recommendations
Generator selection: use xoshiro256++ or xoshiro256**; avoid MT19937, xorshift128+, and LCG variants [1,12]
Seeding: seed once from a high-entropy OS source, initialize via SplitMix64, never pass the seed directly into the state array [1]
Parallel use: use jump functions for non-overlapping subsequences; never seed threads with sequential integers [1]
Floating-point conversion: use (x >> 11) * 0x1.0p-53 [1]
Store the seed with every result set for reproducibility
For non-xoshiro generators, discard the first max(1000, state-size-in-bits) outputs to escape zeroland [10]
Validate with TestU01 SmallCrush routinely, BigCrush before high-stakes use [6]
Version-control the generator implementation
10. Conclusion
The research literature on pseudorandom number generation has advanced substantially since the Mersenne Twister. Blackman and Vigna [1,2] establish, through formal proofs and empirical testing, that the xoshiro256 family achieves faster generation, stronger statistical properties, formal equidistribution guarantees, and correct parallel-use support via jump functions, relative to MT19937.
The HWD test [2] demonstrates that BigCrush passage alone is not sufficient evidence of quality: dSFMT, WELL512, and Tiny Mersenne Twister pass BigCrush while exhibiting detectable Hamming-weight dependencies. xoshiro256++ and xoshiro256** pass both tests to petabyte scale with no known failures [1,2]. The seeding, parallel-use, and floating-point conversion problems are each independent, well-documented sources of error requiring separate attention alongside generator selection.
References
Blackman, D. and Vigna, S. (2021). Scrambled Linear Pseudorandom Number Generators. ACM Transactions on Mathematical Software, 47(4), Article 36.
Blackman, D. and Vigna, S. (2022). A New Test for Hamming-Weight Dependencies. arXiv:2108.13061v2.
Sothanaphan, N. (2017). Determinants of block matrices with noncommuting blocks. Linear Algebra and its Applications, 512, 202–218.
Compagner, A. (1991). The hierarchy of correlations in random binary sequences. Journal of Statistical Physics, 63(5-6), 883–896.
Matsumoto, M. and Nishimura, T. (1998). Mersenne Twister. ACM Transactions on Modeling and Computer Simulation, 8(1), 3–30.
L'Ecuyer, P. and Simard, R. (2007). TestU01. ACM Transactions on Mathematical Software, 33(4), Article 22.
O'Neill, M.E. (2014). PCG. Harvey Mudd College Technical Report HMC-CS-2014-0905.
Turan, M.S., Barker, E., Kelsey, J., McKay, K.A., Baish, M.L. and Boyle, M. (2018). NIST Special Publication 800-90B.
L'Ecuyer, P. (1996). Maximally equidistributed combined Tausworthe generators. Mathematics of Computation, 64(213), 203–213.
Panneton, F., L'Ecuyer, P. and Matsumoto, M. (2006). Improved long-period generators based on linear recurrences modulo 2. ACM Transactions on Mathematical Software, 32(1), 1–16.
Steele, G.L., Lea, D. and Flood, C.H. (2014). Fast splittable pseudorandom number generators. Proceedings of OOPSLA 2014, 453–472.
Vigna, S. (2019). It is high time we let go of the Mersenne Twister. arXiv:1910.06437.
Ackley, D. (2013). NMCS4ALL: Random number generators. New Mexico CS for All. https://www.youtube.com/watch?v=_tN2ev3hO14