A New Test for Hamming-Weight Dependencies

1 INTRODUCTION

Pseudorandom number generators (PRNGs) are algorithms that generate a seemingly random output using a deterministic algorithm. A w-bit PRNG is defined by a state space S, a transition (or next-state) computable function \( \tau :S\rightarrow S \), and a computable output function \( \varphi :S\rightarrow \lbrace \,0,1\,\rbrace ^w \) that maps the state space into w-bit words. One then considers an initial state, or seed \( \sigma \in S \), and computes the sequence of w-bit outputs \( \begin{equation*} \varphi (\sigma), \varphi (\tau (\sigma)), \varphi \bigl (\tau ^2(\sigma)\bigr), \varphi \bigl (\tau ^3(\sigma)\bigr), ... \end{equation*} \) The outputs can be used to generate reals in the unit interval—for example, multiplying them by \( 2^{-w} \). Knuth [8] discusses PRNGs at length.¹

A classic example is given by multiplicative congruential generators, which are defined by a prime modulus \( \mu \) and a multiplier \( \alpha \). Then, \( S=\mathbf {Z}/\mu \mathbf {Z} \), \( \tau :x\mapsto \alpha x \), and the output function is given by the binary representation of x (one tries to choose \( \mu \) close to \( 2^w \)). Another well-known example is the class of \( \mathbf {F}_2 \)-linear generators [10], in which S is a vector of \( \mathbf {F}_2^n \) (i.e., n bits) and \( \tau \) is an \( \mathbf {F}_2 \)-linear transformation on S; however, usually, the transformation can be expressed by basic \( \mathbf {F}_2 \)-linear operations on words, such as rotations, shift, and XORs, rather than in matrix form. The output function might pick a word of w bits from the state: for example, \( n=kw \) for some k and then the state can be represented by k w-bit words, and the output function can just choose one of those. In some generators, moreover, the output function is not \( \mathbf {F}_2 \)-linear.

Several theoretical properties help in the design of PRNGs: however, once designed, a PRNG is submitted to a set of statistical tests, which try to discover some statistical bias in the output of the generator. The tests compute, using the output of the generator, statistics whose distribution is known (at least approximately) under the assumption that the output of the generator is random. Then, by applying the (complementary) cumulative distribution function to the statistics, one obtains a p-value, which should be neither too close to zero nor too close to one (see Knuth [8] for a complete introduction to the statistical testing of PRNGs).

The Hamming weight of a w-bit word x is the number of ones in its binary representation. Tests for Hamming-weight dependencies try to discover some statistical bias in the Hamming weight of the output of the generator. In particular, such tests do not depend on the numerical values of the outputs: indeed, sorting the bits of each examined block (e.g., first all zeroes and then all ones) would not modify the results of the test (albeit the values would now be very small).

Since the number of ones in a random w-bit word has a binomial distribution with w trials and probability of success \( 1/2 \), in the most trivial instance one examines m consecutive outputs \( {x}_0 \), \( {x}_1 \), \( \dots \, \), \( {x}_{m-1} \) and checks that the average of their Hamming weights has the correct distribution (which will be quickly approximated very well by a normal distribution as m grows [5]). Tests may also try to detect dependencies: for example, one can consider (overlapping) pairs of consecutive outputs, and check that the associated pairs of Hamming weights have the expected distribution [11]. Matsumoto and Nishimura [17] have introduced a theoretical figure of merit that can predict after how many samples a \( \mathbf {F}_2 \)-linear generator will fail a specific Hamming-weight test. The NIST statistical test suite for PRNGs [20] contains tests based on Hamming-weight dependencies as well.

In this article, we introduce a new test for Hamming-weight dependencies that improves significantly over the state of the art. We find bias in some old and some new generators for which tests of this type from TestU01 [12], a well-known framework for statistical testing of PRNGs, were unable to find bias even using a large amount of output.

All code used in this article is available under the GNU General Public License.² Code for reproducing the results of this work has been permanently stored on the Zenodo platform.³

2 MOTIVATION

It is known since the early days of \( \mathbf {F}_2 \)-linear generators that sparse transition matrices induce some form of dependency on the Hamming weight of the state. Since the output is computed starting from the state, these dependencies might induce Hamming-weight dependencies in the output as well. For example, if the state has very low Hamming weight—that is, very few ones—one might need a few iterations (or more than a million, in the case of the Mersenne Twister with 19,937 bits of state [16]) before the state contains ones and zeroes approximately in the same amount. However, this is a minor problem because for generators with, say, at least 128 bits of state, the probability of passing through such states is negligible.

