A Report of a Significant Error On a Frequently Used
Pseudo Random Number Generator

arXiv:1408.1900v1 [cs.MS] 7 Aug 2014

Ayşe Ferhan Yeşil, M. Cemal Yalabik
Department of Physics, Bilkent University, 06800 Ankara, Turkey

Abstract
Emergence of stochastic simulations as an extensively used computational
tool for scientific purposes intensified the need for more accurate ways of generating sufficiently long sequences of uncorrelated random numbers. Even
though several different methods have been proposed for this end, deterministic algorithms known as pseudo-random number generators (PRNGs)
emerged to be the most widely used tool as a replicable, portable and easy
to use method to generate such random number sequences. Here, we introduce a simple Poisson process whose simulation gives systematic errors when
the very commonly used random number generator of the GNU C Library
(Glibc) is utilised. The PRNG of Glibc is an additive lagged Fibonacci generator, the family of such PRNGs are accepted as relatively safe among other
PRNGs. The systematic errors indicate complex correlation relations among
random numbers which requires a further explanation.
Keywords: pseudo-random numbers, pseudo-random number generator,
PACS: 02.70.-c, 05.10.-a
1. Introduction
Random numbers are at the core of stochastic simulations for scientific
purposes. Nearly in all of these cases, but especially in Monte Carlo simulations, long sequences of uncorrelated, uniformly distributed random numbers
are required. Several different methods are used to generate such random
number sequences, but as a replicable, portable and easy to use method, deterministic algorithms, namely pseudo-random number generators (PRNGs),
are preferred extensively. However, to be confident of such a number sequence
being “sufficiently” random, i.e. any combination of random numbers within
Preprint submitted to Journal Of Computational Physics

August 15, 2014

the sequence is uncorrelated, is not straightforward since there is no complete
and universal statistical test of analysis. In other words, even though one
cannot be confident that a sequence is totally random by any analysis known
so far, it can still be decided that the sequence is prone to correlations in the
case of empirical or theoretical tests indicating in that direction [1]. As a result, correlation tests for PRNGs have become necessary tools for physicists
[3, 5, 6, 4, 2, 12, 7, 9, 8, 10, 11]. Exemplifying the significance of these tests,
it is important to note that most of the tests are direct results of PRNGs that
go wrong. One of the first such paradigmatic examples, dated back to 1993,
is encountered on the simulation of 2-D Ising Model which shows that results
of various algorithms are different from the exact result and they also differ
from algorithm to algorithm depending on the PRNGs used. This analysis
showed that Generalized Feedback Shift Register (GFSR) PRNGs (a type
of lagged Fibonacci generators (LFGs) which utilise exclusive-or operations)
with smaller lags produce dramatically wrong results, up to hundred times of
the standard deviation, when used with Wolff algorithm [3]. The depth first
nature of the algorithms are suspected to interfere with the intrinsic triple
correlations of the GFSRs. However, Schmid et al showed that even though
the algorithm is randomly updating the system one spin at a time, if the
model has triple structure in itself such as Blume-Capel model, use of random numbers produced by GFSR in the simulations lead to dramatic errors
[4]. On the other hand, GFSRs are not the only faulty PRNGs. LFGs with
addition and subtraction operations also cause systematic errors if they have
small lags [5]. The default random number generator of GlibC, the library of
the commonly used C-compiler unix based environments, is an LFG with addition operation LFG(31,3,+). It’s lag is 31 which is generally considered as
small. Grassberger reported of correlations in Lagged Fibonacci Generators
with correlation lengths proportional to their lags. His method resembles
Vattulainen’s n-block method[10]. Grassberger proposed that it might be related to the triple correlations intrinsic to the random number sequence produced by any LFG. He also showed while simulating 3-D self-avoiding walks
on a cubic lattice, that even an LFG with lag 250, LFG(250,103,+), gives
poor results. This led to the suspicion that LFG(55,24,+) and LFG(17,5,+)
may operate poorly, and was demonstrated to be so [6]. In light of these,
the default PRNG of GlibC should be considered as unreliable since it has a
small lag as well.
In this paper, a simple test model, which shows strong correlations when
utilised with the default random number generator provided with the C com2

