**Generator 5:** We use the same generator as in Generator 5, but we take only one bit from each value retuned by `xorshift64s()` (in our experiments we took the 15th least significant bit of the outcome).

**Generator 6:** In this example we use random bits produced by the well known generator called Mersenne Twister [35]. We used the standard version implemented in 2002 by Takuji Nishimura and Makoto Matsumoto, [46] (we used all 32 bits from each value returned by the generator).

**Generator 7:** In this example we use a permuted congruential generator (PCG), which is described as having good statistical performance despite small and fast code [45]. We used the following code borrowed from the PCG library on https://www.pcg-random.org.

```c
typedef struct { uint64_t state; uint64_t inc; } pcg32_random_t;
uint32_t pcg32_random_r(pcg32_random_t* rng)
{
    uint64_t oldstate = rng->state;
    // Advance internal state
    rng->state = oldstate * 6364136223846793005ULL + (rng->inc|1);
    // Calculate output function (XSH RR), uses old state for max ILP
    uint32_t xorshifted = ((oldstate >> 18u) ^ oldstate) >> 27u;
    uint32_t rot = oldstate >> 59u;
    return (xorshifted >> rot) | (xorshifted << ((-rot) & 31));
}
```

*Observed results:* In our experiments, generators 1, 2, 5, 6, 7 pass the standard tests of Dieharder; also the spectral tests do not reveal for them any abnormality. Generator 3 and 4 fail the battery of tests of Dieharder. It is interesting that Generator 3 produces random graphs with a spectral radius pretty similar to the properties of “truly random” graphs (we performed experiments, in particular, for graphs with 256 vertices and average degree 16). On the other hand, the random graphs (with the same number of vertices and average degrees) obtained with Generator 4 have dramatically different spectral properties.

In our experiments, each graph was represented by a rather short sequence of bits (typically by a block of 4096 bytes). We believe that spectral properties of much larger graphs deserve more careful research in future.

## 6.12 Concluding remarks

Our analysis of the devices is quite superficial: for a more serious one we would need to analyze the PCB, the firmware, make extensive testing in different environmental conditions, etc. However, some observations could be made.

* There are significant improvements since the time when Marsaglia developed *dieharder* for randomness testing (1990s): only one of the devices tested by us failed immediately our health tests (while all three hardware devices tested by Marsaglia did fail miserably, even after correction for the dos-unix conversion problem).

* There are several kinds of threats. First, one may be afraid of hardware instability in different conditions. For that, one should avoid devices that use the environment as a source of randomness (as it is done by gniibe/neug device) since they require extensive testing in different environmental conditions to be sure that enough randomness is available. The devices that use electromagnetic noise are also unreliable in this respect.
* One can be afraid of a trojan horse installed by the device maker or on the way. Then one should avoid device with microprocessors and prefer to use device with easy to check or open architecture (like bitbabbler, or, better, the infinite noise device). The ultimate solution would be to make the device in-house using one of the standard approaches or the open-source hardware data; the bit shift devices or the noise sources with some off-the-shelf ADC or even a commercial audio card could be a reasonable solution.

* One can use inadvertent deficiencies in the device due to bad or instable components. Then the ability to get the raw data for systematic check on the go is important (so you can see if, say, a Zener diode changes its parameters suddenly). For example, from this viewpoint TrueRNGpro defice is much better than TrueRNGv3 one.

* The most paranoid people may be afraid of security holes in widely used cryptographical primitives like hash functions and other tools prescribed by the standard (recall that even accepting unproven assumptions from complexity theory we have no mathematical reasons to trust these standard tools). Then simple and well-understood tools like xor of several sources are preferable even if they require more raw random bits.

* On the other hand, if we care about small statistical regularities and do not worry much about malicious adversaries, the cryptographic primitive are usually more or less enough, and they still provide minimal defense against most intruders even if the statistical properties of the primary source degrade. They also allow a trade-off between speed and security since we may get much more than *N* pseudo-random bits from *N* random bits and still be practically indistinguishable from true random source. The resulting generator can be viewed as pseudo-random generator with random seed and continuous reseeding from hardware random source.

* Note that it is easy to combine protective measures against different kind of threats just by xor-ing the outputs of several generators that are resistant to different attacks: it is enough to have one correctly working generator to guarantee the correct distribution at the output. This assumes that we do not doubt the independence of the sources — but for the most practical cases having several usb devices bought independently from different manufacturers is probably enough (assuming that our application is able to xor the output sequences correctly).

* There are more exotic approaches; for example, some people use lava lamps as sources of randomness: https://qz.com/1642628/cloudflare-uses-lava-lamps-to-generate-a-crucial-resource/ (Figure 6.37) — this is definitely good for advertisements and provides some level of “security through obscurity”. One can also note that a cheap noisy video camera (preferably with minimal post-processing) has a lot of internal noise that is available via a high bandwith channel.

* Of course, one can also use existing high-end solutions to get a high-speed stream of random bits, or some certification, or some additional feeling of security, or a chance to use it for advertisements. Still it is questionable whether this is really much more than a “security theater”.

**Warning:** In this analysis we did not consider an important aspect: how much additional programming / system administration / hardware installation work is needed for this or that generator. This depends on the specific nature of the application (e.g., the operating system in use, where the generators should be installed, etc.), and can make some approaches practically unusable.
# Chapter 7

## Comments on the software used

We start (Section 7.1) with observations made when analyzing the existing code (mainly from `dieharder` suite since it is widely used). Then we describe some software tools that we developed for our experiments. Section 7.2 describes the tools related to the Kolmogorov-Smirnov two-sample test. Section 7.3 describes the utility that includes robust version of more than 15 tests; Section 7.5 describes some graphic tools (very primitive, but useful for quick assessment). Finally, Section 7.6 describes the implementation of a simple two-source extractor with provably good properties based on finite fields (discussed in Section 5.5).

### 7.1 Dieharder: some notes and comments

Since Dieharder [9] utility is often used for practical testing of hardware random generators, it makes sense to look at it more closely and note some peculiarities. Our remarks refer to the version 3.31.1 (which seems to be the current one in March 2021).

* Flag `-T` is mentioned in the code (`dieharder/parsecl.c`, lines 97–100)

```c
if(i == TCNT) {
    fprintf(stderr, "Invalid -T %s option.\n", optarg);
    exit(1);
}
```

It seems that `-T` options is not mentioned in the man page. This message appears when the name given for the `-D` option does not exist in the table of names (so `-T` should be probably replaced by `-D`). The same `-T` is mentioned in the comments in `dieharder/dieharder.h`, lines 18–19:

* `is` set in `set_globals` to a default to be used if `-T 0` is selected.

* `tflag` is otherwise accumulated from a series of `-T` FLAG calls, where

* Too many random bits read in the bitstream test for overlapping (main) mode, `libdieharder/diehard_bitstream.c`, lines 73–78:

* for non-overlapping samples we need $(2^{21})*5/8 = 1310720$ uints, but

* for luck we add one as we'd hate to run out. For overlapping samples,

* we need $2^{21}$ BITS or $2^{18} = 262144$ uints, again plus one to be sure

* we don't run out.

*/

#define BS_OVERLAP 262146
In fact, for overlapping samples we need $2^{21} + 19$ bits, i.e., about $2^{21}/32 = 32$-bits integers (uint). This is $2^{21-5} = 2^{16} = 65536$ (+1 for the +19 additional bits). So the last line should be

```c
#define BS_OVERLAP 65537
```

(Of course, a good generator should pass the test even if some bits are ignored, but in this way the test uses four times more bits than needed, and this requires more time for hardware generator testing.)

* In the same test, there is an error in counting the strings that do not appear as 20-bit substrings of the generator. Namely, line 172 is

```c
w[w20]++;
```

it increments the counter of the number of appearances of bit string w20. There is a small problem: w is defined as `*char` in line 113:

```c
w = (char *)malloc(M*sizeof(char));
```

so the counting happens modulo 256. In the rare case where some string appears 256, 512, etc. times (it should be very rare, since the expected number of occurrences is only 2, still theoretically possible), it will be counted as missing. The obvious correction is to make the counter `bool` (we do not need the count value, only the yes/no information), replacing the line 172 by `w[w20]= 1;` this solves the problem.

* The dieharder bitstream test ignores 8 first bits of the stream both in overlap and non-overlap modes (for different reasons). For the non-overlap mode the bitstream is different, since the memory is considered as a sequence of bytes, not 32-integers, as it is done in the overlap mode, and byte ordering is different. The lines 120-12 in `diehard_bitstream.c` are actually needed only for the overlap case, but increase `i` by 1 anyway, so in the non-overlap mode the first byte is ignored.

* The dieharder bitstream test does not react properly to the `-L` option (always uses overlapped sampling), because line 52 of `diehard_bitstream.c` says so (and local variable is used instead of the global one).

* With small number of samples, `p`-value is sometimes reported as `1.00000` and test fails (this should not happen).

```bash
$ dieharder -g 1 -p 2 -d 1
```

```bash
#============================================================#
```

```bash
# dieharder version 3.31.1 Copyright 2003 Robert G. Brown #
```

```bash
#============================================================#
```

```
rng_name |rands/second| Seed |
cmrg| 4.63e+07 |3941174114|
```

```bash
#============================================================#
```

```
test_name |ntup| tsamples |psamples| p-value |Assessment
```