Yet what we witness very clearly in the case of almost-all-zeros state might be true in general: states with few ones might lead to states with few ones, or due to XOR operations, states with many ones might lead to states with few ones. This kind of dependency is more difficult to detect.

Here we consider as motivation a few generators: xorshift128+ [25] is the stock generator of most Javascript implementations in common browsers; the SFMT (SIMD-Friendly Mersenne Twister) [21] is a recent improvement on the classic Mersenne Twister using SIMD instructions, and we will use the version with 607 bits of state; the dSFMT [22] is another version of the Mersenne Twister that natively generate doubles, and we will use the version with 521 bits of state; WELL is a family of generators with excellent equidistribution properties [19], and we will use the version with 512 bits of state; and finally, we consider a new \( \mathbf {F}_2 \)-linear transformation, xoroshiro, designed by the authors, and the associated generators xoroshiro128+ and xoroshiro1024+ [1].⁴ All of these generators have quite sparse transition matrices (WELL512 having the densest matrix), and one would expect some kind of Hamming-weight dependency to appear.

To check whether this is true, we can test for such dependencies using TestU01 [12], a well-known framework for testing generators, which implements tests related to Hamming weights from L’Ecuyer and Panneton [11] and Rukhin et al. [20] (please refer to the TestU01 guide for a detailed description of the tests). Table 1 shows the basic parameters of the nine tests we performed. The parameters were inspired by the author choices in the BigCrush test suite [12], but instead of analyzing \( 10^{9} \) or fewer outputs, as happens in BigCrush, we analyze up to \( 10^{13} \) 32-bit values (e.g., 40 TB of data), hoping that examining a much larger amount of data might help in finding bias in the generators. Besides the parameter inspired by BigCrush, following a suggestion from a referee, we also tried a number of power-of-two values for the L parameter of HammingIndep and HammingCorr.

Some of the generators have a 64-bit output, but TestU01 expects generators to return 32-bit integers, so we tested both the lower 32 bits, the upper 32 bits, and the upper and lower bits interleaved. (We do not discard any bits, as it is possible in TestU01.) In the case of the dSFMT, there is a specific call to generate a 32-bit value.

The disappointing results are reported in Table 2: despite the very sparse nature of the transition matrices of these generators, and the very large amount of data, only problems with the SFMT (already using \( 10^9 \) values), xorshift128+ (using \( 10^{10} \) values), and xoroshiro128+ (using \( 10^{13} \) values) are reported.⁵ All other p-values at \( 10^{13} \) are within the range \( [0.01..0.99] \).⁶

3 TESTING HAMMING-WEIGHT DEPENDENCIES

In this section, we introduce a new test for Hamming-weight dependencies that will find bias in all generators from Table 2. For the generators whose bias was detected by TestU01, the test will be able to obtain similar or better p-values using an order of magnitude fewer data.

Let us denote with \( \nu x \) [9] the Hamming weight of a w-bit word x—that is, the number of ones in its binary representation. We would like to examine the output of a generator and find medium-range dependencies in the Hamming weights of its outputs.

For example, the current output might have an average weight (close to \( w/2 \)) with higher-than-expected probability if three outputs ago we have seen a word with average weight, or the current output might have a non-average weight (high or low) with higher-than-expected probability depending on whether four outputs ago we have seen average weight and five outputs ago we have seen non-average weight.

We will consider sequences of w-bit values (w even and not too small, say \( \ge 16 \)) extracted from the output of the generator. Usually, w will be the entire output of the generator, but it is possible to run the test on a subset of bits, break the generator output into smaller pieces fed sequentially to the test, or glue (a subset of) the generator output into larger pieces.

The basic idea of the test is that of generating a vector whose coordinates should appear to be drawn from independent random variables with a standard normal distribution, given that the original sequence was random; apply a unitary transformation, obtaining a transformed vector; and derive a p-value using the fact that the coordinates of the transformed vector should still appear to be drawn from independent random variables with a standard normal distribution [24], given that the original sequence was random. The transform will be designed in such a way to make dependencies as those we described emerge more clearly.

First of all, we first must define a window we will be working on: thus, we fix a parameter k and consider overlapping k-tuples of consecutive w-bit values (ideally, the number of bits of state should be less than kw). We will write \( \langle x_0, x_1, ..., x_{k-1}\rangle \) for such a generic k-tuple.

