On testing pseudorandom generators via statistical tests based on the arcsine law

TL;DR

This paper introduces a new statistical test based on the arcsine law for evaluating pseudorandom number generators, including an error analysis using a Berry-Essen type inequality to ensure test reliability.

Contribution

It proposes a second-level statistical test for pseudorandom generators based on the arcsine law, with a detailed error approximation framework.

Findings

The test effectively detects weaknesses in pseudorandom generators.

The Berry-Essen inequality provides accurate error bounds for the test.

The method enhances confidence in the statistical evaluation of randomness.

Abstract

Testing the quality of pseudorandom number generators is an important issue. Security requirements become more and more demanding, weaknesses in this matter are simply not acceptable. There is a need for an in-depth analysis of statistical tests -- one has to be sure that rejecting/accepting a generator as good is not a result of errors in computations or approximations. In this paper we propose a second level statistical test based on the arcsine law for random walks. We provide a Berry-Essen type inequality for approximating the arcsine distribution, what allows us to perform a detailed error analysis of the proposed test.

Tables

Table 1. Table 1 : Results of the ASIN test for several generators with parameters m = 10000 , n = 2 34 , s = 40 formulae-sequence 𝑚 10000 formulae-sequence 𝑛 superscript 2 34 𝑠 40 m=10000,n=2^{34},s=40 . Bolded values indicate that we reject ℋ 0 subscript ℋ 0 \mathcal{H}_{0} at the significance level α = 0.0001 𝛼 0.0001 \alpha=0.0001 .

	$d_{tv}^{asin}$	$p_{{\chi }^2}$
MS Visual C++	0.2093	0.0000
GNU C	0.0255	0.2389
Minstd 48271	0.2089	0.0000
MT19937-64	0.0252	0.2523

Table 2. Table 2: Results ( p 𝑝 p -values) of T a ​ s ​ i ​ n superscript 𝑇 𝑎 𝑠 𝑖 𝑛 T^{asin} and swalk_RandomWalk1 statistics from TestU01 for n = 2 26 𝑛 superscript 2 26 n=2^{26} (for parameters [ TestU01 notation]: N = 1 , n = 10000 , r = 0 , s = 32 , L ​ 0 = L ​ 1 = 2 26 formulae-sequence 𝑁 1 formulae-sequence 𝑛 10000 formulae-sequence 𝑟 0 formulae-sequence 𝑠 32 𝐿 0 𝐿 1 superscript 2 26 N=1,n=10000,r=0,s=32,L0=L1=2^{26} ). Bolded values indicate that we reject ℋ 0 subscript ℋ 0 \mathcal{H}_{0} at the significance level α = 0.0001 𝛼 0.0001 \alpha=0.0001 .

PRNG\Test	$p_{{\chi }^2}$	Statistic H	Statistic M	Statistic J	Statistic R	Statistic C
MS Visual C++	0.0148	0.2700	0.0900	0.6300	0.4200	0.8000
GNU C	0.4731	0.1600	0.9800	0.1100	0.1900	0.4900
Minstd 48271	0.0115	0.0090	0.1400	0.4900	0.0700	0.0044
MT19937-64	0.2548	0.0800	0.1000	0.4200	0.9700	0.3500
${\text{Flawed}}_{\mathrm{MT19937}-64,26,1/66}$	0.0000	0.0000	0.2200	0.0000	0.3900	0.3800

Table 3. Table 3 : Tables present which tests of the BigCrush were not passed. The fist column t ​ n ​ o 𝑡 𝑛 𝑜 tno is the test number, for the p − v ​ a ​ l ​ u ​ e 𝑝 𝑣 𝑎 𝑙 𝑢 𝑒 p-value , ϵ italic-ϵ \epsilon is value that is ϵ < 1.0 ​ e − 15 italic-ϵ 1.0 𝑒 15 \epsilon<1.0e-15 .

tno	test name	parameters	p-value
74	RandomWalk1 R	$L=50,r=0$	$6.1e-4$
80	LinearComp	$r=0$	$1-ϵ$
81	LinearComp	$r=29$	$1-ϵ$
Mersenne Twister

Equations

$$X_i=2B_i-1,S_0=0,S_k=\sum _{i=1}^k{X}_i,k=1,\dots$$

$$\begin{array}[]{llll}&\displaystyle\mathbb{P}\left(\liminf_{n\to\infty}{S_{n}\over\sqrt{2n\log\log n}}=-1\right)&=1,\\ &\displaystyle\mathbb{P}\left(\limsup_{n\to\infty}{S_{n}\over\sqrt{2n\log\log n}}=+1\right)&=1.\end{array}$$

$$D_k=\text{mathbbm}1\left(S_k>0\vee S_{k-1}>0\right),k=1,2,\dots ,$$

$$\begin{array}{rl}\mathbb{P}\left(L_n\le x\cdot n\right)& \underset{n\to \infty }{\to }\frac{1}{\pi }{\int }_0^x\frac{dt}{\sqrt{t(1-t)}}=\frac{2}{\pi }\arcsin \sqrt{x}.\end{array}$$

$$S_n^{asin}=\frac{1}{n}\sum _{k=1}^n{D}_k\in [0,1]$$

$$\mathbb{P}\left(S_n^{asin}\le x|{\mathcal{H}}_0\right)\approx \frac{1}{\pi }{\int }_0^x\frac{dt}{\sqrt{t(1-t)}}=\frac{2}{\pi }\arcsin (\sqrt{x})$$

$$p_j=\mathbb{P}\left(S_n^{asin}>S_{n,j}^{asin}|{\mathcal{H}}_0\right)=1-\frac{2}{\pi }\arcsin \left(\sqrt{S_{n,j}^{asin}}\right),j=1,\dots ,m$$

$$\begin{array}{rl}{P}_1& =\left[0,\frac{1}{2s}\right),\\ P_i& =\left[\frac{2i-3}{2s},\frac{2i-1}{2s}\right),2\le i\le s,\\ P_{s+1}& =\left[1-\frac{1}{2s},1\right].\end{array}$$

$${\mu }_m\left(P_i\right)$$

$${\nu }_m\left(P_i\right)$$

$$T^{asin}=\sum _{i=1}^{s+1}\frac{(O_i-E_i{)}^2}{E_i}=m\cdot \sum _{\mathcal{A}\in {\mathcal{P}}_s}\frac{({\mu }_m(\mathcal{A})-{\nu }_m(\mathcal{A}){)}^2}{{\mu }_m(\mathcal{A})}.$$

$$p_{{\chi }^2}=\mathbb{P}\left(X>T^{asin}\right),$$

$$d_{tv}^{asin}=\frac{1}{2}\sum _{\mathcal{A}\in {\mathcal{P}}_s}|{\mu }_m(\mathcal{A})-{\nu }_m(\mathcal{A})|\in [0,1],$$

$$\mathbb{P}\left(2d_{tv}^{asin}>\epsilon |{\mathcal{H}}_0\right)\le 3exp\left(-\frac{m{\epsilon }^2}{25}\right).$$

$${\mu }_m^{\prime }\left(P_i\right)$$

$${\nu }_m^{\prime }\left(P_i\right)$$

$$T^{asin}=m\cdot \sum _{\mathcal{A}\in {\mathcal{P}}_s}\frac{({\mu }_m^{\prime }(\mathcal{A})-{\nu }_m^{\prime }(\mathcal{A}){)}^2}{{\mu }_m^{\prime }(\mathcal{A})},d_{tv}^{asin}=\frac{1}{2}\sum _{\mathcal{A}\in {\mathcal{P}}_s}|{\mu }_m^{\prime }(\mathcal{A})-{\nu }_m^{\prime }(\mathcal{A})|.$$

$$\underset{x\in [0,1]}{\sup }|F_n(x)-F^{asin}(x)|\le \frac{C}{n},C=\frac{4}{3\pi }\left(2-\frac{3}{2s}\right){\left(\frac{4s^2}{2s-1}\right)}^{\frac{3}{2}}.$$

$$\left|\mathbb{P}\left(S_n^{asin}\in (a,b)\right)-\frac{1}{\pi }{\int }_a^b\frac{dt}{\sqrt{t(1-t)}}\right|\le \frac{C}{n}.$$

$$p_{2k,2\mathfrak{n}}=\left(\genfrac{}{}{0ex}{}{2k}{k}\right)\left(\genfrac{}{}{0ex}{}{2(\mathfrak{n}-k)}{\mathfrak{n}-k}\right)2^{-2\mathfrak{n}}.$$

$$\mathfrak{n}!=\sqrt{2\pi \mathfrak{n}}{\left(\frac{\mathfrak{n}}{e}\right)}^{\mathfrak{n}}\exp \left(\frac{{\theta }_{\mathfrak{n}}}{12\mathfrak{n}}\right).$$