```bash
#============================================================#
```

```
diehard_operm5| 0| 1000000| 2|1.00000000| FAILED
```

The source of the error is that (unfortunately!) the computation of Kolmogorov-Smirnov statistics in `kstest.c` is incorrect (e.g., one should never divide by `count+1`), so for small sample length `n` it produces impossible values (like `p = 0.05` for `n = 2`) and therefore the `ks_test` function returns 1, interpreted as a failure.
* The call `dieharder -g 1 -S 1 -d 1` produces different results when run several times, though `-S 1` should fix the seed. The reason is that in `choose rng.c` the seed is given to the generator before measuring the speed (line 246), and reseeding happens when measuring the speed (lines 251-253) if `-D 8192` (speed measuring) is on, so the result depends on the display options (a rather confusing behavior);

```
\texttt{dieharder -d 1 -g 1 -p 10 -S 1 -s 0 -D pvalues -D 8192}
```

produces different values too (though it seems that strategy prevents reseeding and the seed is fixed). This can be corrected if we reseed again after timing tests in `choose rng.c`. However, `rgb_timing()` changes the global variable seed. To avoid this, one could make variable seed local (as unsigned `long int`) in the `rgb_timing.c`, function `rgb_timing()`.

* The implementation of Marsaglia `opso` test [31] the code (`diehard_opso.c`) does not correspond the description both in the paper and in the code itself. The description suggests to consider a sequence of $n = 2^{21} + 1$ integers, consider them as 32-bit words, and keep only 10 bits in some fixed position. In this way we get a word in 1024-letter alphabet having length $2^{21} + 1$. Then we count how many of the $1024 \times 1024$ two-letter words in this alphabet do not appear among $2^{21}$ factors of length 2 of this $n$-bit strings. (This is a random variable whose distribution is assumed to be very close to a normal one, and this is checked.)

At the same time, code says (here $t$ is an integer variable that is incread by 1 in the for-loop):

```c
if(t%2 == 0) {
    j0 = gsl rng_get(rng);
    k0 = gsl rng_get(rng);
    j = j0 & 0x03ff;
    k = k0 & 0x03ff;
} else {
    j = (j0 >> 10) & 0x03ff;
    k = (k0 >> 10) & 0x03ff;
}
w[j][k] = 1;
```

We see that instead of overlapping pairs (as the name of the test says) here two pairs of 10-bit strings are taken from different (disjoint) positions in a pair of two 32-bit integers obtained from a random number generator ($j0, k0$). These two pairs are registered in the bit array $w$, and then a new pair of integers is obtained. So the code makes something completely different. (In fact, the `dieharder` man page characterizes this test as "suspicious".)

* A similar problem exists with another test from Marsaglia, namely, the `oqso` test [31]. According to the description (reproduced in the `dieharder` source), this test uses 5-bit groups (instead of 10-bit groups), gets a word of length $2^{21} + 3$ in 32-letter alphabet, and then considers $2^{21}$ factors of this word of length 4, and counts the 4-letter words that do not appear as factors. This is done for every of $28 = 32 - 5 + 1$ positions, so we get 28 counters.

The code in `diehard_oqso.c` does something completely different:

```c
for(t=0;t<test[0]->tsamples;t++){
    if(t%6 == 0) {
        i0 = gsl rng_get(rng);
        j0 = gsl rng_get(rng);
```
```c
k0 = gsl_rng_get(rng);
l0 = gsl_rng_get(rng);
boffset = 0;
}
/*
 * Get four "letters" (indices into w)
*/
i = (i0 >> boffset) & 0x01f;
j = (j0 >> boffset) & 0x01f;
k = (k0 >> boffset) & 0x01f;
l = (l0 >> boffset) & 0x01f;

w[i][j][k][l]=1;
boffset+=5;
}
```

Instead of taking (overlapping) factors, it makes 6 quadruples from disjoint bits in four 32-bit integers. (There is no reason to think that the resulting test is worse than the original oqso test, but it does not match the description in the same file.)

* In the runs test (that goes back to Knuth), the code in `dieharder` (file `diehard_runs.c`, lines 132-142) the first and the last elements of the tested block of integers are compared (with no apparent reason, Knuth's description and `runtest.c` did not do anything like that).

* The overlapping sum test is described in [29] as follows. We take the sum of *m* independent random variables uniformly distributed in [0, 1] (obtained from *m* unsigned integers). Compute their sum; its distribution is approximately normal. Then replace one of the number by a fresh random real in [0, 1]; we get another (heavily correlated) random variable whose distribution is close to normal. Doing this *m* − 1 times, we use 2*mal random reals and produce *m* overlapping sums. Each of the sums is approximately normal, and they are correlated. Marsaglia suggests to apply a linear transformation that will convert these (approximately) normally distributed variables to (approximately) independent ones, convert them to (approximately) uniformly distributed independent variables and to use Kolmogorov–Smirnov one-sample test to check this hypothesis (of uniformly null distribution). The original test used *m* = 100, so from every 199 unsigned integer we got one *p*-value. Then the distribution of 100 these *p*-values was assumed to be uniform and again Kolmogorov–Smirnov test was applied to them to get a secondary *p*-value. Finally, ten secondary *p*-values were submitted to one more (“tertiary”) Kolmogorov–Smirnov test.

Obviously this procedure is based on several layers of approximations, and no attempts to provide an error bound were made. (Also the derivation of the suggested linear transformation was not given, just some formulas in file `die.c/cdosum.c` from [29].)

The `dieharder` code (file `diehard_sum.c`) complains about these problems saying that “unfortunately, Marsaglia’s description is pretty opaque”. Still the code from `cdosum.c` for the linear transformation is reproduced, but with a change. Here is a fragment from `diehard_sum.c`:

```c
x[0] = y[0]/sqrt(1.0*m);
x[1] = -x[0]*(m-1)/sqrt(2.0*m - 1.0) + y[1]*sqrt(m/(2.0*m - 1.0));
[...]

for(t=2;t<m;t++){
```
[...]
a = 2.0*m + 1.0 - t;
b = 2.0*a - 2.0;
/* x[t] = y[0]/sqrt(a*b) - y[t-1]*sqrt((a-1.0)/(b+2.0)) + y[t]*sqrt(a/b); */
x[t] = y[t-2]/sqrt(a*b) - y[t-1]*sqrt((a-1.0)/(b+2.0)) + y[t]*sqrt(a/b);

The commented line is taken from Marsaglia's code; however, the first term is changed (for no apparent reason and with no explanation).

It is noted in the beginning of `diehard_sum.c` that "At this point I think there is rock solid evidence that this test is completely useless in every sense of the word. It is broken, and it is so broken that there is no point in trying to fix it. If I crank up tsamples to where one can actually see the asymptotic distribution (which should be uniform) it is visibly non-uniform and indeed the final kstest converges to a non-zero pvalue of 0.09702690 for all rings [random number generators] tested, which hardly seems useful. If one runs the test on only 100 samples (as Marsaglia's code did in both the original fortran or the newer C version) but then increases the number of runs of the test from the default 100, it is easy to make it fail for gold standard generators. The test just doesn't work. It cannot be used to identify a generator that fails the null hypothesis. Don't use it."

We tried Marsaglia's recipe to get (presumably) independent uniformly distributed variables and its version from `dieharder` by computing 20000 overlapping sums from three generators. (Figure 7.1). As we see, these six curves are quite close to each other and differ significantly from the diagonal (which would appear if the distribution is truly uniform). So it seems that indeed this test uses approximation assumptions that are too rough (both in the original form and in the changed one). Its two-sample version, of course, is still valid (as it happens for every test function), only the sensitivity could be a problem.

* File `dieharder_rank_32x32.c` starts with comments

```
/* This is the Diehard BINARY RANK 31x31 test, rewritten from the description in tests.txt on George Marsaglia's diehard site. */
/* This is the BINARY RANK TEST for 31x31 matrices. The leftmost ::
/* 31 bits of 31 random integers from the test sequence are used ::
/* to form a 31x31 binary matrix over the field {0,1}. The rank ::
/* is determined. That rank can be from 0 to 31, but ranks< 28 ::
/* are rare, and their counts are pooled with those for rank 28. ::
/* Ranks are found for 40,000 such random matrices and a chisqua-::
/* re test is performed on counts for ranks 31,30,29 and <=28. ::
```

Still it tests ranks of 32 × 32 matrices (according to the file name and `dieharder` messages).

* For the 6 × 8 rank test the `dieharder` comments look even more confusing. The description says

```
/* This is the Diehard BINARY RANK 6x8 test, rewritten from the description in tests.txt on George Marsaglia's diehard site. */
/* This is the BINARY RANK TEST for 6x8 matrices. From each of ::
/* six random 32-bit integers from the generator under test, a ::
/* specified byte is chosen, and the resulting six bytes form a ::
```
<figure><img src="image_48.png" /><figcaption>Figure 7.1: Six empirical distributions of presumably uniform samples obtaind using three different generators. Three use the original Marsaglia’s transformation recipe, and the other three incorporate dieharder’s changes).</figcaption></figure>
* 6x8 binary matrix whose rank is determined. That rank can be ::
* from 0 to 6, but ranks 0,1,2,3 are rare; their counts are ::
* pooled with those for rank 4. Ranks are found for 100,000 ::
* random matrices, and a chi-square test is performed on ::
* counts for ranks 6,5 and <=4. ::

