A test for Hamming-weight dependencies in pseudorandom number generators

This crate implements the test for Hamming-weight dependencies described in “A New Test for Hamming-Weight Dependencies”, by David Blackman and Sebastiano Vigna. The test finds bias induced by dependencies among the Hamming weights (the number of ones) of the outputs of a pseudorandom number generator (PRNG), in particular for generators based on 𝐅₂-linear transformations such as the dSFMT, xoroshiro128+, and WELL512.

The crate is a port of the original C implementation: for the same generator and parameters, the result of the test is the same. Moreover, the crate can use multiple cores to generate and examine the output of the generator in parallel, again without changing the result of the test.

The test

The test examines a sequence of w-bit words extracted from the output of the generator: usually w is the size of the output, but it is possible to split each output into two halves, or to glue consecutive outputs into larger words. Each word is classified by its Hamming weight as average or extremal: the trit (base-3 digit) of a word is 1 if its weight is one of the 2ℓ + 1 central weights, whose overall probability is close to 1/2 (ℓ = 1 for w = 32, ℓ = 2 for w = 64, ℓ = 3 for w = 128), and otherwise 0 or 2, depending on whether the weight is smaller or larger.

The trits of k consecutive words form a signature, and for each of the 3ᵏ signatures the test computes the average Hamming weight of the words that follow the signature in the sequence. For a random sequence, these averages, suitably normalized, look like independent standard normal values. The vector of normalized averages is then multiplied by the k-th Kronecker power of a 3 × 3 orthonormal matrix, in the same vein as the Walsh–Hadamard transform: since the transform is unitary, for a random sequence the transformed values still look like independent standard normal values. However, each transformed value combines the averages of many signatures: a zero trit in its index means “don't care about that previous word”, a one compares words with low and high weight, and a two compares extremal and average words. Values with few nonzero trits thus make simple dependencies emerge, such as “the weight of a word depends on the weight of the word three steps before”.

The transformed values are divided into C categories by the number of nonzero trits of their index (all values with at least C nonzero trits are in the last category; by default, C = ⌊k / 2⌋ + 1). The test computes the p-value of the most extreme value of each category, corrected for the size of the category, and the final p-value is the smallest category p-value, corrected for the number of categories.

For each category, the test reports also the signature of the most extreme value, written with the most recent word on the right: when a generator fails, the faulty signature provides insight into its structure. For example, xorshift1024 fails with signature 2000000000000001, as each output depends on the first and last word of its 16-word state.

There is also a transitional variant of the test (option -t), which xors the sequence of words, seen as a stream of bits, with itself shifted by one bit, and then examines the result: in this case, the test looks for dependencies between bit transitions.

Running the test

The output of the generator is examined in batches, whose size depends on w and k (it can be reduced with --max-batch-size). The test stops after a given number of bytes, or runs indefinitely if no number is specified. With --progress (or when running indefinitely), p-values are reported after about 10⁸, 1.25 · 10⁸, 1.5 · 10⁸, 1.75 · 10⁸, 2 · 10⁸, 2.5 · 10⁸, 3 · 10⁸, 4 · 10⁸, 5 · 10⁸, 6 · 10⁸, 7 · 10⁸, 8.5 · 10⁸ bytes, and so on, multiplying by ten each time; the test stops as soon as a p-value is below the threshold set with --low-pv. In the paper, tests were performed with w = 32 or w = 64 and k between 8 and 19, stopping after a petabyte of data or at a p-value below 10⁻²⁰.

The test needs 28 · 3ᵏ bytes of memory, which become about 1.2 GB for k = 16 and 33 GB for k = 19. To improve locality, words are counted in packed 32-bit counters, which are moved into 64-bit counters at the end of each batch. Batch sizes are chosen so that the probability that a packed counter overflows is below 10⁻¹⁰⁰: if this happens, the test reports a p-value of 10⁻¹⁰⁰.

Results

The following table, from the page of the test, reports for some 𝐅₂-linear generators the amount of output, in bytes, after which the test yields a p-value below 10⁻²⁰, and the faulty signature. Ranges (→) appear when several variants were tested (e.g., several parameters of the Mersenne Twister): a missing right extreme means that some instances did not fail the test within a petabyte.

Since the generators of this crate are seeded differently (see below), the amount of output and the faulty signatures might be slightly different.

Parallel generation

With option -P, the iterations of the test are split into contiguous ranges, one for each thread: jump-capable generators jump to the start of their range, whereas the others reach it with a sequential pre-scan. Each thread counts the words of its range in its own copy of the packed counters, starting a few words before its range so to have the same signature as a sequential run; at the end of each batch, the copies are added together, and since packed counters are updated by addition modulo 2³², the result is exactly the one of a sequential run, even in case of overflow. Thus, the output of a parallel run is identical to that of a sequential run.

Each thread needs its own copy of the packed counters, so a parallel run with t threads uses at least 4 · (t − 1) · 3ᵏ additional bytes of memory (small values of k use a few more copies to make work units longer). The number of threads is governed by the environment variable RAYON_NUM_THREADS (by default, the number of cores).