$$p_{2k,2\mathfrak{n}}=\frac{1}{\pi \sqrt{k(\mathfrak{n}-k)}}\exp \left(\frac{{\theta }_{2k}-4{\theta }_k}{24k}+\frac{{\theta }_{2(\mathfrak{n}-k)}-4{\theta }_{\mathfrak{n}-k}}{24(\mathfrak{n}-k)}\right).$$

$$\frac{p_{2k,2\mathfrak{n}}}{d_{k,\mathfrak{n}}}\le \exp \left(\frac{1}{24k}+\frac{1}{24(\mathfrak{n}-k)}\right)=\exp \left(\frac{\mathfrak{n}}{24k(\mathfrak{n}-k)}\right)$$

$$\frac{p_{2k,2\mathfrak{n}}}{d_{k,\mathfrak{n}}}\ge \exp \left(\frac{-4}{24k}+\frac{-4}{24(\mathfrak{n}-k)}\right)=\exp \left(-\frac{\mathfrak{n}}{6k(\mathfrak{n}-k)}\right).$$

$$\begin{array}[]{llllll}p_{2k,2\mathfrak{n}}-d_{k,\mathfrak{n}}&\leq&\displaystyle d_{k,\mathfrak{n}}\left(\exp\left(\frac{\mathfrak{n}}{24k(\mathfrak{n}-k)}\right)-1\right)&\leq&\displaystyle d_{k,\mathfrak{n}}\frac{\mathfrak{n}}{12k(\mathfrak{n}-k)},\\[10.0pt] d_{k,\mathfrak{n}}-p_{2k,2\mathfrak{n}}&\leq&\displaystyle d_{k,\mathfrak{n}}\left(1-\exp\left(-\frac{\mathfrak{n}}{6k(\mathfrak{n}-k)}\right)\right)&\leq&\displaystyle d_{k,\mathfrak{n}}\frac{\mathfrak{n}}{6k(\mathfrak{n}-k)},\\ \end{array}$$

$$|p_{2k,2\mathfrak{n}}-d_{k,\mathfrak{n}}|\le d_{k,\mathfrak{n}}\frac{\mathfrak{n}}{6k(\mathfrak{n}-k)}=\frac{\mathfrak{n}}{6\pi {\left(k(\mathfrak{n}-k)\right)}^{\frac{3}{2}}}.$$

$$|p_{2k,2\mathfrak{n}}-d_{k,\mathfrak{n}}|\le \frac{\mathfrak{n}}{6\pi {\left(\delta \mathfrak{n}(\mathfrak{n}-\delta \mathfrak{n})\right)}^{\frac{3}{2}}}=\frac{1}{6\pi {\mathfrak{n}}^2{\left(\delta (1-\delta )\right)}^{\frac{3}{2}}}.$$

$$\left|\sum _{a\le \frac{k}{\mathfrak{n}}\le b}{p}_{2k,2\mathfrak{n}}-\sum _{a\le \frac{k}{\mathfrak{n}}\le b}{d}_{k,\mathfrak{n}}\right|\le \sum _{a\le \frac{k}{\mathfrak{n}}\le b}\left|p_{2k,2n}-d_{k,\mathfrak{n}}\right|\le \sum _{a\le \frac{k}{\mathfrak{n}}\le b}\frac{1}{6\pi {\mathfrak{n}}^2{\left(\delta (1-\delta )\right)}^{\frac{3}{2}}}$$

$$=\frac{\lceil b\mathfrak{n}-a\mathfrak{n}\rceil }{6\pi {\mathfrak{n}}^2{\left(\delta (1-\delta )\right)}^{\frac{3}{2}}}\le \frac{b-a}{3\pi \mathfrak{n}{\left(\delta (1-\delta )\right)}^{\frac{3}{2}}}\le \frac{1}{3\pi \mathfrak{n}{\left(\delta (1-\delta )\right)}^{\frac{3}{2}}}=:\eta .$$

$$\left|{\int }_a^{b}f(x)dx-\sum _{a\le \frac{k}{\mathfrak{n}}\le b}\frac{1}{\mathfrak{n}}f(x_k)\right|\le \sum _{a\le \frac{k}{\mathfrak{n}}\le b}\frac{1}{\mathfrak{n}}(M_i-m_i)=\sum _{a\le \frac{k}{\mathfrak{n}}\le b}\frac{1}{{\mathfrak{n}}^2}|f^{\prime }({\xi }_i)|\le \sum _{a\le \frac{k}{\mathfrak{n}}\le b}\frac{1}{{\mathfrak{n}}^2}\underset{a\le x\le b}{\sup }|f^{\prime }(x)|$$

$$=\frac{\lceil b\mathfrak{n}-a\mathfrak{n}\rceil }{{\mathfrak{n}}^2}\underset{a\le x\le b}{\sup }|f^{\prime }(x)|\le \frac{2(b-a)}{\mathfrak{n}}\underset{a\le x\le b}{\sup }|f^{\prime }(x)|.$$

$$\left|{\int }_a^{b}f(x)dx-\sum _{a\le \frac{k}{\mathfrak{n}}\le b}{d}_{k,\mathfrak{n}}\right|$$

Full text

††journal: Journal of Computational and Applied Mathematics

1 Introduction

Random numbers are key ingredients in various applications, e.g., in cryptography (e.g., for generating cryptographic keys) or in simulations (e.g., in Monte Carlo methods), just to mention a few. No algorithm can produce truly random numbers. Instead, pseudorandom number generators (PRNGs) are used. These are deterministic algorithms producing numbers which we expect to resemble truly random ones in some sense. There are two classes of tests used to evaluate PRNGs, theoretical and statistical ones. Theoretical tests examine the intrinsic structure of a given generator, the sequence does not necessarily need to be generated. Two classical examples are the lattice test [1] and the spectral test described in [2] (Section 3.3.4). See also [3] for a description of some standard tests from this class. This category of tests is very specific to each family of generators e.g., some are designed only for linear congruential generators. On the other hand, the second class of tests – empirical tests – are conducted on a sequence generated by a PRNG and require no knowledge of how it was produced. The main goal of these tests is to check if the sequence of numbers $\mathbf{U}=(U_{1},U_{2},\ldots,U_{n})$ (or bits, depending on the actual implementation) produced by a PRNG has properties similar to those of a sequence generated truly at random. These tests try to find statistical evidence against the null hypothesis $\mathcal{H}_{0}$ stating that the sequence is a sample from independent random variables with uniform distribution. Any function of a finite number of uniformly distributed random variables, whose (sometimes approximate) distribution under hypothesis $\mathcal{H}_{0}$ is known, can be used as a statistical test. Due to the popularity and significance of the problem, a variety of testing procedures have been developed in recent years. Such statistical tests aim at detecting various deviations in generated sequences, what allows for revealing flawed PRNGs producing predictable output. Some of the procedures encompass classical tools from statistics like the Kolmogorov-Smirnov test or the Pearson’s chi-squared test, which are used for comparing the theoretical and empirical distributions of appropriate statistics calculated for a PRNG’s output. It is also possible to adapt tests of normality like the Anderson-Darling or Shapiro-Wilk tests for appropriately transformed pseudorandom sequences. These methods exploit the properties of sequences of i.i.d. random variables. Based on the original sequence $\mathbf{U}$ returned by the examined PRNG we are able to obtain realizations of random variables with known theoretical distributions. Some examples of probabilistic laws used in practice in this kind of tests can be found e.g., in [2]. They include such procedures like the gap test, the permutation test and the coupon collector’s test, just to name a few (see [2] for a more detailed treatment). These methods have also the advantage that they implicitly test the independence of the generator’s output. The main issue with such methods is that a single statistical test looks only at some specific property that holds for sequences of truly random numbers. Hence, for practical purposes bundles of diverse tests are created. Such a test bundle consists of a series of individual procedures based on various stochastic laws from probability theory. A PRNG is then considered as good if the pseudorandom sequences it produces pass all tests in a given bundle. Note that formally it proves nothing, but it increases the confidence in the simulation results. Thus, they are actually tests for non-randomness, as pointed out in [4]. Some examples of such test suites are Marsaglia’s Diehard Battery of Tests of Randomness from 1995, Dieharder developed by Brown et al. (see [5]), TestU01 implemented by L’Ecuyer and Simard (see [6, 7]) and NIST Test Suite [8]. The last one, designed by the National Institute of Standard and Technology, is currently considered as one of the state of the art test bundles. It is often used for the preparation of many formal certifications or approvals.

A result of a single statistical test is typically given in the form of a $p$-value, which, informally speaking, represents the probability that a perfect PRNG would produce “less random” sequence than the sequence being tested w.r.t. the used statistic. We then reject $\mathcal{H}_{0}$ if $p<\alpha$, where $\alpha$ is the significance level (usually $\alpha=0.01$) and accept $\mathcal{H}_{0}$ if $p\geq\alpha$. Such an approach is usually called one level or first level test. Although the interpretation of a single $p$-value has a clear statistical explanation, it is not quite obvious how to interpret the results of a test bundle, i.e., of multiple tests. Under $\mathcal{H}_{0}$ the distribution of $p$-values is uniform. However, in a test bundle several different tests are applied to the same output of a PRNG, hence the results are usually correlated. The documentation of the NIST Test Suite includes some clues on how to interpret the results of their bundle (Section 4.2 in [9]), but in the introduction it is frankly stated: “It is up to the tester to determine the correct interpretation of the test results”.