It is not clear how the words "specified byte" should be understood<fnref n="1" /> The comments in the code say

/*

* We generate 6 random rmax_bits-bit integers and put a
* randomly chosen byte into the LEFTMOST byte position
* of the row/slot of mtx.
*/

But there are traces of the “random choice” of a byte, it seems that just the 8 first bits are
selected from each 32-bit integer. (There are also strange things in the code. For example, six
rows and eight columns with integer values are reserved for the matrix:

```c
mtx = (uint **)malloc(6*sizeof(uint *));
for(i=0;i<6;i++){
    mtx[i] = (uint *)malloc(8*sizeof(uint));
}
```

which is strange since the rank computations in rank.c use individual bits in unsigned integers (and only mtx[i][0] are used later in the code).

* In the birthday tests (`cdbday.c` in `diehard` and `diehard_birthdays.c`) there are several peculiarities in the code. The description (given by Marsaglia are reproduced by Brown) says that *m* days are randomly chosen in a year of *n* days, then the intervals between these “birthdays” are considered and the “number of values that occur more than once in that list” is computed.<fnref n="2" /> This distribution is claimed to be asymptotically close to Poisson distribution with parameter *m*³/4*n*, and a χ<fnref n="2" />-tests is used to measure this closeness. According to Marsaglia’s description, “A sample of 200 *j*'s is taken and a chi-square goodness of fit test provides a *p*-values”. (In fact, Marsaglia’s code uses 500 values of *j*, not 200. Also, Marsaglia uses values *m* = 2¹⁰ and *n* = 2²⁴, both in the comments and the code in `cdbday.c`, but Brown, reproducing Marsaglia’s comments, cites another value *m* = 2⁹, along with *n* = 2²⁴. So the Poisson parameter becomes 2 instead of 16.

Marsaglia gets random 24-bit words (birthdays) from 32-bit random integers (that are tested) using bits 1-24, then using bits 2-25 (from the same integers), etc., thus getting 9 different *p*-values from the *chi*<fnref n="2" />-test. Then “the nine *p*-values provide a sample for a KSTEST”. This process is based on several unfounded assumptions. First, the approximation by Poisson distribution is only an approximation (no derivation is given, not to speak about error bounds). Then χ<fnref n="2" />-test needs to be applied, and the standard practice is to combine the values with small probabilities (rule of thumb says the we should combine tail bin in such a way that the

<footnote n="1">¹It seems that in the original *diehard* code the χ²-test was applied for each of 25 position separately (which is reasonable), but then Kolmogorov-Smirnov test was applied to the 25 results given by χ²-tests, and this is a bad idea, since they are not independent.</footnote>

<footnote n="2">²Strangely, both Marsaglia and Brown consider the interval that crosses the year boundary in a special way: only the part in the new year is counted.</footnote>
expected value for a bin is at least 5). Marsaglia does this according to the rule of thumb, while Brown considers only the right tail and just ignores that bins with expectation less than 5 instead of combining them). As Brown writes in *diehard_bIRTHDAYS.c*, if the number of repetitions (j) is bigger than the allowed maximal value (called kmax), “we simply ignore — it is a bad idea to bundle all the points from the tail into the last bin, as a Poisson distribution can have a lot of points out in that tail!”. (This argument looks very strange, and we cannot explain the reasoning.)

There are two other points to be mentioned. One of them is the interpretation of the words “Let $j$ be the number of values that occur more than once in that list” (of intervals between birthdays). Should we literally count just the values or, if there are three intervals of the same length, we could count them as two repetitions, so $j$ is the difference between the cardinalities of a multiset of intervals and set of intervals? Marsaglia uses the second interpretation, and Brown uses the first one (“We count the number of interval values that occur more than once in the list. Presumably that means that even if an interval occurs 3 or 4 times, it counts only once!”).

The second point is the overlapping sampling done by Marsaglia (bits 1–24, then 2–25, etc.). Brown writes in the comments that “I’m running it on all 32 24-bit strings implicit in the variables. We’ll do this by rotating each variable by one bit position in between a simple run of the test. A full run will therefore be 32 simple (rotated) runs on bits 1–24”. But this is not what happens in the code, it seems<fnref n="3" /> that the 24-bit portions are taken from the input bit stream without any overlap (which is perfectly OK and makes the different $p$-values independent, but is not what is said in the description).

To illustrate how these different assumptions are reflected in the data, we show two empirical distributions graphs. They are done for $n = 2^{24}$, $m = 2^{10}$ (as in Marsaglia’s code). We considered all the intervals (including the one that crosses the year boundary) in the same way, used the Poisson distribution with recommended parameter $m^3/(4n)$ and $chi^2$-test with combined bins in the both tails (until the expectation is at least 5). The size of sample that was compared with Poisson distribution was 500 (as in Marsaglia’s code); for this value bins 8…26 are counted (all numbers smaller than 8 are combined with 8, and all numbers greater than 26 are combined with 26). As one can see from the graphs (Figure 7.2), the assumption of uniform distribution of the $\chi^2$-test $p$-value seems very improbable (for both ways of counting the repetitions), so the validity of the test in the original form is questionable.

## 7.2 Kolmogorov-Smirnov 2-sample tools

As we have mentioned in Section 3.1, one can use Kolmogorov-Smirnov criterion for two samples to get a robust randomness test. Assume (null hypothesis) that $a_1, ..., a_m$ and $b_1, ..., b_n$ are real numbers independently sampled from some distribution (in our example: real functions of blocks of random bits). Let us sort the array $a_1, ..., a_m, b_1, ..., b_n$ and list all the values in the increasing order. Then we forget the values, and remember only where they come from (writing letter a or b). We get a sequence of $n+m$ letters that consists of $m$ letters a and $n$ letters b. The assumption (null hypothesis) implies that all strings of this type could appear with equal probabilities.

Therefore, if we have some subset of the set $S(m, n)$ of all strings of this type that has small cardinality (compared with the total number of all strings of this type, i.e., $\binom{n+m}{n} = (n+m)!/n!m!$), this set can be used as a test set.

<footnote n="3">³It is hard to say it for sure: the code for cutting the random stream in the blocks of requested size, bit.c, has 1441 lines. Still this is what is claimed in the comments.</footnote>
<figure><img src="image_49.png" /><figcaption>Figure 7.2: Empirical distributions for birthday test: 2¹⁰ birthdays, year length 2²⁴, distribution of the number of repeated intervals is compared with Poisson distribution (λ = 16) using samples of size 500. The p-values obtained by χ² test are accumulated from several thousands of trials and their empirical distribution function is shown. Two ways of counting repetitions are used: according to dieharder approach (interval length that repeated k > 2 times is counted as one repetition, left) and to the diehard approach (length that repeated k > 2 times is counted as k − 1 repetitions, right).</figcaption></figure>

Kolmogorov–Smirnov test considers a function on S(m, n) called maximal deviation. It is defined as follows: for each prefix of a string that consists of m letters a and n letters b we count the numbers m' and n' of letters a and b in this prefix, and then compute the deviation as

$$ \left| \frac{m'}{m} - \frac{n'}{n} \right| . $$

One can represent our string as a path from (0, 0) to (1, 1) where each letter a means going right by adding (1/m, 0), and each letter b corresponds to going up by adding (0, 1). Then the deviation is proportional to the maximal distance between the path and the diagonal (Figure 7.3).

For a given deviation value d, one can compute Kolmogorov–Smirnov p-value, i.e., the fraction of strings from S(m, n) that have the same or bigger deviation. The file ks2.mp does this using GNU multi-precision library using standard dynamic programming. The computation is exact (and later the answer is converted to double floats). The function

double ks2bar(int n0, int n1, long deviation_times_n0_n1)

gets two integers (number of both letters), and the maximal deviation multiplied by the product of these two integers (to make it integer) and returns the p-value as described above.