Now we need to classify outputs as “average” or “extremal” with respect to their Hamming weight. We thus consider an integer parameter \( \ell \le w/2 \) and the map \( \begin{equation*} x\stackrel{d}{\mapsto }{\left\lbrace \begin{array}{ll}0& \text{if $\nu x\lt w/2 -\ell $,}\\ 1& \text{if $w/2 -\ell \le \nu x\le w/2 + \ell $,}\\ 2& \text{if $\nu x\gt w/2 + \ell $.} \end{array}\right.} \end{equation*} \) In other words, we compute the Hamming weight of x and categorize x in three classes: left tail (before the \( 2\ell +1 \) most frequent weights), central (the \( 2\ell +1 \) central, most frequent weights), and right tail (after the \( 2\ell +1 \) most frequent weights). The standard choice for \( \ell \) is the integer such that the overall probability of the \( 2\ell +1 \) most frequent weights is closest to \( 1/2 \). For example, for \( w=32 \) we have \( \ell =1 \), whereas for \( w=64 \) we have \( \ell =2 \).

We thus get from the k-tuple \( \langle x_0, x_1, ..., x_{k-1}\rangle \) a signature \( \langle d(x_0), d(x_1), ..., d(x_{k-1})\rangle \) of ktrits (base-3 digits), which we will identify with its value as a number in base 3: \( \begin{equation*} \sum _{i=0}^{k-1}d(x_i)3^{k-1-i}. \end{equation*} \)

Now, given a sequence of m w-bit values, for each signature s we compute the average number of ones in the word appearing after a k-tuple with signature s in the sequence. More precisely, a subsequence of the form \( \langle x_0, x_1, ..., x_k\rangle \) contributes \( \nu x_k \) to the average associated with the signature \( \langle d(x_0), d(x_1), ..., d(x_{k-1})\rangle \).

This bookkeeping can be easily performed using \( 3^k \) integer variables while streaming the generator output. For a large m, this procedure yields \( 3^k \) values with approximately normal distribution [5],⁷ which we normalize to a standard normal distribution; we denote the resulting row vector with \( \mathbf {v}=\langle {v}_0 \), \( {v}_1 \), \( \dots \, \), \( {v}_{3^k-1}\rangle \).⁸\( ^{,} \)⁹

We now apply to \( \mathbf {v} \) a Walsh–Hadamard-like transform, multiplying \( \mathbf {v} \) by the k-th Kronecker power¹⁰ \( T_k \) of the unitary base matrix (1) \( \begin{equation} M=\left(\begin{matrix}\frac{1}{\sqrt 3} & \phantom{-}\frac{1}{\sqrt 2} & \phantom{-}\frac{1}{\sqrt 6}\\ \frac{1}{\sqrt 3} & \phantom{-}0 & -\frac{2}{\sqrt 6}\\ \frac{1}{\sqrt 3} & -\frac{1}{\sqrt 2} & \phantom{-}\frac{1}{\sqrt 6}\\ \end{matrix}\right). \end{equation} \) Assuming that \( T_k \) is indexed using sequences of trits as numerals in base 3, the transform can be implemented recursively in the same vein as the fast Walsh–Hadamard transform (or any transform based on Kronecker powers), as if we write \( \mathbf {v}= [\mathbf {v}^0\;\mathbf {v}^1\;\mathbf {v}^2] \), where \( \mathbf {v}^0 \), \( \mathbf {v}^1 \), and \( \mathbf {v}^2 \) are the three subvectors indexed by signatures starting with 0, 1, and 2, respectively, we have by definition \( T_0=1 \) and (2) \( \begin{equation} \begin{split} \mathbf {v} T_k=[\mathbf {v}^0\quad \mathbf {v}^1\quad \mathbf {v}^2]\left(\begin{matrix}\frac{1}{\sqrt 3}T_{k-1} & \phantom{-}\frac{1}{\sqrt 2}T_{k-1} & \phantom{-}\frac{1}{\sqrt 6}T_{k-1}\\ \frac{1}{\sqrt 3}T_{k-1} & \phantom{--}0 & -\frac{2}{\sqrt 6}T_{k-1}\\ \frac{1}{\sqrt 3}T_{k-1} & -\frac{1}{\sqrt 2}T_{k-1} & \phantom{-}\frac{1}{\sqrt 6}T_{k-1}\\ \end{matrix}\right)\\ =\left[\frac{1}{\sqrt 3}(\mathbf {v}^0 +\mathbf {v}^1 +\mathbf {v}^2)T_{k-1} \quad \frac{1}{\sqrt 2}(\mathbf {v}^0 -\mathbf {v}^2)T_{k-1} \quad \frac{1}{\sqrt 6}(\mathbf {v}^0 -2\mathbf {v}^1 +\mathbf {v}^2)T_{k-1} \right]. \end{split} \end{equation} \) A detailed C implementation of \( T_k \) will be described in Section 4.3.