piler of Linux distributions, is introduced. C is a very commonly used programming language and is also very popular among scientific programming
communities. The GNU C Library (Glibc) is used as the C library in the
GNU systems and most systems with the Linux kernel. The PRNG of this
library is constructed using a linear additive feedback method. We demonstrate the existence of relatively strong correlations through the use of a
simple algorithm which is related to a one dimensional diffusive gas problem.
To inform the unspecialised reader we provide here further information
about the pseudo-random number generators: There are basically two common types of PRNGs: Linear Congruential Generators (LCGs) and Lagged
Fibonacci Generators (LFGs) [5]. LCGs are associated with the general
equation
Xi = (AXi−1 + B) mod M
(1)
which can be represented in short as LCG(A, B, M), with a maximum period
of M that can be reached through a suitable choice of A and B. Here A, B
and M are integers, and Xi is the ith generated random number.
LFGs have the general equation:
Xi = Xi−P ⊙ Xi−Q

(2)

where P and Q are the lags, conventionally chosen as P > Q. LFGs can be
represented as LF G(P, Q, ⊙) [5]. The symbol ⊙ represents any binary arithmetic operation such as addition, subtraction, multiplication or exclusive-or
(xor) operation. If this operation is addition or subtraction, the maximum
period that can be reached for b-bit precision is (2P −1)2b−1 . If the arithmetic
operation is multiplication, the maximum period is (2P − 1)2b−3 [5]. Moreover, if the operator is xor the maximum period becomes 2P −1P . The random
number generator of Glibc uses an additive lagged Fibonacci generator, which
is a combination of a LCG and LFG utilizing 32-bit integer operations. The
lags are produced by a linear congruential generator LCG(16807, 0, 231 − 1)
and then by using those lags a LF G(31, 3, +) produces the random number
sequence. (The first 344 members of the sequence are omitted.) Newest versions of Glibc allow initialization of lags to be done by any seed (X0 ), but
here we stick to the default case in which the initial seed is equal to 1. In
order to apply Eqn 2., an initial set of Xi for i < P must be defined. This is
accomplished by
X0 = 1
3

Xi = 16807Xi−1 mod (231 − 1) for 1 ≤ i < 31
Xi = Xi−31 for i = 31, ..., 33.

(3)

The generation part is,
Xi = (Xi−3 + Xi−31 ) mod 232 for i ≥ 34.

(4)

All Xi ’s for 0 ≤ i ≤ 343 are discarded and in the end ith output of the
RNG becomes Ri = Xi+344 ÷ 2 [13]. The integer division by 2 is carried out
to convert the number to a 31 bit number, a leading zero bit is necessary to
obtain a positive number. In order to obtain real pseudo random numbers,
the result is divided by the biggest 31 bit integer:
ri =

Ri
.
−1

231

(5)

This result then provides a pseudo-random real number linearly distributed
between 0 and 1 with 31 bit precession.

2. Results
The model we study is simple: There are N boxes numbered 1 to N in
a chain structure, all boxes are equally likely to receive a ball at a random
time driven by a Poisson process (Fig. 1). Each time when a ball is received,
a random time increment of δt has elapsed, where δt = − log(r)/N and r is
a uniformly distributed random number between 0 and 1. This then defines
a Poisson process in which an event occurs with rate N per unit time. With
this choice of δt, on the average N balls are placed in unit time. The box
that will receive the ball is chosen randomly as well, such that ith box is
chosen if i − 1 < rN < i (Fig. 2). In between these two steps (choosing the
random time and choosing the random box), n samples are chosen from the
random number sequence in vain, i.e., they are discarded from the used sequence. Furthermore, the simulation is stopped at multiples of a time period
T , so that at each time step when the interval crosses a multiple of T . The
time variable is set to that time (that is, no balls are placed to any box) and
then the simulation continues as described in Fig. 1. The pseudocode for the
model can be seen in the Table 1.

4

Figure 1: a) Each box is equally likely to receive a ball at a random time. b) Time
t is incremented by δt at each step. The small dots on the axis indicate n draws in
vain from the PRNG, and the larger dots indicate a ball is placed to a random box
at that point in time. Whenever δt crosses a multiple of T , the process is paused,
i.e., no ball is placed to any box and time is set to that multiple of T .

The process involving the period (T ) dependence is necessary to observe
the correlation. Indeed Fig. 3 shows that the process is sensitive on the period. In our research, we were analysing the effect of a periodic perturbation
on a system, and that is how was the problem discovered. We then tested
the process for time increment derived from a simpler form δt = r/N, and
still found correlations.
It is expected that over a large period of time each site will be equally
likely to be occupied. However, for n = 0, boxes in the middle have a
tendency to be occupied more than the sites near the both ends (Fig. 2).
Similarly, for n = 1 and n = 2, correlations are observed. However if n > 2,
correlations vanish, in other words there are no preferences left among the
boxes. Note that the induced correlations are significant, i.e., they exceed
5

1.01

n=2
1.005

n=0

n=3
1

0.995

n=1

0.99

