\chapter{Inroduction}

Random number generators are widely used for simulations in 
computational science and engineering. Randomness is often present in 
the formulation of the problem, for example random noise or perturbations,
and quantum processes.  In addition, 
many algorithms are probabilistic, for example Monte Carlo simulation 
and stochastic optimization techniques 
such as simulated annealing or genetic algorithms. 
Random number generators are also used in many other applications, from 
slot machines to cryptography. 
 
Since ``random'' numbers are in practice computed using deterministic
algorithms, these are more accurately called pseudo-random number generators. 
In some simulations, the quality of the pseudo-random numbers (how
closely they resemble truly random sequences) is not that important. 
However in many
of the problems for which random numbers are most heavily used, such 
as Monte Carlo simulation, the quality of the random number generator 
is crucial.
This is especially true in large-scale simulations on parallel 
supercomputers, which consume huge quantities of random numbers, and require 
parallel algorithms for random number generation.
 
As noted in a number of previous review articles 
\cite{KNUTH,MARSAGLIA,PARK,ANDERSON,JAMES,LECUREV,BRENTREV},
random number generators provided by computer vendors or recommended
in some papers and computer science texts have often been of poor quality. 
Even generators that perform well in standard statistical tests for
randomness have sometimes proven to be unreliable for certain applications, 
particularly in Monte Carlo simulations
\cite{KALLE,PARISI,HOOGLAND,FILK,MILCHEV,FERRENBERG,GRASSBERGER,CODDINGTON,FINNS
,FINNS95}.
The many problems caused in the past by inadequate random number generators 
on sequential and vector computers 
are likely to be repeated in a new generation of simulations using parallel 
computers, unless parallel random number generators are very carefully 
studied and tested, and the best algorithms made readily available to users.

...........
%=== The aim of this .....

\newpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Requirement for the Sequencial Random Number Generators}
\label{SEQREQ}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Ideally a pseudo-random number generator would produce a stream of numbers that

\begin{enumerate}
\item
are uniformly distributed,
\item
are uncorrelated,
\item
never repeats itself,
\item
satisfy any statical test for randomness,
\item
are reproduceable (for debugging purpose),
\item
are portable (the same on any computer),
\item
can be changed by adjusting an initial ``seed'' value,
\item
can easily be split into many independent subsequence,
\item
can be generated rapidly using limited computer memory.
\end{enumerate}

In practice it is impossible to satisfy all these requirements exactly.
Since a computer uses finite precision arithmetic to store the state of
the generator, after a certain period the state must match that of a previous
iteration, after which the generator will repeat itself. Also, since the 
numbers must be reproduceable, they are not truly random, but generated by a 
deterministic iterative process, and therefore cannot be completely 
uncorrelated. 
 
For practical purposes, we require that the period of
repetition of the sequence be much larger than the number of 
pseudo-random numbers that might be used in any application, and that
the correlations be small enough that they do not noticeably affect the 
outcome of a computation. The first requirement can be determined fairly
easily, knowing the power of the computer to be used. 
The second is extremely difficult to ascertain, and will generally 
be application-dependent.
 
It has often been the case that the speed of random number generators has 
been overly emphasized, to the detriment of the quality of the generators.
A fast generator requires a minimal number of very simple operations, and
it is this simplicity that often leads to problems with the quality of such
generators.
For most applications the speed of a random number generator is not even
a concern, since the amount of time spent generating the random numbers is 
insignificant compared to the rest of the calculation. 
Users who yearn for super-fast random number generators usually have 
applications that spend most of their time generating random numbers,
and require a huge number of them. These types of applications are 
often the ones that are most sensitive to the quality of the generator, 
in which case it would seem prudent to sacrifice a little speed 
for much better randomness properties.
In using a random number generator, 
it's usually better to be slow than sorry.
 
The speed of a random number generator can usually be increased by allowing
the routine to return an array of values, rather than a single value. This
is clearly advantageous for a vector or parallel implementation, however it
is also usually true for a sequential implementation, since it amortizes the
cost of the function call to the random number generator.
 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Tests for Sequential Random Number Generators}