We will denote the transformed vector by \( \mathbf {v}^{\prime }= \mathbf {v} T_k \), and we shall write \( v^{\prime }_i \) for the transformed values. Since \( T_k \) is unitary, the \( v^{\prime }_i \)’s must appear still to be drawn from a standard normal distribution, and we can thus compute p-values for each of them. We combine p-values by dividing the indices of the vector \( \mathbf {v}^{\prime } \) in C categories \( {\mathscr{C}}_1 \), \( {\mathscr{C}}_2 \), \( \dots \, \), \( {\mathscr{C}}_{C} \) using the number of non-zero trits contained in their base-3 representation—that is, the number of non-zero trits in the associated signature: \( \mathscr{C}_j \), \( 1\le j\lt C \), contains indices with j non-zero trits, whereas \( \mathscr{C}_C \) contains all remaining indices, whose base-3 representation has at least C non-zero trits (\( C\le k \); usually, \( C=\lfloor k/2\rfloor + 1 \)). We discard \( v^{\prime }_0 \).

Given a category, say of cardinality c, we thus have a p-value \( p_i \) for each \( v^{\prime }_i \) in the category; we then consider the minimum of the p-values, say \( \bar{p} \), and compute a final category p-value by composing with the cumulative distribution function of the minimum of c independent uniform random variables in the unit interval [3, Equation (2.2.2)], obtaining the category p-value \( 1-(1-\bar{p})^c \). Finally, we take the minimum category p-value over all categories, and apply again the same cumulative distribution function with parameter C, since we are taking the minimum over C categories: this yields the final p-value of the test. Formally, \( \begin{equation*} p = 1 - \left(1 - \min _{1\le j \le C}\left(1 - \left(1 - \min _{i\in \mathscr{C}_j} p_i \right)^{\left|\mathscr{C}_i\right|}\right)\right)^C. \end{equation*} \)

The point of the transform \( T_k \) is that while \( v_i \) represents the (normalized) average number of ones after k previous outputs with density pattern described by the trit representation of i, \( v^{\prime }_i \) represents a combination of the average number of ones after k previous outputs satisfying different constraints: in the end, a unitary transformation is just a change of coordinates.

As a simple example, let us write the index \( \bar{\imath } \) of a transformed value \( v^{\prime }_{\bar{\imath }} \) as a sequence of trits \( {t}_1 \), \( {t}_2 \), \( \dots \, \), \( {t}_{k} \). If the trits are all zero, looking at (2) one can see that we are just computing the normalized sum of all values, which is of little interest: indeed, we discard \( v^{\prime }_0 \).

On the contrary, if a single trit, say in position \( \bar{\jmath } \), is equal to 1, \( v^{\prime }_{\bar{\imath }} \) is given by the sum of all \( v_i \)’s in which the \( \bar{\jmath } \)-th trit of i is 0 (\( \bar{\jmath } \) steps before we have seen few zeros) minus the sum of all \( v_i \)’s in which the \( \bar{\jmath } \)-th trit of i is 2 (\( \bar{\jmath } \) steps before we have seen many zeros): if the Hamming weight of the output depends on the Hamming weight of the output \( \bar{\jmath } \) steps before, the value of \( v^{\prime }_{\bar{\imath }} \) will be biased. Intuitively, if the transition matrix is very sparse, we expect vectors with low or high Hamming weight to be mapped to vectors with the same property.

If instead a single trit in position \( \bar{\jmath } \) is equal to 2, we will detect a kind of bias in which the Hamming weight of the current value depends on whether the Hamming weight of the output \( \bar{\jmath } \) steps before was extremal or average: more precisely, whether the value is larger when the Hamming weight of the output \( \bar{\jmath } \) steps before was average and smaller when the Hamming weight of the output \( \bar{\jmath } \) steps before was extremal (and vice versa). Intuitively, we expect that the shift/rotate-XOR of states with a very small or very large number of ones will have a small number of ones (in the first case, by sparsity, in the second case, by cancellation).

