1 @cindex random number generators
3 The library provides a large collection of random number generators
4 which can be accessed through a uniform interface. Environment
5 variables allow you to select different generators and seeds at runtime,
6 so that you can easily switch between generators without needing to
7 recompile your program. Each instance of a generator keeps track of its
8 own state, allowing the generators to be used in multi-threaded
9 programs. Additional functions are available for transforming uniform
10 random numbers into samples from continuous or discrete probability
11 distributions such as the Gaussian, log-normal or Poisson distributions.
13 These functions are declared in the header file @file{gsl_rng.h}.
15 @comment Need to explain the difference between SERIAL and PARALLEL random
16 @comment number generators here
19 * General comments on random numbers::
20 * The Random Number Generator Interface::
21 * Random number generator initialization::
22 * Sampling from a random number generator::
23 * Auxiliary random number generator functions::
24 * Random number environment variables::
25 * Copying random number generator state::
26 * Reading and writing random number generator state::
27 * Random number generator algorithms::
28 * Unix random number generators::
29 * Other random number generators::
30 * Random Number Generator Performance::
31 * Random Number Generator Examples::
32 * Random Number References and Further Reading::
33 * Random Number Acknowledgements::
36 @node General comments on random numbers
37 @section General comments on random numbers
39 In 1988, Park and Miller wrote a paper entitled ``Random number
40 generators: good ones are hard to find.'' [Commun.@: ACM, 31, 1192--1201].
41 Fortunately, some excellent random number generators are available,
42 though poor ones are still in common use. You may be happy with the
43 system-supplied random number generator on your computer, but you should
44 be aware that as computers get faster, requirements on random number
45 generators increase. Nowadays, a simulation that calls a random number
46 generator millions of times can often finish before you can make it down
47 the hall to the coffee machine and back.
49 A very nice review of random number generators was written by Pierre
50 L'Ecuyer, as Chapter 4 of the book: Handbook on Simulation, Jerry Banks,
51 ed. (Wiley, 1997). The chapter is available in postscript from
52 L'Ecuyer's ftp site (see references). Knuth's volume on Seminumerical
53 Algorithms (originally published in 1968) devotes 170 pages to random
54 number generators, and has recently been updated in its 3rd edition
56 @comment is only now starting to show its age.
58 It is brilliant, a classic. If you don't own it, you should stop reading
59 right now, run to the nearest bookstore, and buy it.
61 A good random number generator will satisfy both theoretical and
62 statistical properties. Theoretical properties are often hard to obtain
63 (they require real math!), but one prefers a random number generator
64 with a long period, low serial correlation, and a tendency @emph{not} to
65 ``fall mainly on the planes.'' Statistical tests are performed with
66 numerical simulations. Generally, a random number generator is used to
67 estimate some quantity for which the theory of probability provides an
68 exact answer. Comparison to this exact answer provides a measure of
71 @node The Random Number Generator Interface
72 @section The Random Number Generator Interface
74 It is important to remember that a random number generator is not a
75 ``real'' function like sine or cosine. Unlike real functions, successive
76 calls to a random number generator yield different return values. Of
77 course that is just what you want for a random number generator, but to
78 achieve this effect, the generator must keep track of some kind of
79 ``state'' variable. Sometimes this state is just an integer (sometimes
80 just the value of the previously generated random number), but often it
81 is more complicated than that and may involve a whole array of numbers,
82 possibly with some indices thrown in. To use the random number
83 generators, you do not need to know the details of what comprises the
84 state, and besides that varies from algorithm to algorithm.
86 The random number generator library uses two special structs,
87 @code{gsl_rng_type} which holds static information about each type of
88 generator and @code{gsl_rng} which describes an instance of a generator
89 created from a given @code{gsl_rng_type}.
91 The functions described in this section are declared in the header file
94 @node Random number generator initialization
95 @section Random number generator initialization
97 @deftypefun {gsl_rng *} gsl_rng_alloc (const gsl_rng_type * @var{T})
98 This function returns a pointer to a newly-created
99 instance of a random number generator of type @var{T}.
100 For example, the following code creates an instance of the Tausworthe
104 gsl_rng * r = gsl_rng_alloc (gsl_rng_taus);
107 If there is insufficient memory to create the generator then the
108 function returns a null pointer and the error handler is invoked with an
109 error code of @code{GSL_ENOMEM}.
111 The generator is automatically initialized with the default seed,
112 @code{gsl_rng_default_seed}. This is zero by default but can be changed
113 either directly or by using the environment variable @code{GSL_RNG_SEED}
114 (@pxref{Random number environment variables}).
116 The details of the available generator types are
117 described later in this chapter.
120 @deftypefun void gsl_rng_set (const gsl_rng * @var{r}, unsigned long int @var{s})
121 This function initializes (or `seeds') the random number generator. If
122 the generator is seeded with the same value of @var{s} on two different
123 runs, the same stream of random numbers will be generated by successive
124 calls to the routines below. If different values of @var{s} are
125 supplied, then the generated streams of random numbers should be
126 completely different. If the seed @var{s} is zero then the standard seed
127 from the original implementation is used instead. For example, the
128 original Fortran source code for the @code{ranlux} generator used a seed
129 of 314159265, and so choosing @var{s} equal to zero reproduces this when
130 using @code{gsl_rng_ranlux}.
133 @deftypefun void gsl_rng_free (gsl_rng * @var{r})
134 This function frees all the memory associated with the generator
138 @node Sampling from a random number generator
139 @section Sampling from a random number generator
141 The following functions return uniformly distributed random numbers,
142 either as integers or double precision floating point numbers. @inlinefns{}
143 To obtain non-uniform distributions @pxref{Random Number Distributions}.
145 @deftypefun {unsigned long int} gsl_rng_get (const gsl_rng * @var{r})
146 This function returns a random integer from the generator @var{r}. The
147 minimum and maximum values depend on the algorithm used, but all
148 integers in the range [@var{min},@var{max}] are equally likely. The
149 values of @var{min} and @var{max} can determined using the auxiliary
150 functions @code{gsl_rng_max (r)} and @code{gsl_rng_min (r)}.
153 @deftypefun double gsl_rng_uniform (const gsl_rng * @var{r})
154 This function returns a double precision floating point number uniformly
155 distributed in the range [0,1). The range includes 0.0 but excludes 1.0.
156 The value is typically obtained by dividing the result of
157 @code{gsl_rng_get(r)} by @code{gsl_rng_max(r) + 1.0} in double
158 precision. Some generators compute this ratio internally so that they
159 can provide floating point numbers with more than 32 bits of randomness
160 (the maximum number of bits that can be portably represented in a single
161 @code{unsigned long int}).
164 @deftypefun double gsl_rng_uniform_pos (const gsl_rng * @var{r})
165 This function returns a positive double precision floating point number
166 uniformly distributed in the range (0,1), excluding both 0.0 and 1.0.
167 The number is obtained by sampling the generator with the algorithm of
168 @code{gsl_rng_uniform} until a non-zero value is obtained. You can use
169 this function if you need to avoid a singularity at 0.0.
172 @deftypefun {unsigned long int} gsl_rng_uniform_int (const gsl_rng * @var{r}, unsigned long int @var{n})
173 This function returns a random integer from 0 to @math{n-1} inclusive
174 by scaling down and/or discarding samples from the generator @var{r}.
175 All integers in the range @math{[0,n-1]} are produced with equal
176 probability. For generators with a non-zero minimum value an offset
177 is applied so that zero is returned with the correct probability.
179 Note that this function is designed for sampling from ranges smaller
180 than the range of the underlying generator. The parameter @var{n}
181 must be less than or equal to the range of the generator @var{r}.
182 If @var{n} is larger than the range of the generator then the function
183 calls the error handler with an error code of @code{GSL_EINVAL} and
186 In particular, this function is not intended for generating the full range of
187 unsigned integer values @c{$[0,2^{32}-1]$}
188 @math{[0,2^32-1]}. Instead
189 choose a generator with the maximal integer range and zero mimimum
190 value, such as @code{gsl_rng_ranlxd1}, @code{gsl_rng_mt19937} or
191 @code{gsl_rng_taus}, and sample it directly using
192 @code{gsl_rng_get}. The range of each generator can be found using
193 the auxiliary functions described in the next section.
196 @node Auxiliary random number generator functions
197 @section Auxiliary random number generator functions
198 The following functions provide information about an existing
199 generator. You should use them in preference to hard-coding the generator
200 parameters into your own code.
202 @deftypefun {const char *} gsl_rng_name (const gsl_rng * @var{r})
203 This function returns a pointer to the name of the generator.
207 printf ("r is a '%s' generator\n",
212 would print something like @code{r is a 'taus' generator}.
215 @deftypefun {unsigned long int} gsl_rng_max (const gsl_rng * @var{r})
216 @code{gsl_rng_max} returns the largest value that @code{gsl_rng_get}
220 @deftypefun {unsigned long int} gsl_rng_min (const gsl_rng * @var{r})
221 @code{gsl_rng_min} returns the smallest value that @code{gsl_rng_get}
222 can return. Usually this value is zero. There are some generators with
223 algorithms that cannot return zero, and for these generators the minimum
227 @deftypefun {void *} gsl_rng_state (const gsl_rng * @var{r})
228 @deftypefunx size_t gsl_rng_size (const gsl_rng * @var{r})
229 These functions return a pointer to the state of generator @var{r} and
230 its size. You can use this information to access the state directly. For
231 example, the following code will write the state of a generator to a
235 void * state = gsl_rng_state (r);
236 size_t n = gsl_rng_size (r);
237 fwrite (state, n, 1, stream);
241 @deftypefun {const gsl_rng_type **} gsl_rng_types_setup (void)
242 This function returns a pointer to an array of all the available
243 generator types, terminated by a null pointer. The function should be
244 called once at the start of the program, if needed. The following code
245 fragment shows how to iterate over the array of generator types to print
246 the names of the available algorithms,
249 const gsl_rng_type **t, **t0;
251 t0 = gsl_rng_types_setup ();
253 printf ("Available generators:\n");
255 for (t = t0; *t != 0; t++)
257 printf ("%s\n", (*t)->name);
262 @node Random number environment variables
263 @section Random number environment variables
265 The library allows you to choose a default generator and seed from the
266 environment variables @code{GSL_RNG_TYPE} and @code{GSL_RNG_SEED} and
267 the function @code{gsl_rng_env_setup}. This makes it easy try out
268 different generators and seeds without having to recompile your program.
270 @deftypefun {const gsl_rng_type *} gsl_rng_env_setup (void)
271 This function reads the environment variables @code{GSL_RNG_TYPE} and
272 @code{GSL_RNG_SEED} and uses their values to set the corresponding
273 library variables @code{gsl_rng_default} and
274 @code{gsl_rng_default_seed}. These global variables are defined as
278 extern const gsl_rng_type *gsl_rng_default
279 extern unsigned long int gsl_rng_default_seed
282 The environment variable @code{GSL_RNG_TYPE} should be the name of a
283 generator, such as @code{taus} or @code{mt19937}. The environment
284 variable @code{GSL_RNG_SEED} should contain the desired seed value. It
285 is converted to an @code{unsigned long int} using the C library function
288 If you don't specify a generator for @code{GSL_RNG_TYPE} then
289 @code{gsl_rng_mt19937} is used as the default. The initial value of
290 @code{gsl_rng_default_seed} is zero.
296 Here is a short program which shows how to create a global
297 generator using the environment variables @code{GSL_RNG_TYPE} and
301 @verbatiminclude examples/rng.c
305 Running the program without any environment variables uses the initial
306 defaults, an @code{mt19937} generator with a seed of 0,
310 @verbatiminclude examples/rng.out
314 By setting the two variables on the command line we can
315 change the default generator and the seed,
318 $ GSL_RNG_TYPE="taus" GSL_RNG_SEED=123 ./a.out
323 first value = 2720986350
326 @node Copying random number generator state
327 @section Copying random number generator state
329 The above methods do not expose the random number `state' which changes
330 from call to call. It is often useful to be able to save and restore
331 the state. To permit these practices, a few somewhat more advanced
332 functions are supplied. These include:
334 @deftypefun int gsl_rng_memcpy (gsl_rng * @var{dest}, const gsl_rng * @var{src})
335 This function copies the random number generator @var{src} into the
336 pre-existing generator @var{dest}, making @var{dest} into an exact copy
337 of @var{src}. The two generators must be of the same type.
340 @deftypefun {gsl_rng *} gsl_rng_clone (const gsl_rng * @var{r})
341 This function returns a pointer to a newly created generator which is an
342 exact copy of the generator @var{r}.
345 @node Reading and writing random number generator state
346 @section Reading and writing random number generator state
348 The library provides functions for reading and writing the random
349 number state to a file as binary data or formatted text.
351 @deftypefun int gsl_rng_fwrite (FILE * @var{stream}, const gsl_rng * @var{r})
352 This function writes the random number state of the random number
353 generator @var{r} to the stream @var{stream} in binary format. The
354 return value is 0 for success and @code{GSL_EFAILED} if there was a
355 problem writing to the file. Since the data is written in the native
356 binary format it may not be portable between different architectures.
359 @deftypefun int gsl_rng_fread (FILE * @var{stream}, gsl_rng * @var{r})
360 This function reads the random number state into the random number
361 generator @var{r} from the open stream @var{stream} in binary format.
362 The random number generator @var{r} must be preinitialized with the
363 correct random number generator type since type information is not
364 saved. The return value is 0 for success and @code{GSL_EFAILED} if
365 there was a problem reading from the file. The data is assumed to
366 have been written in the native binary format on the same
370 @node Random number generator algorithms
371 @section Random number generator algorithms
373 The functions described above make no reference to the actual algorithm
374 used. This is deliberate so that you can switch algorithms without
375 having to change any of your application source code. The library
376 provides a large number of generators of different types, including
377 simulation quality generators, generators provided for compatibility
378 with other libraries and historical generators from the past.
380 The following generators are recommended for use in simulation. They
381 have extremely long periods, low correlation and pass most statistical
382 tests. For the most reliable source of uncorrelated numbers, the
383 second-generation @sc{ranlux} generators have the strongest proof of
386 @deffn {Generator} gsl_rng_mt19937
387 @cindex MT19937 random number generator
388 The MT19937 generator of Makoto Matsumoto and Takuji Nishimura is a
389 variant of the twisted generalized feedback shift-register algorithm,
390 and is known as the ``Mersenne Twister'' generator. It has a Mersenne
394 @math{2^19937 - 1} (about
396 @math{10^6000}) and is
397 equi-distributed in 623 dimensions. It has passed the @sc{diehard}
398 statistical tests. It uses 624 words of state per generator and is
399 comparable in speed to the other generators. The original generator used
400 a default seed of 4357 and choosing @var{s} equal to zero in
401 @code{gsl_rng_set} reproduces this. Later versions switched to 5489
402 as the default seed, you can choose this explicitly via @code{gsl_rng_set}
403 instead if you require it.
405 For more information see,
408 Makoto Matsumoto and Takuji Nishimura, ``Mersenne Twister: A
409 623-dimensionally equidistributed uniform pseudorandom number
410 generator''. @cite{ACM Transactions on Modeling and Computer
411 Simulation}, Vol.@: 8, No.@: 1 (Jan. 1998), Pages 3--30
415 The generator @code{gsl_rng_mt19937} uses the second revision of the
416 seeding procedure published by the two authors above in 2002. The
417 original seeding procedures could cause spurious artifacts for some seed
418 values. They are still available through the alternative generators
419 @code{gsl_rng_mt19937_1999} and @code{gsl_rng_mt19937_1998}.
422 @deffn {Generator} gsl_rng_ranlxs0
423 @deffnx {Generator} gsl_rng_ranlxs1
424 @deffnx {Generator} gsl_rng_ranlxs2
425 @cindex RANLXS random number generator
427 The generator @code{ranlxs0} is a second-generation version of the
428 @sc{ranlux} algorithm of L@"uscher, which produces ``luxury random
429 numbers''. This generator provides single precision output (24 bits) at
430 three luxury levels @code{ranlxs0}, @code{ranlxs1} and @code{ranlxs2},
431 in increasing order of strength.
432 It uses double-precision floating point arithmetic internally and can be
433 significantly faster than the integer version of @code{ranlux},
434 particularly on 64-bit architectures. The period of the generator is
436 @math{10^171}. The algorithm has mathematically proven properties and
437 can provide truly decorrelated numbers at a known level of randomness.
438 The higher luxury levels provide increased decorrelation between samples
439 as an additional safety margin.
442 @deffn {Generator} gsl_rng_ranlxd1
443 @deffnx {Generator} gsl_rng_ranlxd2
444 @cindex RANLXD random number generator
446 These generators produce double precision output (48 bits) from the
447 @sc{ranlxs} generator. The library provides two luxury levels
448 @code{ranlxd1} and @code{ranlxd2}, in increasing order of strength.
452 @deffn {Generator} gsl_rng_ranlux
453 @deffnx {Generator} gsl_rng_ranlux389
455 @cindex RANLUX random number generator
456 The @code{ranlux} generator is an implementation of the original
457 algorithm developed by L@"uscher. It uses a
458 lagged-fibonacci-with-skipping algorithm to produce ``luxury random
459 numbers''. It is a 24-bit generator, originally designed for
460 single-precision IEEE floating point numbers. This implementation is
461 based on integer arithmetic, while the second-generation versions
462 @sc{ranlxs} and @sc{ranlxd} described above provide floating-point
463 implementations which will be faster on many platforms.
464 The period of the generator is about @c{$10^{171}$}
465 @math{10^171}. The algorithm has mathematically proven properties and
466 it can provide truly decorrelated numbers at a known level of
467 randomness. The default level of decorrelation recommended by L@"uscher
468 is provided by @code{gsl_rng_ranlux}, while @code{gsl_rng_ranlux389}
469 gives the highest level of randomness, with all 24 bits decorrelated.
470 Both types of generator use 24 words of state per generator.
472 For more information see,
475 M. L@"uscher, ``A portable high-quality random number generator for
476 lattice field theory calculations'', @cite{Computer Physics
477 Communications}, 79 (1994) 100--110.
479 F. James, ``RANLUX: A Fortran implementation of the high-quality
480 pseudo-random number generator of L@"uscher'', @cite{Computer Physics
481 Communications}, 79 (1994) 111--114
486 @deffn {Generator} gsl_rng_cmrg
487 @cindex CMRG, combined multiple recursive random number generator
488 This is a combined multiple recursive generator by L'Ecuyer.
493 z_n = (x_n - y_n) \,\hbox{mod}\, m_1
500 z_n = (x_n - y_n) mod m_1
505 where the two underlying generators @math{x_n} and @math{y_n} are,
510 x_n & = (a_1 x_{n-1} + a_2 x_{n-2} + a_3 x_{n-3}) \,\hbox{mod}\, m_1 \cr
511 y_n & = (b_1 y_{n-1} + b_2 y_{n-2} + b_3 y_{n-3}) \,\hbox{mod}\, m_2
519 x_n = (a_1 x_@{n-1@} + a_2 x_@{n-2@} + a_3 x_@{n-3@}) mod m_1
520 y_n = (b_1 y_@{n-1@} + b_2 y_@{n-2@} + b_3 y_@{n-3@}) mod m_2
528 @math{a_3 = -183326},
531 @math{b_3 = -539608},
533 @c{$m_1 = 2^{31} - 1 = 2147483647$}
534 @math{m_1 = 2^31 - 1 = 2147483647}
536 @c{$m_2 = 2145483479$}
537 @math{m_2 = 2145483479}.
539 The period of this generator is
540 @c{$\hbox{lcm}(m_1^3-1, m_2^3-1)$}
541 @math{lcm(m_1^3-1, m_2^3-1)},
542 which is approximately
547 @math{10^56}). It uses
548 6 words of state per generator. For more information see,
552 P. L'Ecuyer, ``Combined Multiple Recursive Random Number
553 Generators'', @cite{Operations Research}, 44, 5 (1996), 816--822.
557 @deffn {Generator} gsl_rng_mrg
558 @cindex MRG, multiple recursive random number generator
559 This is a fifth-order multiple recursive generator by L'Ecuyer, Blouin
560 and Coutre. Its sequence is,
564 x_n = (a_1 x_{n-1} + a_5 x_{n-5}) \,\hbox{mod}\, m
571 x_n = (a_1 x_@{n-1@} + a_5 x_@{n-5@}) mod m
577 @math{a_1 = 107374182},
578 @math{a_2 = a_3 = a_4 = 0},
584 The period of this generator is about
586 @math{10^46}. It uses 5 words
587 of state per generator. More information can be found in the following
591 P. L'Ecuyer, F. Blouin, and R. Coutre, ``A search for good multiple
592 recursive random number generators'', @cite{ACM Transactions on Modeling and
593 Computer Simulation} 3, 87--98 (1993).
597 @deffn {Generator} gsl_rng_taus
598 @deffnx {Generator} gsl_rng_taus2
599 @cindex Tausworthe random number generator
600 This is a maximally equidistributed combined Tausworthe generator by
601 L'Ecuyer. The sequence is,
605 x_n = (s^1_n \oplus s^2_n \oplus s^3_n)
612 x_n = (s1_n ^^ s2_n ^^ s3_n)
622 s^1_{n+1} &= (((s^1_n \& 4294967294)\ll 12) \oplus (((s^1_n\ll 13) \oplus s^1_n)\gg 19)) \cr
623 s^2_{n+1} &= (((s^2_n \& 4294967288)\ll 4) \oplus (((s^2_n\ll 2) \oplus s^2_n)\gg 25)) \cr
624 s^3_{n+1} &= (((s^3_n \& 4294967280)\ll 17) \oplus (((s^3_n\ll 3) \oplus s^3_n)\gg 11))
632 s1_@{n+1@} = (((s1_n&4294967294)<<12)^^(((s1_n<<13)^^s1_n)>>19))
633 s2_@{n+1@} = (((s2_n&4294967288)<< 4)^^(((s2_n<< 2)^^s2_n)>>25))
634 s3_@{n+1@} = (((s3_n&4294967280)<<17)^^(((s3_n<< 3)^^s3_n)>>11))
641 @math{2^32}. In the formulas above
644 denotes ``exclusive-or''. Note that the algorithm relies on the properties
645 of 32-bit unsigned integers and has been implemented using a bitmask
646 of @code{0xFFFFFFFF} to make it work on 64 bit machines.
648 The period of this generator is @c{$2^{88}$}
651 @math{10^26}). It uses 3 words of state per generator. For more
656 P. L'Ecuyer, ``Maximally Equidistributed Combined Tausworthe
657 Generators'', @cite{Mathematics of Computation}, 65, 213 (1996), 203--213.
661 The generator @code{gsl_rng_taus2} uses the same algorithm as
662 @code{gsl_rng_taus} but with an improved seeding procedure described in
667 P. L'Ecuyer, ``Tables of Maximally Equidistributed Combined LFSR
668 Generators'', @cite{Mathematics of Computation}, 68, 225 (1999), 261--269
672 The generator @code{gsl_rng_taus2} should now be used in preference to
676 @deffn {Generator} gsl_rng_gfsr4
677 @cindex Four-tap Generalized Feedback Shift Register
678 The @code{gfsr4} generator is like a lagged-fibonacci generator, and
679 produces each number as an @code{xor}'d sum of four previous values.
683 r_n = r_{n-A} \oplus r_{n-B} \oplus r_{n-C} \oplus r_{n-D}
690 r_n = r_@{n-A@} ^^ r_@{n-B@} ^^ r_@{n-C@} ^^ r_@{n-D@}
694 Ziff (ref below) notes that ``it is now widely known'' that two-tap
695 registers (such as R250, which is described below)
696 have serious flaws, the most obvious one being the three-point
697 correlation that comes from the definition of the generator. Nice
698 mathematical properties can be derived for GFSR's, and numerics bears
699 out the claim that 4-tap GFSR's with appropriately chosen offsets are as
700 random as can be measured, using the author's test.
702 This implementation uses the values suggested the example on p392 of
703 Ziff's article: @math{A=471}, @math{B=1586}, @math{C=6988}, @math{D=9689}.
706 If the offsets are appropriately chosen (such as the one ones in this
707 implementation), then the sequence is said to be maximal; that means
708 that the period is @math{2^D - 1}, where @math{D} is the longest lag.
709 (It is one less than @math{2^D} because it is not permitted to have all
710 zeros in the @code{ra[]} array.) For this implementation with
711 @math{D=9689} that works out to about @c{$10^{2917}$}
714 Note that the implementation of this generator using a 32-bit
715 integer amounts to 32 parallel implementations of one-bit
716 generators. One consequence of this is that the period of this
717 32-bit generator is the same as for the one-bit generator.
718 Moreover, this independence means that all 32-bit patterns are
719 equally likely, and in particular that 0 is an allowed random
720 value. (We are grateful to Heiko Bauke for clarifying for us these
721 properties of GFSR random number generators.)
723 For more information see,
726 Robert M. Ziff, ``Four-tap shift-register-sequence random-number
727 generators'', @cite{Computers in Physics}, 12(4), Jul/Aug
732 @node Unix random number generators
733 @section Unix random number generators
735 The standard Unix random number generators @code{rand}, @code{random}
736 and @code{rand48} are provided as part of GSL. Although these
737 generators are widely available individually often they aren't all
738 available on the same platform. This makes it difficult to write
739 portable code using them and so we have included the complete set of
740 Unix generators in GSL for convenience. Note that these generators
741 don't produce high-quality randomness and aren't suitable for work
742 requiring accurate statistics. However, if you won't be measuring
743 statistical quantities and just want to introduce some variation into
744 your program then these generators are quite acceptable.
746 @cindex rand, BSD random number generator
747 @cindex Unix random number generators, rand
748 @cindex Unix random number generators, rand48
750 @deffn {Generator} gsl_rng_rand
751 @cindex BSD random number generator
752 This is the BSD @code{rand} generator. Its sequence is
756 x_{n+1} = (a x_n + c) \,\hbox{mod}\, m
763 x_@{n+1@} = (a x_n + c) mod m
769 @math{a = 1103515245},
773 The seed specifies the initial value,
774 @math{x_1}. The period of this
777 @math{2^31}, and it uses 1 word of storage per
781 @deffn {Generator} gsl_rng_random_bsd
782 @deffnx {Generator} gsl_rng_random_libc5
783 @deffnx {Generator} gsl_rng_random_glibc2
784 These generators implement the @code{random} family of functions, a
785 set of linear feedback shift register generators originally used in BSD
786 Unix. There are several versions of @code{random} in use today: the
787 original BSD version (e.g. on SunOS4), a libc5 version (found on
788 older GNU/Linux systems) and a glibc2 version. Each version uses a
789 different seeding procedure, and thus produces different sequences.
791 The original BSD routines accepted a variable length buffer for the
792 generator state, with longer buffers providing higher-quality
793 randomness. The @code{random} function implemented algorithms for
794 buffer lengths of 8, 32, 64, 128 and 256 bytes, and the algorithm with
795 the largest length that would fit into the user-supplied buffer was
796 used. To support these algorithms additional generators are available
797 with the following names,
803 gsl_rng_random128_bsd
804 gsl_rng_random256_bsd
808 where the numeric suffix indicates the buffer length. The original BSD
809 @code{random} function used a 128-byte default buffer and so
810 @code{gsl_rng_random_bsd} has been made equivalent to
811 @code{gsl_rng_random128_bsd}. Corresponding versions of the @code{libc5}
812 and @code{glibc2} generators are also available, with the names
813 @code{gsl_rng_random8_libc5}, @code{gsl_rng_random8_glibc2}, etc.
816 @deffn {Generator} gsl_rng_rand48
817 @cindex rand48 random number generator
818 This is the Unix @code{rand48} generator. Its sequence is
822 x_{n+1} = (a x_n + c) \,\hbox{mod}\, m
829 x_@{n+1@} = (a x_n + c) mod m
834 defined on 48-bit unsigned integers with
835 @math{a = 25214903917},
839 The seed specifies the upper 32 bits of the initial value, @math{x_1},
840 with the lower 16 bits set to @code{0x330E}. The function
841 @code{gsl_rng_get} returns the upper 32 bits from each term of the
842 sequence. This does not have a direct parallel in the original
843 @code{rand48} functions, but forcing the result to type @code{long int}
844 reproduces the output of @code{mrand48}. The function
845 @code{gsl_rng_uniform} uses the full 48 bits of internal state to return
846 the double precision number @math{x_n/m}, which is equivalent to the
847 function @code{drand48}. Note that some versions of the GNU C Library
848 contained a bug in @code{mrand48} function which caused it to produce
849 different results (only the lower 16-bits of the return value were set).
852 @node Other random number generators
853 @section Other random number generators
855 The generators in this section are provided for compatibility with
856 existing libraries. If you are converting an existing program to use GSL
857 then you can select these generators to check your new implementation
858 against the original one, using the same random number generator. After
859 verifying that your new program reproduces the original results you can
860 then switch to a higher-quality generator.
862 Note that most of the generators in this section are based on single
863 linear congruence relations, which are the least sophisticated type of
864 generator. In particular, linear congruences have poor properties when
865 used with a non-prime modulus, as several of these routines do (e.g.
866 with a power of two modulus,
871 leads to periodicity in the least significant bits of each number,
872 with only the higher bits having any randomness. Thus if you want to
873 produce a random bitstream it is best to avoid using the least
876 @deffn {Generator} gsl_rng_ranf
877 @cindex RANF random number generator
878 @cindex CRAY random number generator, RANF
879 This is the CRAY random number generator @code{RANF}. Its sequence is
883 x_{n+1} = (a x_n) \,\hbox{mod}\, m
890 x_@{n+1@} = (a x_n) mod m
895 defined on 48-bit unsigned integers with @math{a = 44485709377909} and
897 @math{m = 2^48}. The seed specifies the lower
898 32 bits of the initial value,
899 @math{x_1}, with the lowest bit set to
900 prevent the seed taking an even value. The upper 16 bits of
902 are set to 0. A consequence of this procedure is that the pairs of seeds
903 2 and 3, 4 and 5, etc produce the same sequences.
905 The generator compatible with the CRAY MATHLIB routine RANF. It
906 produces double precision floating point numbers which should be
907 identical to those from the original RANF.
909 There is a subtlety in the implementation of the seeding. The initial
910 state is reversed through one step, by multiplying by the modular
911 inverse of @math{a} mod @math{m}. This is done for compatibility with
912 the original CRAY implementation.
914 Note that you can only seed the generator with integers up to
916 @math{2^32}, while the original CRAY implementation uses
917 non-portable wide integers which can cover all
919 @math{2^48} states of the generator.
921 The function @code{gsl_rng_get} returns the upper 32 bits from each term
922 of the sequence. The function @code{gsl_rng_uniform} uses the full 48
923 bits to return the double precision number @math{x_n/m}.
925 The period of this generator is @c{$2^{46}$}
929 @deffn {Generator} gsl_rng_ranmar
930 @cindex RANMAR random number generator
931 This is the RANMAR lagged-fibonacci generator of Marsaglia, Zaman and
932 Tsang. It is a 24-bit generator, originally designed for
933 single-precision IEEE floating point numbers. It was included in the
934 CERNLIB high-energy physics library.
937 @deffn {Generator} gsl_rng_r250
938 @cindex shift-register random number generator
939 @cindex R250 shift-register random number generator
940 This is the shift-register generator of Kirkpatrick and Stoll. The
941 sequence is based on the recurrence
945 x_n = x_{n-103} \oplus x_{n-250}
952 x_n = x_@{n-103@} ^^ x_@{n-250@}
959 @math{^^} denotes ``exclusive-or'', defined on
960 32-bit words. The period of this generator is about @c{$2^{250}$}
962 uses 250 words of state per generator.
964 For more information see,
967 S. Kirkpatrick and E. Stoll, ``A very fast shift-register sequence random
968 number generator'', @cite{Journal of Computational Physics}, 40, 517--526
973 @deffn {Generator} gsl_rng_tt800
974 @cindex TT800 random number generator
975 This is an earlier version of the twisted generalized feedback
976 shift-register generator, and has been superseded by the development of
977 MT19937. However, it is still an acceptable generator in its own
978 right. It has a period of
980 @math{2^800} and uses 33 words of storage
983 For more information see,
986 Makoto Matsumoto and Yoshiharu Kurita, ``Twisted GFSR Generators
987 II'', @cite{ACM Transactions on Modelling and Computer Simulation},
988 Vol.@: 4, No.@: 3, 1994, pages 254--266.
992 @comment The following generators are included only for historical reasons, so
993 @comment that you can reproduce results from old programs which might have used
994 @comment them. These generators should not be used for real simulations since
995 @comment they have poor statistical properties by modern standards.
997 @deffn {Generator} gsl_rng_vax
998 @cindex VAX random number generator
999 This is the VAX generator @code{MTH$RANDOM}. Its sequence is,
1003 x_{n+1} = (a x_n + c) \,\hbox{mod}\, m
1010 x_@{n+1@} = (a x_n + c) mod m
1016 @math{a = 69069}, @math{c = 1} and
1018 @math{m = 2^32}. The seed specifies the initial value,
1020 period of this generator is
1022 @math{2^32} and it uses 1 word of storage per
1026 @deffn {Generator} gsl_rng_transputer
1027 This is the random number generator from the INMOS Transputer
1028 Development system. Its sequence is,
1032 x_{n+1} = (a x_n) \,\hbox{mod}\, m
1039 x_@{n+1@} = (a x_n) mod m
1044 with @math{a = 1664525} and
1047 The seed specifies the initial value,
1052 @deffn {Generator} gsl_rng_randu
1053 @cindex RANDU random number generator
1054 This is the IBM @code{RANDU} generator. Its sequence is
1058 x_{n+1} = (a x_n) \,\hbox{mod}\, m
1065 x_@{n+1@} = (a x_n) mod m
1070 with @math{a = 65539} and
1072 @math{m = 2^31}. The
1073 seed specifies the initial value,
1074 @math{x_1}. The period of this
1077 @math{2^29}. It has become a textbook example of a
1081 @deffn {Generator} gsl_rng_minstd
1082 @cindex RANMAR random number generator
1083 This is Park and Miller's ``minimal standard'' @sc{minstd} generator, a
1084 simple linear congruence which takes care to avoid the major pitfalls of
1085 such algorithms. Its sequence is,
1089 x_{n+1} = (a x_n) \,\hbox{mod}\, m
1096 x_@{n+1@} = (a x_n) mod m
1101 with @math{a = 16807} and
1102 @c{$m = 2^{31} - 1 = 2147483647$}
1103 @math{m = 2^31 - 1 = 2147483647}.
1104 The seed specifies the initial value,
1106 @math{x_1}. The period of this
1111 This generator is used in the IMSL Library (subroutine RNUN) and in
1112 MATLAB (the RAND function). It is also sometimes known by the acronym
1113 ``GGL'' (I'm not sure what that stands for).
1115 For more information see,
1118 Park and Miller, ``Random Number Generators: Good ones are hard to find'',
1119 @cite{Communications of the ACM}, October 1988, Volume 31, No 10, pages
1124 @deffn {Generator} gsl_rng_uni
1125 @deffnx {Generator} gsl_rng_uni32
1126 This is a reimplementation of the 16-bit SLATEC random number generator
1127 RUNIF. A generalization of the generator to 32 bits is provided by
1128 @code{gsl_rng_uni32}. The original source code is available from NETLIB.
1131 @deffn {Generator} gsl_rng_slatec
1132 This is the SLATEC random number generator RAND. It is ancient. The
1133 original source code is available from NETLIB.
1137 @deffn {Generator} gsl_rng_zuf
1138 This is the ZUFALL lagged Fibonacci series generator of Peterson. Its
1144 t &= u_{n-273} + u_{n-607} \cr
1145 u_n &= t - \hbox{floor}(t)
1153 t = u_@{n-273@} + u_@{n-607@}
1158 The original source code is available from NETLIB. For more information
1162 W. Petersen, ``Lagged Fibonacci Random Number Generators for the NEC
1163 SX-3'', @cite{International Journal of High Speed Computing} (1994).
1167 @deffn {Generator} gsl_rng_knuthran2
1168 This is a second-order multiple recursive generator described by Knuth
1169 in @cite{Seminumerical Algorithms}, 3rd Ed., page 108. Its sequence is,
1173 x_n = (a_1 x_{n-1} + a_2 x_{n-2}) \,\hbox{mod}\, m
1180 x_n = (a_1 x_@{n-1@} + a_2 x_@{n-2@}) mod m
1186 @math{a_1 = 271828183},
1187 @math{a_2 = 314159269},
1190 @math{m = 2^31 - 1}.
1193 @deffn {Generator} gsl_rng_knuthran2002
1194 @deffnx {Generator} gsl_rng_knuthran
1195 This is a second-order multiple recursive generator described by Knuth
1196 in @cite{Seminumerical Algorithms}, 3rd Ed., Section 3.6. Knuth
1197 provides its C code. The updated routine @code{gsl_rng_knuthran2002}
1198 is from the revised 9th printing and corrects some weaknesses in the
1199 earlier version, which is implemented as @code{gsl_rng_knuthran}.
1202 @deffn {Generator} gsl_rng_borosh13
1203 @deffnx {Generator} gsl_rng_fishman18
1204 @deffnx {Generator} gsl_rng_fishman20
1205 @deffnx {Generator} gsl_rng_lecuyer21
1206 @deffnx {Generator} gsl_rng_waterman14
1207 These multiplicative generators are taken from Knuth's
1208 @cite{Seminumerical Algorithms}, 3rd Ed., pages 106--108. Their sequence
1213 x_{n+1} = (a x_n) \,\hbox{mod}\, m
1220 x_@{n+1@} = (a x_n) mod m
1225 where the seed specifies the initial value, @c{$x_1$}
1227 The parameters @math{a} and @math{m} are as follows,
1228 Borosh-Niederreiter:
1229 @math{a = 1812433253}, @c{$m = 2^{32}$}
1232 @math{a = 62089911},
1234 @math{m = 2^31 - 1},
1238 @math{m = 2^31 - 1},
1241 @c{$m = 2^{31}-249$}
1242 @math{m = 2^31 - 249},
1244 @math{a = 1566083941},
1249 @deffn {Generator} gsl_rng_fishman2x
1250 This is the L'Ecuyer--Fishman random number generator. It is taken from
1251 Knuth's @cite{Seminumerical Algorithms}, 3rd Ed., page 108. Its sequence
1256 z_{n+1} = (x_n - y_n) \,\hbox{mod}\, m
1263 z_@{n+1@} = (x_n - y_n) mod m
1268 with @c{$m = 2^{31}-1$}
1269 @math{m = 2^31 - 1}.
1270 @math{x_n} and @math{y_n} are given by the @code{fishman20}
1271 and @code{lecuyer21} algorithms.
1272 The seed specifies the initial value,
1279 @deffn {Generator} gsl_rng_coveyou
1280 This is the Coveyou random number generator. It is taken from Knuth's
1281 @cite{Seminumerical Algorithms}, 3rd Ed., Section 3.2.2. Its sequence
1286 x_{n+1} = (x_n (x_n + 1)) \,\hbox{mod}\, m
1293 x_@{n+1@} = (x_n (x_n + 1)) mod m
1298 with @c{$m = 2^{32}$}
1300 The seed specifies the initial value,
1309 @node Random Number Generator Performance
1310 @section Performance
1313 @comment I made the original plot like this
1314 @comment ./benchmark > tmp; cat tmp | perl -n -e '($n,$s) = split(" ",$_); printf("%17s ",$n); print "-" x ($s/1e5), "\n";'
1317 The following table shows the relative performance of a selection the
1318 available random number generators. The fastest simulation quality
1319 generators are @code{taus}, @code{gfsr4} and @code{mt19937}. The
1320 generators which offer the best mathematically-proven quality are those
1321 based on the @sc{ranlux} algorithm.
1323 @comment The large number of generators based on single linear congruences are
1324 @comment represented by the @code{random} generator below. These generators are
1325 @comment fast but have the lowest statistical quality.
1328 1754 k ints/sec, 870 k doubles/sec, taus
1329 1613 k ints/sec, 855 k doubles/sec, gfsr4
1330 1370 k ints/sec, 769 k doubles/sec, mt19937
1331 565 k ints/sec, 571 k doubles/sec, ranlxs0
1332 400 k ints/sec, 405 k doubles/sec, ranlxs1
1333 490 k ints/sec, 389 k doubles/sec, mrg
1334 407 k ints/sec, 297 k doubles/sec, ranlux
1335 243 k ints/sec, 254 k doubles/sec, ranlxd1
1336 251 k ints/sec, 253 k doubles/sec, ranlxs2
1337 238 k ints/sec, 215 k doubles/sec, cmrg
1338 247 k ints/sec, 198 k doubles/sec, ranlux389
1339 141 k ints/sec, 140 k doubles/sec, ranlxd2
1341 1852 k ints/sec, 935 k doubles/sec, ran3
1342 813 k ints/sec, 575 k doubles/sec, ran0
1343 787 k ints/sec, 476 k doubles/sec, ran1
1344 379 k ints/sec, 292 k doubles/sec, ran2
1347 @node Random Number Generator Examples
1350 The following program demonstrates the use of a random number generator
1351 to produce uniform random numbers in the range [0.0, 1.0),
1354 @verbatiminclude examples/rngunif.c
1358 Here is the output of the program,
1362 @verbatiminclude examples/rngunif.out
1366 The numbers depend on the seed used by the generator. The default seed
1367 can be changed with the @code{GSL_RNG_SEED} environment variable to
1368 produce a different stream of numbers. The generator itself can be
1369 changed using the environment variable @code{GSL_RNG_TYPE}. Here is the
1370 output of the program using a seed value of 123 and the
1371 multiple-recursive generator @code{mrg},
1374 $ GSL_RNG_SEED=123 GSL_RNG_TYPE=mrg ./a.out
1375 @verbatiminclude examples/rngunif.2.out
1378 @node Random Number References and Further Reading
1379 @section References and Further Reading
1381 The subject of random number generation and testing is reviewed
1382 extensively in Knuth's @cite{Seminumerical Algorithms}.
1386 Donald E. Knuth, @cite{The Art of Computer Programming: Seminumerical
1387 Algorithms} (Vol 2, 3rd Ed, 1997), Addison-Wesley, ISBN 0201896842.
1391 Further information is available in the review paper written by Pierre
1395 P. L'Ecuyer, ``Random Number Generation'', Chapter 4 of the
1396 Handbook on Simulation, Jerry Banks Ed., Wiley, 1998, 93--137.
1398 @uref{http://www.iro.umontreal.ca/~lecuyer/papers.html}
1399 in the file @file{handsim.ps}.
1403 The source code for the @sc{diehard} random number generator tests is also
1408 @cite{DIEHARD source code} G. Marsaglia,
1410 @uref{http://stat.fsu.edu/pub/diehard/}
1414 A comprehensive set of random number generator tests is available from
1419 NIST Special Publication 800-22, ``A Statistical Test Suite for the
1420 Validation of Random Number Generators and Pseudo Random Number
1421 Generators for Cryptographic Applications''.
1423 @uref{http://csrc.nist.gov/rng/}
1426 @node Random Number Acknowledgements
1427 @section Acknowledgements
1429 Thanks to Makoto Matsumoto, Takuji Nishimura and Yoshiharu Kurita for
1430 making the source code to their generators (MT19937, MM&TN; TT800,
1431 MM&YK) available under the GNU General Public License. Thanks to Martin
1432 L@"uscher for providing notes and source code for the @sc{ranlxs} and
1433 @sc{ranlxd} generators.
1436 @comment [ LCG(n) := n * 69069 mod (2^32) ]
1437 @comment First 6: [69069, 475559465, 2801775573, 1790562961, 3104832285, 4238970681]
1438 @comment %2^31-1 69069, 475559465, 654291926, 1790562961, 957348638, 2091487034
1440 @comment [q([x1, x2, x3, x4, x5]) := [107374182 mod 2147483647 * x1 + 104480 mod 2147483647 * x5, x1, x2, x3, x4]]
1443 @comment [q1([x1,x2,x3]) := [63308 mod 2147483647 * x2 -183326 mod 2147483647 * x3, x1, x2],
1444 @comment q2([x1,x2,x3]) := [86098 mod 2145483479 * x1 -539608 mod 2145483479 * x3, x1, x2] ]
1445 @comment initial for q1 is [69069, 475559465, 654291926]
1446 @comment initial for q2 is [1790562961, 959348806, 2093487202]
1449 @comment [ b1(x) := rsh(xor(lsh(x, 13), x), 19),
1450 @comment q1(x) := xor(lsh(and(x, 4294967294), 12), b1(x)),
1451 @comment b2(x) := rsh(xor(lsh(x, 2), x), 25),
1452 @comment q2(x) := xor(lsh(and(x, 4294967288), 4), b2(x)),
1453 @comment b3(x) := rsh(xor(lsh(x, 3), x), 11),
1454 @comment q3(x) := xor(lsh(and(x, 4294967280), 17), b3(x)) ]
1455 @comment [s1, s2, s3] = [600098857, 1131373026, 1223067536]
1456 @comment [2948905028, 441213979, 394017882]