\label{SEQTEST}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Over the years many widely-used methods for generating pseudo-random numbers  
have been shown to be inadequate, 
either by theoretical arguments, or empirical tests, or both.
In some cases theoretical arguments can show that there are correlations in
the sequence of numbers, however in many cases the problems only show up in 
empirical tests that statistically compare the results produced by the 
random number generators with results expected from a truly random sequence
of numbers. Many standard tests of this kind are available 
\cite{KNUTH,MARSAGLIA,DUDEWICZ}. % Others??
 
In addition to standard statistical tests, it is useful to apply
application-specific tests that are more relevant to some of the various
applications for which random numbers are used. 
As with the statistical tests, these tests generally 
compare the results obtained using a pseudo-random number 
generator with known exact results that would occur if the numbers 
were truly random. 
Tests of this kind include Monte Carlo simulation \cite{BINDER,GOULD} 
of exactly solvable systems such as the two dimensional 
Ising model \cite{HOOGLAND,KIRKPATRICK,FERRENBERG,CODDINGTON,FINNS,FINNS95}, 
simulations of percolation models~\cite{ZIFF},
and random walks~\cite{GRASSBERGER,ZIFF,FINNS95}. 
Generators that pass standard statistical tests have sometimes been found
to fail these application-specific tests. It is therefore important to
use as wide a variety of empirical tests as possible.
Any application can in principle be used to test random number generators,
by comparing results obtained with two different 
generators \cite{KALLE,PARISI,HOOGLAND}.
 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
\section{Random Number Generators}
\label{RNG}
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
In this section we introduce the most common algorithms used in 
random number generators -- the linear congruential, lagged Fibonacci,
and shift register generators, and techniques which combine two or more
of these generators. 

%More detailed information can be found in a number of other review
%articles \cite{KNUTH,MARSAGLIA,PARK,ANDERSON,JAMES,LECUREV,BRENTREV}. 
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Linear Congruential Generators}
\label{LCG}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
Probably the most commonly-used random number generators are  
linear congruential generators (LCGs) \cite{LEHMER,KNUTH,PARK,JAMES}. 
The standard C and Unix generators {\tt RAND} (32-bit precision),
{\tt DRAND48} and {\tt RANF} (48-bit precision) are of this type. 
LCGs produce a sequence $X_{i}$ of random integers using the relation
%
\begin{equation}
X_{i} = (a * X_{i-1} + c) \quad \bmod M,
\label{eqn:LCG}
\end{equation}
% which we will denote by L($a,c,M$).
where $a$ is known as the multiplier, $M$ is the modulus, and $c$ is
an additive constant that may be set to zero.  The parameters
($a,c,M$) must be chosen carefully to ensure a large period and good
uniformity and randomness properties.  The maximum period is $M$ for
suitably chosen parameters ($M-1$ if $c=0$)
\cite{KNUTH}. 
Standard random number generators return a floating point number $x_{i}$ 
in the interval [0,1), which can be obtained by dividing $X_{i}$
by $M$. 
 
LCGs work very well for most applications but are known to have some
major defects.
The main problem is that the least significant bits of the numbers produced 
are correlated, and 
% as a result %??
a ``scatter-plot'' of ordered tuples {($x_{i}, x_{i+1}$, \ldots)} of
random floating point numbers plotted in the unit hypercube shows regular 
lattice structure \cite{MARSAGLIA,PLANES,LECUYER,NIEDERREITER,BRENTREV}.
This problem becomes worse in higher dimensions, which may affect some 
high-dimensional simulations. 
The problem of lattice structure can be quantified using the Spectral 
Test \cite{KNUTH}.
Examples of LCGs with good parameters that perform well in the Spectral Test 
are given in Refs.~\cite{KNUTH,FISHMAN,CODDINGTON,LECUYER64}.
 