To disclose flaws of PRNGs, a very long sequence is often required. In such situations, the applicability of a statistical test can be limited (depending on the test statistic) by the memory size of the computer. An alternative approach is to use a so-called two level (a term used e.g., in [3]) or second level (a term used e.g., in [4, 10]) test. In this approach we take into account several results from the same test over disjoint sequences generated by a PRNG. We obtain several $p$-values which are uniformly distributed under $\mathcal{H}_{0}$, what is tested by e.g., some goodness-of-fit test (with potentially different level of significance – NIST suggests to use $0.0001$ – obtaining new “final” $p$-value). The authors in [11] observed that this method may be comparable to a first level test in terms of the power of a test (informally speaking, it represents the probability of observing “less random” sequence than the sequence being tested under an alternative hypothesis $\mathcal{H}_{1}$, see [11] for details), but often it produces much more accurate results, as shown in [4]. Roughly speaking, the accuracy is related to the ability, given a non-random PRNG, of recognizing its sequences as non-random (for details see [4]). We will follow this approach.

In the second level approach one has to take under consideration the approximation errors in the computation of a $p$-value. For example, in a first level test one usually calculates a $p$-value of a statistic which – under $\mathcal{H}_{0}$ – is approximately normally distributed. The approximation comes then from the central limit theorem, which lets us substitute the distribution of a given sum with the standard normal distribution. These errors in calculations of individual $p$-values may accumulate, resulting in an error of a $p$-value in a second level test, thus making the test not reliable. Following [4] we say that the second level test is not reliable when, due to errors or approximations in the computation of $p$-values (in the first level), the distribution of $p$-values is not uniform under $\mathcal{H}_{0}$. Fortunately, this approximation error can be bounded using the Berry-Essen inequality and the final error of a second level test can be controlled (see [4, 10] for a detailed example based on the binary matrix rank test). The influence of approximations on the computation of $p$-values in a second level test was also considered in [12, 13]. In this article we present a statistical test based on the arcsine law, in which at some point we approximate a distribution of some random variable with the arcsine distribution. We provide a Berry-Essen type inequality which upper bounds the approximation error, what allows us to control the reliability of our second level test.

An interesting approach for testing PRNGs was presented by Kim et al. in [14]. The concept of their tests is based on the properties of a random walk (the gambler’s ruin algorithm) on the cyclic group $\mathbb{Z}_{n}=\{0,\ldots,n-1\}$ with 0 being an absorbing state – more precisely, on the time till absorption. The authors in [14] propose three different variants of the test. The general idea of the basic procedure is the following. For some fixed $p\in(0,1)$ and $x\in\mathbb{Z}_{n}$, the output $\mathbf{U}=(U_{i})$ of a PRNG is treated as numbers from the unit interval and used to define a random walk starting in $x$ such that if $U_{i}<p$ and the process is in state $s$, then it moves to $(s+1)\mod n$, otherwise it moves to $(s-1)\mod n$. The aim of this test is to compare the theoretical and the empirical distributions of the time to absorption in 0 when starting at $x$. Based on the values of testing statistic, the PRNG is then either accepted or rejected. The authors reported some “hidden defects” in the widely used Mersenne Twister generator. However, one has to be very careful when dealing with randomness. It seems like re-seeding a PRNG with a fixed seed is an error which can lead to wrong conclusions. The criticism was raised by Ekkehard and Grønvik in [15], where the authors also showed that the properly performed tests of Kim et al. [14] do not reveal any defects in the Mersenne Twister PRNG. Recently, the authors in [16] have proposed another gambler’s ruin based procedure for testing PRNGs. In their method they exploited formulas for winning probabilities for arbitrary sequences $p(i)$ and $q(i)$, (i.e., the winning and losing probabilities depend on the current fortune) which are the parameters of the algorithm.

In recent years a novel kind of testing techniques has been introduced for more careful verification of generators. The core idea of this class of methods is based on an observation that the binary sequence $(B_{i})$ produced by a PRNG, after being properly rescaled, can be interpreted as an one-dimensional random walk $(S_{n})_{n\in\mathbb{N}}$ with $S_{k}=\sum_{i=1}^{k}X_{i}$, where $X_{i}=2B_{i}-1$. For random walks defined by truly random binary sequences a wide range of statistics have been considered over the years and a variety of corresponding stochastic laws have been derived (see e.g., [17]). For a good PRNG we may expect that its output will behave like $S_{n}$. Hence, the following idea comes to mind: choose some probabilistic law that holds for truly random bit sequences and compare the theoretical distribution of the corresponding statistic with the empirical distribution calculated for $m$ sequences produced by a given PRNG in $m$ independent experiments. This comparison can be done e.g., by computing the $p$-value of an appropriate test statistic under the null hypothesis that the sequence generated by this PRNG is truly random.

Another concept named statistical distance based testing was suggested in [18]. It relies on calculation of some statistical distances like e.g., total variation distance between the theoretical and empirical distributions for considered characteristics and rejecting a PRNG if the distance exceeds some threshold. We will also follow this approach, indicating the corresponding threshold. In [18] the authors derive their test statistics from the law of iterated logarithm for random walks (the procedure is called the LIL test). The proposed by us procedure uses similar methodology and is based on the arcsine law. We made the code publicly available, see [19]. It includes the arcsine law based as well as the law of iterated logarithm based statistical tests, the implementation of many PRNGs (more than described in this article) including the Flawed generator (see Section 4) and the seeds we used.

Organization of the paper

In the following Section 2 we define a general notion of a PRNG and recall the aforementioned stochastic laws for random walks. The testing method along with the error analysis is described in Section 3. The concise report on experimental results (including the Flawed generator introduced in Section 4) is given in Section 5. In Section 6 we mention other implementations of the tests based on the arcsine law. We conclude in Section 7.

2 Pseudorandom generators and stochastic laws for random walks

2.1 Pseudorandom generators

The intuition behind pseudorandom number generator is clear. However, let us give a strict definition roughly following Asmussen and Glynn [20].

Definition 2.1**.**

A Pseudorandom number generator (PRNG) is a 5-tuple $<E,V,s_{0},f,g>$, where $E$ is a finite state space, $V$ is a set of values, $s_{0}\in E$ is a so-called seed, i.e., an initial state in the sequence $(s_{i})_{i=0}^{\infty}$, a function $f:E\to E$ describes the transition between consecutive states $s_{n}=f(s_{n-1})$ and $g:E\to V$ maps the generator’s state into the output.

Usually $V=(0,1)$ or $V=\{0,1,\ldots,M-1\}$ for some $M\in\mathbb{N}$, the latter one is used throughout the paper. Recall that LCG (linear congruential generator) is a generator which updates its state according to the formula $s_{n}=(as_{n-1}+c)\bmod M$. Thus, it is defined by three integers: a modulus $M$, a multiplier $a$, and an additive constant $c$. In the case $c=0$, the generator is called MCG (multiplicative congruential generator). For a detailed description of some commonly used PRNGs see the surveys [21, 22, 23] or the book [24].

It is clear that both the input and the output of a random number generator can be viewed as a finite sequence of bits. For a PRNG to be considered as good, the output sequences should have some particular property, namely each returned bit has to be generated independently with equal probability of being 0 and 1. We say that the sequence of bits is truly random if it is a realization of a Bernoulli process with success probability $p=\frac{1}{2}$.

Given a PRNG $G$ returning integers from the set $V$, we may obtain a pseudorandom binary sequence with any given length using the following simple procedure. Namely, as long as the bit sequence $s$ is not sufficiently long, generate the next pseudorandom number $a$ and append its binary representation (on $\left\lceil{\log_{2}M}\right\rceil$ bits) to the current content of $s$. In the ideal model with $G$ being truly random number generator, such algorithm produces truly random bit sequences provided that $M$ is a power of 2. Indeed, for $M=2^{k}$ there is one to one correspondence between $k$-bit sequences and the set $V$. Hence, if each number is generated independently with uniform distribution on $V$, then each combination of $k$ bits is equally likely and therefore each bit of the output sequence is independent and equal to 0 or 1 with probability $\frac{1}{2}$.

However, this is not true for $M\neq 2^{k}$. It is easy to observe that in such a case the generator is more likely to output 0s and the generated bits are no longer independent. Thus, rather than simply outputting the bits of $a$, one may instead take $d$ first bits from the binary representation of $\frac{a}{M}$ for some fixed $d$. Such a method has the advantage that it can be easily adopted for an underlying generator returning numbers from the unit interval, what is common for many PRNG implementations.

2.2 Stochastic laws for random walks