More complex trit patterns detect more complex dependencies: the most interesting patterns, however, usually are those with few non-zero trits, as a zero trit acts as a “don’t care about that previous output”: this property is immediate from (2), as the first \( 3^{k-1} \) output values, which correspond to a “don’t care” value in the first position, are obtained by applying recursively \( T_{k-1} \) over the renormalized sum of \( \mathbf {v}_0 \), \( \mathbf {v}_1 \), and \( \mathbf {v}_2 \), thus combining the values associated with signatures identical but for the first trit.

This is also why we assemble p-values by categories: by putting indices with a higher chance of giving low p-values in small categories, the test becomes more sensitive.

4 IMPLEMENTATION DETAILS

We will now discuss some implementation details. To be able to perform our test in the petabyte range, it must be engineered carefully: in particular, the main loop enumerating the output of the generator and computing the values \( v_i \) must be as fast as possible. Counting the number of ones in a word can be performed using single-clock specialized instructions in modern CPUs. The \( v_i \)’s are stored in an array of \( 3^k \) elements indexed by the value of a signature as a numeral in base 3, as required by the recursive implementation of \( T_k \) (see Section 4.3). One can very easily keep track of the current trit signature value by using the update rule \( s\leftarrow \lfloor s / 3\rfloor + t\cdot 3^{k-1} \), where t is the next trit.

We can replace the division with the fixed-point computation \( \lfloor (\lceil 2^{32}/3\rceil s) / 2^{32}\rfloor \) (this strategy works up to \( k=19 \) using 64-bit integers), so by precomputing \( \lceil 2^{32}/3\rceil \) and \( 3^{k-1} \), the costly operations in the update of s can be reduced to two independent multiplications.

4.3 Implementing the Transform T_k

In Figure 1, we show an in-place, recursive C implementation of the transform \( T_k \) defined in Section 3. The code is similar to analogous code for the Walsh–Hadamard transform or similar transforms based on Kronecker powers.

Fig. 1.

View Figure

The code assumes that the \( 3^k \)-dimensional vector \( \mathbf {v} \) is represented in the array v. The value associated with each signature is stored in a position equal to the signature (considered, as usual, as a base-3 numeral). In particular, the first \( 3^{k-1} \) values correspond to signatures of the form 0s, the following \( 3^{k-1} \) values to signatures of the form 1s, and the last \( 3^{k-1} \) values to signatures of the form 2s. The function transform() must be invoked on v with the additional parameter sig set to \( 3^{k-1} \).

We first note that if \( k=1, \) the function will just execute once the body of the for loop, resulting in the in-place multiplication of the 3-dimensional vector \( \mathbf {v} \) by the base matrix M, as expected.

In the general case, the code scans the three subarrays v, p1, and p2, of length \( 3^{k-1} \), which as discussed previously correspond to signatures starting with 0, 1, and 2, respectively. With the notation of (2), these subarrays correspond to the subvectors \( \mathbf {v}^0 \), \( \mathbf {v}^1 \), and \( \mathbf {v}^2 \), respectively, and it is immediate that the three subvectors appearing in the final result of (2) are computed in place by the for loop. After that computation, the recursion applies by induction the transform with one dimension less to each of the three subarrays in place. We conclude that transform() implements correctly in place the transform \( T_k \).

5 CONCLUSION

We have described a new test for Hamming-weight dependencies based on a unitary transform. Properly implemented, the test is very powerful: for example, in a matter of hours, it finds bias in the dSFMT with 521 bits of state and in xoroshiro1024+; it can even find bias in WELL512, even though its transition matrix is much denser, and in the Tiny Mersenne Twister. For these generators, no bias was previously known beyond linearity tests. In particular, the Hamming-weight tests in TestU01 [12], a state-of-the-art testing framework, are unable to find any bias in several generators of Table 3, whereas all those generators fail our test.

Our test is very effective on \( \mathbf {F}_2 \)-linear generators with relatively sparse transitions matrices, particularly when wk is not smaller than the number of bits of state of the generator. In practice, the best results are obtained on generators with less than a few thousand bits of state.

Similar to linearity tests, a failure in our test is an indication of lesser randomness, but in general the impact will depend on the application. We do not expect dependencies like those in xorshift128+ or the SFMT to be pernicious, but they highlight a weakness of the associated \( \mathbf {F}_2 \)-linear transformation.

ACKNOWLEDGMENTS

The authors would like to thank Jeffrey Hunter for useful pointers to the literature about mean passage times in Markov chains and Pierre L’Ecuyer for a number of suggestions that improved the quality of the presentation.