Another problem is that many commonly-used LCGs 
(including {\tt DRAND48} and {\tt RANF}) 
use a modulus $M$ that is a power of 2, 
since it is fast and convenient to implement this on a computer.
However this approach produces highly correlated low-order bits
\cite{KNUTH,ANDERSON,BRENTREV}, 
as well as long-range correlations for intervals that are a power of 2 
\cite{KALLE,FILK,PERCUS,GROTHE,DEMATTEIS}.
This can cause problems for certain types of simulations, for example 
when using a hypercubic grid with a size that is a power of 2, 
or for applications that expect the lower-order bits to be random.
To avoid these problems, it is best to use a modulus that is prime rather
than a power of 2, however for many applications requiring   
single-precision floating point pseudo-random numbers (for which the low-order
bits are irrelevant), {\tt DRAND48} and {\tt RANF} are quite adequate, and 
are often useful for checking the results obtained using other generators.
 
The effects of these regularities can be decreased by increasing
the precision of the generators \cite{MARSAGLIA}, 
for example by using 64-bit rather than 32-bit numbers.
However there are practical limits to this approach -- if $M$ is greater
than machine precision, then much slower multi-precision arithmetic must 
be used, so in practice the precision cannot be made arbitrarily large.  
This means that 48-bit and 64-bit LCGs can be quite slow on 32-bit computers,
however many high-performance computers now use 64-bit processors, which 
should become standard in the near future. In that case, the speed of these
generators is not a problem for most applications.
 
Note that for 32-bit integers the period of these generators 
is at most $2^{32}$, or of order $10^9$. 
On current processors capable of $10^8$  
floating point operations per second, this period can be exhausted 
in seconds, so higher precision (48-bit or more) generators should 
be used. For long simulations on high-performance computers, 
even 48-bit generators may have too small a period.
 
Despite their known problems, large precision LCGs with well-chosen
parameters appear to work very well for all known applications, at least
on sequential computers. However we will see later that there are problems
with implementations of LCGs on parallel computers if the modulus is a 
power of 2.
 
Multiple recursive generators (MRG) \cite{KNUTH,JAMES,LECUREV,LECUYER64}, 
generalize LCGs by using a recurrence of the form
%
\begin{equation}
X_{i} = (a_{1} X_{i-1} + a_{2} X_{i-2} + \ldots + a_{k} X_{i-k} + b) 
\quad \bmod M
\label{eqn:MRG}
\end{equation}
%
for a given value of $k$. Choosing $k>1$ will increase the time taken to
generate each number, but will greatly improve the period and randomness
properties of the generator \cite{LECUREV,LECUYER64}. 
Some practical implementations have been provided for $k=2$
\cite{LECUYER64}. 
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Lagged Fibonnaci Generators}
\label{LFG}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
Lagged Fibonacci generators (LFGs) \cite{KNUTH,JAMES,MARSAGLIA} 
are becoming increasingly popular, 
since they offer a simple method for obtaining very long periods,
and they can be very fast.
The standard C and Unix generator {\tt RANDOM} is of this type. 
The sequence is defined by
%
\begin{equation}
X_{i} = X_{i-p} \odot X_{i-q} \\
\label{eqn:LFG}
\end{equation}
%
which we will denote by F($p,q,\odot$), where $p$ and $q$ are the lags,
$p > q$, and $\odot$ is any binary arithmetic operation, such as 
addition or subtraction modulo $M$, multiplication modulo $M$, or
the bitwise exclusive OR function (XOR).
The arithmetic operations are done modulo any large integer value $M$,
or modulo 1 if the $X$'s are represented as floating point numbers in the 
interval [0,1), as can be done if the operation is addition or subtraction.
Multiplication must be done on the set of odd integers.
The modulus function is not required when using XOR.
This method requires storing the $p$ previous values in the sequence
in an array called a {\it lag table}. 
 