Let $(B_{i})_{i\geq 0}$ be a Bernoulli process with a parameter $p\in(0,1)$, i.e., a sequence of independent random variables with identical distribution $P(B_{1}=1)=1-P(B_{1}=0)=p$. A good PRNG should behave like a generator of Bernoulli process with $p=1/2$ (what we assume from now on). It will be, however, more convenient to consider the following transformed process

[TABLE]

The sequence $X_{i}$ is $\{-1,+1\}$-valued, the process $(S_{n})_{n\in\mathbb{N}}$ is called a random walk.

The law of iterated logarithm

Of course $|S_{n}|\leq n$. However, large values of $|S_{n}|$ occur with small probability and the values of $S_{n}$ are in practice in a much narrower range than $[-n,n]$. The weak and the strong law of large numbers imply that ${S_{n}\over n}\stackrel{{\scriptstyle P}}{{\to}}0,\quad\textrm{and even }{S_{n}\over n}\stackrel{{\scriptstyle a.s.}}{{\to}}0,$ where $\stackrel{{\scriptstyle P}}{{\to}}$ denotes the convergence in probability and $\stackrel{{\scriptstyle a.s.}}{{\to}}$ denotes the almost sure convergence. Thus, the deviations of $S_{n}$ from 0 grow much slower than linearly. On the other hand the central limit theorem states that ${S_{n}\over n}\stackrel{{\scriptstyle D}}{{\to}}\mathcal{N}(0,1)$ (where $\stackrel{{\scriptstyle D}}{{\to}}$ denotes the convergence in distribution), what is in some sense a lower bound on fluctuations of $S_{n}$ – they will leave the interval $[-\sqrt{n},\sqrt{n}]$ since we have $\limsup_{n\to\infty}{S_{n}\over\sqrt{n}}=\infty$ (implied by 0-1 Kolmogorov’s Law, see e.g., Theorem 5.1 in [25]). It turns out that the fluctuations can be estimated more exactly.

Theorem 2.2** (The law of iterated logarithm, [26], cf. also Chapter VIII.5 in [17]).**

For a random walk $S_{n}$ we have

[TABLE]

Thus, to normalize $S_{n}$ dividing by $n$ is too strong and dividing by $\sqrt{n}$ is too weak. The fluctuations of $S_{n}$ from 0 grow proportionally to $\sqrt{2n\log\log n}$.

To depict the law of iterated logarithm, we took $500$ output sequences $\mathcal{B}^{1},\ldots,\mathcal{B}^{500}$ from the Mersenne Twister MT19937 generator, each initialized with a random seed taken from http://www.random.org, where each output $\mathcal{B}^{j}=(B^{j}_{1},\ldots,B^{j}_{n}),B^{j}_{i}\in\{0,1\},j=1,\ldots,500,i=1,\ldots,n$ was of length $n=2^{30}$. In Figure 1 we presented these 500 trajectories $(k,S^{j}_{k}),j=1,\ldots,500,k=0,\ldots,n$, where $S^{j}_{k}=\sum_{i=1}^{k}(2B_{i}^{j}-1)$. Each trajectory is depicted by a single polyline. The darker the image the higher the density of trajectories. We can see that $\pm\sqrt{2n\log\log n}$ roughly corresponds to the fluctuations of $S_{n}$. However, few trajectories after around billion steps are still outside $[-\sqrt{2n\log\log n},\sqrt{2n\log\log n}]$. The law of iterated logarithm tells us that for appropriately large $n$ the trajectories will not leave $[-\sqrt{2n\log\log n},\sqrt{2n\log\log n}]$ with probability $1$, what is not the case in Figure 1. It means that $n$ must be much larger than $2^{30}$.

One could think that the following is a good test for randomness: fix some number, say $100$, and classify the considered PRNG as good if the difference between the number of ones and zeros never exceeds $100$. The large difference may suggest that zeros and ones have different probabilities of occurrence. However, the law of iterated logarithm tells us that this reasoning is wrong. Indeed, we should expect some fluctuations and the absence of them means that a PRNG does not produce bits which can be considered random. This property of random walks was used by the authors in [18] for designing a novel method of testing random number generators.

There is yet another interesting property. Define $S_{n}^{lil}={S_{n}\over\sqrt{2n\log\log n}}$. The law of iterated logarithm implies that $S_{n}^{lil}$ does not converge pointwise to any constant. However, it converges to 0 in probability. Let us fix some small $\varepsilon>0$. For almost all $n$, with an arbitrary high probability $p<1$ the process $S_{n}^{lil}$ will not leave $(-\varepsilon,\varepsilon)$. On the other hand, this tells us that the process will be outside this interval infinitely many times. This apparent contradiction shows how can our intuition be unreliable on phenomena taking place at infinity.

The arcsine law

The observations described previously imply that averaging every $S_{n}$, it will spend half of its time above the $x$-axis and half of its time below. However, the typical situation is counter-intuitive (at first glance): typically the random walk will either spend most of its time above or most of its time below the $x$-axis. This is expressed in the Theorem 2.3 below (for reference see e.g., [17]). Before we formulate the theorem, let us first introduce some notations. For a sequence $X_{1},X_{2},\ldots$, as defined in (1), let

[TABLE]

where $\mathbbm{1}(\cdot)$ is the indicator function. $D_{k}$ is equal to 1 if the number of ones exceeds the number of zeros either at step $k$ or at step $k-1$, and 0 otherwise (in a case of ties, i.e., $S_{k}=0$, we look at the previous step letting $D_{k}=D_{k-1}$). In other words, $D_{k}=1$ corresponds to the situation in which the line segment of the trajectory of the random walk between steps $k-1$ and $k$ is above the $x$-axis.

Theorem 2.3** (The arcsine law).**

Let $(B_{i})_{i\geq 0}$ be a Bernoulli process. Define $X_{i}=2B_{i}-1$ and $L_{n}=\sum_{k=1}^{n}D_{k}$ ($D_{k}$ is given in (2)). For $x\in(0,1)$ we have

[TABLE]

The probability $\mathbb{P}\left(L_{n}\leq x\cdot n\right)$ is the chance that the random walk was above the $x$-axis for at most $x$ fraction of the time. The limiting distribution is called the arcsine distribution. Its density function is given by $f^{asin}(t)={\frac{1}{\pi}\sqrt{t(1-t)}}$ and the cumulative distribution function (cdf) is $F^{asin}(t)=\frac{2}{\pi}\arcsin\sqrt{t}$. The shape of the pdf $f^{asin}(t)$ clearly indicates that the fractions of time spent above and below [math]-axis are more likely to be unequal than close to each other.

3 Testing PRNGs based on the arcsine law

In this Section we will show how to exploit the theoretical properties of random walks from the preceding discussion to design a practical routine for testing PRNGs. We describe our approach based on the arcsine law which we employ for experimental evaluation of several commonly used generators (the results are presented in Section 5). We also perform an error analysis of the proposed testing procedure, providing corresponding bounds on the approximation errors. Finally, we make some remarks on the reliability of our second level test.

3.1 The arcsine law based testing

The general idea of tests is the following. Take a sequence of bits generated by PRNG, rescale them as in (1) and compare the empirical distribution of

[TABLE]

(a fraction of time instants at which ones prevail zeros) with its theoretical distribution assuming that truly random numbers were generated. In terms of hypothesis testing: given the null hypothesis $\mathcal{H}_{0}$ that the bits in the sequence were generated independently and uniformly at random (vs. $\mathcal{H}_{A}$: that the sequence was not randomly generated), the distribution of $S^{asin}_{n}$ follows the arcsine law (Theorem 2.3), i.e., we can conclude that for large $n$ we have

[TABLE]

(we will be more specific on “$\approx$” in Section 3.2). We follow the second level testing approach (cf. [4, 10]), i.e., we take into account several results from the same test over different sequences. To test a PRNG we generate $m$ sequences of length $n$ each, thus obtaining $m$ realizations of the variable $S^{asin}_{n}$. Denoting by $S^{asin}_{n,j}$ the value of $j$-th simulation’s result (we call them a basic tests), we then calculate the corresponding $p$-values

[TABLE]

Under $\mathcal{H}_{0}$ the distribution of $p_{j},j=1,\ldots,m$, should be uniform on $[0,1]$. We fix some partition of $[0,1]$ and count the number of $p$-values within each interval. In our tests we will use an $(s+1)$-element partition $\mathcal{P}_{s}=\{P_{1},\ldots,P_{s+1}\}$, where

[TABLE]

Now we define the measures $\mu_{m}$ (the uniform measure on $\mathcal{P}_{s}$), $\nu_{n}$ (the empirical measure on $\mathcal{P}_{s}$), $E_{i}$ (the expected number of $p$-values within $P_{i}$) and $O_{i}$ (the number of observed $p$-values within $P_{i}$). For $1\leq i\leq s+1$ let

[TABLE]

We perform the Pearson’s goodness-of-fit test, which uses the following test statistic

[TABLE]