0.985
0

2

4

6

8

10

12

14

16

18

20

Figure 2: n = 0, 1, 2 and 3, T = 0.25, N = 20. n is the samples drawn in vain from

the PRNG, T is the period and N is the number of boxes. Dependence on n is
related to the correlations among consecutive numbers.

1% in our example. Also note that, simulations are done for 109 trials, which
is well below the period of the PRNG.
Note that, the abscissa of the graphs in Figures 2 and 3 correspond to the
size of the random number used for selecting a box. Different shapes of
graphs in these plots correspond to a non-uniform distribution for the size
of the random number. Dependences on n and T are then related to the
correlations among the random numbers.
Since we wanted a measure of the systematic error in the probability
distribution P , we looked
P at the first coefficient of the Fourier expansion of
this distribution: c = N
j=1 exp(iπj/N)(P (j) − c0 ) where c0 is the average
value of P . Figure 4 shows magnitude of c, which represents the size of the
systematic error as a function of the period.

6

Pseudocode for the process
1. SET N to lattice size
2. SET period to T
3. SET totalT ime to 0
4. SET n to number of draws in vain
5. SET a to 1
6. SET occupation of k th site o(k) to 0 for 1 ≤ k ≤ N
7. DRAW a random number r
8. SET incrementalT ime to −log(r)/N
9. IF totalT ime plus incrementalT ime is smaller than a × T THEN
10.
INCREASE totalT ime → totalT ime + incrementalT ime
11.
DRAW a random number r
12.
DRAW n random numbers in vain
13.
CHOOSE k th box as k = ⌈rN ⌉
14.
INCREMENT o(k) → o(k) + 1
15. ELSE
16.
INCREMENT a → a + 1
17.
SET totalT ime as a × T
18. REPEAT 6 − 17 steps for 1 billion times

Table 1: The pseudocode which is neccesary to produce the results in this paper.

3. Conclusions
Given the common usage of C as a scientific programming language, we
report of correlations in the default random number generator of Glibc. Even
though it can pass some tests, PRNG of Glibc shows significant correlations
under this very simple test. Systematic errors of this magnitude will cause
problems in simulations which demand random numbers with low correlations, such as Monte Carlo analysis. Scientific community should beware of
accepting default random number generators even for well established, commonly used programming languages as C. We could not identify the precise
mechanism which leads to this behavior. A mathematical analysis of this
mechanism represents an interesting challenge for further analysis.
4. Acknowledgements
This work was supported by Turkish Academy of Sciences.

7

1.015

T=0.3

1.01

T=0.2

1.005
T=0.4
1

T=0.5

0.995

0.99
T=0.1
0.985
0

5

10

15

20

25

Figure 3: n = 0, T = 0.1, 0.2, 0.3, 0.4 and 0.5. For different values of periods we see

different shapes of graphs, however for higher periods correlations cease to exist.

5. References
References
[1] I. Vattulainen, arXiv:cond-mat/9411062
[2] R. M. Ziff, Comput. Phys. 12, 385 (1998)
[3] A. M. Ferrenberg, D. P. Landau and Y. J. Wong, Phys. Rev. Lett. 69,
3382 (1992).
[4] F. Schmid and N. Wilding, Int. J. Mod. Phys. C 6, 781 (1995)
[5] P. D. Coddington, Northeast Parallel Architecture Center, Paper 14
(1994).
[6] P. Grassberger, Physics Letters A, 181 (1993)
[7] S. Mertens and H. Bauke, Phys. Rev. E 69, 055702(R) (2004)

8

0.07

0.06

0.05

c

0.04

0.03

0.02

0.01

0
0

0.5

1

1.5

2

2.5

τ

Figure 4: For different values of period the magnitude of first Fourier coefficient c.

It represents the size of the systematic error as a function of the period.

[8] D. E. Knuth, The Art of Computer Programming, 3rd ed. (AddisonWesley, New York, 1998), Vol. 2
[9] H. G. Katzgraber, arXiv:1005.4117
[10] I. Vattulainen, T. Ala-Nissila and K. Kankaala, Phys. Rev. E 52, 3205
(1995)
[11] I. Vattulainen, K. Kankaala, J. Saarinena, T. Ala-Nissilaa, Comp. Phys.
Comm., 86(1995) 209-226
[12] I. Vattulainen, T. Ala-Nissila and K. Kankaala, Phys. Rev. Lett. 73,
2513 (1994)
[13] This webpage also contains code to generate the random number sequence discussed in this paper. Peter Selinger: The GLIBC Pseudorandom Number Generator. N.p., n.d. Web. 24 Feb. 2014

9