The code for computing the same value used in R statistical package (http://ftp.stjude.org/pub/software/JUMPm/labeled/R-3.1.0/src/library/stats/src/ks.c; the adapted version is in ksmirnov.c) uses the same technique but works with floating point numbers instead of multi-precision library, and includes also an asymptotic formula. One can compare the answer using the utility ks2compare:

```
$ ks2compare 1000 700 50000
Two samples two-side exact Kolmogorov-Smirnov p-value:
p-value for two samples (1000+700) and deviation 0.071429 [=50000/1000*700] is
0.028360425691996
For reference:
```
<figure><img src="image_50.png" /><figcaption>Figure 7.3: A path obtained while comparing two samples (produced by `ks2graph` utility). The deviation from the diagonal is quite visible, and this is confirmed by the numbers: characters appear 500 and 199 times, the deviation is 0.156, and the Kolmogorov-Smirnov *p*-value is 0.00173.</figcaption></figure>

asymptotic formula: 0.029940690
floating point iterations: 0.028360426

The words “two-side” indicate that we consider the deviation in both directions (with absolute values).

The experiments show that for reasonable lengths (several thousands) all the approaches give the answer fast (so it is not a bottleneck) and the floating point iterative computation gives small error. So the asymptotic formula is not really needed (this is good news since there is no easy error estimate) and one may use the computation with floating point numbers and check the answer with an exact computation if needed.

The utility `ks2val` takes a string (made of two different characters) from the standard input and computes the deviation and Kolmogorov-Smirnov *p*-value (using the exact computation):

```
$ echo 0100010111011 | ks2val
Profile:
  length: 13
  character [0] appears [6] times
  character [1] appears [7] times
Maximal deviation: 0.547619
KSvalue: 0.212121212121212
```

The utility `ks2vala` does the same thing as `ks2val` using floating point (long double) computations.

The utility `ks2draw` does the same thing as `ks2val` and also draws the path corresponding to the input string (see Figure 7.3 for example).

Note that for Kolmogorov-Smirnov comparison of two samples only the relative ordering between their elements is important. However, to understand better the nature of the difference between distributions it could be useful to draw the actual distribution functions for several samples. The utility `distribs` takes several files (with the names indicated in the command line); each file contains lines with real numbers considered as a sample. The utility draws (on the same graph) all the distribution functions for given samples (see Figure 7.4).

To obtain data for Kolmogorov-Smirnov analysis, one can use python script `ks2files.py`. It gets two file names from the command line. Each line in both files should start with a floating point
<figure><img src="image_51.png" /><figcaption>Figure 7.4: Utility *distrib* shows distribution functions for several samples obtained by computing spectral radii for random graphs obtained from Moonbase OneRNG generator in different regimes. The difference is clearly visible (and can be confirmed by computing p-values).</figcaption></figure>

real (that is used as an element of the sample); the rest of the line is used for tie breaking (no two lines should be exactly the same). The program sends to the standard output the sequence of characters 1 and 2 that indicate the origin of the data point after merging and sorting (whether it comes from the first or the second file).

Finally, Kolmogorov-Smirnov test can be used to test independence (as described in Section 3.5). This is done by utility *ks2indep*. This utility is applied to two files (the names are given in the command line). Each file is considered as a sequence of bytes, divided into groups of fixed length (in bytes, determined by option -b). Then each group of bytes is converted into a signed integer, interpreting them in the LSB (least significant byte first, option -l) or MSB (most significant byte first, option -m) way. These two sequence of integers (*n*₀, *n*₁, ...) and *m*₀, *m*₁, ...) are used in the following way: if *n*ᵢ satisfies [resp. violates] some condition (compiled-in list of condition is used; one of them is chosen by -c option; now two conditions are implemented: being positive [0] or being even [1]), the corresponding *m*ᵢ is sent into the first or the second subsequence). This is done until a given number of pairs (option -n, default 1000) are classified; then two subsequences are subjected to Kolmogorov-Smirnov test, the p-value is printed and the procedure is again applied to the rest of the files.

This procedure is a valid statistical test for the null hypothesis: both sequences are obtained by independent sampling from two independent sources. The distributions of the sources may be different, but both independence between different integers from the same source and independence between sources are important. Also we assume that there are no ties (if they happen, a warning is printed using *stderr*).

To reduce the independence between integers in one source one could skip several reading after each pair of integers used as described above. The skipped integers in the second source are also used for resolving ties. The number of skipped integers is given as -s option parameter (default 0).

A special regime is allowed with only one file (option -i). Then both integers are taken from the same file (first *n*₀, then *m*₀, then *n*₁, then *m*₁, etc.) and integers are skipped only after *m*₀, *m*₁, ... (skipped integers again are used for tie resolving).
Finally, the numbers of *p*-values in ten buckets [*0, 0.1*), ..., [*0.9, 1*) (ten integers) are printed (if -v verbose option is used).

## 7.3 Robust tests utility

We implemented several tests using the general scheme described in Chapter 3. There is a general driver (rtest.c, utility rtest) that calls different test functions depending on the options. It uses the following options:

* -t test_num the number of test function used (see the list of all test functions below);

* -f test_file specifies the name of the file that is tested. The file is understood as a sequence of 32-bit unsigned integers (four consecutive bytes form an integer) or a sequence of bytes (the interpretation is determined by a test function used);

* -e etalon_file specifies the name of the file that is used for comparison. If no -x option is used, then the test function is applied to data from test_file and etalon_file; in the robust mode (with -x) the test function is applied to data from test_file and to the xor of (new) data from test_file and data from etalon_file;

* -x requests the robust mode (see above for -e option);

* -p num_sample specifies the number of times the test function is applied to (fresh) data from tested file;

* -q num_sample specifies the number of times the test function is applied to etalon data;

* -d dimension indicates the number of values returned by each call to the test function (different test function have different dimensions); this number is the number of *p*-values returned by each Kolmogorov-Smirnov two-sample test;

* -n number determines how many bytes are consumed by each call of the test function; the exact interpretation of this number depends on the test function; some tests ignore this value and decide themselves how many bytes they need; for some other tests only some values of this number are allowed;

* -o directory_name asks to create files with the list of values that were submitted to the Kolmogorov-Smirnov two-sample test, and tells the name of the directory where these files should be created; the directory should not exist before the call. Every Kolmogorov-Smirnov value obtained by the test is accompanied by two files in this directory (that show both samples); one can understand the behavior of the test better when looking at these files.

* -r requires the test to be repeated; if a positive number is provided as an argument, it is interpreted as the number of repetitions; zero argument means that the test should be repeated as long as there are enough data in the input files.

* -k requires the exact computation with multi-precision integers for Kolmogorov-Smirnov two-sample *p*-value (instead of long double computation performed by default); can be useful for very large -p or -q options (exceeding 5000).

The program runs as long as there are data in the input files. After one test is performed and several Kolmogorov-Smirnov values are produced (the number of these values is the dimension of the test), the testing continues with the rest of the bits in both files.
Implementation remarks: the main driver file is rtest.c; the definitions for the test functions are given in test_func.c and additional files (see the macro TESTS_C in the Makefile). To add a new test function, one should put its code in a .c file, put its forward definition in test_func.h and add it to the list of all test functions in test_func.c.

We used hashes in a different way to simplify the code structure: the function returns not only some values (the number of returned values is the dimension), but also returns 64 bits (two 32-bit integers) that are used to resolve the ties: if the test function returns two identical values, these hash values are compared for sorting. Normally 64 bits should be enough to resolve the ties (if this does not happen, this means that the file we are testing is not random enough, so it creates a failure by itself).

The tools to compute the Kolmogorov-Smirnov p-values (for two samples; some tests also use calls to Kolmogorov-Smirnov distribution function for one sample) are taken from another directory (see KS_C macro definition in Makefile).

The following tests are currently implemented (listed according to their test numbers used in -t option; the C name of the corresponding function is shown after the number).

### 0 all_bytes

Dimension 1, ignores -n, counts the number of bytes that one needs to read to see all 256 possible bytes (coupon collector problem) if this is not too large (does not exceed some maximal number that is returned otherwise). Uses about 1.5 Kbyte per test value (this number should be multiplied by -p plus -q option parameters for the full test).

### 1 all_16

A similar test to the previous one but for 16-bit integers (pairs of bytes), so it requires much more data (1.5-2 Mbytes per test value).

### 3 sts_serial

The test depends on two parameters: the number of bits (n) analyzed and the maximal length of substrings analyzed (m). It produces an array of 3m values, i.e., 3 values for each integer from 1 to m. For each integer in this range there are two values called p₁ and p₂ in the NIST document; we report also one auxiliary “raw” value. The value of n is 32 times -n option; the value of m is taken as -d option divided by 3 (integer part). Requires n bits (n/8 bytes, 4 times -n option), the reasonable application is possible when n ≫ 2^m (otherwise the frequencies of factors do not make much sense).

Let X be the sequence of n bits to be analyzed. For each word w of length at most m we count the number count[w] of occurrences of w in X (considered as a cyclic word); the sum of all count[w] for all words w of given length equals n. The we compute $\psi^2[k] = (2^k/n) \sum_w (count[w]^2) - n$ where the sum is taken over all strings w of length k.

For our purposes the values $\psi^2[k]$ (raw values mentioned above) can be used directly, but NIST recommends to consider the first and second differences

$$ \Delta\psi[k] = \psi^2[k] - \psi^2[k-1], \\ \Delta_2\psi[k] = \psi^2[k] - 2\psi^2[k-1] + \psi^2[k-2]. $$

The dieharder code suggests to let $\psi^2[0] = 0$, so $\delta\psi[k]$ is defined for $k \ge 1$ and $\delta\psi_2[k]$ is defined for $k \ge 2$. (Non-existing values are filled by zeros in the test output).

These numbers are converted to p-values according to approximate distributions,

$$ p_1 = igamc(2^{k-2}, \Delta\psi[k]/2), \\ p_2 = igamc(2^{k-3}, \Delta_2\psi[k]/2). $$
### 4 opso

Original description: [31]. (The dieharder implementation deviates from the original description in a rather strange ways, see Section 7.1; we followed the original description). Test consumes $2^{23} + 1$ unsigned 32-bit integers and generates $32 - 10 + 1 = 23$ values according to 23 possible positions of 10-bit substrings in a 32-bit string. For each of 23 positions we get a word of length $n = 2^2 + 1$ in $2^{10}$-letter alphabet. We count 2-letter words that are not factors of this word, and convert this count to a (presumably uniform) $p$-value assuming normal distribution with some parameters (mean about 141909.33 and $\sigma$ about 290.46). The $-n$ parameter should be 2097153, the $-d$ parameter should be 23. Uses about 8 Mbytes per test value.

### 5 oqso

