rng.texi

开放gsl矩阵运算

@cindex RANLXD random number generator

These generators produce double precision output (48 bits) from the
@sc{ranlxs} generator.  The library provides two luxury levels
@code{ranlxd1} and @code{ranlxd2}.
@end deffn


@deffn {Generator} gsl_rng_ranlux
@deffnx {Generator} gsl_rng_ranlux389

@cindex RANLUX random number generator
The @code{ranlux} generator is an implementation of the original
algorithm developed by L@"uscher.  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
@sc{ranlxs} and @sc{ranlxd} described above provide floating-point
implementations which will be faster on many platforms.
The period of the generator is about @c{$10^{171}$} 
@math{10^171}.  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@"uscher
is provided by @code{gsl_rng_ranlux}, while @code{gsl_rng_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,
@itemize @asis
@item
M. L@"uscher, "A portable high-quality random number generator for
lattice field theory calculations", @cite{Computer Physics
Communications}, 79 (1994) 100-110.
@item
F. James, "RANLUX: A Fortran implementation of the high-quality
pseudo-random number generator of L@"uscher", @cite{Computer Physics
Communications}, 79 (1994) 111-114
@end itemize
@end deffn


@deffn {Generator} gsl_rng_cmrg
@cindex CMRG, combined multiple recursive random number generator
This is a combined multiple recursive generator by L'Ecuyer. 
Its sequence is,

@tex
\beforedisplay
$$
z_n = (x_n - y_n) \,\hbox{mod}\, m_1
$$
\afterdisplay
@end tex
@ifinfo
@example
z_n = (x_n - y_n) mod m_1
@end example
@end ifinfo
@noindent
where the two underlying generators @math{x_n} and @math{y_n} are,

@tex
\beforedisplay
$$
\eqalign{ 
x_n & = (a_1 x_{n-1} + a_2 x_{n-2} + a_3 x_{n-3}) \,\hbox{mod}\, m_1 \cr
y_n & = (b_1 y_{n-1} + b_2 y_{n-2} + b_3 y_{n-3}) \,\hbox{mod}\, m_2
}
$$
\afterdisplay
@end tex
@ifinfo
@example
x_n = (a_1 x_@{n-1@} + a_2 x_@{n-2@} + a_3 x_@{n-3@}) mod m_1
y_n = (b_1 y_@{n-1@} + b_2 y_@{n-2@} + b_3 y_@{n-3@}) mod m_2
@end example
@end ifinfo

@noindent
with coefficients 
@math{a_1 = 0}, 
@math{a_2 = 63308}, 
@math{a_3 = -183326},
@math{b_1 = 86098}, 
@math{b_2 = 0},
@math{b_3 = -539608},
and moduli 
@c{$m_1 = 2^{31} - 1 = 2147483647$} 
@math{m_1 = 2^31 - 1 = 2147483647}
and 
@c{$m_2 = 2145483479$}
@math{m_2 = 2145483479}.

The period of this generator is 
@c{$2^{205}$}
@math{2^205} 
(about 
@c{$10^{61}$}
@math{10^61}).  It uses
6 words of state per generator.  For more information see,

@itemize @asis
@item
P. L'Ecuyer, "Combined Multiple Recursive Random Number
Generators," @cite{Operations Research}, 44, 5 (1996), 816--822.
@end itemize
@end deffn

@deffn {Generator} gsl_rng_mrg
@cindex MRG, multiple recursive random number generator
This is a fifth-order multiple recursive generator by L'Ecuyer, Blouin
and Coutre.  Its sequence is,

@tex
\beforedisplay
$$
x_n = (a_1 x_{n-1} + a_5 x_{n-5}) \,\hbox{mod}\, m
$$
\afterdisplay
@end tex
@ifinfo
@example
x_n = (a_1 x_@{n-1@} + a_5 x_@{n-5@}) mod m
@end example
@end ifinfo

@noindent
with 
@math{a_1 = 107374182}, 
@math{a_2 = a_3 = a_4 = 0}, 
@math{a_5 = 104480}
and 
@c{$m = 2^{31}-1$}
@math{m = 2^31 - 1}.

The period of this generator is about 
@c{$10^{46}$}
@math{10^46}.  It uses 5 words
of state per generator.  More information can be found in the following
paper,
@itemize @asis
@item
P. L'Ecuyer, F. Blouin, and R. Coutre, "A search for good multiple
recursive random number generators", @cite{ACM Transactions on Modeling and
Computer Simulation} 3, 87-98 (1993).
@end itemize
@end deffn

@deffn {Generator} gsl_rng_taus
@cindex Tausworthe random number generator
This is a maximally equidistributed combined Tausworthe generator by
L'Ecuyer.  The sequence is,

