|
Packit |
67cb25 |
.. index:: random number generators
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
************************
|
|
Packit |
67cb25 |
Random Number Generation
|
|
Packit |
67cb25 |
************************
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. include:: include.rst
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The library provides a large collection of random number generators
|
|
Packit |
67cb25 |
which can be accessed through a uniform interface. Environment
|
|
Packit |
67cb25 |
variables allow you to select different generators and seeds at runtime,
|
|
Packit |
67cb25 |
so that you can easily switch between generators without needing to
|
|
Packit |
67cb25 |
recompile your program. Each instance of a generator keeps track of its
|
|
Packit |
67cb25 |
own state, allowing the generators to be used in multi-threaded
|
|
Packit |
67cb25 |
programs. Additional functions are available for transforming uniform
|
|
Packit |
67cb25 |
random numbers into samples from continuous or discrete probability
|
|
Packit |
67cb25 |
distributions such as the Gaussian, log-normal or Poisson distributions.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These functions are declared in the header file :file:`gsl_rng.h`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. Need to explain the difference between SERIAL and PARALLEL random
|
|
Packit |
67cb25 |
.. number generators here
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
General comments on random numbers
|
|
Packit |
67cb25 |
==================================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
In 1988, Park and Miller wrote a paper entitled "Random number
|
|
Packit |
67cb25 |
generators: good ones are hard to find." [Commun.: ACM, 31, 1192--1201].
|
|
Packit |
67cb25 |
Fortunately, some excellent random number generators are available,
|
|
Packit |
67cb25 |
though poor ones are still in common use. You may be happy with the
|
|
Packit |
67cb25 |
system-supplied random number generator on your computer, but you should
|
|
Packit |
67cb25 |
be aware that as computers get faster, requirements on random number
|
|
Packit |
67cb25 |
generators increase. Nowadays, a simulation that calls a random number
|
|
Packit |
67cb25 |
generator millions of times can often finish before you can make it down
|
|
Packit |
67cb25 |
the hall to the coffee machine and back.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A very nice review of random number generators was written by Pierre
|
|
Packit |
67cb25 |
L'Ecuyer, as Chapter 4 of the book: Handbook on Simulation, Jerry Banks,
|
|
Packit |
67cb25 |
ed. (Wiley, 1997). The chapter is available in postscript from
|
|
Packit |
67cb25 |
L'Ecuyer's ftp site (see references). Knuth's volume on Seminumerical
|
|
Packit |
67cb25 |
Algorithms (originally published in 1968) devotes 170 pages to random
|
|
Packit |
67cb25 |
number generators, and has recently been updated in its 3rd edition
|
|
Packit |
67cb25 |
(1997).
|
|
Packit |
67cb25 |
It is brilliant, a classic. If you don't own it, you should stop reading
|
|
Packit |
67cb25 |
right now, run to the nearest bookstore, and buy it.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A good random number generator will satisfy both theoretical and
|
|
Packit |
67cb25 |
statistical properties. Theoretical properties are often hard to obtain
|
|
Packit |
67cb25 |
(they require real math!), but one prefers a random number generator
|
|
Packit |
67cb25 |
with a long period, low serial correlation, and a tendency **not** to
|
|
Packit |
67cb25 |
"fall mainly on the planes." Statistical tests are performed with
|
|
Packit |
67cb25 |
numerical simulations. Generally, a random number generator is used to
|
|
Packit |
67cb25 |
estimate some quantity for which the theory of probability provides an
|
|
Packit |
67cb25 |
exact answer. Comparison to this exact answer provides a measure of
|
|
Packit |
67cb25 |
"randomness".
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The Random Number Generator Interface
|
|
Packit |
67cb25 |
=====================================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
It is important to remember that a random number generator is not a
|
|
Packit |
67cb25 |
"real" function like sine or cosine. Unlike real functions, successive
|
|
Packit |
67cb25 |
calls to a random number generator yield different return values. Of
|
|
Packit |
67cb25 |
course that is just what you want for a random number generator, but to
|
|
Packit |
67cb25 |
achieve this effect, the generator must keep track of some kind of
|
|
Packit |
67cb25 |
"state" variable. Sometimes this state is just an integer (sometimes
|
|
Packit |
67cb25 |
just the value of the previously generated random number), but often it
|
|
Packit |
67cb25 |
is more complicated than that and may involve a whole array of numbers,
|
|
Packit |
67cb25 |
possibly with some indices thrown in. To use the random number
|
|
Packit |
67cb25 |
generators, you do not need to know the details of what comprises the
|
|
Packit |
67cb25 |
state, and besides that varies from algorithm to algorithm.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. type:: gsl_rng_type
|
|
Packit |
67cb25 |
gsl_rng
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The random number generator library uses two special structs,
|
|
Packit |
67cb25 |
:type:`gsl_rng_type` which holds static information about each type of
|
|
Packit |
67cb25 |
generator and :type:`gsl_rng` which describes an instance of a generator
|
|
Packit |
67cb25 |
created from a given :type:`gsl_rng_type`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The functions described in this section are declared in the header file
|
|
Packit |
67cb25 |
:file:`gsl_rng.h`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Random number generator initialization
|
|
Packit |
67cb25 |
======================================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: gsl_rng * gsl_rng_alloc (const gsl_rng_type * T)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function returns a pointer to a newly-created
|
|
Packit |
67cb25 |
instance of a random number generator of type :data:`T`.
|
|
Packit |
67cb25 |
For example, the following code creates an instance of the Tausworthe
|
|
Packit |
67cb25 |
generator::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_rng * r = gsl_rng_alloc (gsl_rng_taus);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
If there is insufficient memory to create the generator then the
|
|
Packit |
67cb25 |
function returns a null pointer and the error handler is invoked with an
|
|
Packit |
67cb25 |
error code of :macro:`GSL_ENOMEM`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The generator is automatically initialized with the default seed,
|
|
Packit |
67cb25 |
:data:`gsl_rng_default_seed`. This is zero by default but can be changed
|
|
Packit |
67cb25 |
either directly or by using the environment variable :macro:`GSL_RNG_SEED`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The details of the available generator types are
|
|
Packit |
67cb25 |
described later in this chapter.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: void gsl_rng_set (const gsl_rng * r, unsigned long int s)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function initializes (or "seeds") the random number generator. If
|
|
Packit |
67cb25 |
the generator is seeded with the same value of :data:`s` on two different
|
|
Packit |
67cb25 |
runs, the same stream of random numbers will be generated by successive
|
|
Packit |
67cb25 |
calls to the routines below. If different values of :math:`s \geq 1`
|
|
Packit |
67cb25 |
are supplied, then the generated streams of random
|
|
Packit |
67cb25 |
numbers should be completely different. If the seed :data:`s` is zero
|
|
Packit |
67cb25 |
then the standard seed from the original implementation is used
|
|
Packit |
67cb25 |
instead. For example, the original Fortran source code for the
|
|
Packit |
67cb25 |
:code:`ranlux` generator used a seed of 314159265, and so choosing
|
|
Packit |
67cb25 |
:data:`s` equal to zero reproduces this when using
|
|
Packit |
67cb25 |
:data:`gsl_rng_ranlux`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
When using multiple seeds with the same generator, choose seed values
|
|
Packit |
67cb25 |
greater than zero to avoid collisions with the default setting.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Note that the most generators only accept 32-bit seeds, with higher
|
|
Packit |
67cb25 |
values being reduced modulo :math:`2^{32}`.
|
|
Packit |
67cb25 |
For generators with smaller ranges the maximum seed value will typically be lower.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: void gsl_rng_free (gsl_rng * r)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function frees all the memory associated with the generator
|
|
Packit |
67cb25 |
:data:`r`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Sampling from a random number generator
|
|
Packit |
67cb25 |
=======================================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The following functions return uniformly distributed random numbers,
|
|
Packit |
67cb25 |
either as integers or double precision floating point numbers. |inlinefns|
|
|
Packit |
67cb25 |
To obtain non-uniform distributions, see :ref:`chap_random-number-distributions`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: unsigned long int gsl_rng_get (const gsl_rng * r)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function returns a random integer from the generator :data:`r`. The
|
|
Packit |
67cb25 |
minimum and maximum values depend on the algorithm used, but all
|
|
Packit |
67cb25 |
integers in the range [:data:`min`, :data:`max`] are equally likely. The
|
|
Packit |
67cb25 |
values of :data:`min` and :data:`max` can be determined using the auxiliary
|
|
Packit |
67cb25 |
functions :func:`gsl_rng_max` and :func:`gsl_rng_min`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: double gsl_rng_uniform (const gsl_rng * r)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function returns a double precision floating point number uniformly
|
|
Packit |
67cb25 |
distributed in the range [0,1). The range includes 0.0 but excludes 1.0.
|
|
Packit |
67cb25 |
The value is typically obtained by dividing the result of
|
|
Packit |
67cb25 |
:code:`gsl_rng_get(r)` by :code:`gsl_rng_max(r) + 1.0` in double
|
|
Packit |
67cb25 |
precision. Some generators compute this ratio internally so that they
|
|
Packit |
67cb25 |
can provide floating point numbers with more than 32 bits of randomness
|
|
Packit |
67cb25 |
(the maximum number of bits that can be portably represented in a single
|
|
Packit |
67cb25 |
:code:`unsigned long int`).
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: double gsl_rng_uniform_pos (const gsl_rng * r)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function returns a positive double precision floating point number
|
|
Packit |
67cb25 |
uniformly distributed in the range (0,1), excluding both 0.0 and 1.0.
|
|
Packit |
67cb25 |
The number is obtained by sampling the generator with the algorithm of
|
|
Packit |
67cb25 |
:func:`gsl_rng_uniform` until a non-zero value is obtained. You can use
|
|
Packit |
67cb25 |
this function if you need to avoid a singularity at 0.0.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: unsigned long int gsl_rng_uniform_int (const gsl_rng * r, unsigned long int n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function returns a random integer from 0 to :math:`n-1` inclusive
|
|
Packit |
67cb25 |
by scaling down and/or discarding samples from the generator :data:`r`.
|
|
Packit |
67cb25 |
All integers in the range :math:`[0,n-1]` are produced with equal
|
|
Packit |
67cb25 |
probability. For generators with a non-zero minimum value an offset
|
|
Packit |
67cb25 |
is applied so that zero is returned with the correct probability.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Note that this function is designed for sampling from ranges smaller
|
|
Packit |
67cb25 |
than the range of the underlying generator. The parameter :data:`n`
|
|
Packit |
67cb25 |
must be less than or equal to the range of the generator :data:`r`.
|
|
Packit |
67cb25 |
If :data:`n` is larger than the range of the generator then the function
|
|
Packit |
67cb25 |
calls the error handler with an error code of :macro:`GSL_EINVAL` and
|
|
Packit |
67cb25 |
returns zero.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
In particular, this function is not intended for generating the full range of
|
|
Packit |
67cb25 |
unsigned integer values :math:`[0,2^{32}-1]`.
|
|
Packit |
67cb25 |
Instead choose a generator with the maximal integer range and zero minimum
|
|
Packit |
67cb25 |
value, such as :data:`gsl_rng_ranlxd1`, :data:`gsl_rng_mt19937` or
|
|
Packit |
67cb25 |
:data:`gsl_rng_taus`, and sample it directly using
|
|
Packit |
67cb25 |
:func:`gsl_rng_get`. The range of each generator can be found using
|
|
Packit |
67cb25 |
the auxiliary functions described in the next section.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Auxiliary random number generator functions
|
|
Packit |
67cb25 |
===========================================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The following functions provide information about an existing
|
|
Packit |
67cb25 |
generator. You should use them in preference to hard-coding the generator
|
|
Packit |
67cb25 |
parameters into your own code.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: const char * gsl_rng_name (const gsl_rng * r)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function returns a pointer to the name of the generator.
|
|
Packit |
67cb25 |
For example::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
printf ("r is a '%s' generator\n", gsl_rng_name (r));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
would print something like::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
r is a 'taus' generator
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: unsigned long int gsl_rng_max (const gsl_rng * r)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function returns the largest value that :func:`gsl_rng_get`
|
|
Packit |
67cb25 |
can return.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: unsigned long int gsl_rng_min (const gsl_rng * r)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function returns the smallest value that :func:`gsl_rng_get`
|
|
Packit |
67cb25 |
can return. Usually this value is zero. There are some generators with
|
|
Packit |
67cb25 |
algorithms that cannot return zero, and for these generators the minimum
|
|
Packit |
67cb25 |
value is 1.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: void * gsl_rng_state (const gsl_rng * r)
|
|
Packit |
67cb25 |
size_t gsl_rng_size (const gsl_rng * r)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These functions return a pointer to the state of generator :data:`r` and
|
|
Packit |
67cb25 |
its size. You can use this information to access the state directly. For
|
|
Packit |
67cb25 |
example, the following code will write the state of a generator to a
|
|
Packit |
67cb25 |
stream::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void * state = gsl_rng_state (r);
|
|
Packit |
67cb25 |
size_t n = gsl_rng_size (r);
|
|
Packit |
67cb25 |
fwrite (state, n, 1, stream);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: const gsl_rng_type ** gsl_rng_types_setup (void)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function returns a pointer to an array of all the available
|
|
Packit |
67cb25 |
generator types, terminated by a null pointer. The function should be
|
|
Packit |
67cb25 |
called once at the start of the program, if needed. The following code
|
|
Packit |
67cb25 |
fragment shows how to iterate over the array of generator types to print
|
|
Packit |
67cb25 |
the names of the available algorithms::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const gsl_rng_type **t, **t0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
t0 = gsl_rng_types_setup ();
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
printf ("Available generators:\n");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (t = t0; *t != 0; t++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
printf ("%s\n", (*t)->name);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Random number environment variables
|
|
Packit |
67cb25 |
===================================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The library allows you to choose a default generator and seed from the
|
|
Packit |
67cb25 |
environment variables :macro:`GSL_RNG_TYPE` and :macro:`GSL_RNG_SEED` and
|
|
Packit |
67cb25 |
the function :func:`gsl_rng_env_setup`. This makes it easy try out
|
|
Packit |
67cb25 |
different generators and seeds without having to recompile your program.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. macro:: GSL_RNG_TYPE
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This environment variable specifies the default random number generator.
|
|
Packit |
67cb25 |
It should be the name of a generator, such as :code:`taus` or :code:`mt19937`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. macro:: GSL_RNG_SEED
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This environment variable specifies the default seed for the random
|
|
Packit |
67cb25 |
number generator
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_rng_default
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This global library variable specifies the default random number generator,
|
|
Packit |
67cb25 |
and can be initialized from :macro:`GSL_RNG_TYPE` using :func:`gsl_rng_env_setup`.
|
|
Packit |
67cb25 |
It is defined as follows::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
extern const gsl_rng_type *gsl_rng_default
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_rng_default_seed
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This global library variable specifies the seed for the default random number generator,
|
|
Packit |
67cb25 |
and can be initialized from :macro:`GSL_RNG_SEED` using :func:`gsl_rng_env_setup`.
|
|
Packit |
67cb25 |
It is set to zero by default and is defined as follows::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
extern unsigned long int gsl_rng_default_seed
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: const gsl_rng_type * gsl_rng_env_setup (void)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function reads the environment variables :macro:`GSL_RNG_TYPE` and
|
|
Packit |
67cb25 |
:macro:`GSL_RNG_SEED` and uses their values to set the corresponding
|
|
Packit |
67cb25 |
library variables :data:`gsl_rng_default` and
|
|
Packit |
67cb25 |
:data:`gsl_rng_default_seed`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The value of :macro:`GSL_RNG_SEED` is converted to an :code:`unsigned long int`
|
|
Packit |
67cb25 |
using the C library function :func:`strtoul`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
If you don't specify a generator for :macro:`GSL_RNG_TYPE` then
|
|
Packit |
67cb25 |
:data:`gsl_rng_mt19937` is used as the default. The initial value of
|
|
Packit |
67cb25 |
:data:`gsl_rng_default_seed` is zero.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Here is a short program which shows how to create a global
|
|
Packit |
67cb25 |
generator using the environment variables :macro:`GSL_RNG_TYPE` and
|
|
Packit |
67cb25 |
:macro:`GSL_RNG_SEED`,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. include:: examples/rng.c
|
|
Packit |
67cb25 |
:code:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Running the program without any environment variables uses the initial
|
|
Packit |
67cb25 |
defaults, an :code:`mt19937` generator with a seed of 0,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. include:: examples/rng.txt
|
|
Packit |
67cb25 |
:code:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
By setting the two variables on the command line we can
|
|
Packit |
67cb25 |
change the default generator and the seed::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
$ GSL_RNG_TYPE="taus" GSL_RNG_SEED=123 ./a.out
|
|
Packit |
67cb25 |
GSL_RNG_TYPE=taus
|
|
Packit |
67cb25 |
GSL_RNG_SEED=123
|
|
Packit |
67cb25 |
generator type: taus
|
|
Packit |
67cb25 |
seed = 123
|
|
Packit |
67cb25 |
first value = 2720986350
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Copying random number generator state
|
|
Packit |
67cb25 |
=====================================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The above methods do not expose the random number state which changes
|
|
Packit |
67cb25 |
from call to call. It is often useful to be able to save and restore
|
|
Packit |
67cb25 |
the state. To permit these practices, a few somewhat more advanced
|
|
Packit |
67cb25 |
functions are supplied. These include:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_rng_memcpy (gsl_rng * dest, const gsl_rng * src)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function copies the random number generator :data:`src` into the
|
|
Packit |
67cb25 |
pre-existing generator :data:`dest`, making :data:`dest` into an exact copy
|
|
Packit |
67cb25 |
of :data:`src`. The two generators must be of the same type.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: gsl_rng * gsl_rng_clone (const gsl_rng * r)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function returns a pointer to a newly created generator which is an
|
|
Packit |
67cb25 |
exact copy of the generator :data:`r`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Reading and writing random number generator state
|
|
Packit |
67cb25 |
=================================================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The library provides functions for reading and writing the random
|
|
Packit |
67cb25 |
number state to a file as binary data.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_rng_fwrite (FILE * stream, const gsl_rng * r)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function writes the random number state of the random number
|
|
Packit |
67cb25 |
generator :data:`r` to the stream :data:`stream` in binary format. The
|
|
Packit |
67cb25 |
return value is 0 for success and :macro:`GSL_EFAILED` if there was a
|
|
Packit |
67cb25 |
problem writing to the file. Since the data is written in the native
|
|
Packit |
67cb25 |
binary format it may not be portable between different architectures.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. function:: int gsl_rng_fread (FILE * stream, gsl_rng * r)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This function reads the random number state into the random number
|
|
Packit |
67cb25 |
generator :data:`r` from the open stream :data:`stream` in binary format.
|
|
Packit |
67cb25 |
The random number generator :data:`r` must be preinitialized with the
|
|
Packit |
67cb25 |
correct random number generator type since type information is not
|
|
Packit |
67cb25 |
saved. The return value is 0 for success and :macro:`GSL_EFAILED` if
|
|
Packit |
67cb25 |
there was a problem reading from the file. The data is assumed to
|
|
Packit |
67cb25 |
have been written in the native binary format on the same
|
|
Packit |
67cb25 |
architecture.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Random number generator algorithms
|
|
Packit |
67cb25 |
==================================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The functions described above make no reference to the actual algorithm
|
|
Packit |
67cb25 |
used. This is deliberate so that you can switch algorithms without
|
|
Packit |
67cb25 |
having to change any of your application source code. The library
|
|
Packit |
67cb25 |
provides a large number of generators of different types, including
|
|
Packit |
67cb25 |
simulation quality generators, generators provided for compatibility
|
|
Packit |
67cb25 |
with other libraries and historical generators from the past.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The following generators are recommended for use in simulation. They
|
|
Packit |
67cb25 |
have extremely long periods, low correlation and pass most statistical
|
|
Packit |
67cb25 |
tests. For the most reliable source of uncorrelated numbers, the
|
|
Packit |
67cb25 |
second-generation RANLUX generators have the strongest proof of
|
|
Packit |
67cb25 |
randomness.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: MT19937 random number generator
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_rng_mt19937
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The MT19937 generator of Makoto Matsumoto and Takuji Nishimura is a
|
|
Packit |
67cb25 |
variant of the twisted generalized feedback shift-register algorithm,
|
|
Packit |
67cb25 |
and is known as the "Mersenne Twister" generator. It has a Mersenne
|
|
Packit |
67cb25 |
prime period of :math:`2^{19937} - 1`
|
|
Packit |
67cb25 |
(about :math:`10^{6000}`) and is
|
|
Packit |
67cb25 |
equi-distributed in 623 dimensions. It has passed the DIEHARD
|
|
Packit |
67cb25 |
statistical tests. It uses 624 words of state per generator and is
|
|
Packit |
67cb25 |
comparable in speed to the other generators. The original generator used
|
|
Packit |
67cb25 |
a default seed of 4357 and choosing :data:`s` equal to zero in
|
|
Packit |
67cb25 |
:func:`gsl_rng_set` reproduces this. Later versions switched to 5489
|
|
Packit |
67cb25 |
as the default seed, you can choose this explicitly via :func:`gsl_rng_set`
|
|
Packit |
67cb25 |
instead if you require it.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
For more information see,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* Makoto Matsumoto and Takuji Nishimura, "Mersenne Twister: A
|
|
Packit |
67cb25 |
623-dimensionally equidistributed uniform pseudorandom number
|
|
Packit |
67cb25 |
generator". ACM Transactions on Modeling and Computer
|
|
Packit |
67cb25 |
Simulation, Vol.: 8, No.: 1 (Jan. 1998), Pages 3--30
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The generator :data:`gsl_rng_mt19937` uses the second revision of the
|
|
Packit |
67cb25 |
seeding procedure published by the two authors above in 2002. The
|
|
Packit |
67cb25 |
original seeding procedures could cause spurious artifacts for some seed
|
|
Packit |
67cb25 |
values. They are still available through the alternative generators
|
|
Packit |
67cb25 |
:data:`gsl_rng_mt19937_1999` and :data:`gsl_rng_mt19937_1998`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: RANLXS random number generator
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_rng_ranlxs0
|
|
Packit |
67cb25 |
gsl_rng_ranlxs1
|
|
Packit |
67cb25 |
gsl_rng_ranlxs2
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The generator :code:`ranlxs0` is a second-generation version of the
|
|
Packit |
67cb25 |
RANLUX algorithm of Luscher, which produces "luxury random
|
|
Packit |
67cb25 |
numbers". This generator provides single precision output (24 bits) at
|
|
Packit |
67cb25 |
three luxury levels :code:`ranlxs0`, :code:`ranlxs1` and :code:`ranlxs2`,
|
|
Packit |
67cb25 |
in increasing order of strength.
|
|
Packit |
67cb25 |
It uses double-precision floating point arithmetic internally and can be
|
|
Packit |
67cb25 |
significantly faster than the integer version of :code:`ranlux`,
|
|
Packit |
67cb25 |
particularly on 64-bit architectures. The period of the generator is
|
|
Packit |
67cb25 |
about :math:`10^{171}`.
|
|
Packit |
67cb25 |
The algorithm has mathematically proven properties and
|
|
Packit |
67cb25 |
can provide truly decorrelated numbers at a known level of randomness.
|
|
Packit |
67cb25 |
The higher luxury levels provide increased decorrelation between samples
|
|
Packit |
67cb25 |
as an additional safety margin.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Note that the range of allowed seeds for this generator is :math:`[0,2^{31}-1]`.
|
|
Packit |
67cb25 |
Higher seed values are wrapped modulo :math:`2^{31}`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: RANLXD random number generator
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_rng_ranlxd1
|
|
Packit |
67cb25 |
gsl_rng_ranlxd2
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These generators produce double precision output (48 bits) from the
|
|
Packit |
67cb25 |
RANLXS generator. The library provides two luxury levels
|
|
Packit |
67cb25 |
:code:`ranlxd1` and :code:`ranlxd2`, in increasing order of strength.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: RANLUX random number generator
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_rng_ranlux
|
|
Packit |
67cb25 |
gsl_rng_ranlux389
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The :code:`ranlux` generator is an implementation of the original
|
|
Packit |
67cb25 |
algorithm developed by Luscher. It uses a
|
|
Packit |
67cb25 |
lagged-fibonacci-with-skipping algorithm to produce "luxury random
|
|
Packit |
67cb25 |
numbers". It is a 24-bit generator, originally designed for
|
|
Packit |
67cb25 |
single-precision IEEE floating point numbers. This implementation is
|
|
Packit |
67cb25 |
based on integer arithmetic, while the second-generation versions
|
|
Packit |
67cb25 |
RANLXS and RANLXD described above provide floating-point
|
|
Packit |
67cb25 |
implementations which will be faster on many platforms.
|
|
Packit |
67cb25 |
The period of the generator is about :math:`10^{171}`.
|
|
Packit |
67cb25 |
The algorithm has mathematically proven properties and
|
|
Packit |
67cb25 |
it can provide truly decorrelated numbers at a known level of
|
|
Packit |
67cb25 |
randomness. The default level of decorrelation recommended by Luscher
|
|
Packit |
67cb25 |
is provided by :data:`gsl_rng_ranlux`, while :data:`gsl_rng_ranlux389`
|
|
Packit |
67cb25 |
gives the highest level of randomness, with all 24 bits decorrelated.
|
|
Packit |
67cb25 |
Both types of generator use 24 words of state per generator.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
For more information see,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* M. Luscher, "A portable high-quality random number generator for
|
|
Packit |
67cb25 |
lattice field theory calculations", Computer Physics
|
|
Packit |
67cb25 |
Communications, 79 (1994) 100--110.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* F. James, "RANLUX: A Fortran implementation of the high-quality
|
|
Packit |
67cb25 |
pseudo-random number generator of Luscher", Computer Physics
|
|
Packit |
67cb25 |
Communications, 79 (1994) 111--114
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: CMRG, combined multiple recursive random number generator
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_rng_cmrg
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This is a combined multiple recursive generator by L'Ecuyer.
|
|
Packit |
67cb25 |
Its sequence is,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: z_n = (x_n - y_n) \mod m_1
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where the two underlying generators :math:`x_n` and :math:`y_n` are,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
x_n & = (a_1 x_{n-1} + a_2 x_{n-2} + a_3 x_{n-3}) \mod m_1 \\
|
|
Packit |
67cb25 |
y_n & = (b_1 y_{n-1} + b_2 y_{n-2} + b_3 y_{n-3}) \mod m_2
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
x_n = (a_1 x_{n-1} + a_2 x_{n-2} + a_3 x_{n-3}) mod m_1
|
|
Packit |
67cb25 |
y_n = (b_1 y_{n-1} + b_2 y_{n-2} + b_3 y_{n-3}) mod m_2
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
with coefficients
|
|
Packit |
67cb25 |
:math:`a_1 = 0`,
|
|
Packit |
67cb25 |
:math:`a_2 = 63308`,
|
|
Packit |
67cb25 |
:math:`a_3 = -183326`,
|
|
Packit |
67cb25 |
:math:`b_1 = 86098`,
|
|
Packit |
67cb25 |
:math:`b_2 = 0`,
|
|
Packit |
67cb25 |
:math:`b_3 = -539608`,
|
|
Packit |
67cb25 |
and moduli
|
|
Packit |
67cb25 |
:math:`m_1 = 2^{31} - 1 = 2147483647`
|
|
Packit |
67cb25 |
and
|
|
Packit |
67cb25 |
:math:`m_2 = 2145483479`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The period of this generator is
|
|
Packit |
67cb25 |
:math:`\hbox{lcm}(m_1^3-1, m_2^3-1)`,
|
|
Packit |
67cb25 |
which is approximately
|
|
Packit |
67cb25 |
:math:`2^{185}`
|
|
Packit |
67cb25 |
(about :math:`10^{56}`).
|
|
Packit |
67cb25 |
It uses 6 words of state per generator. For more information see,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* P. L'Ecuyer, "Combined Multiple Recursive Random Number
|
|
Packit |
67cb25 |
Generators", Operations Research, 44, 5 (1996), 816--822.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: MRG, multiple recursive random number generator
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_rng_mrg
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This is a fifth-order multiple recursive generator by L'Ecuyer, Blouin
|
|
Packit |
67cb25 |
and Coutre. Its sequence is,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: x_n = (a_1 x_{n-1} + a_5 x_{n-5}) \mod m
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
with
|
|
Packit |
67cb25 |
:math:`a_1 = 107374182`,
|
|
Packit |
67cb25 |
:math:`a_2 = a_3 = a_4 = 0`,
|
|
Packit |
67cb25 |
:math:`a_5 = 104480`
|
|
Packit |
67cb25 |
and
|
|
Packit |
67cb25 |
:math:`m = 2^{31}-1`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The period of this generator is about
|
|
Packit |
67cb25 |
:math:`10^{46}`.
|
|
Packit |
67cb25 |
It uses 5 words
|
|
Packit |
67cb25 |
of state per generator. More information can be found in the following
|
|
Packit |
67cb25 |
paper,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* P. L'Ecuyer, F. Blouin, and R. Coutre, "A search for good multiple
|
|
Packit |
67cb25 |
recursive random number generators", ACM Transactions on Modeling and
|
|
Packit |
67cb25 |
Computer Simulation 3, 87--98 (1993).
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: Tausworthe random number generator
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_rng_taus
|
|
Packit |
67cb25 |
gsl_rng_taus2
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This is a maximally equidistributed combined Tausworthe generator by
|
|
Packit |
67cb25 |
L'Ecuyer. The sequence is,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: x_n = (s^1_n \oplus s^2_n \oplus s^3_n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
x_n = (s1_n ^^ s2_n ^^ s3_n)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s^1_{n+1} &= (((s^1_n \& 4294967294)\ll 12) \oplus (((s^1_n\ll 13) \oplus s^1_n)\gg 19)) \\
|
|
Packit |
67cb25 |
s^2_{n+1} &= (((s^2_n \& 4294967288)\ll 4) \oplus (((s^2_n\ll 2) \oplus s^2_n)\gg 25)) \\
|
|
Packit |
67cb25 |
s^3_{n+1} &= (((s^3_n \& 4294967280)\ll 17) \oplus (((s^3_n\ll 3) \oplus s^3_n)\gg 11))
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
s1_{n+1} = (((s1_n&4294967294)<<12)^^(((s1_n<<13)^^s1_n)>>19))
|
|
Packit |
67cb25 |
s2_{n+1} = (((s2_n&4294967288)<< 4)^^(((s2_n<< 2)^^s2_n)>>25))
|
|
Packit |
67cb25 |
s3_{n+1} = (((s3_n&4294967280)<<17)^^(((s3_n<< 3)^^s3_n)>>11))
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
computed modulo
|
|
Packit |
67cb25 |
:math:`2^{32}`.
|
|
Packit |
67cb25 |
In the formulas above
|
|
Packit |
67cb25 |
:math:`\oplus`
|
|
Packit |
67cb25 |
denotes *exclusive-or*. Note that the algorithm relies on the properties
|
|
Packit |
67cb25 |
of 32-bit unsigned integers and has been implemented using a bitmask
|
|
Packit |
67cb25 |
of :code:`0xFFFFFFFF` to make it work on 64 bit machines.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The period of this generator is :math:`2^{88}`
|
|
Packit |
67cb25 |
(about :math:`10^{26}`).
|
|
Packit |
67cb25 |
It uses 3 words of state per generator. For more
|
|
Packit |
67cb25 |
information see,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* P. L'Ecuyer, "Maximally Equidistributed Combined Tausworthe
|
|
Packit |
67cb25 |
Generators", Mathematics of Computation, 65, 213 (1996), 203--213.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The generator :data:`gsl_rng_taus2` uses the same algorithm as
|
|
Packit |
67cb25 |
:data:`gsl_rng_taus` but with an improved seeding procedure described in
|
|
Packit |
67cb25 |
the paper,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* P. L'Ecuyer, "Tables of Maximally Equidistributed Combined LFSR
|
|
Packit |
67cb25 |
Generators", Mathematics of Computation, 68, 225 (1999), 261--269
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The generator :data:`gsl_rng_taus2` should now be used in preference to
|
|
Packit |
67cb25 |
:data:`gsl_rng_taus`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: Four-tap Generalized Feedback Shift Register
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_rng_gfsr4
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The :code:`gfsr4` generator is like a lagged-fibonacci generator, and
|
|
Packit |
67cb25 |
produces each number as an :code:`xor`'d sum of four previous values.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: r_n = r_{n-A} \oplus r_{n-B} \oplus r_{n-C} \oplus r_{n-D}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
r_n = r_{n-A} ^^ r_{n-B} ^^ r_{n-C} ^^ r_{n-D}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Ziff (ref below) notes that "it is now widely known" that two-tap
|
|
Packit |
67cb25 |
registers (such as R250, which is described below)
|
|
Packit |
67cb25 |
have serious flaws, the most obvious one being the three-point
|
|
Packit |
67cb25 |
correlation that comes from the definition of the generator. Nice
|
|
Packit |
67cb25 |
mathematical properties can be derived for GFSR's, and numerics bears
|
|
Packit |
67cb25 |
out the claim that 4-tap GFSR's with appropriately chosen offsets are as
|
|
Packit |
67cb25 |
random as can be measured, using the author's test.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This implementation uses the values suggested the example on p392 of
|
|
Packit |
67cb25 |
Ziff's article: :math:`A=471`, :math:`B=1586`, :math:`C=6988`, :math:`D=9689`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
If the offsets are appropriately chosen (such as the one ones in this
|
|
Packit |
67cb25 |
implementation), then the sequence is said to be maximal; that means
|
|
Packit |
67cb25 |
that the period is :math:`2^D - 1`, where :math:`D` is the longest lag.
|
|
Packit |
67cb25 |
(It is one less than :math:`2^D` because it is not permitted to have all
|
|
Packit |
67cb25 |
zeros in the :code:`ra[]` array.) For this implementation with
|
|
Packit |
67cb25 |
:math:`D=9689` that works out to about :math:`10^{2917}`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Note that the implementation of this generator using a 32-bit
|
|
Packit |
67cb25 |
integer amounts to 32 parallel implementations of one-bit
|
|
Packit |
67cb25 |
generators. One consequence of this is that the period of this
|
|
Packit |
67cb25 |
32-bit generator is the same as for the one-bit generator.
|
|
Packit |
67cb25 |
Moreover, this independence means that all 32-bit patterns are
|
|
Packit |
67cb25 |
equally likely, and in particular that 0 is an allowed random
|
|
Packit |
67cb25 |
value. (We are grateful to Heiko Bauke for clarifying for us these
|
|
Packit |
67cb25 |
properties of GFSR random number generators.)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
For more information see,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* Robert M. Ziff, "Four-tap shift-register-sequence random-number
|
|
Packit |
67cb25 |
generators", Computers in Physics, 12(4), Jul/Aug
|
|
Packit |
67cb25 |
1998, pp 385--392.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Unix random number generators
|
|
Packit |
67cb25 |
=============================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The standard Unix random number generators :code:`rand`, :code:`random`
|
|
Packit |
67cb25 |
and :code:`rand48` are provided as part of GSL. Although these
|
|
Packit |
67cb25 |
generators are widely available individually often they aren't all
|
|
Packit |
67cb25 |
available on the same platform. This makes it difficult to write
|
|
Packit |
67cb25 |
portable code using them and so we have included the complete set of
|
|
Packit |
67cb25 |
Unix generators in GSL for convenience. Note that these generators
|
|
Packit |
67cb25 |
don't produce high-quality randomness and aren't suitable for work
|
|
Packit |
67cb25 |
requiring accurate statistics. However, if you won't be measuring
|
|
Packit |
67cb25 |
statistical quantities and just want to introduce some variation into
|
|
Packit |
67cb25 |
your program then these generators are quite acceptable.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: rand, BSD random number generator
|
|
Packit |
67cb25 |
single: Unix random number generators, rand
|
|
Packit |
67cb25 |
single: Unix random number generators, rand48
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: BSD random number generator
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_rng_rand
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This is the BSD :code:`rand` generator. Its sequence is
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: x_{n+1} = (a x_n + c) \mod m
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
with
|
|
Packit |
67cb25 |
:math:`a = 1103515245`,
|
|
Packit |
67cb25 |
:math:`c = 12345` and
|
|
Packit |
67cb25 |
:math:`m = 2^{31}`.
|
|
Packit |
67cb25 |
The seed specifies the initial value,
|
|
Packit |
67cb25 |
:math:`x_1`. The period of this
|
|
Packit |
67cb25 |
generator is
|
|
Packit |
67cb25 |
:math:`2^{31}`,
|
|
Packit |
67cb25 |
and it uses 1 word of storage per generator.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_rng_random_bsd
|
|
Packit |
67cb25 |
gsl_rng_random_libc5
|
|
Packit |
67cb25 |
gsl_rng_random_glibc2
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These generators implement the :code:`random` family of functions, a
|
|
Packit |
67cb25 |
set of linear feedback shift register generators originally used in BSD
|
|
Packit |
67cb25 |
Unix. There are several versions of :code:`random` in use today: the
|
|
Packit |
67cb25 |
original BSD version (e.g. on SunOS4), a libc5 version (found on
|
|
Packit |
67cb25 |
older GNU/Linux systems) and a glibc2 version. Each version uses a
|
|
Packit |
67cb25 |
different seeding procedure, and thus produces different sequences.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The original BSD routines accepted a variable length buffer for the
|
|
Packit |
67cb25 |
generator state, with longer buffers providing higher-quality
|
|
Packit |
67cb25 |
randomness. The :code:`random` function implemented algorithms for
|
|
Packit |
67cb25 |
buffer lengths of 8, 32, 64, 128 and 256 bytes, and the algorithm with
|
|
Packit |
67cb25 |
the largest length that would fit into the user-supplied buffer was
|
|
Packit |
67cb25 |
used. To support these algorithms additional generators are available
|
|
Packit |
67cb25 |
with the following names::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_rng_random8_bsd
|
|
Packit |
67cb25 |
gsl_rng_random32_bsd
|
|
Packit |
67cb25 |
gsl_rng_random64_bsd
|
|
Packit |
67cb25 |
gsl_rng_random128_bsd
|
|
Packit |
67cb25 |
gsl_rng_random256_bsd
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where the numeric suffix indicates the buffer length. The original BSD
|
|
Packit |
67cb25 |
:code:`random` function used a 128-byte default buffer and so
|
|
Packit |
67cb25 |
:data:`gsl_rng_random_bsd` has been made equivalent to
|
|
Packit |
67cb25 |
:data:`gsl_rng_random128_bsd`. Corresponding versions of the :code:`libc5`
|
|
Packit |
67cb25 |
and :code:`glibc2` generators are also available, with the names
|
|
Packit |
67cb25 |
:data:`gsl_rng_random8_libc5`, :data:`gsl_rng_random8_glibc2`, etc.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: rand48 random number generator
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_rng_rand48
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This is the Unix :code:`rand48` generator. Its sequence is
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: x_{n+1} = (a x_n + c) \mod m
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
defined on 48-bit unsigned integers with
|
|
Packit |
67cb25 |
:math:`a = 25214903917`,
|
|
Packit |
67cb25 |
:math:`c = 11` and
|
|
Packit |
67cb25 |
:math:`m = 2^{48}`.
|
|
Packit |
67cb25 |
The seed specifies the upper 32 bits of the initial value, :math:`x_1`,
|
|
Packit |
67cb25 |
with the lower 16 bits set to :code:`0x330E`. The function
|
|
Packit |
67cb25 |
:func:`gsl_rng_get` returns the upper 32 bits from each term of the
|
|
Packit |
67cb25 |
sequence. This does not have a direct parallel in the original
|
|
Packit |
67cb25 |
:code:`rand48` functions, but forcing the result to type :code:`long int`
|
|
Packit |
67cb25 |
reproduces the output of :code:`mrand48`. The function
|
|
Packit |
67cb25 |
:func:`gsl_rng_uniform` uses the full 48 bits of internal state to return
|
|
Packit |
67cb25 |
the double precision number :math:`x_n/m`, which is equivalent to the
|
|
Packit |
67cb25 |
function :code:`drand48`. Note that some versions of the GNU C Library
|
|
Packit |
67cb25 |
contained a bug in :code:`mrand48` function which caused it to produce
|
|
Packit |
67cb25 |
different results (only the lower 16-bits of the return value were set).
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Other random number generators
|
|
Packit |
67cb25 |
==============================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The generators in this section are provided for compatibility with
|
|
Packit |
67cb25 |
existing libraries. If you are converting an existing program to use GSL
|
|
Packit |
67cb25 |
then you can select these generators to check your new implementation
|
|
Packit |
67cb25 |
against the original one, using the same random number generator. After
|
|
Packit |
67cb25 |
verifying that your new program reproduces the original results you can
|
|
Packit |
67cb25 |
then switch to a higher-quality generator.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Note that most of the generators in this section are based on single
|
|
Packit |
67cb25 |
linear congruence relations, which are the least sophisticated type of
|
|
Packit |
67cb25 |
generator. In particular, linear congruences have poor properties when
|
|
Packit |
67cb25 |
used with a non-prime modulus, as several of these routines do (e.g.
|
|
Packit |
67cb25 |
with a power of two modulus,
|
|
Packit |
67cb25 |
:math:`2^{31}` or
|
|
Packit |
67cb25 |
:math:`2^{32}`).
|
|
Packit |
67cb25 |
This leads to periodicity in the least significant bits of each number,
|
|
Packit |
67cb25 |
with only the higher bits having any randomness. Thus if you want to
|
|
Packit |
67cb25 |
produce a random bitstream it is best to avoid using the least
|
|
Packit |
67cb25 |
significant bits.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: RANF random number generator
|
|
Packit |
67cb25 |
single: CRAY random number generator, RANF
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_rng_ranf
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This is the CRAY random number generator :code:`RANF`. Its sequence is
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: x_{n+1} = (a x_n) \mod m
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
defined on 48-bit unsigned integers with :math:`a = 44485709377909` and
|
|
Packit |
67cb25 |
:math:`m = 2^{48}`.
|
|
Packit |
67cb25 |
The seed specifies the lower 32 bits of the initial value,
|
|
Packit |
67cb25 |
:math:`x_1`, with the lowest bit set to
|
|
Packit |
67cb25 |
prevent the seed taking an even value. The upper 16 bits of
|
|
Packit |
67cb25 |
:math:`x_1`
|
|
Packit |
67cb25 |
are set to 0. A consequence of this procedure is that the pairs of seeds
|
|
Packit |
67cb25 |
2 and 3, 4 and 5, etc.: produce the same sequences.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The generator compatible with the CRAY MATHLIB routine RANF. It
|
|
Packit |
67cb25 |
produces double precision floating point numbers which should be
|
|
Packit |
67cb25 |
identical to those from the original RANF.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
There is a subtlety in the implementation of the seeding. The initial
|
|
Packit |
67cb25 |
state is reversed through one step, by multiplying by the modular
|
|
Packit |
67cb25 |
inverse of :math:`a` mod :math:`m`. This is done for compatibility with
|
|
Packit |
67cb25 |
the original CRAY implementation.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Note that you can only seed the generator with integers up to
|
|
Packit |
67cb25 |
:math:`2^{32}`,
|
|
Packit |
67cb25 |
while the original CRAY implementation uses
|
|
Packit |
67cb25 |
non-portable wide integers which can cover all
|
|
Packit |
67cb25 |
:math:`2^{48}`
|
|
Packit |
67cb25 |
states of the generator.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The function :func:`gsl_rng_get` returns the upper 32 bits from each term
|
|
Packit |
67cb25 |
of the sequence. The function :func:`gsl_rng_uniform` uses the full 48
|
|
Packit |
67cb25 |
bits to return the double precision number :math:`x_n/m`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The period of this generator is :math:`2^{46}`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: RANMAR random number generator
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_rng_ranmar
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This is the RANMAR lagged-fibonacci generator of Marsaglia, Zaman and
|
|
Packit |
67cb25 |
Tsang. It is a 24-bit generator, originally designed for
|
|
Packit |
67cb25 |
single-precision IEEE floating point numbers. It was included in the
|
|
Packit |
67cb25 |
CERNLIB high-energy physics library.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index::
|
|
Packit |
67cb25 |
single: shift-register random number generator
|
|
Packit |
67cb25 |
single: R250 shift-register random number generator
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_rng_r250
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This is the shift-register generator of Kirkpatrick and Stoll. The
|
|
Packit |
67cb25 |
sequence is based on the recurrence
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: x_n = x_{n-103} \oplus x_{n-250}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
x_n = x_{n-103} ^^ x_{n-250}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where
|
|
Packit |
67cb25 |
:math:`\oplus`
|
|
Packit |
67cb25 |
denotes *exclusive-or*, defined on
|
|
Packit |
67cb25 |
32-bit words. The period of this generator is about :math:`2^{250}` and it
|
|
Packit |
67cb25 |
uses 250 words of state per generator.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
For more information see,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* S. Kirkpatrick and E. Stoll, "A very fast shift-register sequence random
|
|
Packit |
67cb25 |
number generator", Journal of Computational Physics, 40, 517--526
|
|
Packit |
67cb25 |
(1981)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: TT800 random number generator
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_rng_tt800
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This is an earlier version of the twisted generalized feedback
|
|
Packit |
67cb25 |
shift-register generator, and has been superseded by the development of
|
|
Packit |
67cb25 |
MT19937. However, it is still an acceptable generator in its own
|
|
Packit |
67cb25 |
right. It has a period of
|
|
Packit |
67cb25 |
:math:`2^{800}`
|
|
Packit |
67cb25 |
and uses 33 words of storage per generator.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
For more information see,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* Makoto Matsumoto and Yoshiharu Kurita, "Twisted GFSR Generators
|
|
Packit |
67cb25 |
II", ACM Transactions on Modelling and Computer Simulation,
|
|
Packit |
67cb25 |
Vol.: 4, No.: 3, 1994, pages 254--266.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. The following generators are included only for historical reasons, so
|
|
Packit |
67cb25 |
.. that you can reproduce results from old programs which might have used
|
|
Packit |
67cb25 |
.. them. These generators should not be used for real simulations since
|
|
Packit |
67cb25 |
.. they have poor statistical properties by modern standards.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: VAX random number generator
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_rng_vax
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This is the VAX generator :code:`MTH$RANDOM`. Its sequence is,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: x_{n+1} = (a x_n + c) \mod m
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
with
|
|
Packit |
67cb25 |
:math:`a = 69069`, :math:`c = 1` and
|
|
Packit |
67cb25 |
:math:`m = 2^{32}`.
|
|
Packit |
67cb25 |
The seed specifies the initial value,
|
|
Packit |
67cb25 |
:math:`x_1`. The
|
|
Packit |
67cb25 |
period of this generator is
|
|
Packit |
67cb25 |
:math:`2^{32}`
|
|
Packit |
67cb25 |
and it uses 1 word of storage per
|
|
Packit |
67cb25 |
generator.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_rng_transputer
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This is the random number generator from the INMOS Transputer
|
|
Packit |
67cb25 |
Development system. Its sequence is,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: x_{n+1} = (a x_n) \mod m
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
with :math:`a = 1664525` and
|
|
Packit |
67cb25 |
:math:`m = 2^{32}`.
|
|
Packit |
67cb25 |
The seed specifies the initial value,
|
|
Packit |
67cb25 |
:math:`x_1`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: RANDU random number generator
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_rng_randu
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This is the IBM :code:`RANDU` generator. Its sequence is
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: x_{n+1} = (a x_n) \mod m
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
with :math:`a = 65539` and
|
|
Packit |
67cb25 |
:math:`m = 2^{31}`. The
|
|
Packit |
67cb25 |
seed specifies the initial value,
|
|
Packit |
67cb25 |
:math:`x_1`. The period of this
|
|
Packit |
67cb25 |
generator was only
|
|
Packit |
67cb25 |
:math:`2^{29}`.
|
|
Packit |
67cb25 |
It has become a textbook example of a poor generator.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. index:: RANMAR random number generator
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_rng_minstd
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This is Park and Miller's "minimal standard" MINSTD generator, a
|
|
Packit |
67cb25 |
simple linear congruence which takes care to avoid the major pitfalls of
|
|
Packit |
67cb25 |
such algorithms. Its sequence is,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: x_{n+1} = (a x_n) \mod m
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
with :math:`a = 16807` and
|
|
Packit |
67cb25 |
:math:`m = 2^{31} - 1 = 2147483647`.
|
|
Packit |
67cb25 |
The seed specifies the initial value,
|
|
Packit |
67cb25 |
:math:`x_1`. The period of this
|
|
Packit |
67cb25 |
generator is about
|
|
Packit |
67cb25 |
:math:`2^{31}`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This generator was used in the IMSL Library (subroutine RNUN) and in
|
|
Packit |
67cb25 |
MATLAB (the RAND function) in the past. It is also sometimes known by
|
|
Packit |
67cb25 |
the acronym "GGL" (I'm not sure what that stands for).
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
For more information see,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* Park and Miller, "Random Number Generators: Good ones are hard to find",
|
|
Packit |
67cb25 |
Communications of the ACM, October 1988, Volume 31, No 10, pages
|
|
Packit |
67cb25 |
1192--1201.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_rng_uni
|
|
Packit |
67cb25 |
gsl_rng_uni32
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This is a reimplementation of the 16-bit SLATEC random number generator
|
|
Packit |
67cb25 |
RUNIF. A generalization of the generator to 32 bits is provided by
|
|
Packit |
67cb25 |
:data:`gsl_rng_uni32`. The original source code is available from NETLIB.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_rng_slatec
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This is the SLATEC random number generator RAND. It is ancient. The
|
|
Packit |
67cb25 |
original source code is available from NETLIB.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_rng_zuf
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This is the ZUFALL lagged Fibonacci series generator of Peterson. Its
|
|
Packit |
67cb25 |
sequence is,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: not texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
t &= u_{n-273} + u_{n-607} \\
|
|
Packit |
67cb25 |
u_n &= t - \hbox{floor}(t)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. only:: texinfo
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
t = u_{n-273} + u_{n-607}
|
|
Packit |
67cb25 |
u_n = t - floor(t)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The original source code is available from NETLIB. For more information
|
|
Packit |
67cb25 |
see,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* W. Petersen, "Lagged Fibonacci Random Number Generators for the NEC
|
|
Packit |
67cb25 |
SX-3", International Journal of High Speed Computing (1994).
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_rng_knuthran2
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This is a second-order multiple recursive generator described by Knuth
|
|
Packit |
67cb25 |
in Seminumerical Algorithms, 3rd Ed., page 108. Its sequence is,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: x_n = (a_1 x_{n-1} + a_2 x_{n-2}) \mod m
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
with
|
|
Packit |
67cb25 |
:math:`a_1 = 271828183`,
|
|
Packit |
67cb25 |
:math:`a_2 = 314159269`,
|
|
Packit |
67cb25 |
and
|
|
Packit |
67cb25 |
:math:`m = 2^{31}-1`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_rng_knuthran2002
|
|
Packit |
67cb25 |
gsl_rng_knuthran
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This is a second-order multiple recursive generator described by Knuth
|
|
Packit |
67cb25 |
in Seminumerical Algorithms, 3rd Ed., Section 3.6. Knuth
|
|
Packit |
67cb25 |
provides its C code. The updated routine :data:`gsl_rng_knuthran2002`
|
|
Packit |
67cb25 |
is from the revised 9th printing and corrects some weaknesses in the
|
|
Packit |
67cb25 |
earlier version, which is implemented as :data:`gsl_rng_knuthran`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_rng_borosh13
|
|
Packit |
67cb25 |
gsl_rng_fishman18
|
|
Packit |
67cb25 |
gsl_rng_fishman20
|
|
Packit |
67cb25 |
gsl_rng_lecuyer21
|
|
Packit |
67cb25 |
gsl_rng_waterman14
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
These multiplicative generators are taken from Knuth's
|
|
Packit |
67cb25 |
Seminumerical Algorithms, 3rd Ed., pages 106--108. Their sequence
|
|
Packit |
67cb25 |
is,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: x_{n+1} = (a x_n) \mod m
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
where the seed specifies the initial value,
|
|
Packit |
67cb25 |
:math:`x_1`.
|
|
Packit |
67cb25 |
The parameters :math:`a` and :math:`m` are as follows,
|
|
Packit |
67cb25 |
Borosh-Niederreiter:
|
|
Packit |
67cb25 |
:math:`a = 1812433253`, :math:`m = 2^{32}`,
|
|
Packit |
67cb25 |
Fishman18:
|
|
Packit |
67cb25 |
:math:`a = 62089911`,
|
|
Packit |
67cb25 |
:math:`m = 2^{31}-1`,
|
|
Packit |
67cb25 |
Fishman20:
|
|
Packit |
67cb25 |
:math:`a = 48271`,
|
|
Packit |
67cb25 |
:math:`m = 2^{31}-1`,
|
|
Packit |
67cb25 |
L'Ecuyer:
|
|
Packit |
67cb25 |
:math:`a = 40692`,
|
|
Packit |
67cb25 |
:math:`m = 2^{31}-249`,
|
|
Packit |
67cb25 |
Waterman:
|
|
Packit |
67cb25 |
:math:`a = 1566083941`,
|
|
Packit |
67cb25 |
:math:`m = 2^{32}`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_rng_fishman2x
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This is the L'Ecuyer--Fishman random number generator. It is taken from
|
|
Packit |
67cb25 |
Knuth's Seminumerical Algorithms, 3rd Ed., page 108. Its sequence
|
|
Packit |
67cb25 |
is,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: z_{n+1} = (x_n - y_n) \mod m
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
with :math:`m = 2^{31}-1`.
|
|
Packit |
67cb25 |
:math:`x_n` and :math:`y_n` are given by the :code:`fishman20`
|
|
Packit |
67cb25 |
and :code:`lecuyer21` algorithms.
|
|
Packit |
67cb25 |
The seed specifies the initial value,
|
|
Packit |
67cb25 |
:math:`x_1`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. var:: gsl_rng_coveyou
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
This is the Coveyou random number generator. It is taken from Knuth's
|
|
Packit |
67cb25 |
Seminumerical Algorithms, 3rd Ed., Section 3.2.2. Its sequence
|
|
Packit |
67cb25 |
is,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. math:: x_{n+1} = (x_n (x_n + 1)) \mod m
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
with :math:`m = 2^{32}`.
|
|
Packit |
67cb25 |
The seed specifies the initial value,
|
|
Packit |
67cb25 |
:math:`x_1`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Performance
|
|
Packit |
67cb25 |
===========
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. I made the original plot like this
|
|
Packit |
67cb25 |
.. ./benchmark > tmp; cat tmp | perl -n -e '($n,$s) = split(" ",$_); printf("%17s ",$n); print "-" x ($s/1e5), "\n";'
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. The large number of generators based on single linear congruences are
|
|
Packit |
67cb25 |
.. represented by the :code:`random` generator below. These generators are
|
|
Packit |
67cb25 |
.. fast but have the lowest statistical quality.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The following table shows the relative performance of a selection the
|
|
Packit |
67cb25 |
available random number generators. The fastest simulation quality
|
|
Packit |
67cb25 |
generators are :code:`taus`, :code:`gfsr4` and :code:`mt19937`. The
|
|
Packit |
67cb25 |
generators which offer the best mathematically-proven quality are those
|
|
Packit |
67cb25 |
based on the RANLUX algorithm::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
1754 k ints/sec, 870 k doubles/sec, taus
|
|
Packit |
67cb25 |
1613 k ints/sec, 855 k doubles/sec, gfsr4
|
|
Packit |
67cb25 |
1370 k ints/sec, 769 k doubles/sec, mt19937
|
|
Packit |
67cb25 |
565 k ints/sec, 571 k doubles/sec, ranlxs0
|
|
Packit |
67cb25 |
400 k ints/sec, 405 k doubles/sec, ranlxs1
|
|
Packit |
67cb25 |
490 k ints/sec, 389 k doubles/sec, mrg
|
|
Packit |
67cb25 |
407 k ints/sec, 297 k doubles/sec, ranlux
|
|
Packit |
67cb25 |
243 k ints/sec, 254 k doubles/sec, ranlxd1
|
|
Packit |
67cb25 |
251 k ints/sec, 253 k doubles/sec, ranlxs2
|
|
Packit |
67cb25 |
238 k ints/sec, 215 k doubles/sec, cmrg
|
|
Packit |
67cb25 |
247 k ints/sec, 198 k doubles/sec, ranlux389
|
|
Packit |
67cb25 |
141 k ints/sec, 140 k doubles/sec, ranlxd2
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Examples
|
|
Packit |
67cb25 |
========
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The following program demonstrates the use of a random number generator
|
|
Packit |
67cb25 |
to produce uniform random numbers in the range [0.0, 1.0),
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. include:: examples/rngunif.c
|
|
Packit |
67cb25 |
:code:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Here is the output of the program,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. include:: examples/rngunif.txt
|
|
Packit |
67cb25 |
:code:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The numbers depend on the seed used by the generator. The default seed
|
|
Packit |
67cb25 |
can be changed with the :macro:`GSL_RNG_SEED` environment variable to
|
|
Packit |
67cb25 |
produce a different stream of numbers. The generator itself can be
|
|
Packit |
67cb25 |
changed using the environment variable :macro:`GSL_RNG_TYPE`. Here is the
|
|
Packit |
67cb25 |
output of the program using a seed value of 123 and the
|
|
Packit |
67cb25 |
multiple-recursive generator :code:`mrg`::
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
$ GSL_RNG_SEED=123 GSL_RNG_TYPE=mrg ./a.out
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. include:: examples/rngunif2.txt
|
|
Packit |
67cb25 |
:code:
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
References and Further Reading
|
|
Packit |
67cb25 |
==============================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The subject of random number generation and testing is reviewed
|
|
Packit |
67cb25 |
extensively in Knuth's *Seminumerical Algorithms*.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* Donald E. Knuth, The Art of Computer Programming: Seminumerical
|
|
Packit |
67cb25 |
Algorithms (Vol 2, 3rd Ed, 1997), Addison-Wesley, ISBN 0201896842.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Further information is available in the review paper written by Pierre
|
|
Packit |
67cb25 |
L'Ecuyer,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* P. L'Ecuyer, "Random Number Generation", Chapter 4 of the
|
|
Packit |
67cb25 |
Handbook on Simulation, Jerry Banks Ed., Wiley, 1998, 93--137.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* http://www.iro.umontreal.ca/~lecuyer/papers.html in the file :file:`handsim.ps`.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The source code for the DIEHARD random number generator tests is also
|
|
Packit |
67cb25 |
available online,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* DIEHARD source code, G. Marsaglia, http://stat.fsu.edu/pub/diehard/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
A comprehensive set of random number generator tests is available from
|
|
Packit |
67cb25 |
NIST,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* NIST Special Publication 800-22, "A Statistical Test Suite for the
|
|
Packit |
67cb25 |
Validation of Random Number Generators and Pseudo Random Number
|
|
Packit |
67cb25 |
Generators for Cryptographic Applications".
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
* http://csrc.nist.gov/rng/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Acknowledgements
|
|
Packit |
67cb25 |
================
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Thanks to Makoto Matsumoto, Takuji Nishimura and Yoshiharu Kurita for
|
|
Packit |
67cb25 |
making the source code to their generators (MT19937, MM&TN; TT800,
|
|
Packit |
67cb25 |
MM&YK) available under the GNU General Public License. Thanks to Martin
|
|
Packit |
67cb25 |
Luscher for providing notes and source code for the RANLXS and
|
|
Packit |
67cb25 |
RANLXD generators.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. lcg
|
|
Packit |
67cb25 |
.. [ LCG(n) := n * 69069 mod (2^32) ]
|
|
Packit |
67cb25 |
.. First 6: [69069, 475559465, 2801775573, 1790562961, 3104832285, 4238970681]
|
|
Packit |
67cb25 |
.. %2^31-1 69069, 475559465, 654291926, 1790562961, 957348638, 2091487034
|
|
Packit |
67cb25 |
.. mrg
|
|
Packit |
67cb25 |
.. [q([x1, x2, x3, x4, x5]) := [107374182 mod 2147483647 * x1 + 104480 mod 2147483647 * x5, x1, x2, x3, x4]]
|
|
Packit |
67cb25 |
..
|
|
Packit |
67cb25 |
.. cmrg
|
|
Packit |
67cb25 |
.. [q1([x1,x2,x3]) := [63308 mod 2147483647 * x2 -183326 mod 2147483647 * x3, x1, x2],
|
|
Packit |
67cb25 |
.. q2([x1,x2,x3]) := [86098 mod 2145483479 * x1 -539608 mod 2145483479 * x3, x1, x2] ]
|
|
Packit |
67cb25 |
.. initial for q1 is [69069, 475559465, 654291926]
|
|
Packit |
67cb25 |
.. initial for q2 is [1790562961, 959348806, 2093487202]
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
.. tausworthe
|
|
Packit |
67cb25 |
.. [ b1(x) := rsh(xor(lsh(x, 13), x), 19),
|
|
Packit |
67cb25 |
.. q1(x) := xor(lsh(and(x, 4294967294), 12), b1(x)),
|
|
Packit |
67cb25 |
.. b2(x) := rsh(xor(lsh(x, 2), x), 25),
|
|
Packit |
67cb25 |
.. q2(x) := xor(lsh(and(x, 4294967288), 4), b2(x)),
|
|
Packit |
67cb25 |
.. b3(x) := rsh(xor(lsh(x, 3), x), 11),
|
|
Packit |
67cb25 |
.. q3(x) := xor(lsh(and(x, 4294967280), 17), b3(x)) ]
|
|
Packit |
67cb25 |
.. [s1, s2, s3] = [600098857, 1131373026, 1223067536]
|
|
Packit |
67cb25 |
.. [2948905028, 441213979, 394017882]
|