The test is similar to opso and is described in the same paper. (Again the dieharder implementation does something rather different.) The difference is that the alphabet contains $2^5$ letters (factors of length 5 at different positions in the input integers), and we count factors of length 4 that are missing. Consumes $2^{23} + 3$ integers and produces 28 values. Approximation mean and variance are slightly different (141909.60 and 294.656). The $-n$ parameter should be 2097155, the $-d$ parameter should be 28. Uses about 8 Mbytes per test value.

### 6 bytedistribs

A reimplementation of a test from dieharder, file dab_bytedistribs.c. From every three 32-bit numbers from the generator (dieharder allows also generators of different word size, but we assume the size is 32) one extracts 9 bytes (three from each integer: 8 least significant bits, 8 most significant bits and 8 bits in the middle). This is repeated $n$ times (specified by $-n$ parameter), so we get 9 distributions on $\{0, 1\}^8$ corresponding to 9 different byte positions. Then we apply $\chi^2$-square test to the combined distribution on $9 \cdot 256$ objects with $ndf = 9 \cdot 255$ degrees of freedom (since we have 255 degrees of freedom for each of 9 positions). One resulting value is returned (test dimension is 1). More precisely, each of $256 \cdot 9$ counters has expected value $n/256$ and actual value $x$; we compute the sum $S$ of $(x-e)^2/e$ for all counters and then compute the approximate $p$-value as gsl_sfgamma_inc_Q(ndf/2, S/2). Uses 12 times $[-n]$ bytes; the $-n$ parameter should be much bigger than 256 (otherwise the distributions do not make much sense).

### 7 knuth_runs

Test described in [25, p. 65] (the dieharder deviates from it since it considers the first element in a special way, but probably the results do not differ much). For a sequence of (arbitrary) length $n$ we count the number of parts in its minimal splitting into non-decreasing subsequences (runs). Then we consider 6-vector consisting of numbers $x_1, \dots, x_6$ of runs of lengths 1, 2, 3, 4, 5, $6-\infty$. Then (as explained by Knuth [25]) one computes

$$V = \frac{1}{n} \sum_{1 \le i,j \le 6} a_{ij}(x_i - nb_i)(x_j - nb_j),$$

where coefficients $a_{ij}$ and $b_i$ are given by tables provided by Knuth. Test dimension: 2 (for non-increasing and non-decreasing runs). Uses $[-n]$ option integers $(4 \times [-n]$ bytes).

### 8 osums

This is the “overlapping sum” test described in Marsaglia [29] and then included in dieharder [9] with come changes. Still after some investigation Robert Brown, the author of dieharder, came to the conclusion that this test should not be used (see Section 7.1). Still (like any other test function) it provides a correct test if used in two-sample scheme, so we included it. The test takes any number $m > 1$ as size ($-n$ parameter) and uses $2m - 1$ input 32-bit unsigned
integers (converted to [0, 1] double reals) to produce *m* overlapping sums (using *m* integers in each sum). Then some linear transformation (suggested by Marsaglia) is applied to get *m* presumably independent standard normal variables, they are converted to presumably uniform random variables in [0, 1] and Kolmogorov-Smirnov one-sample test is used. The resulting *p*-value is the only test output (dimension is 1). As we have discussed, it seems that for large *m* there is a significant deviation from the assumed distributions, and Kolmogorov-Smirnov one-sample tests give small values. So it may be wise to stick to the value *m* = 100 as used in original Marsaglia *diehard* test. Uses about 8 × [-n] bytes

9 ent_8_16

This simple test of dimension 4 takes arbitrary number *n* of 32-bit integer (specified in -n option), and returns four numbers (has dimension 4). It considers distributions on bytes and 16-bit integers. For each of these distributions the entropy is computed in a standard way, and also χ²-test is applied (with 255 and 65535 degrees of freedom for 8 and 16 bit versions). The four numbers returned are (entropy for 8 bits, χ²-test *p*-value for 8 bits, entropy for 16 bits, χ²-test *p*-value for 16 bits). Note that the values of entropy are less useful for the subsequent analysis since they are usually very close to 8 and 16, and measure the deviation from the uniform distribution in a similar way. Uses about 8 × [-n] bytes; so -n parameter should be much bigger than 64 (for bytes) and 32768 (for 16-bit integers).

10 fftest

This is a Fourier transform (spectral) test from NIST, as described in [37]. This description includes corrections from citekim-umeno-hasegawa. In this test *n* integers from the generator are converted into 32*n* bit sequence; this sequence is interpreted as array of ±1, and discrete Fourier transform is applied. We assume that *n* is a power of 2, so the simple algorithm for Fourier transform can be used (taken from gs1, GNU scientific library). In this way we get *n* complex numbers (since we start with real numbers, only *n*/2 first number are used). We see how many of them have absolute value that exceeds some threshold (namely, √ln(1/0.05)*n*), the description in NIST document claims that expected value of this number is 0.05*n*/2, and the distribution is close to binomial distribution for *n*/2 independent experiments with 0.05 probability. No derivation for this approximation is given (not to speak about error bounds). It seems that this approximation is quite rough. Figure 7.5 shows the experimental distributions for comparison of some random number generator and its xor with reference generator. One can see that both curves are quite close to each other and quite far from the diagonal (that should appear for uniform distributions). The jumps happen since there are only a few values of the test statistics (the number of values exceeding the thresholds). Let us stress again that the approximation errors do not destroy the validity of test results when using two-sample robust tests.

Test dimension is 1, uses 4 × [-n] bytes; the -n parameter should be a power of 2.

11 rank32x32

This and the following tests (both are implemented in test_rank.c) check ranks of some matrices. This test has dimension 1. It forms [-n option] matrices 32 × 32 (each uses 32 integers). Ranks of these matrices are computed and we count the frequencies of ranks 32, 31, 30, [29 and smaller]. (This grouping is reasonable since matrices with ranks less than 29 are rare.) Then the χ²-statistic is computed and converted to a (presumably) uniform distribution. The same computations are done in *dieharder* (while *diehard* used 31 × 31 matrices). Uses 128 × [-n] bytes.

12 rank6x8
<figure><img src="image_52.png" /><figcaption>Figure 7.5: Left: the distributions of presumably uniform random variables obtained by Fourier transform test, command `rtest -x -f a10bis/raw -e ../bitbabbler/raw_4.sources -p 1000 -q 1000 -n 256 -t 10 -d 1 -o outfft`. Right: a visible difference between test and etalon samples for some bad generator.</figcaption></figure>

This test has dimension 25 (since there are 32 - 8 + 1 = 25 possible position of a 8-bit factor in a 32-bit word). It takes 6 × [-n] integers; for each group of 6 integers it forms 25 matrices 6 × 8 corresponding to 25 positions. Ranks are computed, and for each position a distribution of matrices according to ranks (6, 5, [4 or less]) is analyzed via $\chi^2$ statistics and transformed into a presumably uniform variable (separately for each position). This follows diehard and dieharder with some corrections, see the discussion in Section 7.1. Uses 24 × [-n] bytes.

### 13 bitstream_o

This and the next tests (overlapped and non-overlapped version of a bitstream test from `diehard` and `dieharder`) are implemented in file `test_bitstream.c`. This test considers the input stream of bits (32 integers are concatenated as written from left to right; the bits are also written from left to right). The test uses $2^{21} + 19$ bits (and discards unused bits in the last integer). It considers all 20-bit factors in this long string and counts the missing ones, i.e., 20-bit strings that do not appear as a factor. This number is transformed into a presumably uniformly distributed variable (following `dieharder`) using empirically observed normal distribution with mean 141909 and $\sigma = 428$. The dimension of the test in 1, the `-n` option parameter should be 2097171 (number of bits used), or about 256 Kbytes per test value.

### 14 bitstream_n

Similar to the preceding test, but non-overlapping 20-factors are used (the input bit stream is cut into 20-bit pieces). Uses $2^{21} * 20 = 41943040$ bits, and this number should be the value of `-n` parameter; dimension is 1. Since the words are now independent, the distribution is different, and the `dieharder` uses the same mean 141909 with smaller $\sigma = 290$. About 5 Mbytes per test value

### 15 lz_split

Lempel-Ziv test as it was described in the original version of [37]. The sequence of bits is split into words according to the following Lempel-Ziv rule: the next word is chosen as the minimal prefix of the non-processed yet part that is not a prefix of already chosen word. The algorithm counts the number of full words produced when analyzing 32 × [-n parameter]
<figure><img src="image_53.png" /></figure>

Figure 7.6: Empirical distributions. Each experiment involved 3000 Lempel–Ziv decompositions of a bit sequence of length $10^6$; several experiments are shown with different colors. According to the original version of the test, the empirical distribution functions should be close to a diagonal line. One can see both effects mentioned in [23]: first, the number of different values is small, and this creates the “ladder” structure; second, one can also suspect systematic deviations from the straight line (added for comparison).

input bits. The NIST standard suggest some approximation for the resulting distribution in case of $10^6$ bits (-n 31250), so in this case we convert the values to presumably uniform distribution using this approximation. As we have mentioned in Section 2.3.3, the problems with this approximation were found in [23]. They can be illustrated by a distribution obtained when applying this test (with NIST parameters) may times to several hardware generators (Figure 7.6). Still it can be safely used in the two-sample robust version.

## 16 birthdays

