GSL Shell has moved
The project, its documentation and the release packages are now hosted on Github. Please find here the user manual.
This chapter describes routines for generating random numbers. GSL provides a large collection of random number generators that can be accessed through a uniform interface. You can choose very easily the generator you want to use and even switch between generators without needing to change your commands.
Each instance of a generator keeps track of its own state so that you can have several different generators and they will not interfere with each other. Additional functions are available for transforming uniform random numbers into samples from continuous or discrete probability distributions such as the Gaussian, log-normal or Poisson distributions.
In 1988, Park and Miller wrote a paper entitled “Random number generators: good ones are hard to find.” [park-miller]. Fortunately, some excellent random number generators are available, though poor ones are still in common use. You may be happy with the system-supplied random number generator on your computer, but you should be aware that as computers get faster, requirements on random number generators increase. Nowadays, a simulation that calls a random number generator millions of times can often finish before you can make it down the hall to the coffee machine and back.
A very nice review of random number generators was written by Pierre L’Ecuyer, as Chapter 4 of the book: Handbook on Simulation, Jerry Banks, ed. (Wiley, 1997). The chapter is available in postscript from L’Ecuyer’s ftp site (see references). Knuth’s volume on Seminumerical Algorithms (originally published in 1968) devotes 170 pages to random number generators, and has recently been updated in its 3rd edition (1997). It is brilliant, a classic. If you don’t own it, you should stop reading right now, run to the nearest bookstore, and buy it.
A good random number generator will satisfy both theoretical and statistical properties. Theoretical properties are often hard to obtain (they require real math!), but one prefers a random number generator with a long period, low serial correlation, and a tendency not to “fall mainly on the planes.” Statistical tests are performed with numerical simulations. Generally, a random number generator is used to estimate some quantity for which the theory of probability provides an exact answer. Comparison to this exact answer provides a measure of “randomness”.
In order to generate random numbers you need a “Random Number Generator” or rng for short. A random number generator is an object that can be used to generate number whenever is needed. Of course that is just what you want for a random number generator, but to achieve this effect, the generator must keep track of some kind of “state” variable. Sometimes this state is just an integer (sometimes just the value of the previously generated random number), but often it is more complicated than that and may involve a whole array of numbers, possibly with some indices thrown in. To use the random number generators, you do not need to know the details of what comprises the state, and besides that varies from algorithm to algorithm.
So, if you want to generate random numbers the first step is to use the rng.new() function,
This function returns a “random number generator” object of the specified type name. If you do not specify a particular generator the default taus2 generator will be used. See the paragraph on algorithms for a description of the different algorithms available.
This function returns a real number uniformly distributed in the range [0,1). The range includes 0.0 but excludes 1.0. Some generators compute this ratio internally so that they can provide floating point numbers with more than 32 bits of randomness.
This function returns a random integer from 0 to n-1 inclusive by scaling down and/or discarding samples from the generator R. All integers in the range [0,n-1] are produced with equal probability.
This method set the seed of the generator to the given integer value.
Return an array with all the list of all the supported generator type.
So here a simple example that prints 20 integer numbers between 0 and 999:
r = rng.new() -- we create a random number generator
for k=1,20 do
local n = r:getint(1000) -- we obtain a random integer < 1000
print(n)
end
and here another example that creates a 5x5 matrix and initialize it with random numbers like before:
r = rng.new() -- we create a random number generator
m = new(5, 5, |i,j| r:getint(1000)) -- create the matrix
The following generators are recommended for use in simulation. They have extremely long periods, low correlation and pass most statistical tests. For the most reliable source of uncorrelated numbers, the second-generation RANLUX generators have the strongest proof of randomness.
The MT19937 generator of Makoto Matsumoto and Takuji Nishimura is a variant of the twisted generalized feedback shift-register algorithm, and is known as the “Mersenne Twister” generator. It has a Mersenne prime period of 219937 - 1 (about 106000 ) and is equi-distributed in 623 dimensions. It has passed the DIEHARD statistical tests. It uses 624 words of state per generator and is comparable in speed to the other generators. The original generator used a default seed of 4357 and choosing S equal to zero in ‘set’ reproduces this. Later versions switched to 5489 as the default seed, you can choose this explicitly via ‘set’ instead if you require it.
For more information see,
The generator ‘mt19937’ uses the second revision of the seeding procedure published by the two authors above in 2002. The original seeding procedures could cause spurious artifacts for some seed values. They are still available through the alternative generators ‘mt19937_1999’ and ‘mt19937_1998’.
The ‘ranlux’ generator is an implementation of the original algorithm developed by Lüscher. It uses a lagged-fibonacci-with-skipping algorithm to produce “luxury random numbers”. It is a 24-bit generator, originally designed for single-precision IEEE floating point numbers. This implementation is based on integer arithmetic, while the second-generation versions RANLXS and RANLXD described above provide floating-point implementations which will be faster on many platforms. The period of the generator is about 10171 . The algorithm has mathematically proven properties and it can provide truly decorrelated numbers at a known level of randomness. The default level of decorrelation recommended by Lüscher is provided by ‘ranlux’, while ‘ranlux389’ gives the highest level of randomness, with all 24 bits decorrelated. Both types of generator use 24 words of state per generator.
For more information see,
This is a combined multiple recursive generator by L’Ecuyer. Its sequence is,
zn = (xn - yn ) mod m1
where the two underlying generators xn and yn are,
xn = (a1 xn-1 + a2 xn-2 + a3 xn-3 ) mod m1
yn = (b1 yn-1 + b2 yn-2 + b3 yn-3 ) mod m2
with coefficients a1 = 0, a2 = 63308, a3 = -183326, b1 = 86098, b2 = 0, b3 = -539608, and moduli m1 = 231 - 1 = 2147483647 and m2 = 2145483479.
The period of this generator is lcm(m13 -1, m23 -1), which is approximately 2185 (about 1056 ). It uses 6 words of state per generator. For more information see,
This is a fifth-order multiple recursive generator by L’Ecuyer, Blouin and Coutre. Its sequence is,
xn = (a1 xn-1 + a5 xn-5 ) mod m
with a1 = 107374182, a2 = a3 = a4 = 0, a5 = 104480 and m = 231 - 1.
The period of this generator is about 1046 . It uses 5 words of state per generator. More information can be found in the following paper,
This is a maximally equidistributed combined Tausworthe generator by L’Ecuyer. The sequence is,
xn = (s1n ^^ s2n ^^ s3n)
where,
s1n+1 = (((s1n & 4294967294)<<12)^^(((s1n <<13)^^s1n)>>19))
s2n+1 = (((s2n & 4294967288)<< 4)^^(((s2n << 2)^^s2n)>>25))
s3n+1 = (((s3n & 4294967280)<<17)^^(((s3n << 3)^^s3n)>>11))
computed modulo 232 . In the formulas above ^^ denotes “exclusive-or”. Note that the algorithm relies on the properties of 32-bit unsigned integers and has been implemented using a bitmask of ‘0xFFFFFFFF’ to make it work on 64 bit machines.
The period of this generator is 288 (about 1026 ). It uses 3 words of state per generator. For more information see,
The generator ‘taus2’ uses the same algorithm as ‘taus’ but with an improved seeding procedure described in the paper,
The generator ‘taus2’ should now be used in preference to ‘taus’.
The ‘gfsr4’ generator is like a lagged-fibonacci generator, and produces each number as an ‘xor’‘d sum of four previous values.
rn = rn-A ^^ rn-B ^^ rn-C ^^ rn-D
Ziff (ref below) notes that “it is now widely known” that two-tap registers (such as R250, which is described below) have serious flaws, the most obvious one being the three-point correlation that comes from the definition of the generator. Nice mathematical properties can be derived for GFSR’s, and numerics bears out the claim that 4-tap GFSR’s with appropriately chosen offsets are as random as can be measured, using the author’s test.
This implementation uses the values suggested the example on p392 of Ziff’s article: A=471, B=1586, C=6988, D=9689.
If the offsets are appropriately chosen (such as the one ones in this implementation), then the sequence is said to be maximal; that means that the period is 2D - 1, where D is the longest lag. (It is one less than 2D because it is not permitted to have all zeros in the ‘ra[]’ array.) For this implementation with D=9689 that works out to about 102917 .
Note that the implementation of this generator using a 32-bit integer amounts to 32 parallel implementations of one-bit generators. One consequence of this is that the period of this 32-bit generator is the same as for the one-bit generator. Moreover, this independence means that all 32-bit patterns are equally likely, and in particular that 0 is an allowed random value. (We are grateful to Heiko Bauke for clarifying for us these properties of GFSR random number generators.)
For more information see,
| [park-miller] | Commun. ACM, 31, 1192-1201 |