4.

Generating random numbers

Central to any MC simulation are the random numbers. Hence it is important to have a good
source of random numbers available for the simulations. Getting ’good’ random numbers is in fact
not quite as easy as many people think it is, so we will spend quite some time on this topic.

Random numbers can be divided into (at least) three categories: true random numbers,
pseudorandom numbers and quasirandom numbers.
• True random : no way to predict what the next random number is
• Pseudorandom: a sequence of numbers which is algorithmically produced and not really
random
– It can be fully repeated if the initial conditions and algorithm are known
• Quasirandom numbers: act as random numbers in some sort of simulations, but are
well-ordered in some other types.
More precise explanations for these categories are given later in this section.

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 1

 4.1. Non-automated means

Although these are almost never used anymore, it is good to keep in mind that there do exist
non-automated means to generate random numbers.
The most obvious one is to just write one manually.

Although that has obvious problems for generating large sets of numbers, it is widely used for
selecting the seed number (see below) for electronic random number generators.

Another one, which was historically used to some extent, is to select numbers from some number
sequence, e.g. the phone book or the decimals of π .
• The former method is highly nonadvisable, as there obviously can be strong non-random
features in phone numbers.
• The latter is not so bad, since the decimals of π are not supposed to have correlations

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 2

One could for instance get rough random numbers between 0 and 1 by selecting always 4 numbers
at a time from π and divide these by 10000:
3.141592653589793238462643383279502884197169399375105
These non-automated random numbers are essentially pseudorandom.

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 3

 4.2. Mechanically generated random numbers

[Knuth, Seminumerical algorithms]

As technology evolved, mechanical devices were devised which produced random numbers.
• The machines used in lottery or roulette
• Also specifically designed machines which could generate thousands of numbers have
been built. In 1955 the RAND corporation actually published a book with 1 million
random numbers. (undoubtedly a strong contender for the title of most boring book ever
published).
Since the advent of electronic computers, such devices have become obsolete, except in applications
(such as the lotto games) where it is important that even a layman can easily check that there is
no scam involved.

Mechanical random numbers can be true random numbers, such as the ones generated in a lotto
machine or a roulette, (unless maybe if the casino is run by the mob...). Numbers generated by a
mechanical calculating device may also be pseudorandom.

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 4

 4.3. Electronic hardware-generated random numbers