As with the LCGs, it is important that the parameters ($p,q,M$) be carefully
chosen in order to provide good randomness properties and the largest
possible period.
If $M$ is taken to be $2^b$ (i.e. the $X$'s have $b$-bit precision),
the maximal period is obtained if 
% and only if %?? - but what about Marsaglia's paper about arrays??
the lags are the exponents of a primitive polynomial \cite{KNUTH,BRENT}. 
In that case the period is
$2^{p}-1$ for XOR, 
$(2^{p}-1)2^{b-1}$ for addition and subtraction,  
and $(2^{p}-1)2^{b-3}$ for multiplication 
\cite{BRENT,BRENTREV,KNUTH,MARSAGLIA,TSAY}. 
Tables of suitable lags are available in the literature
\cite{KNUTH,BRENT,HERINGA,ZIERLER,TSAY}. 
% Should we include a table of values here???
An advantage of this generator is that
the period can be made arbitrarily large by just increasing the lag $p$.
This also improves the randomness properties \cite{MARSAGLIA,CODDINGTON}, 
since smaller lags mean higher correlations between the numbers in the 
sequence. The low-order bits can have particularly poor randomness properties
if small lags are used.
 
Empirical tests have shown that the randomness properties of LFGs are
best when multiplication is used, with addition (or subtraction) being
next best, and XOR being by far the worst \cite{MARSAGLIA,CODDINGTON,FINNS,FINNS
95}. 
This is intuitively reasonable in that multiplication mixes the bits in two 
numbers much more than addition, which is in turn much better than XOR.  
 
LFGs using addition (or subtraction) are the most popular because they are 
very simple and very fast. 
All the computation can be done in floating point, which 
avoids a costly integer to floating point conversion, and large periods can 
be obtained on 32-bit processors without having to use slow multi-precision 
arithmetic. Each pseudo-random 
number can be generated with just a single floating point addition and a
modulus operation. 
 
Great care needs be taken in choosing the lags for this type of generator.
Many implementations use (or recommend in the documentation) lags that are
much too small to give adequate randomness properties, even though
it has been known for many years that additive LFGs fail some standard
statistical tests for very small lags (such as $p = 17$) \cite{MARSAGLIA},
and that increasing the lag improves the randomness properties of the 
generator.
More recently, it was shown that even lags on the order of hundreds can 
give incorrect results in a number of tests based on common applications 
such as Monte Carlo simulation, 
percolation and random walks \cite{GRASSBERGER,ZIFF,CODDINGTON,FINNS,FINNS95}.
Unfortunately this information seems not to have percolated through the
mathematical and computational science community, and it is still extremely 
rare to see a lag of greater than 1000 recommended for an additive LFG.
We recommend using ($p,q$) of at least (1279,1063), 
and preferably much larger values.
% see the table and refs ???  give a value that is around 10,000 ?
The standard Unix generator {\tt RANDOM} is an additive LFG with a default 
lag of 31, which is much too small, however it is possible to initialize it 
with a larger lag.
 
Setting a minimum acceptable lag is obviously a moving target, since computers
are continually becoming faster, allowing for more stringent randomness tests 
that use more random numbers. 
The recommendations given here are based on the results of the current set of 
application-specific tests mentioned above, using currently available 
computers. These values may not be adequate for future applications that use
Teraflop computers. In any simulation, the largest feasible lag should be 
used, and the results should always be checked using a different generator.
 
Small lags were necessary in the past because of limited computer memory,
however on current high-performance computers (and even personal computers), 
the additional memory 
requirement of a lag table with a few thousand entries can easily be handled.
However we will see later that memory constraints may be a problem
in implementing this algorithm on parallel computers for some applications. 
The choice of the lag may affect the speed of the generator, depending
on the type of computer used. For example, if vector processors are used, 
a larger lag may improve performance, since the vector lengths are longer. 
If a scalar processor with limited cache memory is used, having a 
very large lag may cause cache misses and reduce the performance. 
 
Multiplicative LFGs have seen little use, which is somewhat surprising 
considering their excellent randomness properties and extremely long period. 
Although slower than additive LFGs, 
they are just as fast as 32-bit LCGs, and much faster than LCGs that require 
multiple-precision arithmetic. 
Multiplicative LFGs can also be used with much smaller lags than for
additive LFGs. Many different tests are failed by additive LCGs with small
lags (less than 100) 
\cite{MARSAGLIA,FERRENBERG,GRASSBERGER,CODDINGTON,FINNS,FINNS95}, 
however no currently published results in any test show failure of a
multiplicative LFG for a lag as small as 17. However we recommend using a 
lag of at least 100, to ensure a large period and better randomness
properties.
 
One of the obstacles in implementing multiplicative LFGs is handling the
possible overflow of the multiplication. In most languages a portable 
implementation would require multiple precision arithmetic, which would be 
quite slow. However if the C programming language is used, and the $X$'s 
are specified to be of type {\tt unsigned int}, then the language 
specification for the multiplication of two {\tt unsigned int} 
values in C guarantees that the result will be correct modulo $M=2^b$
if $b$ is the word size (the number of bits in $X$), without having
to worry about overflow or the use of multiple precision arithmetic. 
Many Fortran compilers allow the calling of subroutines written in C.
 
The randomness properties of LFGs can be improved (without sacrificing 
too much in speed) by using multiple lags (or ``taps'') 
\cite{GOLOMB,BRENT,ZIFF,GRASSBERGER,CODDINGTON,FINNS95}, 
i.e. by combining three or more previous elements of the sequence, 
rather than two. 
This type of generator has not yet been widely used or studied, 
however it seems likely that a 3- or 4-lag additive LFG would 
be a very fast and effective random number generator.
 
% Add something here about how this can be done?? How lags are chosen??
% Give tables of suitable lags??
 
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Shift Register Generators}
\label{SHIFTREG}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
Shift register (or Tausworthe) 
generators \cite{KNUTH,MARSAGLIA,TAUSWORTHE,KIRKPATRICK,GOLOMB}
are generally used in a form where they can be considered as a special 
case of a lagged Fibonacci generator using XOR.
XOR gives by far the worst randomness
properties of any operation for an LFG \cite{KNUTH,MARSAGLIA,CODDINGTON},
so these generators are not recommended.  
 
