Review: Generators, uniform floating-point draws, and empirical testing¶
Citations
- Blackman, D., & Vigna, S. (2021). Scrambled linear pseudorandom number generators. ACM Transactions on Mathematical Software, 47(4), 1-32. DOI
- Goualard, F. (2022). Drawing random floating-point numbers from an interval. ACM Transactions on Modeling and Computer Simulation, 32(3), 1-24. DOI
- L'Ecuyer, P., & Simard, R. (2007). TestU01: A C library for empirical testing of random number generators. ACM Transactions on Mathematical Software, 33(4), 1-40. DOI
Abstracts¶
Blackman & Vigna (2021). \( \mathbb{F}_2 \)-linear pseudorandom number generators are very popular due to their high speed, to the ease with which generators with a sizable state space can be created, and to their provable theoretical properties. However, they suffer from linear artifacts that show as failures in linearity-related statistical tests such as the binary-rank and the linear-complexity test. In this article, we give two new contributions. First, we introduce two new \( \mathbb{F}_2 \)-linear transformations that have been handcrafted to have good statistical properties and at the same time to be programmable very efficiently on superscalar processors, or even directly in hardware. Then, we describe some scramblers, that is, nonlinear functions applied to the state array that reduce or delete the linear artifacts, and propose combinations of linear transformations and scramblers that give extremely fast pseudorandom number generators of high quality. A novelty in our approach is that we use ideas from the theory of filtered linear-feedback shift registers to prove some properties of our scramblers, rather than relying purely on heuristics. In the end, we provide simple, extremely fast generators that use a few hundred bits of memory, have provable properties, and pass strong statistical tests.
Goualard (2022). Drawing a floating-point number uniformly at random from an interval \( [a, b) \) is usually performed by a location-scale transformation of some floating-point number drawn uniformly from \( [0, 1) \). Due to the weak properties of floating-point arithmetic, such a transformation cannot ensure respect of the bounds, uniformity or spatial equidistributivity. We investigate and quantify precisely these shortcomings while reviewing the actual implementations of the method in major programming languages and libraries, and we propose a simple algorithm to avoid these shortcomings without compromising performances.
L'Ecuyer & Simard (2007). We introduce TestU01, a software library implemented in the ANSI C language, and offering a collection of utilities for the empirical statistical testing of uniform random number generators (RNGs). It provides general implementations of the classical statistical tests for RNGs, as well as several others tests proposed in the literature, and some original ones. Predefined tests suites for sequences of uniform random numbers over the interval (0, 1) and for bit sequences are available. Tools are also offered to perform systematic studies of the interaction between a specific test and the structure of the point sets produced by a given family of RNGs. [...] Besides introducing TestU01, the article provides a survey and a classification of statistical tests for RNGs. It also applies batteries of tests to a long list of widely used RNGs.
Hifuku runs the same survey on two engines. The host engine draws from
numpy.random.Generator, and the device engine draws from xoroshiro128+ in
numba.cuda.random. Three questions follow from that arrangement, and these
three papers answer one each: what the device generator is and which of its bits
are sound (Blackman & Vigna), how many distinct floats a draw should produce and
how to produce them (Goualard), and how to decide whether a generator matches the
distribution it claims (L'Ecuyer & Simard).
Linear engines and scramblers¶
Blackman & Vigna start from a tension in \( \mathbb{F}_2 \)-linear generators.
Linearity is what makes the useful properties provable: full period, large state
in few operations, and a connection to linear-feedback shift registers that
carries an existing theory. Linearity is also detectable. The binary-rank test
and the linear-complexity test "are failed by all linear generators," and
TestU01 implements them under the names MatrixRank and LinearComp.
A scrambler is a nonlinear mapping from the state of the linear engine to the
output word. Section 4 describes four of them, and divides them into two classes.
The sum of two words of state (+, section 4.1) and multiplication of one word
by an odd constant (*, section 4.2) are weak. Sum-rotate-sum (++, section
4.3) and multiply-rotate-multiply (**, section 4.4) are strong, because the
rotation carries the high bits of the first operation down into the low
positions.
The weakness of + follows from the algebra of the operation:
the lowest bit output by the + scrambler is just a xor of bits following the same linear recurrence, and thus follows, in turn, the same linear recurrence. [...] As we consider higher bits, there is still a linear recurrence describing their behavior, but it becomes quickly of such a high linear complexity to become undetectable.
The introduction states the same point in terms of the output:
while empirical tests usually do not show linear artifacts anymore, the lower bits are unchanged or just slightly modified by such operations. Thus, those bits in isolation (or combined with a sufficiently small number of good bits) will fail linearity tests.
The remedy the authors give is to take the high bits. Discussing xoshiro256+,
they write that "its lowest bits have low linear complexity, but since one needs
just the upper 53 bits, the resulting floating-point values have no linear bias."
Section 5.3 repeats the advice for the whole + family, including xoroshiro128+,
which the paper offers as the choice "if space is an issue."
Table 11 puts numbers on which bits are weak. For xoroshiro128+ the estimated linear complexity of the five lowest bits is 128, 8256, 349632, 11017632, and 275584032, a factor of 25 or more at every step. The text reads the table this way: "while an accurate linear-complexity test might catch the fourth lowest bit of xoroshiro128+, the degree raises quickly to the point the linearity is undetectable." The weakness is therefore concentrated in roughly the low four bits.
The paper also records one residual defect. Under the authors' own Hamming-weight dependency test, xoroshiro128+ fails, but "bias can be detected only after 5 TB of data, which makes it unlikely to affect applications in any way."
Seeding and separating streams¶
Two points bear on running one generator across many parallel walkers.
Seeding must not use the generator being seeded. The authors recommend SplitMix, "as research has shown that initialization must be performed with a generator radically different in nature from the one initialized to avoid correlation on similar seeds."
Separation is available by construction. All the generators in the paper "have jump functions that make it possible to move ahead quickly by any number of next-state steps." Where random starting points give only a probabilistic guarantee of non-overlap, "one can also use jumping to guarantee the absence of overlap."
The testing protocol¶
The paper's own empirical work sets a threshold worth reusing. The authors run BigCrush from several seeds and "consider a test failed if its p-value is outside of the interval [0.001..0.999]." A failure that appears for every seed is called systematic, and only systematic failures are reported.
How many floats a draw should produce¶
Goualard's section 2.2 gives the counting argument for the unit interval. Let \( p \) be the width of the significand. Then
is the largest set of equally spaced floats in \( [0, 1) \), and the paper states the consequence directly:
considering a set with less than \( 2^p \) values is wasteful (we could draw more); however, allowing more than \( 2^p \) values destroys uniformity.
Of the two procedures the survey finds in real libraries, the first is to "generate a random integer in the interval \( [0, 2^k - 1] \) and divide it by \( 2^k \)," and for that procedure "the best choice for \( k \) is \( p \)." The second procedure builds a float in \( [1, 2) \) by writing a random fraction into the significand and subtracting one; it reaches only \( 2^{p-1} \) values, half of what the format allows.
The rest of the paper concerns the location-scale transformation \( a + (b - a) x \) used to move a draw from \( [0, 1) \) to an arbitrary interval, and it quantifies three separate failures of that step: the right bound can be returned, the reachable values are not equally spaced, and they are not equally likely. One measurement gives the flavor. Mapping the \( 2^{24} \) values of \( S_0^1 \) into \( [0.25, 1) \) in binary32 reaches 13,981,013 of the 16,777,216 floats available there, about 83 percent, with a non-zero variance in the spacing between reachable values.
Testing a generator against its claimed distribution¶
L'Ecuyer & Simard state the hypothesis under test. For a generator of real numbers, \( \mathcal{H}_0^A \) is that the successive outputs \( u_0, u_1, \ldots \) "are independent random variables from the uniform distribution over the interval (0, 1), that is, i.i.d. \( U(0,1) \)." For a bit generator, \( \mathcal{H}_0^B \) is that the output bits are independent and equally likely. The two are linked: any prespecified sequence of bits taken from the \( u_i \) must be independent under \( \mathcal{H}_0^A \), so bit tests test the real-valued hypothesis indirectly.
Neither hypothesis can hold exactly for an algorithmic generator. The vectors \( (u_0, \ldots, u_{t-1}) \) take values only in the finite set \( \Psi_t \) of \( t \)-dimensional vectors the generator can produce, whose size is bounded by the number of admissible seeds. The bit hypothesis fails for the same reason as soon as the sequence length exceeds the number of bits in the state. The design goal is therefore that \( \Psi_t \) be evenly distributed over the unit cube, and the empirical goal is more modest than proof:
no universal test or battery of tests can guarantee, when passed, that a given generator is fully reliable for all kinds of simulations. But even if statistical tests can never prove that an RNG is foolproof, they can certainly improve our confidence in it. [...] The difference between the good and bad RNGs [...] is that the bad ones fail very simple tests whereas the good ones fail only very complicated tests that are hard to figure out or impractical to run.
The paper argues against fixing a rejection region in advance. The textbook procedure suits a fixed, usually small sample; when testing a generator "the sample sizes are huge and can usually be increased at will," so the authors report the p-value
and read it on a scale. A p-value below about \( 10^{-10} \) is a clear failure. A p-value that is not close to 0 or 1 means the test detected nothing. A p-value in between, 0.002 for example, "is suspicious but does not clearly indicate rejection," and the response is to replicate the test on disjoint output from the same generator "until either failure becomes obvious or suspicion disappears." Multiplicity is the reason for the caution: across many tests, p-values below 0.01 or above 0.99 "are often obtained by chance even if the RNG behaves correctly with respect to these tests (such values should normally appear approximately 2% of the time)."
The Kolmogorov-Smirnov statistics appear in the paper as equations 1 to 3, used with Anderson-Darling in the two-level procedure that compares \( N \) replicated first-order p-values against the uniform distribution:
Section 6.3 defines the predefined batteries. SmallCrush runs in seconds and is meant "to detect gross defects in generators or errors in their implementation." Crush consumes about \( 2^{35} \) numbers over 96 tests and 144 statistics. BigCrush consumes about \( 2^{38} \) numbers over 106 tests and 160 statistics. The authors make no claim that the selected tests are independent or exhaustive.
Relevance to Hifuku¶
The conversion in _uniform_float¶
chain_tree.py converts a draw to a float32 with
where _UNIT_FLOAT32_SCALE is \( 2^{-24} \). The two papers each fix one half of
that expression. Goualard fixes the count: binary32 has \( p = 24 \), so the
output set is \( \{k \cdot 2^{-24} \mid 0 \leq k < 2^{24}\} \), which is the
largest equally spaced set the format supports. Blackman & Vigna fix which bits:
the + scrambler leaves the low bits weak, so the 24 kept bits are the top ones,
40 bits above the low four that Table 11 identifies.
The host side follows the same rule at its own precision. numpy.random.Generator
produces float64 draws that are multiples of \( 2^{-53} \) with a maximum below
one, which is \( k = p = 53 \), the choice Goualard recommends.
The library conversion xoroshiro128p_uniform_float32 instead keeps 53 bits and
narrows the result to float32. That admits \( 2^{53} \) values where the format
supports \( 2^{24} \), the case Goualard identifies, and the arithmetic
consequence is a return of exactly 1.0 for every draw at or above
\( 2^{64} - 2^{39} \), a fraction \( 2^{-25} \) of the range, or about one draw
in \( 3.4 \times 10^7 \). All four claims in this paragraph and the one before it
are checked in
goualard2022drawing-unit-float32.py.
Separating the walkers¶
create_xoroshiro128p_states seeds state 0 with SplitMix64 and gives each later
state the previous one advanced \( 2^{64} \) steps by the jump function. Both
halves of that are the paper's recommendations: SplitMix for initialization
because it is unlike the generator it seeds, and jumping rather than random
starting points because jumping guarantees the absence of overlap instead of
making it improbable. The separation between walkers is therefore a property of
the seeding and not of luck, which matters because walkers that shared draws
would share trajectories.
A parameter difference worth recording¶
The numba implementation uses the rotation and shift amounts (55, 14, 36). Table 2 of Blackman & Vigna gives (24, 16, 37) for the xoroshiro128 engine. The paper's tabulated test results are for its own parameters, so they describe the numba generator only as far as the two parameter sets behave alike, which the paper does not address. This is a reason to measure the library's output rather than to infer its quality from the paper.
What the test suite asserts¶
tests/test_random_draws.py tests \( \mathcal{H}_0^A \) for the device uniform,
tests it for the device normal against the standard normal, and compares each
against the host engine with a two-sample Kolmogorov-Smirnov test. The threshold
is \( 10^{-3} \), the same bound Blackman & Vigna use in their BigCrush protocol.
Two limits of that suite follow from L'Ecuyer & Simard, and both are deliberate.
The assertion is one-sided, pvalue > alpha, where the published protocol
rejects outside \( [0.001, 0.999] \). A sample that is uniform to an implausible
degree passes. The upper tail detects excessive uniformity, which is a property
of the generators themselves rather than of the conversion arithmetic these tests
are aimed at.
The seeds are fixed, so a suspicious p-value is reproducible rather than random, and the remedy the paper prescribes is replication on disjoint output. A single low p-value from these tests is a reason to rerun with different seeds before treating it as a defect, in the same way that Blackman & Vigna count only systematic failures.
Neither engine is under test as a generator. Two hundred thousand draws against Crush's \( 2^{35} \) is a coarse instrument, and the generators have already been tested at that scale in the literature. What these tests cover is the part that is local to Hifuku: the conversion from a draw to a float, the seeding of the walkers, and the claim that the two engines sample the same distributions. The normal draws make the last point plainly, because the device uses Box-Muller over two device uniforms and the host uses numpy's ziggurat. The algorithms differ, the streams differ, and the distribution is what has to agree.
The Hamming-weight caveat in scale¶
Blackman & Vigna detect Hamming-weight dependency in xoroshiro128+ after 5 TB of output. A survey draws a handful of values per variation, so a run of a few hundred thousand variations consumes megabytes of 64-bit output. The detection threshold is five to six orders of magnitude above any survey budget the code accepts.
In-document navigation¶
| Paper | Section | Content |
|---|---|---|
| Blackman & Vigna | 1, Introduction | Why linear generators fail MatrixRank and LinearComp |
| Blackman & Vigna | 4, Scramblers | The *, **, +, and ++ scramblers, and the weak low bits of + |
| Blackman & Vigna | 5, Table 1 | Systematic BigCrush failures and the Hamming-weight results |
| Blackman & Vigna | 5.3, Choosing a Generator | xoroshiro128+ for floating-point output; SplitMix seeding |
| Blackman & Vigna | 7, Equidistribution | Equidistribution of +-scrambled generators |
| Blackman & Vigna | 9, Table 11 | Linear complexity of the five lowest bits |
| Goualard | 2.2 | The set \( S_0^1 \), the \( 2^p \) counting argument, and the two procedures |
| Goualard | 3 | Survey of 15 languages and libraries |
| Goualard | 4 | Error analysis of the location-scale transformation |
| Goualard | 5 | The \( \gamma \)-section algorithm |
| L'Ecuyer & Simard | 2 | \( \mathcal{H}_0^A \), \( \mathcal{H}_0^B \), and the point set \( \Psi_t \) |
| L'Ecuyer & Simard | 3 | p-values rather than rejection regions; the two-level procedure |
| L'Ecuyer & Simard | 4, 5 | The tests, classified by what they detect |
| L'Ecuyer & Simard | 6.3 | SmallCrush, Crush, BigCrush, and the bit-sequence batteries |
| L'Ecuyer & Simard | 7, Table I | Battery results for widely used generators |