@tex
\beforedisplay
$$
x_n = (s^1_n \oplus s^2_n \oplus s^3_n) 
$$
\afterdisplay
@end tex
@ifinfo
@example
x_n = (s1_n ^^ s2_n ^^ s3_n) 
@end example
@end ifinfo

@noindent
where,

@tex
\beforedisplay
$$
\eqalign{
s^1_{n+1} &= (((s^1_n \& 4294967294)\ll 12) \oplus (((s^1_n\ll 13) \oplus s^1_n)\gg 19)) \cr
s^2_{n+1} &= (((s^2_n \& 4294967288)\ll 4) \oplus (((s^2_n\ll 2) \oplus s^2_n)\gg 25)) \cr
s^3_{n+1} &= (((s^3_n \& 4294967280)\ll 17) \oplus (((s^3_n\ll 3) \oplus s^3_n)\gg 11))
}
$$
\afterdisplay
@end tex
@ifinfo
@example
s1_@{n+1@} = (((s1_n&4294967294)<<12)^^(((s1_n<<13)^^s1_n)>>19))
s2_@{n+1@} = (((s2_n&4294967288)<< 4)^^(((s2_n<< 2)^^s2_n)>>25))
s3_@{n+1@} = (((s3_n&4294967280)<<17)^^(((s3_n<< 3)^^s3_n)>>11))
@end example
@end ifinfo

@noindent
computed modulo 
@c{$2^{32}$}
@math{2^32}.  In the formulas above 
@c{$\oplus$}
@math{^^}
denotes ``exclusive-or''.  Note that the algorithm relies on the properties
of 32-bit unsigned integers and has been implemented using a bitmask
of @code{0xFFFFFFFF} to make it work on 64 bit machines.

The period of this generator is @c{$2^{88}$}
@math{2^88} (about
@c{$10^{26}$}
@math{10^26}).  It uses 3 words of state per generator.  For more
information see,

@itemize @asis
@item
P. L'Ecuyer, "Maximally Equidistributed Combined Tausworthe
Generators", @cite{Mathematics of Computation}, 65, 213 (1996), 203--213.
@end itemize
@end deffn

@deffn {Generator} gsl_rng_gfsr4
@cindex Four-tap Generalized Feedback Shift Register
The @code{gfsr4} generator is like a lagged-fibonacci generator, and 
produces each number as an @code{xor}'d sum of four previous values.

@tex
\beforedisplay
$$
r_n = r_{n-A} \oplus r_{n-B} \oplus r_{n-C} \oplus r_{n-D}
$$
\afterdisplay
@end tex
@ifinfo
@example
r_n = r_@{n-A@} ^^ r_@{n-B@} ^^ r_@{n-C@} ^^ r_@{n-D@}
@end example
@end ifinfo

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 the example on p392 of
Ziff's article: @math{A=471}, @math{B=1586}, @math{C=6988}, @math{D=9689}.

If the offsets are appropriately chosen (such the one ones in
this implementation), then the sequence is said to be maximal.
I'm not sure what that means, but I would guess that means all
states are part of the same cycle, which would mean that the
period for this generator is astronomical; it is
@c{$(2^K)^D \approx 10^{93334}$}
@math{(2^K)^D \approx 10^@{93334@}}
where @math{K=32} is the number of bits in the word, and D is the longest
lag.  This would also mean that any one random number could 
easily be zero; ie 
@c{$0 \le r < 2^{32}$}
@math{0 <= r < 2^32}.

Ziff doesn't say so, but it seems to me that the bits are
completely independent here, so one could use this as an efficient
bit generator; each number supplying 32 random bits.  The quality of the
generated bits depends on the underlying seeding procedure, which 
may need to be improved in some circumstances.

For more information see,
@itemize @asis
@item
Robert M. Ziff, "Four-tap shift-register-sequence random-number 
generators", @cite{Computers in Physics}, 12(4), Jul/Aug
1998, pp 385-392.
@end itemize
@end deffn

@node Unix random number generators
@section Unix random number generators

The standard Unix random number generators @code{rand}, @code{random}
and @code{rand48} are provided as part of GSL. Although these
generators are widely available individually often they aren't all
available on the same platform.  This makes it difficult to write
portable code using them and so we have included the complete set of
Unix generators in GSL for convenience.  Note that these generators
don't produce high-quality randomness and aren't suitable for work
requiring accurate statistics.  However, if you won't be measuring
statistical quantities and just want to introduce some variation into
your program then these generators are quite acceptable.

@cindex BSD random number generator, rand
@cindex Unix random number generators, rand
@cindex Unix random number generators, rand48

@deffn {Generator} gsl_rng_rand
@cindex BSD random number generator
This is the BSD @code{rand()} generator.  Its sequence is