Despite their serious drawbacks, shift register generators have been very 
popular in the past, mainly because they were comparatively fast.
However on modern processors addition is not markedly slower than XOR,
so there is little reason to use these generators.
 
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Combined Generators}
\label{COMBINED}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
Combining two different generators has been shown (both theoretically and
empirically) to produce an improved quality generator in many 
circumstances \cite{MARSAGLIA,LECUYER,LECUREV,DENG}.
 
Based on an algorithm introduced by Wichmann and Hill \cite{WICHMANN},   
L'Ecuyer \cite{LECUYER} has shown how to additively combine two different 
32-bit LCGs to produce a generator that passes all known statistical tests 
and has a long period of around $10^{18}$, 
thus overcoming the major drawbacks of standard 32-bit LCGs. 
This has been implemented in a program known as {\tt RANECU} 
\cite{LECUYER,JAMES}. 
Combining two LCGs in this way is effectively a more efficient way of
implementing an LCG with a much larger modulus \cite{TEZUKA}. 
Recently L'Ecuyer {\it et al.} \cite{LECUYER64} have implemented
combined 48-bit and 64-bit LCGs and MRGs, with even larger periods
and better randomness properties. 
A combined 32-bit LCG is substantially slower than a standard 32-bit LCG, 
although it is more appropriate to compare it to a 64-bit LCG (which has
the same period), in which case the performance is similar, at least on
a 32-bit machine where multiple-precision arithmetic is required for the
64-bit LCG. 
 
Other proposed combined generators include algorithms combining 
an LFG with an LCG \cite{MARSAGLIA},
or an LFG with a simple Weyl (or arithmetic sequence) generator, which 
is the basis for the {\tt RANMAR} generator \cite{JAMES,RANMAR}
commonly used in computational physics applications.
The addition of the Weyl generator greatly improves the
randomness properties over the single additive LFG, but {\tt RANMAR}
still fails some Monte Carlo tests \cite{CODDINGTON}, since the lag used  
($p=97$) is much too small. If you are using this generator, you should
greatly increase the lags, to at least (1279,1063).
 
 
 
 