All the 𝐅₂-linear generators of the paper provide arbitrary jumps using their characteristic polynomial (see the f2 module).

Differences from the C implementation

• Parameters are run-time options rather than compile-time macros: -w (HWD_BITS), --prng-bits (HWD_PRNG_BITS), -k (HWD_DIM), -c (HWD_NUMCATS), and --max-batch-size (MAX_BATCH_SIZE). The generator is selected at compile time using Cargo features.

• Generators return 64-bit values: 32-bit generators return their output in the upper 32 bits, and must be tested with --prng-bits 32. With 64-bit outputs, -w 32 splits each output into two words, upper half first.

• Generators are initialized from a 64-bit seed (option -S), and generators with more than 64 bits of state fill it using SplitMix64: the faulty signatures and the amount of data needed to find bias are thus slightly different from those reported in the paper.

• p-values are printed with full precision, and the reports contain the exact number of bytes examined, but not the current date.

Usage

The generator to test is selected at compilation time using Cargo features: besides incr, a counter useful only to test the crate, the available generators are those of the paper, xorshift128, xorshift128plus, xorshift1024, xorshift1024plus, xoroshiro128, xoroshiro128plus, xoroshiro1024, xoroshiro1024plus, and the 32-bit generator well512a. For example,

cargo run -r -F xoroshiro128plus -- -P --progress --low-pv=1e-20 1e15
Generator: xoroshiro128+ (24, 16, 37)
Seed: 0x0000000000000000
Running a test for Hamming-weight dependencies with k = 8 and 5 categories using 10 parallel generators (jump-ahead) on 64-bit words, the full outputs, analyzing bits (batches of 2311072 words, 0.002 GiB RAM)
Examining 1000000000000000 bytes, stopping if p < 1e-20
mix3 extreme = 2.57485 (sig = 00020000) weight 1 (16), p-value = 0.14893176321045537
mix3 extreme = 2.95929 (sig = 00001020) weight 2 (112), p-value = 0.2924060708331487
mix3 extreme = 2.97046 (sig = 01002200) weight 3 (448), p-value = 0.736610976840799
mix3 extreme = 2.96745 (sig = 21200100) weight 4 (1120), p-value = 0.9655482751525518
mix3 extreme = 4.32354 (sig = 22120101) weight >=5 (4864), p-value = 0.0719628674515154
bits per word = 64 (analyzing bits); min category p-value = 0.0719628674515154

processed 110931456 bytes in 0.057 seconds (1.9524 GB/s, 7.0286 TB/h)

p = 0.3116223400948752
------

[...]

mix3 extreme = 2.21810 (sig = 00002000) weight 1 (16), p-value = 0.3498212859734305
mix3 extreme = 11.28875 (sig = 00000012) weight 2 (112), p-value = 1.6703369397814162e-27
mix3 extreme = 3.43838 (sig = 20000210) weight 3 (448), p-value = 0.23067832880149697
mix3 extreme = 3.21723 (sig = 20221000) weight 4 (1120), p-value = 0.7655706409615056
mix3 extreme = 3.37608 (sig = 21101022) weight >=5 (4864), p-value = 0.9720588981323016
bits per word = 64 (analyzing bits); min category p-value = 1.6703369397814162e-27

processed 10000008544000 bytes in 435.037 seconds (22.9866 GB/s, 82.7517 TB/h)

p = 8.351684698907081e-27
Test completed in 435.04 seconds

will test xoroshiro128+ with w = 64 and k = 8 on an Apple M1 Max with 10 cores, finding bias after 10¹³ bytes in about seven minutes, with the same faulty signature reported in the paper (00000012). The test processes about 23 GB/s; the C implementation, which is sequential, processes about 3.3 GB/s on the same hardware, and would need about 50 minutes.

The default parameters (w = 64, k = 8) are those used for 64-bit generators with 128 bits of state; generators with a larger state need a larger k. For example, to replicate the tests of the paper on xorshift1024 and on the 32-bit generator WELL512a, use

cargo run -r -F xorshift1024 -- -P -k 16 --progress --low-pv=1e-20 1e15
cargo run -r -F well512a -- -P -w 32 --prng-bits 32 -k 16 --progress --low-pv=1e-20 1e15

The repository configures Cargo to compile for the native CPU (-C target-cpu=native): on x86-64, for example, otherwise population counts do not use the POPCNT instruction. Since cargo install does not use this configuration, install the test with, for example,

RUSTFLAGS="-C target-cpu=native" cargo install hwd -F xoroshiro128plus

Adding your own generator

To add a new generator, add a feature in Cargo.toml and a corresponding implementation in the prng module. If skipping is possible, you can implement the try_skip method, which must succeed for every offset or for none (see the prng module documentation): for 𝐅₂-linear generators, you just need to implement the LinearGenerator trait and provide the characteristic polynomial of the transition map, as done, for example, for xoroshiro128.

The prng module has the same interface as that of the coll-birth crate.