[http://www.fourmilab.ch/hotbits/]

4.3.1. True random numbers

It might seem an obvious idea to design microprocessors, or parts of them, to be able to generate
random numbers electronically - that is, design an electronic part which delivers a signal which
randomly gets translated to 0’s and 1’s when they are translated to digital form. But this is in
practice very rarely done.

But even an electronic random number generator could have its problems; it is easy to imagine that
minute electronic disturbances from the environment could affect the results produced by it.

For a physicist, an obvious idea to obtain random numbers independent of any reasonably possible
outside disturbance is to use radioactive decay. Since nuclear processes are almost completely

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 5

independent of all every-day electronic interactions, this should (if properly calibrated) deliver true
randomness.

A quick google search reveals that there does indeed exist random numbers based on radioactivity,
see http://www.fourmilab.ch/hotbits/. Here is a sample of 16 truly random bytes:
unsigned char hotBits[16] = 168, 193, 241, 195, 37, 251, 190, 121, 105, 137, 173, 226, 95, 181,
239, 231 ;
But unfortunately downloading random numbers from the internet would be way too slow for most
MC simulations, and putting a radioactive source inside a chip is not a very good idea either.

4.3.2. Pseudorandom numbers

Before we discuss the common algorithmic pseudorandom numbers, we mention something which is
something between true and pseudorandom numbers.
4.3.2.1. HAVEGE – HArdware Volatile Entropy Gathering and Expansion
[http://www.irisa.fr/caps/projects/hipsor/HAVEGE.html]

There is a serious effort going on in generating random numbers based on computer hardware. This
is not done in the chip electronics, but by software which reads chip registers.

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 6

 • These are unpredictable in the sense that it is practically impossible afterwards to know
what values they had at the moment of reading.
This is especially important for modern cryptography, which heavily relies on random numbers.
• In case cryptography uses a pseudorandom number generator where the seed could be
somehow figured out, then knowledge of what algorithm is used could be used to help
crack the code. Hence it is important to have absolutely non-repeatable random number
generators.

There are several such approaches, but most are pretty slow, capable of producing only 10’s or 100’s
of bits per second.

In a quite new (2002 on) development, Seznec and Sendrier have formed a generator which uses
several advanced features of modern superscalar processors, such as caches, branch predictors,
pipelines, instruction caches, etc. Every invocation of the operating system modifies thousands of
these binary volatile states. So by reading from them, one can obtain data which will immediately
change afterwards and hence is practically impossible to reproduce.

So for practical purposes they do generate true random numbers, although formally they do not.

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 7

“HAVEGE (HArdware Volatile Entropy Gathering and Expansion) is a user-level software
unpredictable random number generator for general-purpose computers that exploits these
modifications of the internal volatile hardware states as a source of uncertainty.”
“During an initialization phase, the hardware clock cycle counter of the processor is used to gather
part of this entropy: tens of thousands of unpredictable bits can be gathered per operating system
call in average.” In practice, this is done for 16384 entries in a storage array called “Walk”.
After this, HAVEGE works as a pseudorandom number generator which keeps modifying the Walk
array. But additional hardware-generated uncertainty is read in and used to read entries from the
Walk table in a chaotic order. Also, hardware-generated data is used to modify the Walk table
entries. Their pseudorandom-number generator relies on XOR operations (see below). Moreover,
the data is hidden so that even the user of the routine can not read it.

It is even possible to randomly personalize the generator source code during installing, and not only
with numbers but actually even by randomizing bit-level logical and shift operators in the source
code.

On a 1GHz Pentium, the routine can throughput randomness at a rate of 280 Mbit/s. On Pentium
III’s, the uncertainty is gathered from the processor instruction cache, data cache, L2 cache, as well
as the data translation buffer.

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 8

This routine clearly seems very promising for efficient generation of random bits (numbers) on
modern processors.

The downside here is that the implementation can not possibly be fully portable. As of the writing
of this (Jan 2006), the routine works for Sun UltraSparc’s, Pentium II, III, 4, Celeron, Athlon and
Duron, Itanium, and PowerPC G4. A nice list, but by no means complete.
For scientific simulations, as discussed in the next section, it is actually usually desirable to use
predictable random numbers, to enable repeatability. For this kind of application, HAVEGE is
clearly unsuitable.

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 9

 4.4. Algorithmic pseudorandom numbers

[http://www.physics.helsinki.fi/∼vattulai/rngs.html. Numerical Recipes. Knuth. Karimäkis notes]

The by far most common way of obtaining random numbers is by using some algorithm which
produces a seemingly random sequence of numbers. But the numbers are not truly random; if the
call to the routine is repeated with the same input parameters, the numbers will always be the
same. Such numbers are called pseudorandom numbers.

The most common way to implement the generators is to enable initializing them once with an
integer number, called the seed number. For a given seed number, the generator produces always
the same sequence of numbers.
The seed number could be set from the system clock, or just selected manually. Doing the latter is
actually almost always advisable, since this allows one to repeat the simulation identically, i.e. the
simulation is repeatable.
• Science should per definition be repeatable!
• Repeatability is also useful from a practical point of view: can get more details if needed

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 10

The sequence in any pseudorandom generator will (has to) eventually repeat itself. The repeat
interval (period) is an important measure of the quality of the generator.
• Period should be larger than amount of numbers needed!

The repeat interval in the most used decent generators is of the order of the number of integers
that can be represented by a 4-byte integer, 232 ≈ 4 × 109. While this might seem large enough,
it is not necessarily so!
• A simple calculation shows it is quite easy to exceed 4 × 109 random numbers nowadays.

Another important quantity to know for integer generators is the maximum value the routine can
return. This is often (but not always) about the same as the maximum possible value of a (signed
or unsigned) integer or long integer variable, i.e. e.g. 231 − 1 or 232 − 1.
To then get floating point random numbers divided evenly between 0.0 and 1.0 one should simply

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 11

convert the integer to a floating point number and divide it with the maximum m,
rinteger
rfloat = (1)
m
Random numbers divided evenly in an interval are called uniform random numbers.

Finally, it is often important to know exactly what the limits of the generator are. That is, for a
floating point generator, can exactly 0 and 1 (within the numerical precision) be returned or not?
• Is the interval [0, 1], [0, 1[, ]0, 1] or ]0, 1[.
Purely mathematically this does not matter: the probability of hitting exactly the limit is exactly 0.
But on a computer it may matter a lot!. Especially if we work with 4-byte floating point values,
which have a mantissa with only about 7 digit accuracy, the probability of hitting the limit is
actually about 1 in 10 million. Even for 8-byte floating point variables, with a 15-digit mantissa,
the probability of hitting the limit may not be negligible.
• It is easy to think of cases where this matters (imagine e.g. what happens if you take a
logarithm of a random number which can be 0).
• Different applications need different limits, and many times it really does not matter after
all.

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 12

 4.4.1. The basic linear congruential algorithm

I will now attempt to take a pedagogical approach to generating random numbers, presenting one
of the most common approaches step-by-step.
One of the simplest decent, and probably still most used method to generate random numbers is to
use the following equation

Ij+1 = aIj + c (mod m) (2)

Here “mod” us the modulus (remainder of division) operation.

This approach is called the linear congruential algorithm or if c = 0 the multiplicative c. a.

For the sake of example, let us take as the seed number I0 = 4, and select the constants as
a = 7, c = 4 and m = 24 − 1 = 15. Plugging in the numbers gives

aI0 + c = 32

and after taking the modulus we get I1 = 2. Continuing gives the following list of numbers:

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 13

 Step i Ii
0 4
1 2
2 3
3 10
4 14
5 12
6 13
7 5
8 9
9 7
10 8
11 0
12 4
13 2
14 3
15 10
16 14
So we see that I12 = I0, i.e. the period is 12. This is pretty good, considering that obviously the
period cannot exceed m.

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 14

Otherwise the sequence of numbers does indeed look fairly much random. But if we instead use
e.g. I0 = 11 the period becomes only 3:
Step i Ii
0 11
1 6
2 1
3 11
4 6
5 1
6 11
This is of course terrible!

This example illustrates many important things:

1. It is not good to have a generator where the result depends on the seed number.

2. Hence care should be taking in selecting not only the random number generator but also the
constants in it.

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 15

3. Do not fiddle around with the constants, thinking that a “random” change in one of the numbers
would improve on the randomness of the results. Since the constants should be carefully selected,
this is likely to lead to a much worse result (this will be illustrated in the exercises)

For the linear congruential generator there are actually quite simple rules which guarantee that the
period reaches the maximum m − 1, i.e. we obtain all integers in the range [0, m − 1] before
repetition. These are:
• c and m should have no common prime factors
• a − 1 is divisible with all the prime factors of m
• a − 1 is divisible with 4 if m is divisible with 4.
Since the numbers chosen above, a = 7, c = 4 and m = 15, do not fulfill the second criterion, it
is indeed to be expected that the full period is not achieved, and that the period length may depend
on the seed. This is a terrible feature for a generator.

But if we had chosen m = 17, the full period will be achieved for any seed, except 5 (since 7
×5 + 4mod 17 = 5). In general, if a solution in the interval [0, m − 1] exists for the equation
aI + c (mod m) = I there will always be one seed value I which fails in them.

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 16

But these criteria are not enough to obtain an optimal generator. Or in other words they are a
necessary but not a sufficient set of rules. There may still be bad correlations in the numbers, and
much more extensive testing should be carried out to obtain a truly good linear generator. Hence it
is best not to fiddle around with the numbers on your own.

A note on coding generators in practice. There are two basic approaches one can use to code these
generators. The first is to have two separate subroutines, one which sets the seed number, another
which gives the random numbers. The first obviously takes the seed number as an argument,
whereas the latter does not need any argument at all.
This is the approach taken e.g. in the C language. The system-provided random number generator
is given a seed with the function
void srand(unsigned int seed)
and the actual random numbers are generated using the function
int rand(void)
This means that the subroutine srand has to pass the seed number to rand() somehow. But this
can of course be easily implemented e.g. with an external variable, or a module in Fortran90. In C:

unsigned long seed=1;

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 17

int ansistandardrand(void) /* Not recommended for anything */
{
long a=1103515245;
long c=12345;
long div=32768;

seed = seed*a + c;
return (unsigned int) (seed/65536) % div;
}

void ansistandardsrand(unsigned int seed_set)

{
seed=seed_set;
}

Another approach is to require the seed number to always hang along in the function call, and be
treated identically no matter what. In that case one simply sets seed once, then lets it change
value without touching it outside the main routine. E.g.

seed=1;
for(i=0;i<=10000;i++) {

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 18

 ... (code not modifying seed) ...
randomnumber=rand(&seed)
... (code not modifying seed) ...
}

Which approach is better to use will of course depend on the other aspects of the code.
• The first one is slightly simpler in large codes, as one does not have to drag along the
seed number over many places of the code.
• The latter one has the advantage that one could use different random numbers sequences
at the same time easily, by using variables seed1, seed2 and so on

But keeping the possible problems in mind still, the simple equation 2 is often quite good enough for
many applications, especially with a careful choice of the parameters. A very widely used generator
is the “Minimal standard” generator proposed by Park and Miller. It has
5 31
a = 7 = 16807 c = 0(!) m=2 −1

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 19

This is by no means perfect, and indeed in many applications it is horribly bad.

When implementing it one meets a practical problem. Since the modulus factor m is almost 231,
the values returned by the generator can also be close to 231. This means that on a computer one
needs at least 31 bits to describe them. Then if we look at the product aI , one sees that the values
can easily exceed 232.
Why is this a problem? Because many compilers to date can only handle integers with sizes up to
32 bits.
For instance on most 32-bit computers (including Pentiums) the sizes of the C integer variables are
16 bits (data type short) or 32 bits (data type int and long).
So doing the product above directly is not even possible on all machines.
Fortunately there exists a well-established solution by Schrage. It involves an approximate
factorization of m, which we will not go into here. An implementation of this, from Numerical
Recipes in C is:

#define IA 16807
#define IM 2147483647

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 20

#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define MASK 123459876

float ran0(long *idum)

{
long k;
float ans;

*idum ^= MASK;
k=(*idum)/IQ;
*idum=IA*(*idum-k*IQ)-IR*k;
if (*idum < 0) *idum += IM;
ans=AM*(*idum);
*idum ^= MASK;
return ans;
}

(if you want the actual integer seed number idum returned, just remove the multiplication with AM.
)

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 21

This returns a random number in the interval [0, 1[, which is a common convention in the field.
The XOR function (ˆ ) with MASK is just a trick to prevent problems if the routine is called with a
seed of 0.

When using C, it is actually possible to avoid this implementation problem using a dirty trick.
• The C standard actually specifies that if you multiply two 32-bit integers and place the
result into a 32-bit int register, the result returned is the low-order 32-bit of the real, up
to 64-bit result.
• So this is exactly the operation of taking a modulus with 232 for unsigned integers. Hence
once can reduce the entire linear congruential generator into a single line:

unsigned long seed=1;

...
seed=1664525L*seed+1013904223L;
...

The constants are by Knuth and Lewis [ref. in Numerical Recipes].

This is highly non-portable, since it requires that the length of “long” be 32 bits (not true on e.g.

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 22

Alphas). But it should not be much worse than the Park-Miller generator, and it does have the
advantage of being extremely fast (only one multiplication and addition required), so its use might
be justifiable in a temporary application which will never be transferred anywhere.

4.4.2. Problems with the linear congruential algorithm

The problems with the simple approach 2 can be divided into two categories: poor implementations,
and fundamental problems.
4.4.2.1. Implementation problems
Although it might not seem necessary to list ways in which one can screw things up in programming,
there is one important example which should be mentioned.
In the ANSI standard of the C language, the way to implement the system functions (in stdlib.h)
rand() and srand(seed) is not specified. But they do specify an example of a generator, which
is (here I have modified the code to add clarifying variable names).:

unsigned long seed=1;

int ansistandardrand(void) /* Not recommended for anything */
{
long a=1103515245;
long c=12345;

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 23

 long div=32768;

seed = seed*a + c;
return (unsigned int) (seed/65536) % div;
}
void ansistandardsrand(unsigned int seed_set)
{
seed=seed_set;
}

Note that the modulus m is not given explicitly, but really is 232. But then the returned value
is divided by the very low number div=32768, meaning that only 32768 distinct values can be
returned, and that the smallest non-zero value which can be returned is 1/32768. Hence here the
repeat interval is much longer than the number of returned values.
It is obvious that there are many many applications where values of less than 1/32768 can be
returned. So this example should essentially not be used for anything.

Unfortunately this generator, and even worse modifications of this, is in fact implemented on many
compilers (both C and Fortran, don’t know about Java). This leads us to the important rule of
thumb: never use the compiler standard random number generator! And do not think that

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 24

even if you test on one system the compiler generator, that you can use it in many places. Since
good codes should be portable to any system, this is not an acceptable line of thinking.

4.4.2.2. Fundamental problems

Some thought also reveals that there are quite basic fundamental problems associated with the
linear congruential sequence generators.
One is that there is sequential correlation between successive calls. This can be a particular problem
when generating manydimensional data (an example is given in the exercises); the data will appear
as planes in the many-dimensional space.

Related to this is this simple problem: if you consider a very small value of I in the Park-Miller
generator
Ij+1 = 16807Ij (mod 2147483647)
say Ij = 10. Then Ij+1 = 168070, i.e. still much less than the modulus 2147483647. So when I
is divided by the modulus to give a floating-point value, we see that we first get 4.6566128 × 10−9,
then 4.6566128 × 10−8. I.e. a very small value will always be followed by another small value,
whereas for truly random numbers of course it could be followed by any number between 0 and 1 !

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 25

And yet another problem: the sequence can not return the same number twice after each other,
even though in true random numbers even this should be possible.

There are also more subtle problems, such as so called short-range serial correlations. And they also
often can have really terrible problems in a 2D plane; this is illustrated in an exercise.

4.4.3. Developments of the linear congruential algorithm

To overcome these problems, it seems like a pretty obvious idea to ’shuffle’ the numbers somehow.
This should at least solve the problem with successive small numbers and the 2D problem, and
might help with the short-range correlation as well.
This can be simply achieved by having a table which holds a number of random numbers, and
returns one of them in a random sequence. In practice, already a small table is enough to improve
on the results significantly. Numerical Recipes presents a solution which has an array of 32 elements.

#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 26

#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

float ran1(long *idum)

{
int j;
long k;
static long iy=0;
static long iv[NTAB];
float temp;

if (*idum <= 0 || !iy) {

if (-(*idum) < 1) *idum=1;
else *idum = -(*idum);
for (j=NTAB+7;j>=0;j--) {
k=(*idum)/IQ;
*idum=IA*(*idum-k*IQ)-IR*k;
if (*idum < 0) *idum += IM;
if (j < NTAB) iv[j] = *idum;
}

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 27

 iy=iv[0];
}
k=(*idum)/IQ;
*idum=IA*(*idum-k*IQ)-IR*k;
if (*idum < 0) *idum += IM;
j=iy/NDIV;
iy=iv[j];
iv[j] = *idum;
if ((temp=AM*iy) > RNMX) return RNMX;
else return temp;
}

How does this work? The first long if clause initializes the sequence the first time the generator is
called. It first generates 8 numbers which are thrown away, then fills in the array with 32 random
number elements.
After this the ordinary Park-Miller generator follows, except that in the middle we have the operations

j=iy/NDIV;
iy=iv[j]
iv[j] = *idum

Here j first acquires a random value from the iy value, which is the random number generated in

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 28

the previous call. NDIV has a size such that j will be in the range 0-31, as it should. Then the
returned random number is set from the stored number iv[j], after which the random number
generated in this call is written into the array in the position j.
This generator clearly solves many of the problems mentioned above. But even it can not return
the same number twice after each other (in that case the period would be reduced to 32 or less).

4.4.4. Combined linear congruential generators

Proceeding in the development of the congruential generators, one can combine two single generators
to form one with a very much longer period. This can be done by generating two sequences, then
subtracting the result of one of them from the other (subtraction prevents an integer overflow). If
the answer is negative, the number is wrapped to the positive side by adding the modulus m of one
of the generators.
This forms a random-number sequence whose period can be the multiple of the period of the two
generators. With generators similar to the Park-Miller generator with m ∼ 231, one can thus reach
a period of the order of 262 ≈ 4.6×1018. Important in selecting the moduli m1 and m2 is that
the periods they form do not share many common factors. In the following generator the periods
are
m1 − 1 = 2 × 3 × 7 × 631 × 81031 = 2147483562
and
m2 − 1 = 2 × 19 × 31 × 1019 × 1789 = 2147483398

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 29

so they share only the factor of 2, and the period of the combined generator thus becomes
≈ 2.3×1018. Thus at least period exhaustion is practically impossible on present-day computers,
although not necessarily in 10 or 20 years.

#define IM1 2147483563

#define IM2 2147483399
#define AM (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

float ran2(long *idum)

{
int j;

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 30

 long k;
static long idum2=123456789;
static long iy=0;
static long iv[NTAB];
float temp;

if (*idum <= 0) {
if (-(*idum) < 1) *idum=1;
else *idum = -(*idum);
idum2=(*idum);
for (j=NTAB+7;j>=0;j--) {
k=(*idum)/IQ1;
*idum=IA1*(*idum-k*IQ1)-k*IR1;
if (*idum < 0) *idum += IM1;
if (j < NTAB) iv[j] = *idum;
}
iy=iv[0];
}
k=(*idum)/IQ1; ! A
*idum=IA1*(*idum-k*IQ1)-k*IR1; ! A
if (*idum < 0) *idum += IM1; ! A

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 31

 k=idum2/IQ2; ! B
idum2=IA2*(idum2-k*IQ2)-k*IR2; ! B
if (idum2 < 0) idum2 += IM2; ! B
j=iy/NDIV; ! S+G
iy=iv[j]-idum2; ! S+G
iv[j] = *idum; ! S+G
if (iy < 1) iy += IMM1; ! P
if ((temp=AM*iy) > RNMX) return RNMX;
else return temp;
}

What is going on here? The if clause is again initialization. Part A calculates the first random
number using Schrage’s method, part B likewise the second. Part S+G handles the shuffle and
generates the combined random number, using the periodicity clause on line P.
Numerical Recipes places great trust in this algorithm, they even promise $ 1000 to anyone who
demonstrates that it fails in any known application. As far as I know, it has not been reported to fail
so far, despite being subject to quite some testing [www.cs.adelaide.edu.au/users/paulc/papers/sccs-526/sccs-526.ps.gz].
(note that RAN2 in the first edition of numerical recipes is entirely different, and this has been
reported to have problems [http://www.lysator.liu.se/c/num-recipes-in-c.html].)

4.4.5. Generalized feedback shift register algorithm

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 32

[Lewis, Payne, Journal of the ACM 20 (1973) 465]

Another, widely used and independent line of generators are the so called GFSR generators, which
were originally, in 1973, developed to overcome some problems in the simplest linear congruential
algorithms.
In this method, one starts with p random integer numbers ai, i = 0, 1, . . . , p − 1 generated
somehow in advance. Then the new elements k, with k ≥ p, can be generated as

ak = ak−p+q ⊕ ak−p

where p and q are constants, p > q , and ⊕ is the XOR logical operation.
(The XOR logical operation on two bits returns 1 is exactly one of the bits is 1, 0 otherwise.)

The difficult part is generating the initial numbers. In the original paper an older random
number generator, so called Kendall sequence, was used to generate a first number of binary
“1111100011011101010000100101100”, for p = 5. This number is the first number in a so called
W sequence. The second number in the W is then the same but shifted forwards by 6 bits, i.e.
1011001111100011011101010000100’. Then the next number is again shifted forwards, and so on.
This process gives the following numbers (the shift point is indicated by the horizontal line).

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 33

 Index Binary
0 1111100011011101010000100|101100
1 1011001111100011011101010|000100
2 0001001011001111100011011|101010
3 1010100001001011001111100|011011
4 0110111010100001001011001|111100
Note that one more shift would give back the original number, so no more shifts are useful.
Now the first five columns of this binary matrix is used as the originating sequence, i.e. we have
the first 5 numbers as
Element Binary
a0 11010
a1 10001
a2 11011
a3 11100
a4 10011
This example is used for p = 5, q = 2. Let us look at in detail what happens after this. We
have the numbers k = 0, . . . , 4 from above. We now want the number k = 5. Then from the
definition of the algorithm

a5 = a5−5+2 ⊕ a5−5 = a2 ⊕ a0 (3)

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 34

The XOR operation can easily be done manually and gives

a5 = 00001 (4)

Then the next step would give

a6 = a3 ⊕ a1 = 01101 (5)
and so on...
If we actually look at the original numbers in the first table we see these are just the columns of the
number!
In algorithmic form the GFSR generator can be written as

1◦ k=0
2◦ Set a0, . . . , ap−1 to some suitable initial values
3◦ Output ak
4◦ ak = ak+q mod p ⊕ ak
5◦ k = k + 1 mod p
6◦ Goto step 3◦

and as an actual C code (headers excluded)

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 35

main()
{
unsigned int xn[5];
unsigned int p,q;
unsigned int k;

p=5; q=2;

/* GFSR generator of original example in Lewis/Payne paper */

xn[0]=26;
xn[1]=17;
xn[2]=27;
xn[3]=28;
xn[4]=19;

k=0;
while(1) {
printf("%d\n",xn[k]);
xn[k]=xn[ ((k+q)%p) ] ^ xn[k];
k=((k+1)%p);
}

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 36

}

If we run this, we see that the period is 31. This actually follows from the full derivation and
motivation: it can be shown that the period of the GFSR generator is
p
2 −1 (6)

From this we also observe an obvious weakness of this original formulation: to get the initial values
we actually already generated all the 31 independent numbers the method can produce.
• This is, however, not really a general problem: Lewis and Payne showed that more
generally one needs to form a p × p matrix with linearly independent rows, and gave a
recipe for how to do that.

With that approach, the method actually does have some clear advantages:
1. It is extremely fast after initialization: note above that a single loop step, without the print,
only needs three memory references, one XOR and a couple of additions and moduli. All of these
can easily be made in a single bit-level machine language operation (for instance a multiplication
is on the bit-level waaay slower than an addition or logical operation). Hence it is well suited
for situations where very fast generation of random numbers is needed with minimal demands on
hardware, e.g. in small devices without too much computing power like mobile phones.

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 37

2. By increasing p one can rapidly increase the period.
The quality of the GFSR algorithms clearly depends on the choice of p and q . Like the parameters
in the congruential algorithms, these have to be carefully chosen. The smallest ones are quite bad,
but sequences where p is of the order of 1000 or more pass even quite advanced tests [Vattulainen, PRL
73 (1994) 2513].

The basic GFSR does have known weaknesses, the most obvious of which is that the schemes
to initialize the array in a good way are complicated and quite slow. But there are many
further developments of it, e.g. using several XOR operations, and non-linear versions.
[http://www.math.sci.hiroshima-u.ac.jp/∼m-mat/MT/ARTICLES/tgfsr3.pdf]. We will discuss one of those later on.

4.4.5.1. Historical note

The original implementation of the XOR part (not the initialization of the method) is

FUNCTION RAND(M,P,Q,INTSIZ)
C
C M(P)=TABLE OF P PREVIOUS RANDCM NUMBERS.
C P,O=POLYNOMIAL PARAMETERS:X**P+X**Q+1.
C .NOT. OPERATOR IMPLEMENTED IN ARITHMETIC.
C INTSIZ=INTEGER SIZE (BITS) OF HOST MACHINE: E.G.,
C IBM 360, 31; CDC 6000, 48; SRU 1100, 35; HP 2100, 15.
C

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 38

 LOGICAL AA,BB,LCOMPJ,LCOMPK
INTEGER A,B,P,Q,INTSIZ,M(1)
EQUIVALENCE (AA,A),(BB,B),(MCOMPJ,LCOMPJ),(MCOMPK,LCOMPK)
DATA J/0/
N=(2**(INTSIZ-1)-1)*2+1
J=J+l
IF(J.GT.P) J=1
K=J+Q
IF(K.GT.P) K=K-P
MCOMPJ=N-M(J)
MCOMPK=N-M(K)
A=M(K)
B=M(J)
BB=LCOMPJ.AND.AA.OR.LCOMPK.AND.BB
M(J)=B
RAND=FLOAT(M(J))/FLOAT(N)
RETURN
END

which is pretty hard to decipher for anyone not grown up with 1960’s Fortran... Note the extremely
dirty trick of using EQUIVALENCE between logical and integer variables and the construction of
an XOR operation by a combination of AND and OR. At that time Fortran did not have built-in
bit-level logical operations, which lead to the above implementation. The authors simply assumed
that a logical AND and OR also worked as bit-level operators on the whole word size. They then

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 39

tested it on 4 different computers of the day and concluded essentially “it seems to work fine on all
of them”...

4.4.6. Nonlinear generators

[http://random.mat.sbg.ac.at/software/; especially ftp://random.mat.sbg.ac.at/pub/data/weingaThesis.ps]

Many of the most modern and promising generators are based on the idea of using generators which
have nonlinear properties.
4.4.6.1. Nonlinear congruential generators
The basic idea of the inverse congruential generators is to generate numbers using

yn+1 = aȳn + b (mod M )

This at first looks exactly like the equation for the ordinary linear congruential generators. But the
difference is given by the bar on y . This signifies the solution for the equation

cc̄ = 1 (mod M )

where c is given and c̄ is the unknown. So instead of directly calculating the new number from yn,
one first calculates the congruence ȳn.

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 40

Calculating c̄ is actually not quite easy, and won’t be discussed here. The method to do can
be found among the software given on the link above (look for the prng-3.0.tar library and
subroutines prng inverse ? there).
This basic generator is already clearly better than the ordinary linear generators when the constants
constants and modulus are comparable. It has been further developed by replacing the ȳ operation
with other operators.

Another development is simply to add a quadratic term to the linear generator to produce the
quadratic congruential generator,

2
yn+1 = ayn + byn + c (mod M )

Several other developments of generators along this line exist.

4.4.6.2. Nonlinear GFSR generators
Matsumoto and Kurita have recently (1990’s) developed the basic GFSR generator further. The
basic change is called a “twist”, i.e. instead of using the simple GFSR operation of

ak = ak−p+q ⊕ ak−p (7)

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 41

one uses

ak = ak−p+q ⊕ ak−pA (8)

where A is a binary matrix of dimensions w × w where w is the word size (e.g. 32 for an unsigned
4-byte integer). ak−p must here be considered a row vector for the vector×matrix operation to
make sense.
This is a “twisted GFSR” or TGFSR generator. [Matsumoto and Kurita, 1992 ACM transactions on modelling
and Computer Similations 2 (179-194); http://www.math.sci.hiroshima-u.ac.jp/ m-mat/MT/ARTICLES/tgfsr3.pdf]. With suitable
choices of p, q and A this generator can achieve the maximal period of

pw
2 −1 (9)

which is about a factor of 2w more than the period of the basic GFSR. This is a huge improvement,
of course; even for 4-byte integers it is 232 more.

An additional major advantage of the TGFSR is that A can be chosen to make the algorithmic
implementation almost as simple as that of the basic GFSR. That is, the algorithm can be identical
to the GFSR one given above except that step 4 changes. The TGFSR algorithm is

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 42

 1◦ k=0
2◦ Set a0, . . . , ap−1 to some suitable initial values
3◦ Output ak

0 if least significant bit of ak = 0
4◦ ak = ak+q mod p ⊕ shiftrightak ⊕
a if least significant bit ofak = 1
5◦ k = k + 1 mod p
6◦ Goto step 3◦

Here a is a fixed constant binary vector of size w and the shift means a one-bit shift operation.

The TGFSR generator is proven to have several major advantages over the GFSR and pass quite
advanced tests. One obvious practical advantage is that the number of initialization vectors p does
not need to be very large to achieve a long period. E.g. already w = 16, p = 25, q = 11, a =
A875 (last one is hexadecimal) gives a period of 2400 − 1.
• Moreover, also the initialization vector is much less sensitive than in the GFSR, any
initialization vector (except all zero) will in fact do since the algorithm in any case goes
through the full possible set of ways to obtains sequences of p integers from the set of 2w
possible integers.

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 43

A further development of the TGFSR is the Mersenne Twister from 1996-1997. It is based on the
same idea as the TGFSR, the difference is that now the crucial step is

u t
ak = ak−p+q ⊕ (ak−p|ak−p+1)A (10)

where auk−p denotes the upper w − r bits of the vector ak−p, and atk−p+1 the lower r bits of the
vector ak−p+1. That is r is now an extra constant between 0 and w. The alert reader will notice
that if r = 0 this reduces back to the TGFSR.
• The A can still be handled in the same manner as in the TGFSR, and initialization is
equally easy.

With suitable good choices of the parameters this generator can achieve a period of

pw−r
2 −1 (11)

The values recommended by the authors are (w, p, q, r) = (32, 624, 397, 31) and a =
9908B0DF , which gives the period
19937
2 −1
which should be enough for a few years to come... The name “Mersenne” comes from the choice
of w, p and r to give a period which is a so called Mersenne prime number.

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 44

In addition to being proven to have many good features, this generator also has the advantage of
being readily available in source code from the web, see http://www.math.sci.hiroshima-u.ac.jp/∼m-
mat/MT/emt.html
The Mersenne twister has passed a wide range of tests on random numbers, and is despite its fairly
complex nature very efficient computationally; sometimes even faster than the standard C language
generator.

4.4.7. Combined generators

It is of course also possible to combine generators of different type, which when carried out well
should minimize the chances of the algorithm of one generator messing things up. One famous
generator which seems to be generally held in good respect is the RANMAR generator by G.
Marsaglia (source code with built-in documentation ranmar.f available on the course home page).
It has a period of ≈ 2144 which also should be good enough for some time to come...

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 45

 4.5. Tests of generators

4.5.1. Theoretical tests

[Mainly from G+T 12.19, also parts from Knuth]

We have already seen several very basic tests of algorithms, but we will reiterate what they are.
0. The very first thing to consider is that the period of the generator is larger than the number of
random numbers needed in your simulations. Testing for this simply amounts to finding or knowing
the period.

1. Another very basic test for a generator is of course to check that it produces a uniform distribution
between 0 and 1 in 1D. But failing this test is so basic that no sensible bug-free generator does it.

2. Many-dimensional smoothness tests. As an extension of the above, an important test is that the
random-number distribution is flat also in many dimensions. Here actually many generators already
start to have problems, as we shall see in the exercises.
Testing for problems 1 and 2 is simple, just doing the plots is the basic tests in 1, 2 and 3D.
Zooming in on parts of the test region may be useful in case the overall region looks OK.

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 46

3. χ2 test
One basic test for generators which have passed test 1, is to repeat the test but look at the
fluctuations from the average value, when one generates some finite number of N random numbers.
(FIGURE DRAWN DURING LECTURE)
• I.e. if we say collect generate 100 random numbers between 0 and 1 and make statistics
of them in bins of width 0.1, one should not normally get exactly 10 in every bin, but say
8 in one, 13 in the next and so on.
– One can then look at the fluctuations, and find whether they are what is expected
from probability theory.
If we consider generating N random number and placing them in M bins, then it is clear that the
expected value Ei for bin i is Ei = N/M . If the observed value in each bin is yi, then we can
calculate the χ-square statistic of this test with

M
2
X (yi − Ei)2
χ =
i=1
Ei

For large values of M (M > 30) there should be a 50 % probability that χ2 > M − 2/3, and a
50 % probability that χ2 < M − 2/3. Now we can test the generator by calculating χ2 numerous

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 47

times (for different random number sequences of course), and seeing whether it has the average
M − 2/3.

It is also possible to predict the probabilities for different percentage points, and then calculate how
often χ2 exceeds this point. For instance, for M = 50 χ2 should exceed 1.52M at most 1 % of
the time. This may actually be a much stronger test of whether the distribution is what it should
be.
(n.b. Gould-Tobochnik has an error in their equation 12.19.c : they say χ2 ≤ M , when it should
be χ2 ≈ M .)

Finding correlation problems may be difficult. We described the “two small numbers in sequence”
problem above, and mentioned that other exist. In fact, the most complex correlation tests are the
empirical tests mentioned in the next section.
4. Autocorrelation tests
One way to look for short-range correlations in a sequence of random numbers xi is to use the
autocorrelation function

hxi+k xii − hxii2

C(k) =
hxixii − hxii2

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 48

where hxi+k xii is found by forming all possible products xi+k xi for a given k and dividing by the
number of terms in the product. C(k) should become zero when k → ∞.

4.5.2. Empirical tests

The tests described above were of a rather general nature, and a good generator should pass them
all. However, even these tests do almost certainly not prove that a generator is good enough.
A logical continuation of testing is to use tests which are close to the system studied. For instance,
if one simulates the Ising model, a good way to test the generator is to simulate the Ising model
in a case where the accurate answer is known. These tests are called empirical, physical or
application-specific.

The research group of Vattulainen, Ala-Nissilä and Kankaala (formerly at HIP, now at HUT) have
developed several such tests. Here is a quick overview of some of them; codes to test them can be
found in http://www.physics.helsinki.fi/∼vattulai/rngs.html.

1. SN test. The test uses random walkers on a line and calculates the total number of sites visited
by the random walkers (versus the number of jumps made). In the test, N random walkers move
simultaneously without any interaction such that, at any jump attempt, they can make a jump to
the left or to the right with equal probability. After n >> 1 jumps by all random walkers, the

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 49

mean number of sites visited SN,n has an asymptotic form

1/2 x
SN,n ∼ log N n

where the exponent x should be 1/2 according to theory. The value of the exponent x observed
from simulations serves as a measure of correlations.

2. Interface roughening in 1+1 dimensions. In this case the roughening of a 1-dimensional interface
is followed with time. Consider two (independent) 1D random walks, which determine the heights
h1(n) and h2(n) of two interfaces versus the number of jumps made n. The height of the interface
between the two random walkers is then

h1(n) − h2(n),

whose height-height correlation function follows a power law in distance n where the roughening
exponent should be 1/2.

3. Ising model autocorrelation test. In this test, averages of some physical quantities such as
the energy, the susceptibility, and the updated cluster size in the 2D Ising model are calculated.
Additionally, their autocorrelation functions and corresponding integrated autocorrelation values are

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 50

determined. The exact value is known only for the energy; for the other quantities, the test works
by comparing results of different generators.
This is a fairly stringent test; e.g. RAN3 from Numerical Recipes first edition fails this test, as do
several of the GFSR generators.

4.5.3. So, what generator should I use?

To summarize all of this discussion, I give my personal view of what generator to use when.
I can see almost no reason to ever use anything less than the Park-Miller minimal generator. And
since even this has many known problems, it should be used only in cases where the random
numbers are of secondary importance. This can be for instance when only a few hundred random
numbers are needed for e.g. selecting impact points on a 2D surface, or when initializing velocities
of atoms in an MD simulation.
In cases where random numbers are used all through a large simulation run in crucial parts of the
code, I would recommend using something better than the Park-Miller generator.

Since there are no guarantees any generator is good enough for a problem which has not been
studied before, a good strategy would be to choose a few completely different generators and repeat
some of the central simulations with all of these. If no dependence on the choice of generator is

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 51

found, there probably is no problem. The generators chosen could be e.g. the RAN2 from numerical
Recipes second edition, the Mersenne twister and RANMAR.
(And remember that if you find problems with RAN2, you can claim the $1000 reward!)

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 52

 4.6. Generating non-uniform random numbers

[Numerical Recipes, Karimäki lecture notes]

So far we have only discussed generating random numbers in a uniform distribution, but at quite
some length. There is a good reason to this - random numbers in any other distribution are almost
always generated starting from random numbers distributed uniformly between 0 and 1.
But in the natural sciences it is clear that data often comes in many other forms. E.g. the peaks in
x-ray or γ spectra have a Gaussian or Lorentzian shape, the decay activity of a radioactive sample
follows an exponential function, and so on.
There are two basic approaches to generate other distributions than the uniform. The first is the
analytical one, the other the numerical von Neumann rejection method.

4.6.1. Analytical approach (inversion method)

We want to calculate random numbers which are distributed as some arbitrary function f (x). To
be a reasonable probability distribution, the function must have the properties

f (x) > 0 for all x

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 53

and Z ∞
f (x)dx = 1
−∞

Otherwise there are no limits on what f can be.

In the derivation, we will also need the cumulative distribution function (“kertymäfunktio”)

Z x
F (x) = f (t)dt
−∞

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 54

and its inverse function F −1(s),
−1
s = F (x) ⇐⇒ x = F (s)

Let us denote our generator for uniform random numbers Pu(0, 1). We now postulate that to
generate random numbers r distributed as f (x), we should perform the following steps:
1◦ Generate a uniformly distributed number u = Pu(0, 1)
2◦ Calculate x = F −1(u)

To prove this, we will show that the cumulative distribution function F 0(x) of the numbers x is the
function F (x). Since each function f has a unique cumulative function, this is enough as a proof.
Consider the probability that a point x is below the given point r ,
0
F (r) = P (x ≤ r)

This is the cumulative distribution function of the x, but we do not yet know what the function is.
But now, using (2◦ ),
0 −1
F (r) = P (x ≤ r) = P (F (u) ≤ r)

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 55

Now we can apply the function F on both sides of the inequality in the parentheses, and get

0 −1
F (r) = P (F (F (u)) ≤ F (r)) = P (u ≤ F (r))

But because u is just a uniformly distributed number, we have simply P (a ≤ b) = b and hence

0
F (r) = P (u ≤ F (r)) = F (r) q.e.d.

So the algorithm is very simple, but requires that it is possible to calculate the inverse function of
the integral of the function we want. Since all functions are not integrable, using this analytical
approach is not always possible.

Let us illustrate how this works in practice with an example. Say we want to have random numbers
with an exponential decay above 0, i.e.

e−x

x>0
f (x) =
0 otherwise

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 56

Now we first calculate
Z x Z x
−x −x
F (x) = f (t)dt = e =1−e
0 0

Then we check that the distribution is normalized, which is now easy:

−∞
F (∞) = 1 − e =1

(if it were not, we would need to add a normalization factor in front),

and then solve

−x −1
s = F (x) = 1 − e =⇒ x = − log(1 − s) i.e. F (u) = − log(1 − u)

But because u is a random number between 0 and 1, so is 1 − u, and we can reduce this to
−1
F (u) = − log(u)

To test whether this really works, I wrote the following small gawk script (since this is for demo
only, it is excusable to use the system random number generator):

gawk ’BEGIN {

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 57

 # Initialize random number generator from system clock
srand();
# Generate 10000 exponential deviates:
for(i=0;i<10000;i++) {
print -log(rand());
}
exit;
}’

The statistics of this is, from two runs:

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 58

((Hint for those of you using Unix: a single command line is enough to generate simple statistics
like this:
expdev | awk ’{ printf "%.1f\n",int($1*10)*0.1; }’ | sort -n | uniq -c |
awk ’{ print $2,$1 }’ > expdev.stat1
))

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 59

So it really is a nice exponential dependence (since it looks linear on a log scale).
4.6.1.1. Discrete distribution
For a discrete distribution, we can also use the inversion method. Let us say we have points
pi, i = 1, . . . , N which define the probability distribution function for some set of points xi. We
then have to generate the cumulative distribution function as a discrete function Fj , which can be
achieved by summation:
X j
Fj = pi
i=1
for all j = 1, . . . , N . We further have to set F0 = 0 and ensure that FN = 1 to handle the ends
correctly.

Then the generation algorithm becomes

1◦ Generate a uniformly distributed number u = Pu(0, 1)
2◦ Find k such that Fk−1 < u ≤ Fk
which gives integers k whose probability of occurring is proportional to pk .
Step 2◦ is very easy and efficient to do using a binary search.

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 60

Note that in case it is not possible to find an inverse F −1(x) for an analytical function, one
possibility is to tabulate it as a discrete distribution, then use the method described here.
• Interpolation may of course be used to make the answer more accurate

• However, note that if the set of points pi has a region where it is zero in the middle of the
distribution, Fj will be flat in this region. If you use interpolation schemes when doing
step 2◦, you may still create a point in the forbidden region
– may well lead to either your program crashing, or (even worse!) completely unphysical
results.

4.6.2. von Neumann rejection method

There is another way to generate random numbers, which is purely numerical and works for any
finite-valued function in a finite interval, regardless of whether it can be integrated or inverted. It
is called the (von Neumann) rejection method or hit-and-miss method
The idea is straightforward. Consider a function f (x) defined in some finite interval x ∈ [a, b]. It
has to be normalized to give probabilities. Let M be an number which is ≥ f (x) for any x in the
interval:

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 61

Now a rather obvious algorithm to generate a random number in this distribution is
1◦ Generate a uniformly distributed number x = Pu(a, b)
2◦ Generate a uniformly distributed number y = Pu(0, M )
3◦ If y > f (x) this is a miss: return to 1◦
4◦ Otherwise this is a hit: return x

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 62

This way we obtain random numbers x which are distributed according to the given distribution
f (x). Note that y only carries the role of a checking variable and is not returned.

This seems nice and easy. The downside is of course that we do some redundant work: all the
“miss” numbers were generated in vain. The probability to get a hit is

Rb
a
f (x)dx 1
P (hit) = =
M (b − a) M (b − a)

For a function like that plotted in the figure above, this is not too bad. But if the function is highly
peaked, or has a long weak tail, the number of misses will be enormous.

Consider for instance an exponential function e−x, x > 0. If the problem is such that one can
neglect very small values of e−x (large values of x) one can use some cutoff value b >> 1 in the
generation of random numbers. One could then use the hit-and-miss algorithm in the interval [0, b]
to generate random numbers for the exponential function. The normalized function would be

e−x e−x
f (x) = R b = −b
, x ∈ [0, b]
−x
e dx 1 − e
0

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 63

and the probability of misses

1 1
P (miss) = 1 − =1−
M (b − 0) Mb

and since the maximum of f (x) is at 0, we can use

1
M =
1 − e−b

and get
1 − e−b
P (miss) = 1 −
b
If for instance b = 100, we have ≈ 99 % misses, i.e. terrible efficiency. Obviously it would be
much better in this case to use the analytical approach.

4.6.3. Combined analytical-rejection method

Unfortunately for many functions it will not be possible to do the analytical approach. But there
may still be a way to do better than the basic hit-and-miss algorithm.
In case the shape of the function is known (which certainly almost always is the case in 1D), then
maybe we can find a function which is always larger than the one to be generated, but only slightly

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 64

so. Then if we can generate analytically random numbers for the larger function, we can again use
the hit-and-miss method, but with much less misses.

To put this more precisely, say we can find a function g(x) for which a constant A exists such that

Ag(x) ≥ f (x) for all x ∈ [a, b].

It is important to include the constant A here because both g(x) and f (x) are probabilities
normalized to one. For this to be useful, we further have to demand that it is possible to form the
inverse of the cumulative function G−1(x) of g(x).

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 65

Then the algorithm becomes:
1◦ Generate a uniformly distributed number u = Pu(0, 1)
2◦ Generate a number distributed as g(x): x = G−1(u)
3◦ Generate a uniformly distributed number y = Pu(0, Ag(x))
4◦ If y > f (x) this is a miss: return to 1◦
5◦ Otherwise this is a hit: return x

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 66

Note that we do not ever have to generate numbers in the interval [a, b] in this case, so the limits
can be at infinity as well. This can be a major advantage over the pure hit-and-miss method, which
requires a finite cutoff to be used for functions extending to infinity.

4.6.4. Generating Gaussian random numbers

The Gaussian function

1 x−µ 2
“ ”
1
−2
f (x) = √ e σ
2πσ
is of course one of the most common functions in science, and there are numerous applications where
one wants to generate random numbers distributed as f (x). Unfortunately it is not integrable, as
we all know, so using the inversion method directly is not possible.

But as you probably all also remember, the definite integral from −∞ to ∞ of f (x) can
be evaluated by a trick using two dimensions. For simplicity, let’s work with the the Gaussian
distribution in the form centered at 0 with σ = 1,

1 − 1 x2
f (x) = √ e 2
2π

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 67

The integral can be calculated by taking the square of the integral
»Z ∞ –2 Z ∞ Z ∞ Z ∞ Z ∞
−1 x2 1 x2 − 1 y 2
−2 −1 2 2
2 (x +y )
e 2 dx = e e 2 dxdy = e dxdy
−∞ −∞ −∞ −∞ −∞

and switching to polar coordinates r 2 = x2 + y 2, dxdy = rdrdφ:

Z ∞ Z 2π Z ∞ „ «
1
− r 2 1
− r 2 1 2
= e 2 rdrdφ = 2π e 2 d r = 2π
0 0 0 2

The Box-Muller method to generate random numbers with a Gaussian distribution relies on a
similar trick.
In this method, we also consider two Gaussian distributions in 2D. Their joint distribution is
1 − 1 x2 1 − 1 y2 1 − 1 (x2+y2)
f (x, y) = √ e 2 √ e 2 = e 2
2π 2π 2π

Switching again to polar coordinates, and remembering that the surface element transforms as
dxdy = rdrdφ we get the polar density distribution function
1 −1 r2
g(r, φ) = f (x, y)r = re 2
2π

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 68

We can now separate the r and φ contributions:

g(r, φ) = fφ(φ)fr (r)

where
1
fφ(φ) =
2π
−1 2
fr (r) = re 2r

So if we can generate fφ and fr separately, we will also be able to generate the joint 2D distribution
and hence two Gaussian numbers at a time.

Generating fφ is trivial, we just need to form a uniform number and multiply by 2π to get an even
1
distribution in the range [0, 2π], which has the value 2π everywhere. fr can also be handled since
Z r 2 1 r2
−1
2r −2
Fr (r) = re =1−e
0

which can be inverted to give

q
−1
Fr (u) = −2 log(1 − u)

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 69

So we can obtain both r and φ. After this the x and y can be obtained easily using

x = r cos φ
y = r sin φ

So the polar or Box-Muller algorithm becomes

1◦ Generate a uniformly distributed number u1 = Pu(0, 1)
2◦ Calculate φ = 2πu1
3◦ Generate a uniformly distributed number u2 = Pu(0, 1)
◦ √
4 Calculate r = −2 log u2
5◦ Obtain x and y using

x = r cos φ
y = r sin φ

So this gives two numbers at a time. In practice, the subroutine is best written so that steps 1◦ -
5◦ are carried out every second step to get both x and y . The first step returns x and stores y .
The second step only returns y without any calculations at all.

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 70

Algorithmically this method is fine, there is no wasted effort in terms of “misses” or something like
that.
• But in terms of computational efficiency, it leaves much to be hoped for. It requires
√
computation of the functions , log, sin and cos.
– Calculating all of these, and especially the last three ones, is excruciatingly slow
compared to the simple arithmetic operations. So it would be nice to have a more
efficient routine.

It turns out that a variety of a hit-and-miss algorithm is actually almost always faster. The trick is
to avoid having to calculate the sine and cosine explicitly.
Here we consider a simple unit circle:

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 71

and generate a point (v1, v2) inside a unit circle centered at the origin:
1◦ Obtain v1 = Pu(−1, 1) and v2 = Pu(−1, 1) and w = v12 + v22
Then we check whether the point is really inside the circle:
2◦ If w ≥ 1 return to step 1◦
Then from basic geometry we know that the sine and cosine for the φ (which we do not know) of
the point (v1, v2) are
 √
cos φ = v1/ w
√
sin φ = v2/ w

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 72

(Here we have to make sure w can not be exactly 0!). Moreover, w is of the form Pu(0, 1). Now
we can again obtain the desired x and y :
√
3◦ Calculate r = −2 log w
4◦ Calculate  √
x = r cos φ = rv1/ w
√
y = r sin φ = rv2/ w

The advantage here is that by using the latter parts of the equation above, we do not ever have to
calculate the sine and cosine explicitly. This makes this approach faster, even though we do have
to reject
4 − π12 π
= 1 − ≈ 21%
4 4
of all values in the hit-and-miss steps 1 and 2.

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 73

 4.7. Generating random numbers on the surface of a sphere

[Own knowledge and derivation]

I will not discuss generating random numbers for multidimensional cases. The very basic cases are
straightforward extensions of the 1D case: an N -dimensional uniform distribution is just a vector
of N 1-dimensional deviates, and the hit-and-miss approach works in any dimensions. If interested
on generating Gaussian distributions in multidimensions, a good discussion can be found in the MC
lectures of Karimäki (available on the course web page).

But there is one serious pitfall related to generating random deviates in 2 dimensions that it needs
to be mentioned here separately.

The problem is simply to select a random direction in 3 dimensions. That is, in spherical coordinates
any vector can of course be given in the unit system (r, θ, φ). To give a random direction, one then
simply has to generate θ and φ randomly. So the obvious solution would seem to be to generate

θ = πPu(0, 1)
(WRONG)
φ = 2πPu(0, 1)

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 74

This is outright wrong!

Why is this? Consider the unit sphere r = 1. To generate random numbers on the surface of it
(which is equivalent to generating a random direction), we want to have an equal number of points
per surface area everywhere on the sphere. The surface area element dσ is

sin θdθdφ

i.e. the area is not an even function of θ . Hence one has to generate the θ random numbers
distributed as
sin θ

rather than uniformly. Fortunately we know how to do this from section .

Let’s for generality consider an arbitrary angular interval [α, β] lying inside the [0, π] region. Then
the normalization factor is

Z β
fN = sin(θ)dθ = cos α − cos β
α

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 75

and the normalized probability function is

sin θ
f (θ) =
cos α − cos β

The cumulative function for an arbitrary interval [α, β] is

θ
sin θ cos α − cos θ
Z
F (θ) = dθ =
α cos α − cos β cos α − cos β

and the inverse can be solved as follows:

cos α − cos θ
u=
cos α − cos β

=⇒ − cos θ = u(cos α − cos β) − cos α

whence
−1
θ = cos (cos α − u(cos α − cos β))

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 76

For the original case of a random direction anywhere in 3D, α = 0 and β = π and

−1
θ = cos (1 − 2u)

which gives the correct algorithm

θ = cos−1 (1 − 2Pu(0, 1))



φ = 2πPu(0, 1)

If you do not believe me, compare the following two plots. They show 10000 points generated on
the unit sphere, the left one using the wrong way to generate θ , the right one using the correct
way. The rotation is exactly the same in both figures, so they really are comparable.

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 77

The problem is quite obvious in the left figure.

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 78

 4.8. Non-random (quasi-) random distributions

Finally, we discuss pseudo-random numbers which are not random at all. Still, in some cases
using these may be more efficient than using real pseudorandom numbers. This can be understood
intuitively from the following two figures:

The left one shows 30000 random points in the unit box, the right one 300 points. From the left
side we see that there is nothing wrong with the distribution itself (on this scale at least): it fills

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 79

the space. But the right side shows that with low statistics, there are clumps in the data here and
there, and gaps elsewhere. This is as it should be for truly random numbers.

But this brings us to think that in case we want to integrate over the 2D plane, would it not be
more efficient to arrange the numbers so that it is guaranteed that even for low statistics points lie
fairly smoothly everywhere in the interval? This is the idea behind quasirandom numbers.
In fact, it has been shown that the convergence in MC integration behaves as follows:
∝ √1 for true and pseudo-random numbers.
n
1
∝ n at best for quasi-random numbers

4.8.1. Stratified sampling

The simplest imaginable way to achieve this is to use fully stratified sampling.
• This means that we first decide how many points we want to generate, the divide our
sample into exactly this many equal-sized boxes.
• Then we generate exactly one random number inside each box.
This way it is guaranteed that the point distribution is fairly smooth.
This is illustrated in the following two figures: one has 1024 completely random numbers in 2D, the

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 80

other 1024 numbers generated with fully stratified MC using a 32 × 32 grid (you can figure out
yourself which one is which):

Coding this is of course trivial, here is an example:

# Size of integration box

box=1.0;
# Interval size
gridsize=32;
for(i=0;i<gridsize;i++) {

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 81

 for(j=0;j<gridsize;j++) {
x=(i+rand())*(box/gridsize);
y=(j+rand())*(box/gridsize);
print x,y;
# Evaluate function to be integrated here
}
}

This method has a significant practical drawback, however: one has to decide on the number of
points needed in advance, and any intermediate result is worthless since parts of the 2D space have
not been examined at all. One could introduce a mixing scheme which selects the minor boxes in
random order to overcome the latter problem, but this would not solve the former.

A somewhat more flexible solution is to use partially stratified sampling. In here, we also divide
the integration sample into boxes, but then select several points per box. This has the advantage
that we can choose a somewhat smaller number of boxes, then do several loops where in every
loop we select one point per box. This way we can stop the simulation anytime the outermost loop
finishes.
This is probably clearer in code:

# Size of integration box

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 82

 box=1.0;
# Number of intervals to do
ninterval=4;
# Interval size
gridsize=16;
for (interval=0;interval<ninterval;interval++) {
for(i=0;i<gridsize;i++) {
for(j=0;j<gridsize;j++) {
x=(i+rand())*(box/gridsize);
y=(j+rand())*(box/gridsize);
print x,y;
# Evaluate function to be integrated here
}
}
# Intermediate result could be printed and simulation
# stopped here
}

So we could in this example stop the run after 256, 512 or 768 steps in case we see we already have
good enough statistics.
But even this is not very convenient. In large simulations, where the actual evaluation of a function

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 83

can take hours or even days, the runs often die or have to be stopped in the middle, to be restarted
later. For the stratified MC schemes one then would have to devise a restart which knows at which
i and j values to restart. Quite possible to implement, but rather tedious.

4.8.2. Quasi-random numbers

[Karimäki notes p. 9, Numerical Recipes]

Fortunately it is quite possible to implement a number sequence such that it fills space evenly
without having any specific “filling interval” like the stratified schemes. Such random number
sequences are called low-discrepancy numbers or quasi-random numbers. Contrary to true
random numbers, they are designed to be highly correlated, in a way such that they will fill space
fairly evenly.
We give three examples of methods to generate such sequences.
4.8.2.1. Richtmeyer-sequences
We want to have a set of vectors xi in k dimensions

xi = (xi1, xi2, . . . , xik )

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 84

Each vector element is obtained as follows:

p
xij = i Nj (mod 1)

where (mod 1) means we take the decimal part, and Nj is the j :th prime number.
This is simple in concept, but not very efficient if k is large, since one also needs a list of prime
numbers.
4.8.2.2. Van der Corput-sequences
This method is so simple that it can be easily explained, but in the default version it works only in
1D (you can think about yourself whether you could invent an extension to 2D).
The method is as follows:
1◦ Take an integer, and write it in some base (e.g. binary)
2◦ Reverse the digits in the number
3◦ Put a decimal point in front and interpret it as a decimal
Using e.g. the binary system, we get:

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 85

4.8.2.3. Sobol sequences
The Sobol sequence is somewhat more complicated in definition, based on XOR operations and
needing a list of initializing numbers. The code (and explanation of how it works) can be found in
Numerical Recipes 2nd ed. chapter 7.7.
But here is an illustration of how it works:

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 86

We see that there are gaps, but these are always filled in on successive calling on the routine.
The Sobol sequence is not necessarily quite as efficient as an ideal quasi-random number generator.

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 87

The Sobol sequence efficiency for MC integration in n dimensions is
n
O((ln N ) /N
√
whereas the optimal efficiency is 1/N . Still, this is not bad at all compared to the 1/ N efficiency
obtained by basic pseudo-random numbers.
To summarize, if you think quasi-random numbers might work better than pseudorandom numbers
in your application, what should you do? First check the literature on whether someone has
examined this in a problem similar to yours. If not, simply test it. I personally would first go for
partially stratified MC or a Sobol sequence (since the code is readily available). And if the number
of dimensions is low (< 6), I would also test a regular grid (no randomness at all!).

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 88

 4.9. Final remark

Back in the 1940’s, when the very first generators were designed, von Neumann thought that:

Anyone who considers arithmetic methods of

producing random digits
is, of course, in a state if sin.
J. von Neumann

Now, when generators tend to be at least fairly decent and are very widely used, this fortunately is
no longer true!

Basics of Monte Carlo simulations, Kai Nordlund 2006 JJ J  I II × 89