Sparse Serial Tests of Uniformity for Random Number Generators.
Different versions of the  serial test for testing the uniformity and independence of vectors of successive values produced by a (pseudo)random number generator are studied. These tests partition the t-dimensional unit hypercube into k cubic cells of equal volume, generate n points (vectors) in this hypercube, count how many points fall in each cell, and compute a test statistic defined as the sum of values of some univariate function f applied to these k individual counters. Both overlapping and nonoverlapping vectors are considered. For different families of generators, such as linear congruential, Tausworthe, nonlinear inversive, etc., different ways of choosing these functions and of choosing k are compared, and formulas are obtained for the (estimated) sample size required to reject the null hypothesis of independent  uniform random variables, as a function of the period length of the generator. For the classes of alternatives that correspond to linear generators, the most efficient tests turn out to have $k \gg n$ (in contrast to what is usually done or recommended in simulation books) and to use overlapping vectors.
Introduction
. The aim of this paper is to examine certain types of serial
tests for testing the uniformity and independence of the output sequence of general-purpose
uniform random number generators (RNGs) such as those found in software
libraries. These RNGs are supposed to produce \imitations" of mutually independent
random variables uniformly distributed over the interval [0; 1) (i.i.d. U(0; 1), for short).
Testing an RNG whose output sequence is U amounts to testing the null
hypothesis are i.i.d. U(0; 1)."
To approximate this multidimensional uniformity, good RNGs are usually designed
(theoretically) so that the multiset t of all vectors rst
successive output values, from all possible initial seeds, covers the t-dimensional
unit hypercube [0; 1) t very evenly, at least for t up to some t 0 , where t 0 is chosen
somewhere between 5 and 50 or so. When the initial seed is chosen randomly, this
t can be viewed in some sense as the sample space from which points are chosen at
random to approximate the uniform distribution over [0; 1) t . For more background
on the construction of RNGs, see, for example, [13, 17, 21, 35].
For large t, the structure of t is typically hard to analyze theoretically. Moreover,
even for a small t, one would often generate several successive t-dimensional vectors
of the form (u statistical testing then comes into play
because the dependence structure of these vectors is hard to analyze theoretically. An
excessive regularity of t implies that statistical tests should fail when their sample
P. L'Ecuyer and R. Simard, Departement d'Informatique et de Recherche Operationnelle,
Universite de Montreal, C.P. 6128, Succ. Centre-Ville, Montreal, H3C 3J7, Canada. e-
mail: lecuyer@iro.umontreal.ca and simardr@iro.umontreal.ca. S. Wegenkittl, Institute of
Mathematics, University of Salzburg, Hellbrunnerstrasse 34, A-5020 Salzburg, Austria, e-mail:
ste@random.mat.sbg.ac.at This work has been supported by the National Science and Engineering
Research Council of Canada grants # ODGP0110050 and SMF0169893, by FCAR-Quebec grant #
93ER1654, and by the Austrian Science Fund FWF, project no. P11143-MAT. Most of it was performed
while the rst author was visiting Salzburg University and North Carolina State University,
in 1997-98 (thanks to Peter Hellekalek and James R. Wilson).
sizes approach the period length of the generator. But how close to the period length
can one get before trouble begins?
Several goodness-of-t tests for H 0 have been proposed and studied in the past
(see, e.g., [13, 9, 26, 41] and references therein). Statistical tests can never certify for
good an RNG. Dierent types of tests detect dierent types of deciencies and the
more diversied is the available battery of tests, the better.
A simple and widely used test for RNGs is the serial test [1, 6, 8, 13], which
operates as follows. Partition the interval [0; 1) into d equal segments. This determines
a partition of [0; 1) t into cubic cells of equal size. Generate nt random
numbers U construct the points V
and let X j be the number of these points falling into cell j, for
has the multinomial distribution with parameters
1=k). The usual version of the test, as described for example in [6, 13, 14]
among other places, is based on Pearson's chi-square statistic
is the average number of points per cell, and the distribution of X 2
under H 0 is approximated by the chi-square distribution with k 1 degrees of freedom
when   5 (say).
In this paper, we consider test statistics of the general form
where f n;k is a real-valued function which may depend on n and k. We are interested
for instance in the power divergence statistic
real-valued parameter (by we understand the limit as
could also consider - seems unnecessary
in the context of this paper. Note that D . The power divergence statistic
is studied in [39] and other references given there. A more general class is the '-
divergence family, where f n;k (X Other forms of
f n;k that we consider are f n;k (where I denotes the indicator function),
which the corresponding Y is the
number of cells with at least b points, the number of empty cells, and the number of
collisions, respectively.
We are interested not only in the dense case, where  > 1, but also in the sparse
case, where  is small, sometimes much smaller than 1. We also consider (circular)
overlapping versions of these statistics, where U replaced
by V i .
In a slightly modied setup, the constant n is replaced by a Poisson random
variable  with mean n. Then, (X a vector of i.i.d. Poisson random
variables with mean  instead of a multinomial vector, and the distribution of Y
becomes easier to analyze because of this i.i.d. property. For large k and n, however,
the dierence between the two setups is practically negligible, and our experiments
are with
A rst-order test observes the value of Y , say y, and rejects H 0 if the p-value
is much too close to either 0 or 1. The function f is usually chosen so that p too
close to 0 means that the points tend to concentrate in certain cells and avoid the
others, whereas p close to 1 means that they are distributed in the cells with excessive
uniformity. So p can be viewed as a measure of uniformity, and is approximately a
random variable under H 0 if the distribution of Y is approximately continuous.
A second-order (or two-level) test would obtain N \independent" copies of Y ,
say is the theoretical distribution
of Y under H 0 , and compare their empirical distribution to the uniform. Such a
two-level procedure is widely applied when testing RNGs (see [6, 13, 16, 29, 30]).
Its main supporting arguments are that it tests the RNG sequence not only at the
global level but also at a local level (i.e., there could be bad behavior over short
subsequences which \cancels out" over larger subsequences), and that it permits one
to apply certain tests with a larger total sample size (for example, the memory size of
the computer limits the values of n and/or k in the serial test, but the total sample
size can exceed n by taking N > 1). Our extensive empirical investigations indicate
that for a xed total sample size Nn, when testing RNGs, a test with
typically more e-cient than the corresponding test with N > 1. This means that
for typical RNGs, the type of structure found in one (reasonably long) subsequence
is usually found in (practically) all subsequences of the same length. In other words,
when an RNG started from a given seed fails spectacularly a certain test, it usually
fails that test for most admissible seeds.
The common way of applying serial tests to RNGs is to select a few specic
generators and some arbitrarily chosen test parameters, run the tests, and check if
H 0 is rejected or not. Our aim in this paper is to examine in a more systematic
way the interaction between the serial tests and certain families of RNGs. From each
family, we take an RNG with period length near 2 e , chosen on the basis of the usual
theoretical criteria, for integers e ranging from 10 to 40 or so. We then examine,
for dierent ways of choosing k and constructing the points V i , how the p-value of
the test evolves as a function of the sample size n. The typical behavior is that
takes \reasonable" values for a while, say for n up to some threshold n 0 , then
converges to 0 or 1 exponentially fast with n. Our main interest is to examine the
relationship between n 0 and e. We adjust (crudely) a regression model of the form
log
e++ where
and  are two constants and  is a small noise. The result
gives an idea of what size (or period length) of RNG is required, within a given family,
to be safe with respect to these serial tests for the sample sizes that are practically
feasible on current computers. It turns out that for popular families of RNGs such as
the linear congruential, multiple recursive, and shift-register, the most sensitive tests
choose k proportional to 2 e and yield
which means that n 0 is
a few times the square root of the RNG's period length.
The results depend of course on the choice of f in (1.2) and on how d and t are
chosen. For example, for linear congruential generators (LCGs) selected on the basis
of the spectral test [6, 13, 24], the serial test is most sensitive when k  2 e , in which
case
k). These \most e-cient" tests are very sparse (  1). Such large
values of k yield more sensitive tests than the usual ones (for which k  2 e and
5 or so) because the excessive regularity of LCGs really shows up at that level
of partitioning. For k  2 e , the partition eventually becomes so ne that each cell
contains either 0 or 1 point, and the test loses all of its sensitivity.
For xed n, the non-overlapping test is typically slightly more e-cient than the
overlapping one, because it relies on a larger amount of independent information.
However, the dierence is typically almost negligible (see Section 5.3) and the non-overlapping
test requires t times more random numbers. If we x the total number
of U i 's that are used, so the non-overlapping test is based on n points whereas the
overlapping one is based on nt points, for example, then the overlapping test is typically
more e-cient. It is also more costly to compute and its distribution is generally
more complicated. If we compare the two tests for a xed computing budget, the
overlapping one has an advantage when t is large and when the time to generate the
random numbers is an important fraction of the total CPU time to apply the test.
In Section 2, we collect some results on the asymptotic distribution of Y for
the dense case where k is xed and n ! 1, the sparse case where both
and the very sparse case where n=k ! 0. In
Section 3 we do the same for the overlapping setup. In Section 4 we brie
y discuss the
e-ciency of these statistics for certain classes of alternatives. Systematic experiments
with these tests and certain families of RNGs are reported in Section 5. In Section 6,
we apply the tests to a short list of RNGs proposed in the literature or available in
software libraries and widely used. Most of these generators fail miserably. However,
several recently proposed RNGs are robust enough to pass all these tests, at least for
practically feasible sample sizes.
2. Power Divergence Test Statistics for Non-Overlapping Vectors. We
brie
y discuss some choices of f n;k in (1.2) which correspond to previously introduced
tests. We then provide formulas for the exact mean and variance, and limit theorems
for the dense and sparse cases.
2.1. Choices of f n;k . Some choices of f n;k are given in Table 2.1. In each case,
Y is a measure of clustering: It tends to increase when the points are less evenly
distributed between the cells. The well-known Pearson and loglikelihood statistics,
are both special cases of the power divergence, with
respectively [39]. H is related to G 2 via the relation 2). The
statistics N b , W b , and C count the number of cells that contain exactly b points (for
b  0), the number of cells that contain at least b points (for b  1), and the number
of collisions (i.e., the number of times a point falls in a cell that already has a point in
respectively. They are related by N
and
2.2. Mean and Variance. Before looking at the distribution of Y , we give
expressions for computing its exact mean and variance under H 0 .
If the number of points is xed at n, Denoting
one obtains after some algebraic manipulations:
x

Table
Some choices of f n;k and the corresponding statistics.
divergence
loglikelihood
negative entropy
number of cells with exactly b points
number of cells with at least b points
number of empty cells
number of collisions
x
x
x
min(n
x
y
(f(x) )(f(y)
Although containing a lot of summands, these formulas are practical in the sparse
case since for the Y 's dened in Table 2.1, when n and k are large and
small, only the terms for small x and y in the above sums are non-negligible. These
terms converge to 0 exponentially fast as a function of x y, when x
rst two moments of Y are then easy to compute by truncating the sums after a small
number of terms. For example, with 1000, the relative errors on E[H ] and
are less than the sums are stopped at of 1000, and
if the sums are stopped at similar behavior is observed
for the other statistics.
The expressions (2.1) and (2.2) are still valid in the dense case, but for larger ,
more terms need to be considered. Approximations for the mean and variance of D -
when   1, with error terms in o(1=n), are provided in [39], Chapter 5, page 65.
In the Poisson setup, where n is the mean of a Poisson random variable, the X j
are i.i.d. Poisson() and the expressions become
2.3. Limit Theorems. The limiting distribution of D - is a chi-square in the
dense case and a normal in the sparse case. Two-moment-corrected versions of these
6 PIERRE L'ECUYER, RICHARD SIMARD, AND STEFAN WEGENKITTL
results are stated in the next proposition. This means that D (C)
- and D (N)
- in the
proposition have exactly the same mean and variance as their asymptotic distribution
(e.g., 0 and 1 in the normal case). Read and Cressie [39] recommend this type
of standardization, which tends to be closer to the asymptotic distribution than a
standardization by the asymptotic mean and variance. The two-moment corrections
become increasingly important when - gets away from around 1. The mean and
variance of D - can be computed as explained in the previous subsection. Another
possibility would be to correct the distribution itself, e.g., using Edgeworth-type expansions
[39], page 68. This gives extremely complicated expressions, due in part to
the discrete nature of the multinomial distribution, and the gain is small.
Proposition 2.1. For - > 1, the following holds under H 0 .
(i) [Dense case] If k is xed and n !1, in the multinomial setup
convergence in distribution, and
is the chi-square distribution with k 1 degrees of freedom. In the
Poisson setup, D (C)
(ii) [Sparse case] For both the multinomial and Poisson setups, if
is the standard normal distribution.
Proof. For the multinomial setup, part (i) can be found in [39], page 46, whereas
part (ii) follows from Theorem 1 of [11], by noting that all the X j 's here have the same
distribution. The proofs simplify for the Poisson setup, due to the independence. The
p n=k are i.i.d. and asymptotically N(0; 1) in the dense case, so
their sum of squares, which is X 2 , is asymptotically  2 (k).
We now turn to the counting random variables N b , W b , and C. These are not
approximately chi-square in the dense case. In fact, if xed k, it is clear
that xed b. This implies that W b ! k and
random variables are all degenerate.
For the Poisson setup, each X i is Poisson(), so p b
b  0 and N b is BN(k; p b ), a binomial with parameters k and p b . If k is large and p b
is small, N b is thus approximately Poisson with (exact) mean
for b  0:
The next result covers other cases as well.
Proposition 2.2. For the Poisson or the multinomial setup, under H 0 , suppose
that k !1 and n !1, and let 1 ,
positive constants.
one also has C
(ii) For
Proof. In (i), since 0, one has for the Poisson case E[N b+1 ]=E[N b
0. The relative contribution of W b+1 to the sum W
(a sum of correlated Poisson random variables) is then negligible compared
with that of N b , so N b and W b have the same asymptotic distribution (this follows
from Lemma 6.2.2 of [2]). Likewise, under these conditions with 2, C has the
same asymptotic distribution as N 2 , because
e  1. For the multinomial setup, it has been shown (see [2], Section 6.2) that N b
and W b , for b  2, are asymptotically Poisson(kp b ) when  ! 0, the same as for the
Poisson setup. The same argument as for W 2 applies for C, using again their Lemma
6.2.2, and this proves (i). For for the Poisson setup, we saw already that N 0
is asymptotically
For the multinomial case, the same result follows from Theorem 6.D of [2], and this
proves (ii). Part (iii) is obtained by applying Theorem 1 of [11].
The exact distributions of C and N 0 under H 0 , for the multinomial setup, are
given by
where the
are the Stirling numbers of the second kind (see [13], page 71,
where an algorithm is also given to compute all the non-negligible probabilities in
time O(n log n)).
In our implementation of the test based on C, we used the Poisson approximation
for   1=32, the normal approximation for  > 1=32 and n > 2 15 , and the exact
distribution otherwise.
3. Overlapping vectors. For the overlapping case, let X (o)
t;j be the number of
overlapping vectors j. Now, the formulas (2.1)
and (2.2) for the mean and variance, and the limit theorems in Propositions 2.1 and
2.2, no longer stand. The analysis is more di-cult than for the disjoint case because
in general P [X (o)
depends on i and P [X (o)
depends on the pair
(i; in a non-trivial way.
Theoretical results have been available in the overlapping multinomial setup, for
the Pearson statistic in the dense case. Let
and let X 2
be the equivalent of X 2
for the overlapping vectors of dimension t 1:
Consider the statistic ~
Good [8] has shown
that E[X 2
exactly (see his Eq. (5) and top of page 280) and that when
page 284). This setup, usually with
n=k  5 or so, is called the overlapping serial test or the m-tuple test in the literature
and has been used previously to test RNGs (e.g., [1, 29, 30]). The next proposition
generalizes the result of Good to the power divergence statistic in the dense case.
Further generalization is given by Theorem 4.2 of [43].
Proposition 3.1. Let
the power divergence statistic for the t-dimensional overlapping vectors, and dene
~
in the multinomial setup, if - > 1, k is xed,
and n !1, ~
Proof. The result is well-known for 1. Moreover, a Taylor series expansion
of D -;(t) in powers of X (o)
easily shows that D
probability as Theorem A6.1). Therefore, ~
D -;(t) has
the same asymptotic distribution as ~
D 1;(t) and this completes the proof.
For the sparse case, where k; our
simulation experiments support the conjecture that
~
~
The overlapping empty-cells-count test has been discussed in a heuristic way in a
few papers. For calls it the overlapping pairs sparse occupancy
(OPSO) and suggests a few specic parameters, without providing the underlying the-
ory. Marsaglia and Zaman [32] speculate that N 0 should be approximately normally
distributed with mean ke  and variance ke  (1 3e  ). This make sense only if  is
not too large or not too close to zero. We studied empirically this approximation and
found it reasonably accurate only for 2    5 (approximately). The approximation
could certainly be improved by rening the variance formula.
Proposition 2.2 (i) and (ii) should hold in the overlapping case as well. Our
simulation experiments indicate that the Poisson approximation for C is very accurate
for (say)  < 1=32, and already quite good for   1, when n is large.
4. Which Test Statistic and What to Expect?. The LFSR, LCG, and MRG
generators in our lists are constructed so that their point sets t over the entire period
are superuniformly distributed. Thus, we may be afraid, if k is large enough, that
very few cells (if any) contain more than 1 point and that D - , C, N 0 , N b and W b
for b  2 are smaller than expected. In the extreme case where assuming
that the distribution of C under H 0 is approximately Poisson with mean n 2 =(2k),
the left p-value of the collision test is . For a xed
number of cells, this p-value approaches 0 exponentially fast in the square of the
sample size n. For example,
k, and 16
k, respectively. Assuming that k is near the RNG's period length, i.e.,
means that the test starts to fail abruptly when the sample size exceeds
approximately 4 times the square root of the period length. As we shall see, this is
precisely what happens for certain popular classes of generators. If we use the statistic
b instead of C, in the same situation, we have
and the sample size required to obtain a p-value less than a xed (small) constant is
2. In this setup, C and N 2 are equivalent to W 2 , and choosing
e-cient test.
Suppose now that we have the opposite: Too many collisions. One simple model
of this situation is the alternative are i.i.d. uniformly distributed over
boxes, the other k k 1 boxes being always empty." Under H 1 , W b is approximately
Poisson with mean
is large and  1 is small) instead of
Therefore, for a given  0 , and x 0 such that
the power of the test at level  0 is
where x 0 depends on b. When b increases, for a xed  0 , x 0 decreases and  1 decreases
as well if n=k 1 maximizes the power unless n=k 1 is large. In fact
the test can have signicant power only if  1 exceeds a few units (otherwise, with
large probability, one has W not rejected). This means
which can be approximated by O(k (b 1)=b
is reasonably large. Then, is the best choice. If k 1 is small,  1 is maximized
(approximately) by taking
The alternative H 1 just discussed can be generalized as follows: Suppose that the
have a probability larger than 1=k, while the other k k 1 cells have a smaller
probability. H 1 is called a hole (resp., peak , split) alternative if k 1 =k is near 1 (resp.,
near 0, near 1/2). We made extensive numerical experiments regarding the power of
the tests under these alternatives and found the following. Hole alternatives can be
detected only when n=k is reasonably large (dense case), because in the sparse case
one expects several empty cells anyway. The best test statistics to detect them are
those based on the number of empty cells N 0 , and D - with - as small as possible (e.g.,
0). For a peak alternative, the power of D - increases with - as a concave
function, with a rate of increase that typically becomes very small for - larger than
3 or 4 (or higher, if the peak is very narrow). The other test statistics in Table 2.1
are usually not competitive with D 4 (say) under this alternative, except for W b which
comes close when b  n=k 1 (however it is hard to choose the right b because k 1 is
generally unknown). The split alternative with the probability of the k k 1 low-probability
cells equal to 0 is easy to detect and the collision test (using C or W 2 ) is
our recommendation. The power of D - is essentially the same as that of C and W 2 ,
for most -, because E[W 3 ] has a negligible value, which implies that there is almost a
one-to-one correspondence between C, W 2 , and D - . However, with the small n that
su-ces for detection in this situation, E[W 2 ] is small and the distribution of D - is
concentrated on a small number of values, so neither the normal nor the chi-square
is a good approximation of its distribution. Of course, the power of the test would
improve if the high-probability cells were aggregated into a smaller number of cells,
and similarly for the low-probability cells. But to do this, one needs to know where
these cells are a priori .
These observations extend (and agree with) those made previously by several
authors (see [39] and references therein), who already noted that for D - , the power
decreases with - for a hole alternative and increases with - for a peak alternative.
This implies in particular that G 2 and H are better [worse] test statistics than X 2 to
detect a hole [a peak]. In the case of a split alternative for which the cell probabilities
are only slightly perturbed, X 2 is optimal in terms of Pitman's asymptotic e-ciency
is optimal in terms of Bahadur's e-ciency (see [39] for details).
5. Empirical Evaluation for RNG Families.
5.1. Selected Families of RNGs. We now report systematic experiments to
assess the eectiveness of serial tests for detecting the regularities in specic families
of small RNGs. The RNG families that we consider are named LFSR3, GoodLCG,
BadLCG2, MRG2, CombL2, InvExpl. Within each family, we constructed a list of
specic RNG instances, with period lengths near 2 e for (integer) values of e ranging
from 10 to 40. These RNGs are too small to be considered for serious general purpose
softwares, but their study gives good indication about the behavior of larger instances
from the same families. At step n, a generator outputs a number un 2 [0; 1).
The LFSR3s are combined linear feedback shift register (LFSR) (or Tausworthe)
generators with three components, of the form
where  means bitwise exclusive-or, and are constant parameters
selected so that the k j are reasonably close to each other, and the sequence fung
has period length (2 k1 1)(2 k2 1)(2 k3 1) and is maximally equidistributed (see
[19] for the denition and further details about these generators).
The GoodLCGs are linear congruential generators (LCGs), of the form
where m is a prime near 2 e and a is selected so that the period length is m 1 and
so that the LCG has an excellent behavior with respect to the spectral test (i.e., an
excellent lattice structure) in up to at least 8 dimensions. The BadLCG2s have the
same structure, except that their a is chosen so that they have a mediocre lattice
structure in 2 dimensions. More details and the values of a and m can be found in
[24, 26]. The MRG2 are multiple recursive generators of order 2, of the form
period length m 2 1, and excellent lattice structure as for the GoodLCGs [17, 21].
The CombL2s combine two LCGs as proposed in [15]:
so that the combined generator has period length (m 1 1)(m 2 1)=2 and an excellent
lattice structure (see [28] for details about that lattice structure).
InvExpl denotes a family of explicit inversive nonlinear generators of period length
m, dened by
where m is prime and (an) 1 mod
5.2. The Log-p-values. For a given test statistic Y taking value y, let
the log-p-value of the test as
For example, means that the right p-value is between 0.01 and 0.001. For a given
class of RNGs, given Y , t, and a way of choosing k, we apply the test for dierent
values of e and with sample size
e+ , for where the
constant
is chosen so that the test starts to fail at approximately the same value of
for all (or most) e. More specically, we dene ~  (resp.   ) as the smallest values
of  for which the absolute log-p-value satises j'j  2 (resp. j'j  14) for a majority
of values of e. These thresholds are arbitrary.
5.3. Test Results: Examples and Summary. Tables 5.1 and 5.2 give the log-
p-values for the collision test applied to the GoodLCGs and BadLCG2s, respectively,
in dimensions, with . Only the log-
p-values ' outside of the set f1; 0; 1g, which correspond to p-values less than 0:01,
are displayed. The symbols and ! mean '  14 and '  14, respectively. The
columns not shown are mostly blank on the left of the table and lled with arrows on
the right of the table. The small p-values appear with striking regularity, at about
the same  for all e, in each of these tables. This is also true for other values of e not
shown in the table. One has ~

Table

5.1, while ~
in

Table

5.2. The GoodLCGs fail because their structure is too regular (the left p-values
are too small because there are too few collisions), whereas the BadLCG2s have
the opposite behavior (the right p-values are too small because there are too many
collisions; their behavior correspond to the split alternative described in Section 4).

Table

5.3 gives the values of ~
and   for the selected RNG families, for the
collision test in 2 and 4 dimensions. All families, except InvExpl, fail at a sample size
proportional to the square root of the period length . At
1=2 , the left or
right p-value is less than 10 14 most of the time. The BadLCG2s in 2 dimensions are
the rst to fail: They were chosen to be particularly mediocre in 2 dimensions and the
test detects it. Apart from the BadLCG2s, the generators always fail the tests due to
excessive regularity. For the GoodLCGs and LFSR3s, for example, there was never
a cell with more than 2 points in it. For the LFSR3s, we distinguish two cases: One
where d was chosen always odd and one where it was always the smallest power of 2
such that . In the latter case, the number of collisions is always 0, since
no cell contains more than a single point over the entire period of the generator, as a
consequence of the \maximal equidistribution" property of these generators [19]. The
left p-values then behave as described at the beginning of Section 4. The InvExpl
resist the tests until after their period length is exhausted. These generators have
their point set t \random-looking" instead of very evenly distributed. However,
they are much slower than the linear ones.
We applied the power divergence tests with 4, and in most cases
the p-values were very close to those of the collision test. In fact, when no cell count
which we have observed frequently), there is a one-to-one
correspondence between the values of C and of D - for all - > 1. Therefore, all these
statistics should have similar p-values if both E[W 3 ] and the observed value of W 3 are
small (the very sparse situation). For the overlapping versions of the tests, the values
of
, and   are exactly the same as those given in Table 5.3. This means that the

Table
The log-p-values ' for the GoodLCGs with period length   2 e , for the collision test (based on
C), in cells, and sample size . The table entries give the
values of '. The symbols and ! mean '  14 and '  14, respectively. Here, we have ~
and
22 3 11
26 2 6
28 2
overlapping tests are more e-cient than the non-overlapping ones, because they call
the RNG t times less.
We applied the same tests with smaller and larger numbers of cells, such as
found that ~
and   increase when
moves away from 2 e . A typical example: For the GoodLCGs with
6, 5, and 7 for the four choices of k given above, respectively, whereas
. The classical way of applying the serial test for RNG testing uses a large
average number of points per cell (dense case). We applied the test based on X 2 to
the GoodLCGs, with k  n=8, and found empirically
This means that the required sample size now increases as O( 2=3 ) instead of O( 1=2 )
as before; i.e., the dense setup with the chi-square approximation is much less e-cient
than the sparse setup. We observed the same for D - with other values of - and other
values of t, and a similar behavior for other RNG families.
For the results just described, t was xed and d varied with e. We now x
(i.e., we take the rst two bits of each number) and vary the dimension as

Table

5.4 gives the results of the collision test in this setup. Note the change in
for
the GoodLCGs and BadLCG2s: The tests are less sensitive for these large values of
t.
We also experimented with two-level tests, where a test of sample size n is replicated
times independently. For the collision test, we use the test statistic C T , the
total number of collisions over the N replications, which is approximately Poisson
with mean Nn 2 e n=k =(2k) under H 0 . For the power divergence tests, we use as test
statistics the sum of values of D (N)
- and of D (C)
- , which are approximately N(0; N)
and  2 (N(k 1)) under H 0 , respectively. We observed the following: The power

Table
The log-p-values ' for the collision test, with the same setup as in Table 5.1, but for the
BadLCG2 generators. Here, ~
22
26
28

Table
Collision tests for RNG families, in t dimensions, with k  2 e . Recall that ~
(resp.   ) is the
smallest integer  for which j'j  2 (resp. j'j  14) for a majority of values of e, in tests with sample
e+ .
RNG family
LFSR3, d power of
of a test with (N; n) is typically roughly the same as that of the same test at level
one with sample size n
N . Single-level tests thus need a smaller total
sample size than the two-level tests to achieve the same power. On the other hand,
two-level tests are justied when the sample size n is limited by the memory size of
the computer at hand. (For n  k, the counters X j are implemented via a hashing
14 PIERRE L'ECUYER, RICHARD SIMARD, AND STEFAN WEGENKITTL

Table
Collision tests with divisions in each dimension and dimensions.
Generators
~
CombL2
technique, for which the required memory is proportional to n instead of k). Another
way of doing a two-level test with D - is to compute the p-values for the N replicates
and compare their distribution with the uniform via (say) a Kolmogorov-Smirnov or
Anderson-Darling goodness-of-t test. We experimented extensively with this as well
and found no advantage in terms of e-ciency, for all the RNG families that we tried.
6. What about real-life LCGs?. From the results of the preceding section
one can easily predict, conservatively, at which sample size a specic RNG from a
given family will start to fail. We verify this with a few commonly used RNGs, listed
in

Table

6.1. (Of course, this list is far from exhaustive).

Table
List of selected generators.
LCG1. LCG with
LCG2. LCG with
LCG3. LCG with
LCG4. LCG with
LCG5. LCG with
LCG6. LCG with
LCG7. LCG with
LCG8. LCG with
LCG9. LCG with
RLUX. RANLUX with
WEY1. Nested Weyl with
(see [10]).
WEY2. Shued nested Weyl with
(see [10]).
CLCG4. Combined LCG of [25].
CMRG96. Combined MRG in Fig. 1 of [18].
CMRG99. Combined MRG in Fig. 1 of [23].
Generators LCG1 to LCG9 are well-known LCGs, based on the recurrence x
at step i. LCG1 and LCG2 are recommended
by Fishman [7] and a FORTRAN implementation of LCG1 is given by
Fishman [6]. LCG3 is recommended in [14], among others, and is used in the SIMSCRIPT
II.5 and INSIGHT simulation languages. LCG4 is in numerous software
systems, including the IBM and Macintosh operating systems, the Arena and SLAM
II simulation languages (note: the Arena RNG has been replaced by CMRG99 after
we wrote this paper), MATLAB, the IMSL library (which also provides LCG1 and

Table
The log-p-values for the collision test in cells, and sample size
m.
Generator

Table
The log-p-values for the two-level collision test (based on C T ) in
cells, sample size for each replication, and replications.
Generator
LCG5), the Numerical Recipes [38], etc., and is suggested in several books and papers
(e.g., [3, 36, 40]). LCG6 is used in the VAX/VMS operating system and on Convex
computers. LCG5 and LCG9 are the rand and rand48 functions in the standard
libraries of the C programming language [37]. LCG7 is taken from [6] and LCG8 is
used in the CRAY system library. LCG1 to LCG4 have period length 2 31 2, LCG5,
LCG6, AND LCG9 have period length m, and LCG7 and LCG8 have period length
RLUX is the RANLUX generator implemented by James [12], with luxury level
24. At this luxury level, RANLUX is equivalent to the subtract-with-borrow
generator with modulus 43 and proposed in [31] and
used, for example, in MATHEMATICA (according to its documentation). WEY1 is
a generator based on the nested Weyl sequence dened by
(see [10]). WEY2 implements the shued nested Weyl sequence proposed in
dened by
CLCG4, CMRG96, and CMRG99 are the combined LCG of [25], the combined MRG
given in Figure 1 of [18], and the combined MRG given in Figure 1 of [23].

Table

6.2 gives the log-p-values for the collision test in two dimensions, for LCG1
to LCG6, with k  m and m. As expected, suspect values start to appear
at sample size n  4
these LCGs are denitely rejected with n
m.
LCG4 has too many collisions whereas the others have too few. By extrapolation,
LCG7 to LCG9 are expected to start failing with n around 2 26 , which is just a bit
more than what the memory size of our current computer allowed when we wrote this
paper. However, we applied the two-level collision test with
. Here, the total number of collisions C T is approximately Poisson with
mean . The log-p-values are in Table 6.3. With a total
sample size of 32  2 24 , LCG7 and LCG8 fail decisively; they have too few collisions.
We also tried 4, and the collision test with overlapping, and the results were
similar.
We tested the other RNGs (the last 5 in the table) for several values of t ranging
from 2 to 25. RLUX passed all the tests for t  24 but failed spectacularly in 25
dimensions. With the log-p-value for
the collision test is are 239 collisions, while E[CjH 0 ]  166). For a
two-level test with the total number of collisions
was C much more than This result is not
surprising, because for this generator all the points V i in 25 dimensions or more lie
in a family of equidistant hyperplanes that are 1=
3 apart (see [20, 42]). Note that
RANLUX with a larger value of L passes these tests, at least for t  25. WEY1
passed the tests in 2 dimensions, but failed spectacularly for all t  3: The points are
concentrated in a small number of boxes. For example, with
a sample size as small as
('  14). WEY2, CLCG4, CMRG96, and CMRG99 passed all the tests that we tried.
7. Conclusion. We compared several variants of serial tests to detect regularities
in RNGs. We found that the sparse tests perform better than the usual (dense)
ones in this context. The choice of the function f n;k does not seem to matter much.
In particular, collisions count, Pearson, loglikelihood ratio, and other statistics from
the power divergence family perform approximately the same in the sparse case. The
overlapping tests require about the same sample size n as the non-overlapping ones
to reject a generator. They are more e-cient in terms of the quantity of random
numbers that need to be generated.
It is not the purpose of this paper to recommend specic RNGs. For that, we
refer the reader to [22, 23, 27, 33], for example. However, our test results certainly
eliminate many contenders. All LCGs and LFSRs fail these simple serial tests as
soon as the sample size exceeds a few times the square root of their period length,
regardless of the choice of their parameters. Thus, when their period length is less
than 2 50 or so, which is the case for the LCGs still encountered in many popular
software products, they are easy to crack with these tests. These small generators
should no longer be used. Among the generators listed in Table 6.1, only the last four
pass the tests described in this paper, with the sample sizes that we have tried. All
others should certainly be discarded.



--R


Oxford Science Publica- tions
A Guide to Simulation

Inversive congruential pseudorandom numbers: A tutorial


The serial test for sampling numbers and other tests for randomness
A Guide to Chi-Squared Testing
Pseudorandom number generator for massively parallel molecular-dynamics simulations
Asymptotic normality and e-ciency for certain goodness-of- t tests
RANLUX: A Fortran implementation of the high-quality pseudorandom number generator of Luscher
The Art of Computer Programming
Simulation Modeling and Analysis










A random number generator based on the combination of four LCGs
Selection criteria and testing
An object-oriented random-number package with many long streams and substreams
Structural properties for two classes of combined random number generators
Inversive and linear congruential pseudorandom number generators in empirical tests
A current view of random number generators
A new class of random number generators


Asymptotic divergence of estimates of discrete distri- butions

Good ones are hard to
The Standard C Library
Portable random number generators

Thoughts on pseudorandom number generators
Tests for the uniform distribution
On the add-with-carry and subtract-with-borrow random number generators

--TR

--CTR
Makoto Matsumoto , Takuji Nishimura, Sum-discrepancy test on pseudorandom number generators, Mathematics and Computers in Simulation, v.62 n.3-6, p.431-442, 3 March
Peter Hellekalek , Stefan Wegenkittl, Empirical evidence concerning AES, ACM Transactions on Modeling and Computer Simulation (TOMACS), v.13 n.4, p.322-333, October
Pierre L'Ecuyer , Jacinthe Granger-Pich, Combined generators with components from different families, Mathematics and Computers in Simulation, v.62 n.3-6, p.395-404, 3 March
Pierre L'Ecuyer, Software for uniform random number generation: distinguishing the good and the bad, Proceedings of the 33nd conference on Winter simulation, December 09-12, 2001, Arlington, Virginia
random number generation;collision test;goodness-of-fit;m-tuple test;multinomial distribution;OPSO;serial test