A “year” consists of $2^{24}$ days; we choose 1024 random days (“birthdays”) and compute 1024 intervals between birthdays. Some intervals may coincide, and we count the number of coincidences (i.e., the difference between 1024 and the number of different intervals). We repeat this experiment several times (-n option) and compare the resulting distribution with presumed approximation (Poisson distribution with $\lambda = 16$), getting one *p*-value that is considered as the output of the test. More precisely, we use first 24 bits from the 32-bit integers from the input sequence; this gives us one *p*-value, than next 24 bits from the same 32-bit integers, and so on (32 cyclic shifts, as described — but not implemented — in *dieharder*). The test dimension is 32, and it uses 1024 times [-n option] integers ($4096 \times [-n]$ bytes) for each set of 32 values. Note that the test is robust (as any other two-sample test) but its usefulness is not clear (due to the systematic deviation of the observed distribution from the presumed one, see Section 7.1, Figure 7.2).

There are also some shell scripts that run these tests for suitable parameters (adjusted for 1 Mb,
10 Mb, 100 Mb, 1 Gb, 10 Gb files). Some tests take a long time, so it may have sense to delete some tests from these files if they take too much time. Two arguments are required: the name of the tested file and the name of the reference file.

A short reminder to be taken into account when interpreting the test values:

* The *p*-values produced by these files are reliable (not based on any unproven assumptions about the tests and/or approximations<fnref n="4" />). For any given test this means that for every *p* > 0 and for every reference file the probability of the event “the test with random bits in the tested file will produce value at most *p*” is at most *p*. Note also that the tests never “rewind” the file (using the same bits several times): if there is not enough bits in the file, no value is produced.

* There is no claim of uniform distributions for the test output (for a random bits input), so output values close to 1 do not indicate any problems.

* The tests (some individual tests and the shell scripts) produce many numbers at once. These number are not independent. So one should *not* multiply them to get a more impressive number. Moreover, one should apply the Bonferroni correction (that makes the conclusion less impressive): If for a chosen family of tests one of *m* number (where *m* is the total number of outputs, depending on the chosen family) is *p*, the resulting *p*-value is *mp*. For example, if one of the 40 test values is less than 0.01, one should not become suspicious about the generator, since 40 · 0.01 is 0.4, definitely not enough to start worrying.

* On the other hand, if you repeat some test several times (say, using -r option, or just applying it to different runs of the same random bit generator), the output values *are* independent. This can be used to amplify the non-randomness conclusions. For example, if a test (or a group of tests) fixed before the experiment produces 37 values, we made 15 runs with the same generator (and the same or different reference files, fixed in advance), and it turned out that one of 37 values is less than 0.001 for all the experiments, the correctly computed *p*-values is 37 · (0.001)¹⁵. One should distinguish this from the other case: for each of 15 experiments there is some output value less than 0.001, but this value can be produced by different tests in different experiments. In this case the correct computation of *p*-value gives (37 · 0.001)¹⁵.

* The recommended approach when discrediting some generator: make some preliminary experiments to determine which tests look more suspicious. Then try the selected test on the fresh input from the generator and report the *p*-value obtained for this test. But be careful if a small *p*-value appears after several attempts to discredit the generator: a honest researcher should multiply the resulting *p*-value by the number of attempts.

## 7.4 Spectral tests

The spectral tests currently involve a lot of parameters, so they are not yet included in the robust testing suite. Currently we do *not* have examples where they beat all the other tests in the sense that they find some deficiency in a generator that cannot be found by other tests. Still it would be desirable to have more robust tests especially if they use different approaches.

The call

`example -g gap -e err -s size -d deg -o iter -n number -l -u -h] filename`