Under the null hypothesis, $T^{asin}$ has approximately $\chi^{2}(s)$ distribution. We calculate the corresponding $p$-value

[TABLE]

where $X$ has a $\chi^{2}(s)$ distribution. Large values of $T^{asin}$ – and thus small values of $p_{\chi^{2}}$ – let us suspect that a given PRNG is not good. Typically, we reject $\mathcal{H}_{0}$ (i.e., we consider the test failed) if $p_{\chi^{s}}<\alpha,$ where $\alpha$ is a predefined level of significance (for a second level test we use $\alpha=0.0001$, as suggested by NIST). Note that the probability of rejecting $\mathcal{H}_{0}$ when the sequence is generated by a perfect random generator (so-called Type I error) is exactly $\alpha$.

Another approach relies on the statistical distance based testing, which is the technique presented in [18]. We consider the statistic

[TABLE]

i.e., a total variation distance between the theoretical distribution $\mu_{m}$ and the empirical distribution $\nu_{m}$. Similarly, large values of $d_{tv}^{asin}$ indicate that a given PRNG is not good. Concerning Type I error we will make use of the following lemma (see Lemma 3 in [27] or its reformulation, Lemma 1 in [28]).

Lemma 3.4**.**

Assume $\mathcal{H}_{0}$ and consider the partition $\mathcal{P}_{s}$. Then, for all $\varepsilon\geq\sqrt{20(s+1)/m}$ we have

[TABLE]

To summarize, for a given PRNG we generate $m$ sequences of length $n$ each. and we choose $s$ (and thus the partition $\mathcal{P}_{s}$). We then calculate $d^{asin}_{tv}$ and $T^{asin}$ together with its $p_{\chi^{2}}$-value. We specify the thresholds for $p_{\chi^{2}}$-value and $d^{asin}_{tv}$ indicating whether the test failed or not (the details are presented in Section 5). We denote the described procedure as the ASIN test.

Remark. Note that the described procedure for calculating $T^{asin}$ and $d_{tv}^{asin}$ is equivalent to the following one. Instead of calculating $p$-values of $S^{asin}_{n,j}$, we could directly count the number of $S^{asin}_{n,j}$ falling into each interval $P_{i},i=1,\ldots,s+1$ and compare the empirical distribution with the theoretical one. To be more precise, for $1\leq i\leq s+1$ let

[TABLE]

where $\mathbb{P}\left(S^{asin}_{n}\in P_{i}\right)=F^{asin}(b)-F^{asin}(a)$ for $P_{i}=[a,b]$. Then statistics $T^{asin}$ and $d_{tv}^{asin}$ can be rewritten as

[TABLE]

This technique was presented in [18] (for the total variation and few other distances) and this is how our implementation of the ASIN test [19] calculates the statistics.

We could also calculate just one $p$-value of the statistic $S_{n^{\prime}}^{asin}$ for a longer sequence (say, for $n^{\prime}=n\cdot m$) – i.e., perform a first level test. However, as mentioned in Section 1, the second level approach produces more accurate results (roughly speaking, the accuracy is related to the ability, given a non-random PRNG, of recognizing its sequences as non-random, see details in [4]).

It is worth noting that the following approach can be applied when $d^{asin}_{tv}$ or $p_{\chi^{2}}$ are slightly outside the acceptance region (e.g., if $p_{\chi^{2}}\in(10^{-4},10^{-2}$), what suggests rejecting $\mathcal{H}_{0}$, but is not a strong evidence). Namely, double the length of the sequence, take a new output from the PRNG and apply the test again. Repeat the procedure (at most some predefined number of times) until the evidence is strong enough (e.g., $p_{\chi^{2}}<10^{-4}$) or $\mathcal{H}_{0}$ is accepted (e.g., $p_{\chi^{2}}>0.01$). This method, called “automation of statistical tests on randomness”, was proposed and analyzed in [29].

3.2 Error analysis

3.2.1 Bounding errors in approximating $p$-values in basic tests

In this subsection we will show a bound on the approximation error in (3). Recall that $F^{asin}(x)={2\over\pi}\arcsin\sqrt{t}$.

Lemma 3.5**.**

Fix a partition $\mathcal{P}_{s},s\geq 2$ and an even $n\geq 2$. Let $F_{n}$ be the cdf of the empirical distribution of $S^{asin}_{n}$ under $\mathcal{H}_{0}$ (stating that the bits $B_{1},\ldots,B_{n}$ were generated uniformly at random), i.e., $F_{n}(x)=\mathbb{P}\left(S^{asin}_{n}\leq x|\mathcal{H}_{0}\right)$. Then we have

[TABLE]

Proof. We will show that for fixed $a$ and $b$ such that $0\leq a<b\leq 1$ we have

[TABLE]

Let us assume that $n=2\mathfrak{n}$. Let $p_{2k,2\mathfrak{n}}$ denote the probability that during $2k$ steps in the first $2\mathfrak{n}$ steps the random walk was above 0-axis, i.e., $p_{2k,2\mathfrak{n}}=\mathbb{P}\left(L_{2\mathfrak{n}}=2k\right)$. The classical results on a simple random walk state that

[TABLE]

The standard proof of Theorem 2.3 (see, e.g., Chapter XII.8 in [30]) shows that $p_{2k,2\mathfrak{n}}$ converges to $d_{k,\mathfrak{n}}={1\over\pi\sqrt{k(\mathfrak{n}-k)}}$. In the following, we will bound the difference $|p_{2k,2\mathfrak{n}}-d_{k,\mathfrak{n}}|$. We will use a version of Stirling’s formula stating that for each $\mathfrak{n}$ there exists $\theta_{\mathfrak{n}}$, $0<\theta_{\mathfrak{n}}\leq 1$, such that

[TABLE]

Plugging (5) into each factorial appearing in (4) we have

[TABLE]

Thus, we get

[TABLE]

and

[TABLE]

For any $x$ it holds that $1-e^{-x}\leq x$ and for $x\in[0,1.25]$ we have that $e^{x}-1\leq 2x$. Note that $\frac{\mathfrak{n}}{24(\mathfrak{n}-k)}\leq 1.25$, what is equivalent to $\mathfrak{n}\geq{30k^{2}\over 30k-1}$, what holds for any $0<k<\mathfrak{n}$. Hence,

[TABLE]

what implies

[TABLE]

Fix $\delta>0$ and assume furthermore that $\delta\leq\frac{k}{\mathfrak{n}}\leq 1-\delta$. The function $k\mapsto\left(k(\mathfrak{n}-k)\right)^{3/2}$ achieves the minimum value at the endpoints of the considered interval, thus

[TABLE]

We will estimate the approximation error in (3) in two steps. First, take two numbers $a,b$ such that $\delta\leq a<b\leq 1-\delta$. We have

[TABLE]

The second kind of errors in probability estimates given by (3) is caused by approximating the sum by an integral. Let us consider an arbitrary function $f$ differentiable in the interval $(a,b)$. Split $(a,b)$ into subintervals of length $\frac{1}{\mathfrak{n}}$ and let $x_{k}$ be an arbitrary point in the interval containing $\frac{k}{\mathfrak{n}}$. Denote by $M_{k}$ and $m_{k}$ the maximum and the minimum value of $f$ on that interval, respectively. Using the Lagrange’s mean value theorem we obtain

[TABLE]

For $f(x)=\frac{1}{\pi\sqrt{x(1-x)}}$ we have $f^{\prime}(x)=\frac{2x-1}{2\pi(x(1-x))^{3/2}}$ and $\frac{1}{\mathfrak{n}}f\left(\frac{k}{\mathfrak{n}}\right)=d_{k,\mathfrak{\mathfrak{n}}}.$ Hence, in the considered interval $(a,b)\subseteq(\delta,1-\delta)$ we have

[TABLE]

We also have

[TABLE]

Taking $\delta={1\over 2s}$ we obtain

[TABLE]

what justifies the approximation (3) for $\delta\leq a<b\leq 1-\delta$. To complete the analysis we need to investigate the errors “on the boundaries” of a unit interval, i.e., for $(0,\delta)$ (and, by symmetry, for $(1-\delta,\delta)$). We get

[TABLE]

where the last inequality follows directly from the preceding calculations. ∎

Remark. Let $X_{1},X_{2},\ldots$ be zero-average i.i.d. random variables with $E|X_{i}|^{3}<\infty$. Denote $EX_{i}^{2}=\sigma^{2}$. The central limit theorem states that $\mathcal{N}$, a normal random variable $N(0,1)$ (denote its cdf by $\Phi$), is the limiting distribution of $Y_{n}=\sum_{i=1}^{n}{X_{i}\over\sigma\sqrt{n}}$ (denote its cdf by by $F^{Y}_{n}$). It means that for large $n$ we can approximate $Y_{n}$ by $\mathcal{N}$ and the approximation error is bounded by the Berry-Essen inequality

[TABLE]