@tex
\beforedisplay
$$
x_{n+1} = (a x_n + c) \,\hbox{mod}\, m
$$
\afterdisplay
@end tex
@ifinfo
@example
x_@{n+1@} = (a x_n + c) mod m
@end example
@end ifinfo

@noindent
with 
@math{a = 1103515245}, 
@math{c = 12345} and 
@c{$m = 2^{31}$}
@math{m = 2^31}.
The seed specifies the initial value, 
@math{x_1}.  The period of this
generator is 
@c{$2^{31}$}
@math{2^31}, and it uses 1 word of storage per
generator.
@end deffn

@deffn {Generator} gsl_rng_random_bsd
@deffnx {Generator} gsl_rng_random_libc5
@deffnx {Generator} gsl_rng_random_glibc2
These generators implement the @code{random()} family of functions, a
set of linear feedback shift register generators originally used in BSD
Unix.  There are several versions of @code{random()} in use today: the
original BSD version (e.g. on SunOS4), a libc5 version (found on
older GNU/Linux systems) and a glibc2 version.  Each version uses a
different seeding procedure, and thus produces different sequences.

The original BSD routines accepted a variable length buffer for the
generator state, with longer buffers providing higher-quality
randomness.  The @code{random()} function implemented algorithms for
buffer lengths of 8, 32, 64, 128 and 256 bytes, and the algorithm with
the largest length that would fit into the user-supplied buffer was
used.  To support these algorithms additional generators are available
with the following names,

@example
gsl_rng_random8_bsd
gsl_rng_random32_bsd
gsl_rng_random64_bsd
gsl_rng_random128_bsd
gsl_rng_random256_bsd
@end example

@noindent
where the numeric suffix indicates the buffer length.  The original BSD
@code{random} function used a 128-byte default buffer and so
@code{gsl_rng_random_bsd} has been made equivalent to
@code{gsl_rng_random128_bsd}.  Corresponding versions of the @code{libc5}
and @code{glibc2} generators are also available, with the names
@code{gsl_rng_random8_libc5}, @code{gsl_rng_random8_glibc2}, etc.
@end deffn

@deffn {Generator} gsl_rng_rand48
@cindex rand48 random number generator
This is the Unix @code{rand48} generator.  Its sequence is

@tex
\beforedisplay
$$
x_{n+1} = (a x_n + c) \,\hbox{mod}\, m
$$
\afterdisplay
@end tex
@ifinfo
@example
x_@{n+1@} = (a x_n + c) mod m
@end example
@end ifinfo
@noindent
defined on 48-bit unsigned integers with 
@math{a = 25214903917}, 
@math{c = 11} and 
@c{$m = 2^{48}$} 
@math{m = 2^48}. 
The seed specifies the upper 32 bits of the initial value, @math{x_1},
with the lower 16 bits set to @code{0x330E}.  The function
@code{gsl_rng_get} returns the upper 32 bits from each term of the
sequence.  This does not have a direct parallel in the original
@code{rand48} functions, but forcing the result to type @code{long int}
reproduces the output of @code{mrand48}.  The function
@code{gsl_rng_uniform} uses the full 48 bits of internal state to return
the double precision number @math{x_n/m}, which is equivalent to the
function @code{drand48}.  Note that some versions of the GNU C Library
contained a bug in @code{mrand48} function which caused it to produce
different results (only the lower 16-bits of the return value were set).
@end deffn

@node Numerical Recipes generators
@section Numerical Recipes generators
@cindex Random number generators, Numerical recipes
@cindex Numerical recipes, random number generators
@comment
The following generators are provided for compatibility with
@cite{Numerical Recipes}.  Note that the original Numerical Recipes
functions used single precision while we use double precision.  This will
lead to minor discrepancies, but only at the level of single-precision
rounding error.  If necessary you can force the returned values to single
precision by storing them in a @code{volatile float}, which prevents the
value being held in a register with double or extended precision.  Apart
from this difference the underlying algorithms for the integer part of
the generators are the same.

@deffn {Generator} gsl_rng_ran0 
Numerical recipes @code{ran0} implements Park and Miller's @sc{minstd}
algorithm with a modified seeding procedure.
@end deffn

@deffn {Generator} gsl_rng_ran1 
Numerical recipes @code{ran1} implements Park and Miller's @sc{minstd}
algorithm with a 32-element Bayes-Durham shuffle box.
@end deffn

@deffn {Generator} gsl_rng_ran2 
Numerical recipes @code{ran2} implements a L'Ecuyer combined recursive
generator with a 32-element Bayes-Durham shuffle-box.
@end deffn

@deffn {Generator} gsl_rng_ran3 
Numerical recipes @code{ran3} implements Knuth's portable
subtractive generator.