<footnote n="4">⁴The only compromise here happens when we use floating point numbers for computing the Kolmogorov-Smirnov value; looking at the output of `ks2compare` utility, we may hope that this does not create any problems for samples of size up to 5000, but one should be careful before pronouncing the final verdict and use option `-k` that asks for an exact computation.</footnote>
The utility example analyses one file with name filename (so this utility should be combined
with the xoring made separately to get a robust test). It uses bytes from this file to construct a
sequence of graphs. The size of the graph is determined by option parameter -s size; it should be
a power of 2 not exceeding 65536. Each (undirected) graph is constructed by emitting deg random
edges from each vertex. (Note that the resulting undirected graph may not have constant degree,
but the average degree is close to 2 · deg.) Then the algorithm from [43] is used to compute the
second eigenvalue for the resulting graph. This is an iterative algorithm that has some parameters.
It constructs some lower and upper bounds for the second eigenvalue and stops when the difference
between them becomes smaller than gap (-g parameter). Also it uses some internal randomness
(now crudely replaced by some fixed pseudo-random generator<fnref n="5" />. Also the algorithm uses some
upper bound for the number of iterations that can be specified by option -o. The option -n parameter
specifies the maximal number of graphs that should be constructed (if there is enough bits in the
input file filename).

The output of the utility is a sequence of lines. Each of them corresponds to one graph, and the
contents of the lines is determined by the options -l, -u, -h (at least one of them should be present).
The output line contains

The output of the utility is a sequence of lines. Each of them corresponds to one graph, and the
contents of the lines is determined by the options -l, -u, -h (at least one of them should be present).
The output line contains

- the lower bound for the second eigenvalue, if -l is present;

- the upper bound for the second eigenvalue, if -l is present;

- some hash value (useful for tie breaking), if -h is present.

Also an option -v (verbose) can be used to provide more information.

7.5 Other tools

There are few simple tools that visualize data from .wav and raw data files and perform conversion
between formats.

* `readwav [-m|-s|-y] filename.wav` shows data (signal readings interpreted as signed integers) from a .wav file. For a stereo file [-s] it takes values from both channels and draws a pair. The option [-y] does the same but corrects the 1-sample lag in the right channel (that happens in Yamaha MW10 and Behringer 1204USB mixers we used); mono signal has identical readings in both channels after the correction. The [-m] option splits the movo .wav file into pairs of sample values (signed integers) and draws the corresponding points. The number of bits per value should be 8, 16 or 24.

* `getdatawav [-m|-s|-y] source.wav dest1 [dest2]` converts standard .wav file into raw format, deleting the header and splitting stereo data into two separate files. Options -s and -y require two destination files (dest1 is for the left channel, dest2 is for the right channel). Option -m is used for mono files, and just deletes the header. The input file should have resolution 8, 16, or 24. The order of bytes remains unchanged.

* `drawpair -b num_bytes -l|-m -c` draws pairs obtained from two files (in parallel) or from one file (-c; two consecutive data form a pair). The `num_bytes` parameter gives the number of bytes interpreted as a signed integer; `-l/-m` option (one of two should be present) says whether the integers are assumed to appear in least or most significant byte first.

<footnote n="5">⁵Since our xor-trick guaranteed correctness of *p*-values independently of the function used for testing, we decided to use some simple pseudo-random sequence with fixed seed.</footnote>
## 7.6 Two-source extractor tools

We implemented the two-source extractor described in Section 5.5. We used the finite field of size $2^{1024}$ whose elements are represented as polynomials in $\mathbb{F}_2[x]$ considered modulo

$$p(x) = x^{1024} + x^{19} + x^6 + x + 1.$$

This polynomial is suggested as low-weight irreducible polynomial in [51]. “Low weight” means that only a few coefficients are non-zero; this allows us to simplify the division modulo this polynomial. We used the natural basis

$$1, x, x^2, \dots, x^{1023}$$

in the field (the quotient of $\mathbb{F}_2[x]$ modulo $p(x)$); the multiplication by the basis elements can be obtained by repeated multiplications by $x$. And a multiplication by $x$ is just a shift in the array of coefficients, followed by a reduction modulo $p(x)$. This reduction is applied to a polynomial of degree 1024; we have to change coefficients in the degrees involved (20, 7, 2, 1). This simple operations are implemented in file `extract.c`; we produce a bit string of length 128 from two strings of length 1024 (the field elements).

The file `extractor.c` calls this extraction function. Since in our experiments we used `.wav` files where most significant bits (of a 16-bit words from analog-digital converter are often zeros), we decided to use `xor` of two halves of an 16-bit word to get 8 bits that have (hopefully) reasonably high min-entropy. The command

```
extractor file1 file2
```

applies this operation (`xoring` of two bytes in a 16-bit integer) to each of files `file1`, `file2`, applies the extractor function as described, and send the resulting bits to standard output.

In this way from 2048 + 2048 bits from both files we generate 1024 + 1024 bits (by `xoring`), and then generate 128 output bits, so the production rate is 1 : 32 (one output bit for 32 input bits).
# Chapter 8

## Conclusions

Let us summarize some observations made in this report.

* The main goal of this project was to review the current situation in hardware random number generation (devices, tests, standards) taking into account the corresponding theory (algorithmic information theory, pseudorandom generators, extractors). It turns out to be quite unsatisfactory: the basic principles are not always understood by the authors of the standards and developers of tests, as well as the device makers. As a result, tests are unreliable in many ways and testing results and claims made by the manufacturers are often meaningless.

* The reliable test methodology can be developed based on two-sample Kolmogorov–Smirnov tests. The general scheme is easy to implement. We adapted a dozen of often used tests into this framework and showed that this methodology works. Still one should add more tests to it, so that existing tests (and new ones that are developed by different authors) could be used in a reliable way to produce real p-values, not just some numbers with no guarantees.

* There are several physical processes that can be used for random number generators. The devices we tested used Zener diodes (most of them), analog multipliers (Infinite Noise, Bit-babblers) and photon reflection (Quantis). None of them give a bit sequence that can pass the basic tests without some “conditioning”, but all three approaches seem to be suitable for high quality generators after some post-processing, and there are no obvious reasons to prefer expensive “quantum randomness” (whatever it is).

* One can use postprocessing based on functions used in cryptography. There are no theoretical reasons or assumptions that guarantee anything about this type of post-processing, but the currently known tests do not show obvious weaknesses in this approach. Still the possibility of some hidden trapdoor remains.

* Instead, one can use much simpler methods based on xoring several random streams (generated by different sources or in different time) or using other linear transformations of bits. Known theoretical results about that are not really applicable to practical situations, but this does not prevent these methods to work quite well in practice. Simple postprocessing of this type often gives sequences that do not exhibit suspicious behavior with known tests.

* There are also more complicated (and not yet used in practical generators) tools for randomness extraction based on two-source extractors. They are significantly slower than xoring, but still are quite practical and can be easily implemented in software.

* All these approaches can easily provide streams of random bits with rates 100–1000 kilobits per second using easily available devices (few hundred euros) and some additional postprocessing. The choice of the device and postprocessing methods depends on the speed requirements and
the paranoia level. Most paranoid people can achieve rates about 10–100 kilobits per second without using any commercial random bits generator, just by combining zener diodes and audio digitizing equipment.

* One can obtain rates 10–100 megabits per second or higher using several devices in parallel or more expensive (and sometimes less transparent) devices. Faster random bit production would require serious parallelism or much more expensive devices.

## Acknowledgments

The authors are grateful to their colleagues who worked on the RaCAF project (Laurent Bienvenu, Bruno Durand, Fabien Givors, Ruslan Ishkuvatov, Gregory Lafitte) and to the members of the ES-CAPE team (LIRMM, Montpellier), Kolmogorov Seminar (Moscow), Theoretical Computer Science Lab (Moscow, HSE), and all others from whom we learned about randomness, including Vladimir Uspensky, Leonid Levin, Alexander Zvonkin, Nikolay Vereshchagin, Vladimir Vovk, Vladimir Vyugin, Péter Gács, Wolfgang Merkle, Paul Vitányi, Daniil Musatov, Andrei Rumyantsev, Mikhail Andreev, Gleb Novikov, Bruno Bauwens, Konstantin Makarychev, Yury Makarychev, Ilya Razensteyn, Gleb Posobin, Alexey Vinogradov. The work was supported by ANR RaCAF grant ANR-15-CE40-0016-0. We are also grateful to the device makers who kindly agreed to answer some of our questions.
Bibliography

[1] ASA Publication: Ronald L. Wasserstein, Nicole A. Lazar, Editorial: The ASA's Statement on p-Values: Context, Process, and Purpose, *The American Statistician*, 70(2), 129-133 (2016), http://dx.doi.org/10.1080/00031305.2016.1154108

[2] Andrian Belinski, *Method and apparatus of entropy source with multiple hardware random noise sources and continuous self-diagnostic logic*, US9477443B1 patent, https://patents.google.com/patent/US9477443.

[3] Daniel J. Benjamin et al. (72 authors), Redefine statistical significance, *Nature human behaviour*, 2, 6–10 (2018)

[4] Laurent Bienvenu, Péter Gács, Mathieu Hoyrup, Cristobal Rojas, Alexander Shen, Algorithmic tests and randomness with respect to a class of measures, *Proceedings of the Steklov Institute of Mathematics*, 274, 34–89 (2011). See also http://arxiv.org/abs/1103.1529.

[5] Lenore Blum, Manuel Blum, Mike Shub, A Simple Unpredictable Pseudo-Random Number Generator, *SIAM Journal on Computing*, 15(2), 364–383 (1986), https://doi.org/10.1137/0215025

[6] Manuel Blum, Silvio Micali, How to Generate Cryptographically Strong Sequences of Random Bits, *SIAM Journal on Computing*, 13(4), 850–864 (1984), https://doi.org/10.1137/0213053 (preliminary version was presented at FOCS 1982 conference).

[7] Émile Borel, *Le Hasard*, Paris, Librarire Félix Alcan, 1920. 312 pp.

[8] Jan Bouda, Matej Pivoluska, Martin Plesch, Improving the Hadamard extractor, *Theoretical Computer Science*, 459 (2012), 69–76, http://dx.doi.org/10.1016/j.tcs.2012.07.030.

[9] Robert G. Brown, *Die Harder: A Gnu Public License Random Number Generator*, version 3.31.1, http://www.phy.duke.edu/~rgb/General/dieharder.php (2006–2018).

[10] A. Philip David, Steven de Rooij, Glenn Shafer, Alexander Shen, Nikolai Vereshchagin, Vladimir Vovk, Insuring against loss of evidence in game-theoretic probability, *Statistics and Probability Letters*, 81, 157–162 (2011), https://doi.org/10.1016/j.spl.2010.10.013

[11] Robert Davies, *Hardware random number generators*. Presented at 15th Australian Statistics Conference, July 2000, and 51st conference of New Zealand Statistical Association, September 2000, http://robertnz.net/hwrng.htm.

[12] Yevgeniy Dodis, Ariel Elbaz, Roberto Oliveira, Ran Raz, Improved Randomness Extraction from Two Independent Sources. In: *Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques. RANDOM 2004, APPROX 2004*. Lecture Notes in Computer Science, 3122, Springer, Berlin, Heidelberg. https://doi.org/10.1007/978-3-540-27821-4_30.

[13] Daniela Frauchiger, Renato Renner, Matthias Troyer, *True randomness from realistic quantum devices*, arxiv preprint, https://arxiv.org/pdf/1311.4547.pdf
[14] Yuri Gurevich, Grant Olney Passmore, Impugning Randomness, Convincingly, *Studia Logica*, **100**(1-2), 193-222 (April 2012), https://link.springer.com/article/10.1007/s11225-012-9375-1, see also https://arxiv.org/pdf/1601.00665.pdf, https://www.cl.cam.ac.uk/~gp351/Gurevich-Passmore-IRC.pdf.

[15] Yuri Gurevich, Vladimir Vovk. *Test statistics and p-values*. Working paper #16, On-line compression modelling project (new series), http://www.alrw.net/articles/16.pdf. See also https://arxiv.org/pdf/1702.02590.pdf.

[16] Johan Håstad, Russell Impagliazzo, Leonid A. Levin, Michael Luby, A Pseudorandom Generator from any One-way Function, *SIAM Journal on Computing*, **28**(4), 1364-1396 (1999, original conference papers from 1989 and 1990), https://doi.org/10.1137/S0097539793244708.

[17] Miguel Herrero-Collantes, Juan Carlos Garcia-Escartin, Quantum Random Number Generators, *Review of Modern Physics*, **89**, 015004 – Published 22 February 2017, https://journalsaps.org/rmp/pdf/10.1103/RevModPhys.89.015004. See also https://arxiv.org/abs/1604.03304v2.

[18] S. Hoory, N. Linial, A. Wigderson, Expander graphs and their applications, *Bulletin of the American Mathematical Society*, **43**(4), 439-561 (2006)

[19] Darren Hurley-Smith, Julio Hernandez-Castro, *Quam Bene Non Quantum: Bias in a Family of Quantum Random Number Generators*, https://eprint.iacr.org/2017/842.

[20] Darren Hurley-Smith, Julio Hernandez-Castro, *Quantum Leap and Crash: Searching and Finding Bias in Quantum Random Number Generators*, *ACM Transactions on Privacy and Security*, **23**(3). 1-25 (2020), ISSN 2471-2566, see also https://kar.kent.ac.uk/81957/11/Quantum_Leap_TOPS_Submission_FINAL.pdf

[21] David Johnston, *Random Number Generators — Principles and Practices. A Guide for Engineers and Programmers*, Walter de Gruyter, Berlin/Boston, ISBN 978-1-5015-1513-2, https://doi.org/10.1515/9781501506062-201, 436 pp.

[22] Wolfgang Killmann, Werner Schindler, *A proposal for: Functionality classes for random number generators*, Version 2.0, 18 September 2011, AIS20/AIS31, https://www.bsi.bund.de/SharedDocs/Downloads/DE/BSI/Zertifizierung/Interpretationen/AIS_20_Functionality_classes_for_random_number_generators_e.pdf

[23] Song-Yu Kim, Ken Umeno, Akio Hasegava, *Corrections of the NIST Statistical Test Suite for Randomness*, https://eprint.iacr.org/2004/018.pdf.

[24] Alexander Kireev, *On the falsified results of the “referendum” in Sevastopol* [In Russian], blog post in LiveJournal, https://kireev.livejournal.com/1095568.html.

[25] Donald Knuth, *The Art of Computer Programming. Volume 2. Seminumerical Algorithms*. Second edition, Reading et al., Addison-Wesley, 1981. ISBN 0-201-03822-6.

[26] Alexey Kupriyanov, Gauss against Churov: preliminary conclusions [In Russian], *Troitskii variant* newspaper, 08.05.2019, https://trv-science.ru/2018/05/08/gauss-protiv-churovapromezhutochnyj-itog/.

[27] Patrick Lacharme, Post-Processing Functions for a Biased Physical Random Number Generator. In: *Fast Software Encryption. FSE 2008*. Lecture Notes in Computer Science, **5086**. Springer, Berlin, Heidelberg. https://doi.org/10.1007/978-3-540-71039-4_21
[28] George Marsaglia, A Current View of Random Number Generators, *Computer Science and Statistics*, Sixteenth Symposium on the Interface, Elsevier, North-Holland (1985), 3-10.

[29] George Marsaglia, *Random Numbers CDROM including the Diehard Battery of Tests of Randomness*, 1995, was available at http://stat.fsu.edu/pub/diehard/; now (2019) still available as snapshots from https://web.archive.org. Contains the preprint version of [31, 28].

[30] George Marsaglia, Wai Wan Tsang, Some difficult-to-pass tests of randomness, *Journal of Statistical Software*, 7(3), 2002, https://www.jstatsoft.org/article/view/v007i03.

[31] George Marsaglia, Arif Zaman, Monkey Tests for Random Number Generators, *Computers and Mathematics with Applications*, 26(9), 1-10 (November 1993).

[32] George Marsaglia, Xorshift RNGs, *Journal of Statistical Software*, 8, issue 14 (July 2003), http://dx.doi.org/10.18637/jss.v008.i14

[33] George Marsaglia, https://github.com/omec-project/il_trafficgen/blob/master/pktgen/app/xorshift64star.h

[34] Ueli M. Maurer, A Universal Statistical Test for Random Bit Generators, *Journal of Cryptology*, 5(2), 89-105 (1992), https://link.springer.com/article/10.1007/BF00193563.

[35] M. Matsumoto, T. Nishimura, Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator, *ACM Transactions on Modeling and Computer Simulation*, 8 (1): 3-30 (1998)

[36] FIPS PUB 140-2. *Security requirements for cryptographic modules*. Federal information processing standards publication (Supercedes FIPS PUB 140-1, 1994 January 11). Issued May 25, 2001. Information Technology Laboratory, National Institute of Standards and Technology, http://csrc.nist.gov/publications/fips/fips140-2/fips1402.pdf

[37] NIST Special Publication 800-22: Andrew Rukhin, Juan Soto, James Nechvatal, Miles Smid, Elaine Barker, Stefan Leigh, Mark Levenson, Mark Vangel, David Banks, Alan Heckert, James Dray, San Vo, *A Statistical Test Suite for Random and Pseudorandom Number Generators for Cryptographic Applications*, Revised: April 2010, Lawrence E. Bassham III, National Institute of Standards and Technology, Technology Administration, U.S. Department of Commerce (NIST), Revision 1a (2010), 131 pp., https://www.nist.gov/publications/statistical-test-suite-random-and-pseudorandom-number-generators-cryptographic. Previous version seems to be unavailable at this site, but the review of Elaine B. Barker, *ITL Bulletin, December 2000, 3 pp.* is available at https://tsapps.nist.gov/publication/get_pdf.cfm?pub_id=151231. The Lempel-Ziv test, criticised in [23], was there (#10) according to the review; it is missing in the updated version.

[38] NIST Special Publication 800-90A: Elaine Barker, John Kelsey, *Recommendation for Random Number Generation Using Deterministic Random Bit Generators*, National Institute of Standards and Technology, Technology Administration, U.S. Department of Commerce (NIST), June 2015, https://csrc.nist.gov/publications/detail/sp/800-90a/rev-1/final (previous version: January 2012).

[39] NIST Special Publication 800-90B: Meltem Sönmez Turan, Elaine Barker, John Kelsey, Kerry McKay, Mary Baish, Michael Boyle, *Recommendation for the Entropy Sources Used for Random Bit Generation*, National Institute of Standards and Technology, Technology Administration, U.S. Department of Commerce (NIST), January 2018, https://csrc.nist.gov/publications/detail/sp/800-90b/final.
[40] NIST Special Publication 800-90C (Second Draft): Elaine Barker, John Kelsey, *Recommendation for Random Bit Generator (RBG) Constructions*, April 2016, https://csrc.nist.gov/CSRC/media/Publications/sp/800-90c/draft/documents/sp800_90c_second_draft.pdf.

[41] NIST Special Publication 800-107: Quynh Dang, Recommendation for Applications Using Approved Hash Algorithms Revision 1 (August 2012), https://nvlpubs.nist.gov/nistpubs/Legacy/SP/nistspecialpublication800-107r1.pdf.

[42] Gleb Novikov, Randomness deficiencies, *Computability in Europe, 2017: Unveiling Dynamics and Complexity*, Lecture Notes in Computer Science, **10307**, Springer, 338-350, https://doi.org/10.1007/978-3-319-58741-7_32.

[43] Yann Ollivier, *Spectral gap of a graph*, a public domain code for finding the spectral gap of a graph, http://www.yann-ollivier.org/specgraph/specgraph (2004)

[44] Yuval Peres. Iterative von Neumann's procedure for extracting random bits, *The Annals of Statistics*, **1**(1), 590–597 (March 1992), https://www.jstor.org/stable/2242181

[45] Melissa E. O'Neill, Melissa, (5 September 2014). PCG: *A Family of Simple Fast Space-Efficient Statistically Good Algorithms for Random Number Generation* (Technical report). Harvey Mudd College. HMC-CS-2014-0905. See also https://www.pcg-random.org/download.html

[46] Takuji Nishimura and Makoto Matsumoto. A C-program for MT19937, with initialization improved 2002/1/26. http://www.math.sci.hiroshima-u.ac.jp/m-mat/MT/MT2002/CODES/mt19937ar.c

[47] R, a freely available language and environment for statistical computing, see *The Comprehensive R Archive Network*, cran-r.project.prg (January 18, 2020, version R-4.0.3), file R-4.0.3/src/library/stats/src/ks.c, functions psmirnov2x and pkstwo.

[48] RAND corporation, *A Million Random Digits with 100,000 Normal Deviates*, Free press, 1955. (Reissued in 2001 as ISBN 0-8330-3047-7.)

[49] Omer Reingold, Salil Vadhan, Avi Wigderson, *A Note on Extracting Randomness from Santha-Vazirani Sources*, report available from Reingold, https://omereingold.files.wordpress.com/2014/10/svsources.pdf.

[50] Miklos Santha, Umesh V. Vazirani, Generating Quasi-random Sequences from Semi-random sources, *Journal of Computer and System Sciences*, **33**, 75–87 (1986), https://doi.org/10.1016/0022-0000(86)90044-9.

[51] Gadiel Seroussi, *Table of Low-Weight Binary Irreducible Polynomials*, Hewlett-Packard Computer Science Laboratory technical report, HPL-98-135 (August 1998), https://www.hpl.hp.com/techreports/98/HPL-98-135.pdf

[52] Alexander Shen, Around Kolmogorov complexity: Basic Notions and Results, *Measures of Complexity*. Festschrift for Alexey Chervonenkis, Springer (2015), 75–116, https://link.springer.com/chapter/10.1007/978-3-319-21852-6_7, see also https://arxiv.org/abs/1504.04955v1.

[53] Alexander Shen, *Election and statistics: the case of “United Russia”*, 2009–2018, preprint, https://arxiv.org/abs/1204.0307.

[54] Alexander Shen, *Making randomness tests more robust*, HAL archive, 2018, https://hal.archives-ouvertes.fr/hal-01707610.
[55] Alexander Shen, *Randomness Tests: Theory and Practice*. In: Fields of Logic and Computation III, Lecture Notes in Computer Science, 12180, Springer-Verlag, 258–290 https://link.springer.com/chapter/10.1007%2F978-3-030-48006-6_18 (2020). See also: https://hal-lirmm.ccsd.cnrs.fr/lirmm-03065320

[56] Alexander Shen, Vladimir A. Uspensky, Nikolai K. Vereshchagin, *Kolmogorov Complexity and Algorithmic Randomness*, American Mathematical Society (2017), http://www.lirmm.fr/~ashen/kolmbook-eng-scan.pdf.

[57] Richard Simard, Pierre L'Ecuyer, Computing the Two-Sided Kolmogorov-Smirnov Distribution, *Journal of Statistical Software*, 39(11), 18 pp. (2011), available at https://www.jstatsoft.org/index.

[58] Tom Stoppard, *Rosencrantz and Guildenstern Are Dead*, a play (1966), ISBN978-0-8021-3275-8 (1971 edition).

[59] TrueRNG v3 documentation from Ubld.it, http://ubld.it/truerng_v3.

[60] Nikolay Vereshchagin, Alexander Shen, Algorithmic Statistics Revisited, *Measures of Complexity*. Festschrift for Alexey Chervonenkis, Springer (2015), 235–252, https://link.springer.com/chapter/10.1007/978-3-319-21852-6_17, see also https://arxiv.org/abs/1504.04950v2.

[61] Nikolay Vereshchagin, Alexander Shen, Algorithmic statistics: forty years later, *Computability and Complexity*. Essays Dedicated to Rodney G. Downey on the Occasion of His 60th Birthday. Springer (2017), Lecture Notes in Computer Science, **10010**, 669–737, see also https://arxiv.org/abs/1607.08077.

[62] Andrew C. Yao, Theory and application of trapdoor functions, *23rd Annual Symposium on Foundations of Computer Science* (FOCS), 80–91 (1982), http://ieeexplore.ieee.org/document/4568378/.