where $C_{0}$ is a positive constant (in original paper [31] it was shown that $C_{0}\leq 7.59$, in [32] it was shown that $C_{0}\leq 0.4785$). Lemma 3.5 is thus a Berry-Essen type inequality for approximating $S^{asin}_{n}$ by a random variable with cdf $F^{asin}$, tailored to our needs.

3.2.2 Reliability of the results from the second level test

Following [10], we say that a basic test (calculating $S^{asin}_{n,j}$) is not reliable if, due to approximation errors in the computations of $p_{j}$-values, the distribution of $p_{j},j=1,\ldots,m$ for truly random numbers is not uniform. We test the uniformity via $T^{asin}$ and $d^{asin}_{tv}$. Since we compare two continuous distributions, some discretization needs to be applied. In our testing procedure we use a partition $\mathcal{P}_{s}$ for this, splitting the interval $[0,1]$ into $(s+1)$ intervals (i.e., the bins). Lemma 3.5 states that a maximum error in the computation of $p_{j}$ is bounded by ${C\over n}$ (note that $C$ implicitly depends on $s$). It means that a $p_{j}$-value that should belong to a given bin can be found in the neighboring ones only if the distance between $p_{j}$ and one of the endpoints of a given bin is less than ${C\over n}$. Thus, this is also the fraction of $p_{j}$-values that can be found in wrong bins. The maximum propagated deviation is twice the error (since most bins have two neighbors), i.e.,

[TABLE]

Under $\mathcal{H}_{0}$ the distribution of the numbers $p_{j},j=1,\ldots,m$ in the bins $1,\ldots,s+1$ is a multinomial distribution. Indeed, this is equivalent to throwing $m$ balls independently into $s+1$ bins, where the probability of choosing first and last bin is ${1\over 2s}$ and ${1\over s}$ for all remaining bins. The variance of the ratio of number of balls in bin $j\in\{1,s+1\}$ is equal to ${2s-1\over m4s^{2}}$, and for bin $j\in\{2,\ldots,s\}$ is equal to ${s-1\over ms^{2}}$. We have $\sigma>{\sqrt{s-1\over s^{2}m}}$, where $\sigma$ is the expected statistical deviation of the ratio of $p_{j}$-values found in a given bin. We expect that the error in approximating $p_{j}$-values propagates into an additional deviation. If the deviation is smaller than the statistical deviation, i.e., if

[TABLE]

then we say that the second level test is reliable. Note that the reliability of a test imposes a restriction on a relation between the length of a sequence used for each base test (i.e., $n$) and the number of basic tests ($m$). Inequality (6) implies a lower bound on $m$, namely

[TABLE]

4 The Flawed PRNG

In this section we present $\textsf{Flawed}_{rng,N,\tau}$ – a family of PRNGs. The family depends on three parameters: $rng$ (a PRNG, e.g., the Mersenne Twister), $N$ (a small integer parameter, e.g., $N=30$) and $\tau\in[0,1]$. $\textsf{Flawed}_{rng,N,\tau}(seed)$ generates the same output as $rng(seed)$ for a fraction $1-\tau$ of all possible seeds. For the remaining fraction $\tau$ of seeds it outputs bits such that the corresponding walk of the length $n=2^{N}$ spends exactly half of the time ($2^{N-1}$ steps) above zero and exactly half of the time below zero. In the following we will denote $n=2\mathsf{n}=2^{N}$.

4.1 Dyck Paths

To generate walks with the aforementioned property we will use Dyck paths, i.e., walks starting and ending at 0 with the property that for each prefix the number of ones is not smaller than the number of zeros.

Definition 4.6**.**

A sequence of $2\mathsf{n}$ bits $B_{1},\ldots,B_{2\mathsf{n}}$ is called a Dyck path if the corresponding walk $S_{k}$ fulfills $S_{k}=\sum_{i=1}^{k}(2B_{i}-1)\geq 0,k=1,\ldots,2\mathsf{n}-1$ and $S_{2\mathsf{n}}=0$. A set of all Dyck paths of length $2\mathsf{n}$ is denoted by $\mathcal{D}_{2\mathsf{n}}$.

Thus, a Dyck path of length $2\mathsf{n}$ corresponds to a valid grouping of $\mathsf{n}$ pairs of parentheses. We have $|\mathcal{D}_{2\mathsf{n}}|=C_{\mathsf{n}}={1\over\mathsf{n}+1}{2\mathsf{n}\choose\mathsf{n}}$ ($C_{\mathsf{n}}$ is the $\mathsf{n}$-th Catalan number).

4.2 Sampling Dyck Paths

We are interested in generating Dyck paths uniformly at random. To achieve this goal we will use the following three ingredients.

(1) Walk sampling

Let $\mathcal{I}_{2\mathsf{n}+1}^{-1}$ be the set of sequences of bits $B_{1},\ldots,B_{2\mathsf{n}+1}$ such that the corresponding walk $S_{k}$ ends at $-1$, i.e., $S_{2\mathsf{n}+1}=-1$. One can easily sample a sequence $I\in\mathcal{I}^{-1}_{2\mathsf{n}+1}$ uniformly at random – it is enough to make a random permutation of the vector of bits $(0,\ldots,0,1,\ldots,1)$, consisting of $\mathsf{n}+1$ zeros and $\mathsf{n}$ ones.

(2) $f_{\textsf{Dyck}}$ transformation

One can obtain a Dyck path of length $2\mathsf{n}$ from $I\in\mathcal{I}_{2\mathsf{n}+1}^{-1}$ using Algorithm 1.

Observe that $f_{\textsf{Dyck}}$ transforms $I\in\mathcal{I}_{2\mathsf{n}+1}^{-1}$ into a Dyck path. This follows from simple observations:

$I$ has exactly $\mathsf{n}+1$ zeros and $\mathsf{n}$ ones; 2. 2.

since $t=\min(\textrm{argmin}\{S_{k}:k\in\{1,\ldots,2\mathsf{n}+1\}\})$ then $B_{t}=0$ and after $B_{t}$ is removed then $f_{\textsf{Dyck}}(I)$ has exactly $\mathsf{n}$ bits equal to [math] and $\mathsf{n}$ bits equal to $1$; 3. 3.

from the definition of $t$ (which enforces in particular that $B_{t+1}=1$), the walk that corresponds to bits $(B_{t+1},\ldots,B_{2\mathsf{n}+1},B_{1},\ldots,B_{t-1})$ cannot go below [math].

An example of a $f_{\textsf{Dyck}}$ transformation is presented in the Figure 2.

(3) The Cycle Lemma

The correspondence between the set $\mathcal{I}_{2\mathsf{n}+1}^{-1}$ and the set $\mathcal{D}_{2\mathsf{n}}$ is expressed by the Cycle Lemma (see, e.g., [33]).

Lemma 4.7** (The Cycle Lemma).**

For any $I\in\mathcal{I}_{2\mathsf{n}+1}^{-1}$ the path $f_{\textsf{Dyck}}(I)$ is a Dyck path. Moreover, any Dyck path in $\mathcal{D}_{2\mathsf{n}}$ is the image of exactly $2\mathsf{n}+1$ paths in $\mathcal{I}_{2\mathsf{n}+1}^{-1}$.

Thus, to obtain a random sampling of a Dyck path of length $2\mathsf{n}$ one needs to run Algorithm 2

4.3 The Flawed generator

As mentioned at the beginning of this section, the $\textsf{Flawed}_{rng,N,\tau}(seed)$ generator – described as Algorithm 3 – works exactly the same as the underlying generator for a fraction $(1-\tau)$ of seeds (lines 1-2). For the remaining fraction of $\tau$ of seeds (lines (3-8), the “else” branch) the output is generated in the following way.

The first $2^{N-2}$ bits are exactly the same as the first $2^{N-2}$ bits of $rng(seed)$ (line 4). 2. 2.

The next $2^{N-2}$ bits $(z_{2^{N-2}+1},\ldots,z_{2^{N-1}})$ are generated as follows:

1. (a)

a pseudorandom permutation $\pi$ is generated (line 5), 2. (b)

the bit $z_{2^{N-2}+i}$ is set to be equal to $1-z_{\pi(i)}$ (lines 6-8).

As the result, there is the same number of zeros and ones in the first $2^{N-1}$ bits – these bits are denoted as $\mathbf{z}_{L}$ (i.e., the corresponding walk is at zero at step $2^{N-1}$). 3. 3.

The remaining $2^{N-1}$ bits (denoted as $\mathbf{z}_{R}$) in the block are obtained by calling $\mathbf{z}_{R}\leftarrow\textsf{DyckPaths}(N,\mathbf{z}_{L},rng,seed)$ (Algorithm 4). As the result, the whole block $\mathbf{z}_{L}\mathbf{z}_{R}$ of $2^{N}$ output bits (concatenated blocks of $\mathbf{z}_{L}$ and $\mathbf{z}_{R}$) has the property that the corresponding walk spends the same number of steps above and below 0 (the description of DyckPaths algorithm is below).

Example 4.8** (Flawed).**

Let $N=5$ and a $seed$ is such that lines 4-10 are executed. Let the result of line 4 be: $(z_{1},\ldots,z_{8})=(1,1,0,0,0,0,0,1)\leftarrow rng(seed)$ and line 5 returns a permutation:

[TABLE]

Then the bits computed in lines 6-8 are:

[TABLE]

Then bits $z_{\mathbf{L}}=(z_{1},\ldots,z_{8},z_{9},\ldots,z_{16})$ are used as an input to function $\mathsf{DyckPaths}$ in line 10. The example is continued as Example 4.9.

Let us assume that Algorithm 3 (Flawed) has generated bits $\mathbf{z}_{L}=(z_{1},\ldots,z_{2^{N-1}})$ by executing lines $4-8$. Then the corresponding random walk spent $T_{U}(\mathbf{z}_{L})$ time under the $x$-axis and $T_{O}(\mathbf{z}_{L})$ time over the $x$-axis ($T_{U}(\mathbf{z}_{L})=2^{N-1}-T_{O}(\mathbf{z}_{L})$). The goal of the procedure $\textsf{DyckPaths}(N,\mathbf{z}_{L},rnd,seed)$ is to generate bits $\mathbf{z}_{R}=(z_{2^{N-1}+1},\ldots,z_{2^{N}})\leftarrow\textsf{DyckPaths}(N,(z_{1},\ldots,z_{2^{N-1}}),rng,seed)$ in such a way that:

$T_{U}(\mathbf{z}_{R})=T_{O}(\mathbf{z}_{L})$ and

• 2.

$T_{O}(\mathbf{z}_{R})=T_{U}(\mathbf{z}_{L})$.

Then if one concatenates sequences $\mathbf{z}_{L}$ and $\mathbf{z}_{R}$, for the corresponding random walk it holds that $T_{U}(\mathbf{z}_{L}\mathbf{z}_{R})=T_{O}(\mathbf{z}_{L}\mathbf{z}_{R})$.

The following is an informal explanation of the procedure $\textsf{DyckPaths}(N,(z_{1},\ldots,z_{2^{N-1}}),rng)$ (the formal description is provided by Algorithm 4).

The walk corresponding to $\mathbf{z}_{L}$ is given by $S_{0}:=0,S_{k}=\sum_{i=1}^{k}(2z_{i}-1),k=1,\ldots,2^{N-1}$. 2. 2.

The sequence $D_{k}$ is defined as: $D_{k}=\mathbbm{1}(S_{k}>0\vee S_{k-1}>0)$, for $k=1,\ldots,2^{N-1}$. 3. 3.

The set of points where the walk changes its sign is defined as $R:=\{i:D_{i+1}\neq D_{i}\}\cup\{2^{N-1}\},i=1,\ldots,2^{N-1}-1$. 4. 4.

Elements of $R$ are sorted in increasing order (obtaining $(r_{1},\ldots,r_{w})$). 5. 5.

The sequence $\{l_{i}\}$ is defined as: $l_{1}:=1$, and the next “left-ends” as $l_{i}=r_{i-1}+1$ (for $i=2,\ldots,w$). 6. 6.

The set $\{O_{i}\}$ is defined as $O_{i}:=\{l_{i},l_{1}+1,\ldots,r_{i}\}$, for $i=1,\ldots,w$. 7. 7.

Bits $2^{N-1}+1,\ldots,2^{N}$ ($\mathbf{z}_{R}$) are chosen so that the whole walk spends the same number of steps over and under $x$-axis. Dyck’s paths are generated222The definitions of $l_{i}$ and $r_{i}$ imply that $|O_{i}|$ is even, $i=1,\ldots,w$: $DP_{i}=\textsf{sampleDyckPath}(|O_{i}|/2,rng,h(seed,D_{r_{i}},i),D_{r_{i}})$, for $i=1,\ldots,w$, for some hash function $h$. We use here a hash function $h$ to obtain differently sampled Dyck’s paths. This is achieved by re-seeding the generator $rng$ to be dependent on:

• (a)

the current $seed$,

• (b)

a single bit equal to $D_{r_{i}}$ which corresponds to the type of the sequence one wants to get (over or under the $x$-axis),

• (c)

the path number. 8. 8.

A relative ordering of the paths (generated in the previous step) is obtained from a permutation $\rho\leftarrow\textsf{RandPerm}(w,rng(seed))$. 9. 9.

The resulting bits are obtained by concatenating permuted Dyck’s paths.

Example 4.9** (DyckPaths).**

Let input to $\mathsf{DyckPaths}$ be $z_{\mathbf{L}}=(z_{1},\ldots,z_{16})=(1,1,0,0,0,0,0,1,1,1,1,0,1,0,0,1)$. Then $(r_{1},r_{2},r_{3},r_{4})=(4,10,14,16)$ and $(l_{1},l_{2},l_{3},l_{4})=(1,5,11,15)$ and thus $O_{1}=\{1,2,3,4\},O_{2}=\{5,6,7,8,9,10\},O_{3}=\{11,12,13,14\},O_{4}=\{15,16\}$.

Let the output of sampleDyckPath (called in lines 7-9 of DyckPaths) are $DP_{1}=(0,1,0,1)$, $DP_{2}=(1,1,0,1,0,0)$, $DP_{3}=(0,0,1,1)$, $DP_{4}=(1,0)$.

Let

[TABLE]

Then $(z_{17},\ldots,z_{32})=DP_{\rho(1)}DP_{\rho(2)}DP_{\rho(3)}DP_{\rho(4)}=(0,0,1,1,1,1,0,1,0,0,0,1,0,1,1,0)$.

Ten sample trajectories of the Flawed generator (all generated by the Dyck path-based part of Algorithm 3) are depicted in Figure 3 (the instance of $\textsf{Flawed}_{rng,N,\tau}$ was initialized with the following parameters: $N={18}$, rng – the Mersenne Twister).

5 Experimental results

In this section we briefly report our experimental results of testing some widely used PRNGs implemented in standard libraries in various programming languages. We have applied the ASIN test to different generators including the implementations of the standard C/C++ linear congruential generators, the standard generator rand from the GNU C Library, the Mersenne Twister, the Minstd and the CMRG generator. As our last example we show the results of testing the Flawed generator. Flawed is identified by our ASIN test as non-random, whereas it passed many other tests, including all closely related procedures (swalk_RandomWalk1 test from TestU01 with statistics: H, M, J, R, C, see Table 2).

Each considered PRNG was tested by generating $m=10000$ sequences of length $n\in\{2^{26},\,2^{30},\,2^{34}\}$, using the partition $\mathcal{P}^{asin}_{40}$, i.e., $s=40$. For these parameters our second level test is reliable (see Section 3.2.2) – $\sigma$, the expected statistical deviation of the ratio of $p_{j}$ values found in a given bin is greater than $\sqrt{{s-1\over s^{2}m}}=0.0015$, what significantly exceeds the maximum propagated error $\Delta={2C\over n}$, i.e., (6) holds. Note that for $s=40$ the inequality (7) yields:

$\Delta=7.0703\cdot 10^{-8},m\leq 4.8760\cdot 10^{12}$ for $n=2^{34}$,

• 2.

$\Delta=0.00000113,m\leq 1.9047\cdot 10^{10}$ for $n=2^{30}$,

• 3.

$\Delta=0.000018$, $m\leq 7.44027\cdot 10^{7}$ for $n=2^{26}$.

In the experiments we used our custom implementations of tested PRNGs (except the Mersenne Twister). We used 64-bit version of C++11 implementation of the Mersenne Twister, i.e., the class std::mt19937_64, which is, however, known to have some problems [34]. The generators were initialized with random seeds from http://www.random.org [35] and each sequence was generated using different seed.

The results are presented in Table 1. The values indicating that $\mathcal{H}_{0}$ should be rejected (w.r.t. significance level $\alpha=0.0001$, a value suggested by NIST for second level test) are bolded. For $p_{\chi^{2}}$ these are simply the values smaller or equal to $\alpha$. Concerning the values of $d^{asin}_{tv}$, Lemma 3.4 implies that for $\varepsilon\leq\sqrt{20\cdot 41/10000}\leq 0.2862$ we have $\mathbb{P}\left(d_{tv}^{asin}>{\varepsilon/2}\right)\leq 3\exp\left(-400\varepsilon^{2}\right).$ It can be checked that $3exp(-400\varepsilon^{2})\leq 0.0001$ for $\varepsilon\leq 0.1605$, in other words

[TABLE]

i.e., we reject $\mathcal{H}_{0}$ if the value of $d_{tv}^{asin}$ is larger than $0.0802$. Note that for the results in Table 1 either both statistics $d_{tv}^{asin}$ and $p$-value of $T^{asin}$ reject $\mathcal{H}_{0}$ or both accept it.

We have also calculated the swalk_RandomWalk1 statistics from TestU01 for 10000 sequences of length $2^{26}$ of each PRNG. The following parameters for swalk_RandomWalk1 were used: $N=1,n=10000,r=0,s=32,L0=L1=67108864=2^{26}$. The results are given in Table 2 (including the Flawed generator described in Section 4). For each Statistic H, M, J, R and C, the corresponding $p$-values were obtained using the chi-square statistics. For convenience, $p$-values of $T^{asin}$ are also included in the Table 2 (in the column $p_{\chi^{2}}$).

Our ASIN test would reject the MS Visual C++ PRNG and the Minstd with a multiplier 48271 (The Minstd with a multiplier 16807 gave similar results - not reported here) as good PRNGs. Note that this is indicated by both $p_{\chi^{2}}$ and the value of $d_{tv}^{asin}$. We also conducted the experiments for the procedure rand from the standard library in the Borland C/C++ (not included here). The outcomes are very akin to those for a standard PRNG in the MS Visual C++. Note that for $n=2^{26}$ none of the $p$-values calculated by the swalk_RandomWalk1 from TestU01 suggests rejecting the hypothesis that the MS Visual C++ PRNG is good, whereas the statistics H and C indicate that there can be some flaws in the Minstd 48271. It is worth mentioning that the MS Visual C++ PRNG passes the NIST Test Suite [8], as pointed out in [18]. Minstd, despite its weaknesses, became a part the C++11 standard library. It is implemented by the classes std::minstd_rand0 (with the multiplier 16807) and std::minstd_rand (with the multiplier 48271). Concerning the GNU C and the MT19937-64 – as can be seen in both Table 1 and Table 2 – they can be both considered as good. It is worth mentioning that the results for the CMRG generator (not reported here) were similar to those for the MT19937-64.

The open source code of our implementation is publicly available, see [19] (it includes the Flawed PRNG as well as the Law of Iterated Logarithm test from [18]).

Discussion on the influence of the parameter $\tau$ of the Flawed PRNG on the statistic $T^{asin}$

Recall that in the Algorithm 3 the parameter $\tau$ corresponds to a fraction of simulations which are exactly half of the time above and half of the time below $x$-axis, i.e., we have $S_{n,j}^{asin}=0.5$ for $\lfloor\tau m\rfloor$ simulations. Note that the $p_{j}$-value is then also equal to $0.5$. The remaining $\lceil(1-\tau)m\rceil$ simulations come from the $rng$. Let us assume that the $rng$ returns truly random numbers.

Concerning $T^{asin}$ statistic, we have $E_{1}=E_{s+1}={m\over 2s}$ and $E_{i}={m\over s},i=2,\ldots,s$. Set $r:=\lceil{s\over 2}\rceil+1$. For an “ideal“ $rng$ we would have $O_{1}=O_{s+1}={m(1-\tau)\over 2s},O_{k}={m(1-\tau)\over s},k\in\{2,\ldots,s\}\setminus\{r\}$ and $O_{r}={m(1-\tau)\over s}+\tau m$. Thus,

[TABLE]

For the parameters $m=10000,\tau={1\over 66},s=40$ we have $T^{asin}=89.5316$ and the corresponding $p$-value is less than $0.0000118$. It means that even for an $rng$ producing truly random numbers, the ASIN test should recognize the Flawed generator as not good.

5.1 Results of TestU01 for Flawed

We have run several general-purpose tests against the Flawed generator. For SmallCrush all 15 out of 15 tests were passed. For the Mersenne Twister (MT) and the $\textsf{Flawed}_{MT,30,1/66}$ we run BigCrush. Tests for which generators failed are presented in the Table 3.

6 Notes on Takashima’s method for testing PRNGs and the arcsine test implementation from TestU01

The idea of using the arcsine law for developing statistical tests for an empirical evaluation of PRNGs was formerly proposed by Takashima in [36, 37, 38]. In this series of articles, test statistics based on the arcsine law were applied for assessing the randomness of the output of maximum-length linearly recurring sequences ($m$-sequences in short). The experimental results presented there clearly show that the bits produced by this family of PRNGs are biased. Besides revealing the weakness of $m$-sequences, these outcomes have also proved that Takashima’s tests are effective methods, worth applying in practice.

The approach introduced in [36, 37] can be briefly described as follows. After an initialization of a PRNG, a sequence of $2\mathfrak{n}m$ bits is generated and divided into $m$ subsequences of length $n=2\mathfrak{n}$. Then, each subsequence is used for constructing a random walk. For each of these $m$ sample random walks, the value of a test statistic based on the arcsine law is calculated. The investigated statistic, called in [36, 37] the sojourn time – denote them by $t_{n}^{j},j=0,\ldots,m-1$ – is the time spent by a random walk above the $x$-axis. From $m$ realizations of this statistic an empirical distribution of the sojourn time $f_{2i}=|\{j:t^{j}_{n}=2i\}|$, $i=0,\ldots,\mathfrak{n}$ is then derived and compared with its theoretical distribution via a chi-square test. The whole procedure is repeated $\lambda\geq 1$ times, yielding a set of $\chi^{2}$ test statistics’ values $\{\chi^{2}_{k}\},k=0,\ldots,\lambda-1$. The final step of the Takashima’s testing method is to count the number of $\chi_{k}^{2}$ values falling between $90$-th and $95$-th percentile and those bigger than $95$-th percentile of a respective $\chi^{2}$ distribution. These two counts are then the basis for deciding if $\mathcal{H}_{0}$ should be rejected. Note that for $\lambda>1$ this is a third level test, which in general is not reliable, as shown in [3].

The author in [37] considers also a slightly modified variant of the procedure, where the chi-square test is combined with the Kolomogorov-Smirnov test. Another method, presented in [38], exploits the relations between the sojourn time and the last visit time for one-dimensional random walks.

It is worth noting that in our simulations we used binary sequences of length at least $n=2^{26}$. Thus, a direct application of the Takashima’s methods from [36, 37] would require large amount of additional memory to store the values of $f_{2i},i=0,\ldots,\mathfrak{n}$.

As the arcsine law based statistical tests were proven to be useful in detecting flaws of some PRNGs, such procedures were implemented in the TestU01 library (see [7]). This tool, developed by L’Ecuyer and Simard, provides a big variety of functions for empirical examining of PRNGs. One of the test modules, swalk, contains a procedure swalk_RandomWalk1, which calculates a bunch of test statistics for a sample of $m$ random walks constructed from chosen bits of generated binary sequences. Among them, there is the Statistic J, which implements the test based on the arcsine law. This procedure is similar to ours. Namely, $m$ calculated values of the test statistic are grouped according to some partition and their empirical distribution is compared with the theoretical one by means of the chi-square test. The main difference is that in our testing method the partition size $s$ is a parameter chosen by the user, whereas the partition used by swalk_RandomWalk1 is calculated automatically, depending on the tested sequence. Moreover, we provide bounds on approximation errors in the computations of $p$-values (a Berry-Essen type inequality), assuring the reliability of the whole testing procedure.

7 Conclusions

In this paper we analyzed a method for testing PRNGs based on the arcsine law for random walks. Our procedure is a second level statistical test. We also provided a detailed error analysis of the proposed method. The approximation errors in the calculation of $p$-values are bounded by a Berry-Essen type inequality, what allows to control the overall error, assuring the reliability of the test. We evaluate the quality of PRNGs via the chi-square statistics as well as by calculating a statistical distance (the total variation distance) between the empirical distribution of the considered characteristic for generated pseudorandom output and its theoretical distribution for truly random binary sequences.

The experimental results presented in this paper show that our testing procedure can be used for detecting weaknesses in many common PRNGs’ implementations. Likewise the Law of Iterated Logarithm test from [18], the ASIN test has also revealed some flaws and regularities in generated sequences not necessarily being identified by other current state of the art tools like the NIST SP800-22 Testing Suite or TestU01. Thus, these kind of testing techniques seem to be very promising, as they allow also for recognition of different kinds of deviations from those detected by existing tools. Nevertheless, like other statistical tests, the ASIN test is not universal and encompasses only one from an immense range of characteristics of random bit strings and does not capture all known flaws. Therefore, the testing procedures relying on properties of random walks like the ASIN test should be used along with other tests for more careful assessment of pseudorandom generators. This issue is well depicted by the provided example of obviously non-random generator Flawed for which the LIL test has failed to detect its weaknesses, but the ASIN test has turned out to be very sensitive for that kind of deviations. Hence, an important line of further research is to develop another novel tests utilizing various properties of random walks. Such tests, when combined together, should be capable of detecting more hidden dependencies between the consecutive bits in the sequences generated by PRNGs. This could lead to designing more robust test suites for evaluating the quality of random numbers generated by the new implementations of PRNGs as well as those being already in use, especially for cryptographic purposes.

Acknowledgements

We would like to thank the anonymous reviewers whose suggestions and insightful comments helped significantly improve and clarify this manuscript. In particular we thank one of the reviewers for pointing out the article [4] on second level tests.