$Id: proposal.html,v 1.50 2003/04/11 07:26:32 jmaurer Exp $

A Proposal to Add an Extensible Random Number Facility to the Standard Library (Revision 2)

Any one who considers arithmetical methods of producing random digits is, of course, in a state of sin.

John von Neumann, 1951

Revision history

• 2003-04-10: Publication in the Post-Oxford mailing as N1452.
• Fix editing glitch.
• 2002-11-10: Publication in the Post-Santa Cruz mailing as N1398.
• The seed(first, last) interface now needs "unsigned long" values.
• Introduce "variate_generator", adjust distribution interface accordingly.
• Add "add-on packages" discussion.
• All distribution parameters must be defaulted.
• Add "target audience" subsection to "motivation" section.
• Add discussion of manager class.
• Engines are independent of distributions, thus consider respective lifetimes.
• Add "sharing of engines" as a major requirement.
• Add some open issues.
• 2002-10-11: First publication on the C++ committee's library reflector.

I. Motivation

Why is this important? What kinds of problems does it address, and what kinds of programmers, is it intended to support? Is it based on existing practice?

• numerics (simulation, Monte-Carlo integration)
• games (shuffling card decks, non-deterministic enemy behavior)
• testing (generation of test input data for good coverage)
• security (generation of cryptographic keys)

Programmers in all of the above areas have to find ways to generate random numbers. However, the difficulty to find generators that are both efficient and have good quality is often underestimated, and so ad-hoc implementations often fail to meet either or both of these goals.

The C++ standard library includes std::rand, inherited from the C standard library, as the only facility to generate pseudo-random numbers. It is underspecified, because the generation function is not defined, and indeed early C standard library implementations provided surprisingly bad generators. Furthermore, the interface relies on global state, making it difficult or inefficient to provide for correct operation for simultaneous invocations in multi-threaded applications.

There is a lot of existing practice in this area. A multitude of libraries, usually implemented in C or Fortran, is available from the scientific community. Some implement just one random number engine, others seek to provide a full framework. I know of no comprehensive C++ framework for generating random numbers that adheres to the design principles put forth in section III.

Random number generators are appropriate for this TR because they fall into one of the domains (numerics) identified in N1314 as a target for the TR.

Target Audience

• programmers that provide additional engines
• programmers that provide additional distributions
• programmers that provide generic add-on packages
• programmers that need random numbers

II. Impact On the Standard

What does it depend on, and what depends on it? Is it a pure extension, or does it require changes to standard components? Does it require core language changes?

The ISO C99 extension that specify integral types having a given minimum or exact bitwidth (e.g. int32_t) aids in implementing this proposal, however these types (or the equivalent thereof under another name) can be defined with template metaprogramming in standard C++, so these are not strictly necessary.

In case the ISO C99 extensions become part of the TR, section IV should be reviewed whether some requirements could be reformulated with the ISO C99 extensions.

In case a standard reference-counted smart pointer becomes part of the TR, section IV should be reviewed and instances of the smart pointer be added to the acceptable template parameters for a variate_generator.

III. Design Decisions

Why did you choose the specific design that you did? What alternatives did you consider, and what are the tradeoffs? What are the consequences of your choice, for users and implementors? What decisions are left up to implementors? If there are any similar libraries in use, how do their design decisions compare to yours?

• CLHEP (original at http://wwwinfo.cern.ch/asd/lhc++/clhep/index.html, modifications from FermiLab at (anonymous CVS) :pserver:anonymous@zoomcvs.fnal.gov:/usr/people/cvsuser/repository)
• crng 1.1: Random-number generators (RNGs) implemented as Python extension types coded in C (at http://www.sbc.su.se/~per/crng/)
• Swarm 2.1.1 (multi-agent simulation of complex systems), random number package, using a Smalltalk-like programming language (at http://www.santafe.edu/projects/swarm/swarmdocs/set/swarm.random.sgml.reference.html)
• GNU Scientific Library: general scientific computing library implemented in C, comprehensive coverage of random number engines and distributions (at http://sources.redhat.com/gsl)

• Donald E. Knuth, "The Art of Computer Programming Vol. 2"
• William H. Press et al., "Numerical Recipes in C"

A. Overview on Requirements

• allows users to choose in speed / size / quality trade-offs
• has a tight enough specification to get reliable cross-platform results
• allows storage of state on non-volatile media (e.g., in a disk file) to resume computation later
• does not impede sequence "jump-ahead" for parallel computation
• provides a variety of base engines, not just one
• allows the user to write its own base engines and use it with the library-provided distributions
• provides the most popular distributions
• allows the user to write its own distributions and use it with the library-provided engines
• allows sharing of engines by several distributions
• does not prevent implementations with utmost efficiency
• provides both pseudo-random number engines (for simulations etc.) and "true" non-deterministic random numbers (for cryptography)

B. Pseudo-Random vs. Non-Deterministic Random Numbers

In this proposal, a pseudo-random number engine is defined as an initial internal state x(0), a function f that moves from one internal state to the next x(i+1) := f(x(i)), and an output function o that produces the output o(x(i)) of the generator. This is an entirely deterministic process, it is determined by the initial state x(0) and functions f and o only. The initial state x(0) is determined from a seed. Apparent randomness is achieved only because the user has limited perception.

A non-deterministic random-number engine provides a sequence of random numbers x(i) that cannot be foreseen. Examples are certain quantum-level physics experiments, measuring the time difference between radioactive decay of individual atoms or noise of a Zehner diode. Relatively unforeseeable random sources are also (the low bits of) timing between key touches, mouse movements, Ethernet packet arrivals, etc. An estimate for the amount of unforeseeability is the entropy, a concept from information theory. Completely foreseeable sequences (e.g., from pseudo-random number engines) have entropy 0, if all bits are unforeseeable, the entropy is equal to the number of bits in each number.

Pseudo-random number engines are usually much faster than non-deterministic random-number engines, because the latter require I/O to query some randomness device outside of the computer. However, there is a common interface feature subset of both pseudo-random and non-deterministic random-number engines. For example, a non-deterministic random-number engine could be employed to produce random numbers with normal distribution; I believe this to be an unlikely scenario in practice.

Other libraries, including those mentioned above, only provide either pseudo-random numbers, suitable for simulations and games, or non-deterministic random numbers, suitable for cryptographic applications.

C. Separation of Engines and Distributions

This proposal honours this conceptual separation, and provides a class template to merge an arbitrary engine with an arbitrary distribution on top. To this end, this proposal sets up requirements for engines so that each of them can be used to provide uniformly distributed random numbers for any of the distributions. The resulting freedom of combination allows for the utmost re-use.

Engines have usually been analyzed with all mathematical and empirical tools currently available. Nonetheless, those tools show the absence of a particular weakness only, and are not exhaustive. Albeit unlikely, a new kind of test (for example, a use of random numbers in a new kind of simulation or game) could show serious weaknesses in some engines that were not known before.

This proposal attempts to specify the engines precisely; two different implementations, with the same seed, should return the same output sequence. This forces implementations to use the well-researched engines specified hereinafter, and users can have confidence in their quality and the limits thereof.

On the other hand, the specifications for the distributions only define the statistical result, not the precise algorithm to use. This is different from engines, because for distribution algorithms, rigorous proofs of their correctness are available, usually under the precondition that the input random numbers are (truely) uniformly distributed. For example, there are at least a handful of algorithms known to produce normally distributed random numbers from uniformly distributed ones. Which one of these is most efficient depends on at least the relative execution speeds for various transcendental functions, cache and branch prediction behaviour of the CPU, and desired memory use. This proposal therefore leaves the choice of the algorithm to the implementation. It follows that output sequences for the distributions will not be identical across implementations. It is expected that implementations will carefully choose the algorithms for distributions up front, since it is certainly surprising to customers if some distribution produces different numbers from one implementation version to the next.

Other libraries usually provide the same differentiation between engines and distributions. Libraries rarely have a wrapper around both engine and distribution, but it turns out that this can hide some complexities from the authors of distributions, since some facitilies need to be provided only once. A previous version of this proposal had distributions directly exposed to the user, and the distribution type dependent on the engine type. In various discussions, this was considered as too much coupling.

Since other libraries do not aim to provide a portable specification framework, engines are sometimes only described qualitatively without giving the exact parameterization. Also, distributions are given as specific functions or classes, so the quality-of-implementation question which distribution algorithm to employ does not need to be addressed.

D. Templates vs. Virtual Functions

For virtual functions in a class hierarchy, the core language requires a (nearly) exact type match for a function in a derived classes overriding a function in a base class. This seems to be unnecessarily restrictive, because engines can sometimes benefit from using different integral base types. Also, with current compiler technology, virtual functions prevent inlining when a pointer to the base class is used to call a virtual function that is overridden in some derived class. In particular with applications such as simulations that sometimes use millions of pseudo-random numbers per second, losing significant amounts of performance due to missed inlining opportunities appears to not be acceptable.

The CLHEP library bases all its engines on the abstract base class HepRandomEngine. Specific engines derive from this class and override its pure virtual functions. Similarly, all distributions are based on the base class HepRandom. Specific distributions derive from this class, override operator(), and provide a number of specific non-virtual functions.

The GNU Scientific Library, while coded in C, adheres to the principles of object-structuring; all engines can be used with any of the distributions. The technical implementation is by mechanisms similar to virtual functions.

E. Parameterization and Initialization for Engines

Since random-number engines are mathematically designed with computer implementation in mind, parameters are usually integers representable in a machine word, which usually coincides nicely with a C++ built-in type. The parameters could either be given as (compile-time) template arguments or as (run-time) constructor arguments.

Providing parameters as template arguments allows for providing predefined parameterizations as simple "typedef"s. Furthermore, the parameters appear as integral constants, so the compiler can value-check the given constants against the engine's base type. Also, the library implementor can choose different implementations depending on the values of the parameters, without incurring any runtime overhead. For example, there is an efficient method to compute (a*x) mod m, provided that a certain magnitude of m relative to the underlying type is not exceeded. Additionally, the compiler's optimizer can benefit from the constants and potentially produce better code, for example by unrolling loops with fixed loop count.

As an alternative, providing parameters as constructor arguments allows for more flexibility for the library user, for example when experimenting with several parameterizations. Predefined parameterizations can be provided by defining wrapper types which default the constructor parameters.

Other libraries have hard-coded the parameters of their engines and do not allow the user any configuration of them at all. If the user wishes to change the parameters, he has to re-implement the engine's algorithm. In my opinion, this approach unnecessarily restricts re-use.

Regarding initialization, this proposal chooses to provide "deterministic seeding" with the default constructor and the seed function without parameters: Two engines constructed using the default constructor will output the same sequence. In contrast, the CLHEP library's default constructed engines will take a fresh seed from a seed table for each instance. While this approach may be convenient for a certain group of users, it relies on global state and can easily be emulated by appropriately wrapping engines with deterministic seeding.

In addition to the default constructor, all engines provide a constructor and seed function taking an iterator range [it1,it2) pointing to unsigned integral values. An engine initializes its state by successively consuming values from the iterator range, then returning the advanced iterator it1. This approach has the advantage that the user can completely exploit the large state of some engines for initialization. Also, it allows to initialize compound engines in a uniform manner. For example, a compound engine consisting of two simpler engines would initialize the first engine with its [it1,it2). The first engine returns a smaller iterator range that it has not consumed yet. This can be used to initialize the second engine.

The iterator range [it1,it2) is specified to point to unsigned long values. There is no way to determine from a generic user program how the initialization values will be treated and what range of bits must be provided, except by enumerating all engines, e.g. in template specializations. The problem is that a given generator might have differing requirements on the values of the seed range even within one seed call.

For example, imagine a

   xor_combine<lagged_fibonacci<...>, mersenne_twister<...> >

How does a concise programming interface for those increasingly complex and varying requirements on [first, last) look like? I don't know, and I don't want to complicate the specification by inventing something complicated here.

Thus, the specification says for each generator how it uses the seed values, and how many are consumed. Additional features are left to the user.

In a way, this is similar to STL containers: It is intended that the user can exchange iterators to various containers in generic algorithms, but the container itself is not meant to be exchanged, i.e. having a Container template parameter is often not adequate. That is analogous to the random number case: The user can pass an engine around and use its operator() and min and max functions generically. However, the user can't generically query the engine attributes and parameters, simply because most are entirely different in semantics for each engine.

The seed(first, last) interface can serve two purposes:

1. In a generic context, the user can pass several integer values >= 1 for seeding. It is unlikely that the user explores the full state space with the seeds she provides, but she can be reasonably sure that her seeds aren't entirely incorrect. (There is no formal guarantee for that, except that the ability to provide bad seeds usually means the parameterization of the engine is bad, e.g. a non-prime modulus for a linear congruential engine.) For example, if the user wants a seed(uint32_t) on top of seed(first, last), one option is to use a linear_congruential generator that produces the values required for seed(first, last). When the user defines the iterator type for first and last so that it encapsulates the linear_congruential engine in operator++, the user doesn't even need to know beforehand how many values seed(first, last) will need.
2. If the user is in a non-generic context, he knows the specific template type of the engine (probably not the template value-based parameterization, though). The precise specification for seed(first, last) allows to know what values need to be passed in so that a specific initial state is attained, for example to compare one implementation of the engine with another one that uses different seeding.
3. If the user requires both, he needs to inject knowledge into (1) so that he is in the position of (2). One way to inject the knowledge is to use (partial) template specialization to add the knowledge. The specific parameterization of some engine can then be obtained by querying the data members of the engines.

I haven't seen the iterator-based approach to engine initialization in other libraries; most initialization approaches rely on a either a single value or on per-engine specific approaches to initialization.

An alternative approach is to pass a zero-argument function object ("generator") for seeding. It is trivial to implement a generator from a given iterator range, but it is more complicated to implement an iterator range from a generator. Also, the exception object that is specified to be thrown when the iterator range is exhausted could be configured in a user-provided iterator to generator mapping. With this approach, some engines would have three one-argument constructors: One taking a single integer for seeding, one taking a (reference?) to a (templated) generator, and the copy constructor. It appears that the opportunities for ambiguities or choosing the wrong overload are too confusing to the unsuspecting user.

F. Parameterization and Initialization for Distributions

The probability density functions of distributions usually have parameters. These are mapped to constructor parameters, to be set at runtime by the library user according to her requirements. The parameters for a distribution object cannot change after its construction. When constructing the distribution, this allows to pre-compute some data according to the parameters given without risk of inadvertently invalidating them later.

Distributions may implement operator()(T x), for arbitrary type T, to meet special needs, for example a "one-shot" mode where each invocation uses different distribution parameters.

G. Properties as Traits vs. In-Class Constants

H. Which Engines to Include

Engines may be categorized according to the following dimensions.

• integers or floating-point numbers produced (Some engines produce uniformly distributed integers in the range [min,max], however, most distribution functions expect uniformly distributed floating-point numbers in the range [0,1) as the input sequence. The obvious conversion requires a relatively costly integer to floating-point conversion plus a floating-point multiplication by (max-min+1)^-1 for each random number used. To save the multiplication, some engines can directly produce floating-point numbers in the range [0,1) by maintaining the state x(i) in an appropriately normalized form, given a sufficiently good implementation of basic floating-point operations (e.g. IEEE 754).
• quality of random numbers produced (What is the cycle length? Does the engine pass all relevant statistical tests? Up to what dimension are numbers equidistributed?)
• speed of generation (How many and what kind of operations have to be performed to produce one random number, on average?)
• size of state (How may machine words of storage are required to hold the state x(i) of the random engine?)
• option for independent subsequences (Is it possible to move from x(i) to x(i+k) with at most O(log(k)) steps? This allows to efficiently use subsequences x(0)...x(k-1), x(k)...x(2k-1), ..., x(jk)...x((j+1)k-1), ..., for example for parallel computation, where each of the m processors gets assigned the (independent) subsequence starting at x(jk) (0 <= k < m).)

engine	int / float	quality	speed	size of state	subsequences	comments
linear_congruential	int	medium	medium	1 word	yes	cycle length is limited to the maximum value representable in one machine word, passes most statisticial tests with chosen parameters.
mersenne_twister	int	good	fast	624 words	no	long cycles, passes all statistical tests, good equidistribution in high dimensions
subtract_with_carry	both	medium	fast	25 words	no	very long cycles possible, fails some statistical tests. Can be improved with the discard_block compound engine.
discard_block	both	good	slow	base engine + 1 word	no	compound engine that removes correlation provably by throwing away significant chunks of the base engine's sequence, the resulting speed is reduced to 10% to 3% of the base engine's.
xor_combine	int	good	fast	base engines	yes, if one of the base engines	compound engine that XOR-combines the sequences of two base engines

Some engines were considered for inclusion, but left out for the following reasons:

engine	int / float	quality	speed	size of state	subsequences	comments
shuffle_output	int	good	fast	base engine + 100 words	no	compound engine that reorders the base engine's output, little overhead for generation (one multiplication)
lagged_fibonacci	both	medium	fast	up to 80,000 words	no	very long cycles possible, fails birthday spacings test. Same principle of generation as subtract_with_carry, i.e. x(i) = x(i-s) (*) x(i-r), where (*) is either of +, -, xor with or without carry.
inversive_congruential (Hellekalek 1995)	int	good	slow	1 word	no	x(i+1) = a x(i)-1 + c. Good equidistribution in several dimensions. Provides no apparent advantage compared to ranlux; the latter can produce floating-point numbers directly.
additive_combine (L'Ecuyer 1988)	int	good	medium	2 words	yes	Combines two linear congruential generators. Same principle of combination as xor_combine, i.e. z(i) = x(i) (*) y(i), where (*) is one of +, -, xor.
R250 (Kirkpatrick and Stoll)	int	bad	fast	~ 20 words	no	General Feedback Shift Register with two taps: Easily exploitable correlation.
linear_feedback_shift	int	medium	fast	1 word	no	cycle length is limited to the maximum value representable in one machine word, fails some statistical tests, can be improved with the xor_combine compound engine.

The GNU Scientific Library and Swarm have additional engine that are not mentioned in the table below.

Engine	this proposal	CLHEP	crng	GNU Scientific Library	Swarm	Numerical Recipes	Knuth
LCG(231-1, 16807)	minstd_rand0	-	ParkMiller	ran0, minstd	-	ran0	p106, table 1, line 19
LCG(232, a=1664525, c=1013904223)	linear_congruential< ..., 1664525, 1013904223, (1 << 32) >	-	-	-	LCG1gen	-	p106, table 1, line 16
LCG1 + LCG2 + LCG3	-	-	WichmannHill	-	-	-	-
(LCG1 - LCG2 + LCG3 - LCG4) mod m0	-	-	-	-	C4LCGXgen	-	-
LCG(231-1, 16807) with Bays/Durham shuffle	shuffle_output<minstd_rand0, 32> (shuffle_output not in this proposal)	-	-	ran1	PMMLCG1gen	ran1	Algorithm "B"
(LCG(231-85, 40014) + LCG(231-249, 40692)) mod 231-85	ecuyer1988 (additive_combine not in this proposal)	Ranecu	LEcuyer	-	C2LCGXgen	-	-
(LCG(231-85, 40014) with Bays/Durham shuffle + LCG(231-249, 40692)) mod 231-85	additive_combine< shuffle_output< linear_congruential<int, 40014, 0, 2147483563>, 32>, linear_congruential<int, 40692, 0, 2147483399> > (additive_combine and shuffle_output not in this proposal)	-	-	ran2	-	ran2	-
X(i) = (X(i-55) - X(i-33)) mod 109	-	-	-	ran3	~SCGgen	ran3	-
X(i) = (X(i-100) - X(i-37)) mod 230	-	-	-	-	-	-	ran_array
X(i) = (X(i-55) + X(i-24)) mod 232	lagged_fibonacci< ..., 32, 55, 24, ...> (lagged_fibonacci not in this proposal)	-	-	-	ACGgen	-	-
DEShash(i,j)	-	-	-	-	-	ran4	-
MT	mt19937	MTwistEngine	MT19937	mt19937	MT19937gen	-	-
X(i) = (X(i-37) - X(i-24) - carry) mod 232	subtract_with_carry< ..., (1<<32), 37, 24, ...>	-	-	-	SWB1gen	-	-
X(i) = (X(i-43) - X(i-22) - carry) mod 232-5	subtract_with_carry< ..., (1<<32)-5, 43, 22, ...>	-	-	-	PSWBgen	-	-
RCARRY with block discard by Lüscher	discard_block< subtract_with_carry<...>, ...>	RanluxEngine, Ranlux64Engine	Ranlux	ranlx*	-	-	-
Hurd	-	Hurd160, Hurd288	-	-	-	-	-
physical model by Ranshi	-	Ranshi	-	-	-	-	-
return predefined data	-	NonRandom	-	-	-	-	-
RANMAR: z(i) = (z(i-97) - z(i-33)) mod 224; y(i+1) = (y(i)-c) mod 224-3; X(i) = (z(i) - y(i)) mod 224	additive_combine< lagged_fibonacci< (1<<24), 97, 33, ... >, linear_congruential< (1<<24)-3, 1, c, ...> (additive_combine and lagged_fibonacci not in this proposal)	JamesRandom	-	ranmar	-	-	-
Taus88	taus88 = xor_combine ...	-	Taus88	taus, taus2	-	-	-
Taus60	xor_combine< linear_feedback_shift< 31, 13, 12 >, 0, linear_feedback_shift< 29, 2, 4 >, 2, 0> (linear_feedback_shift not in this proposal)	-	-	-	C2TAUSgen	-	-
GFSR, 4-tap	-	-	-	gfsr4	-	-	-
MRG32k3a	-	-	MRG32k3a	-	-	-	-

I. Which Distributions to Include

• Integer uniform
• Floating-point uniform
• Exponential
• Normal
• Gamma
• Poisson
• Binomial
• Geometric
• Bernoulli

Distribution	this proposal	CLHEP	crng	GNU Scientific Library	Swarm	Numerical Recipes	Knuth
uniform (int)	uniform_int	-	-	-	UniformIntegerDist	-	-
uniform (float)	uniform_real	RandFlat	UniformDeviate	flat	UniformDoubleDist	-	uniform
exponential	exponential_distribution	RandExponential	ExponentialDeviate	exponential	ExponentialDist	exponential	exponential
normal	normal_distribution	RandGauss*	NormalDeviate	gaussian	NormalDist	normal (gaussian)	normal
lognormal	-	-	-	lognormal	LogNormalDist	-	-
gamma	gamma_distribution	RandGamma	GammaDeviate	gamma	GammaDist	gamma	gamma
beta	-	-	BetaDeviate	beta	-	-	beta
poisson	poisson_distribution	Poisson	PoissonDeviate	poisson	PoissonDist	poisson	poisson
binomial	binomial_distribution	RandBinomial	BinomialDeviate	binomial	-	binomial	binomial
geometric	geometric_distribution	-	GeometricDeviate	geometric	-	-	geometric
bernoulli	bernoulli_distribution	-	BernoulliDeviate	bernoulli	BernoulliDist	-	-
random bit	-	RandBit	-	-	RandomBitDist	-	-
breit-wigner	-	RandBreitWigner	-	-	-	-	-
chi-square	-	RandChiSquare	-	chisq	-	-	chi-square
landau	-	Landau	-	landau	-	-	-
F	-	-	-	F	-	-	F (variance-ratio)
t	-	-	-	t	-	-	t

J. Taxonomy of Concepts

NumberGenerator ---- UniformRandomNumberGenerator ---- PseudoRandomNumberGenerator
                \--- RandomDistribution

K. Validation

This is considered an important feature for library implementors and serious users to check whether the provided library on the given platform returns the correct numbers. It could be argued that a library implementor should provide a correct implementation of some standard feature in any case.

No other library I have encountered provides explicit validation values in either their specification or their implementation, although some of them claim to be widely portable.

Another design option for validation that was part of early drafts of this proposal is moving the reference number (10000th value in the sequence) from specification space to implementation space, thus providing a validation(x) static member function for each engine that compares the hard-coded 10000th value of the sequence with some user-provided value x presumeably obtained by actually invoking the random-number engine object 10000 times. Due to the template-based design, this amounted to a "val" template value parameter for each engine, and the validation(x) function reduced to the trivial comparison "val == x". Handling validation for floating-point engines required more machinery, because template value parameters cannot be of floating-point type. Also, from a conceptual perspective, it seemed odd to demand a validation decision from the very entitiy which one wanted to validate.

L. Non-Volatile Storage of Engine and Distribution State

Engines and distributions use the usual idioms of operator<< and operator>>. If the user needs additional processing before or after storage on non-volatile media, there is always the option to use a temporary std::stringstream.

Distributions sometimes store values from their associated source of random numbers across calls to their operator(). For example, a common method for generating normally distributed random numbers is to retrieve two uniformly distributed random numbers and compute two normally distributed random numbers out of them. In order to reset the distribution's random number cache to a defined state, each distribution has a reset member function. It should be called on a distribution whenever its associated engine is exchanged or restored.

M. Values vs. References

The wrapper variate_generator can store engines either by value or by reference, explicitly chosen by the template parameters. This allows to use references to a single engine in several variate_generator, but makes it explicit to the user that he is responsible for the management of the lifetime of the engine.

N. Providing the Probability Density Function in Distributions

O. Implementation-defined behaviour

The precise state-holding base data types for the various engines are left implementation-defined, because engines are usually optimized for binary integers with 32 bits of word size. The specification in this proposal cannot foresee whether a 32 bit quantity on the machine is available in C++ as short, int, long, or not at all. It is up to the implementation to decide which data type fits best. The implementation is required to document the choice of data type, so that users can (non-portably) rely on the precise type, for example for further computation. Should the ISO C99 extensions become part of ISO C++, the implementation-defined types could be replaced by e.g. int_least32_t.

The method how to produce non-deterministic random numbers is considered implementation-defined, because it inherently depends on the implementation and possibly even on the runtime environment: Imagine a platform that has operating system support for randomness collection, e.g. from user keystrokes and Ethernet inter-packet arrival timing (Linux /dev/random does this). If, in some installation, access to the operating system functions providing these services has been restricted, the C++ non-deterministic random number engine has been deprived of its randomness. An implementation is required to document how it obtains the non-deterministic random numbers, because only then can users' confidence in them grow. Confidence is of particular concern in the area of cryptography.

The algorithms how to produce the various distributions are specified as implementation-defined, because there is a vast variety of algorithms known for each distribution. Each has a different trade-off in terms of speed, adaptation to recent computer architectures, and memory use. The implementation is required to document its choice so that the user can judge whether it is acceptable quality-wise.

P. Lower and upper bounds on UniformRandomNumberGenerator

Those bounds are not specified to be tight, because for some engines, the bounds depend on the seeds. The seed can be changed during the lifetime of the engine object, while the values returned by min() and max() are invariant. Therefore, min() and max() must return conservative bounds that are independent of the seed.

Q. With or without manager class

Earlier versions of this propsoal made (potentially user-written) distributions directly visible to (some other) user that wants to get random numbers distributed accordingly ("client"), there was no additional management layer between the distribution and the client.

The following additional features could be provided by the management layer:

• The management layer contains an adaptor (to convert the engine's output into the distribution's input) in addition to the engine and the distribution.
• Adaptors and distributions do not store state, but instead, in each invocation, consume an arbitrary number of input values and produce a fixed number of output values. The management layer is responsible for connecting the engine - adaptor - distribution chain, invoking each part when more numbers are needed from the next part of the chain.
• On request, the management layer is responsible for saving and restoring the buffers that exist between engine, adaptor, and distribution.
• On request, the management layer shares engines with another instance of the management layer.

• Retrieve a uniform integer with value either 1 or 2.
• If 1, return a number with normal distribution.
• If 2, return a number with gamma distribution.

The ability to share engines is important. This proposal makes lifetime management issues explicit by requiring pointer or reference types in the template instantiation of variate_generator if reference semantics are desired. Without a management class, providing this features is much more cumbersome and imposes additional burden on the programmers that produce distributions. Alternatively, reference semantics could always be used, but this is an error-prone approach due to the lack of a standard reference-counted smart pointer. I believe it is impossible to add a reference-counted sharing mechanism in a manager-class-free design without giving its precise type. And that would certainly conflict in semantic scope with a smart pointer that will get into the standard eventually.

The management layer allows for a few common features to be factored out, in particular access to the engine and some member typedefs.

Comparison of other differing features between manager and non-manager designs:

• When passing a variate_generator as a an argument to a function, the design in this proposal allows selecting only those function signatures during overload resolution that are intended to be called with a variate_generator. In contrast, misbehaviour is possible without a manager class, similar to iterators in function signatures: template<class BiIter> void f(BiIter it) matches f(5), without regard to the bidirectional iterator requirements. An error then happens when instantiating f. The situation worsens when several competing function templates are available and the wrong one is chosen accidentally.
• With the engine passed into the distribution's constructor, the full type hierarchy of engine (and possibly adaptor) are available to the distribution, allowing to cache information derived from the engine (e.g. its value range) . Also, (partial) specialization of distributions could be written that optimize the interaction with a specific engine and/or adaptor. In this proposal's design, this information is available in the variate_generator template only. However, optimizations by specialization for the engine/adaptor combination are perceived to possibly have high benefit, while those for the (engine+adaptor) / distribution combination are presumed to be much less beneficial.
• Having distribution classes directly exposed to the client easily allows that the user writes a distribution with an additional arbitrary member function declaration, intended to produce random numbers while taking additional parameters into account. In this proposal's design, this is possible by using the generic operator()(T x) forwarding function.

R. Add-on packages

• unique seeding: Make it easy for the user to provide a unique seed for each instance of a pseudo-random number engine. Design idea:

    class unique_seed;
  
    template<class Engine>
    Engine seed(unique_seed&);
  

  The "seed" function retrieves some unique seed from the unique_seed class and then uses the seed(first, last) interface of an engine to implant that unique seed. Specific seeding requirements for some engine can be met by (partial) template specialization.
• runtime-replaceable distributions and associated save/restore functionality: Provide a class hierarchy that invokes distributions by means of a virtual function, thereby allowing for runtime replacement of the actual distribution. Provide a factory function to reconstruct the distribution instance after saving it to some non-volatile media.

S. Adaptors

Z. Open issues

• Some engines require non-negative template arguments, usually bit counts. Should these be given as "int" or "unsigned int"? Using "unsigned int" sometimes adds significant clutter to the presentation. Or "size_t", but this is probably too large a type?

IV. Proposed Text

This subclause defines a facility for generating random numbers.

Random number requirements

In the following table, X denotes a number generator class returning objects of type T, and u is a (possibly const) value of X.

Number generator requirements (in addition to function object)
expression	return type	pre/post-condition	complexity
X::result_type	T	std::numeric_limits<T>::is_specialized is true	compile-time

In the following table, X denotes a uniform random number generator class returning objects of type T, u is a value of X and v is a (possibly const) value of X.

Uniform random number generator requirements (in addition to number generator)
expression	return type	pre/post-condition	complexity
u()	T	-	amortized constant
v.min()	T	Returns some l where l is less than or equal to all values potentially returned by operator(). The return value of this function shall not change during the lifetime of v.	constant
v.max()	T	If std::numeric_limits<T>::is_integer, returns l where l is less than or equal to all values potentially returned by operator(), otherwise, returns l where l is strictly less than all values potentially returned by operator(). In any case, the return value of this function shall not change during the lifetime of v.	constant

In the following table, X denotes a pseudo-random number engine class returning objects of type T, t is a value of T, u is a value of X, v is an lvalue of X, it1 is an lvalue and it2 is a (possibly const) value of an input iterator type It having an unsigned integral value type, x, y are (possibly const) values of X, os is convertible to an lvalue of type std::ostream, and is is convertible to an lvalue of type std::istream.

A pseudo-random number engine x has a state x(i) at any given time. The specification of each pseudo-random number engines defines the size of its state in multiples of the size of its result_type, given as an integral constant expression.

Pseudo-random number engine requirements (in addition to uniform random number generator, CopyConstructible, and Assignable)
expression	return type	pre/post-condition	complexity
X()	-	creates an engine with the same initial state as all other default-constructed engines of type X in the program.	O(size of state)
X(it1, it2)	-	creates an engine with the initial state given by the range [it1,it2). it1 is advanced by the size of state. If the size of the range [it1,it2) is insufficient, leaves it1 == it2 and throws invalid_argument.	O(size of state)
u.seed()	void	post: u == X()	O(size of state)
u.seed(it1, it2)	void	post: If there are sufficiently many values in [it1, it2) to initialize the state of u, then u == X(it1,it2). Otherwise, it1 == it2, throws invalid_argument, and further use of u (except destruction) is undefined until a seed member function has been executed without throwing an exception.	O(size of state)
u()	T	given the state u(i) of the engine, computes u(i+1), sets the state to u(i+1), and returns some output dependent on u(i+1)	amortized constant
x == y	bool	== is an equivalence relation. The current state x(i) of x is equal to the current state y(j) of y.	O(size of state)
x != y	bool	!(x == y)	O(size of state)
os << x	std::ostream&	writes the textual representation of the state x(i) of x to os, with os.fmtflags set to ios_base::dec|ios_base::fixed|ios_base::left and the fill character set to the space character. In the output, adjacent numbers are separated by one or more space characters. post: The os.fmtflags and fill character are unchanged.	O(size of state)
is >> v	std::istream&	sets the state v(i) of v as determined by reading its textual representation from is. post: The is.fmtflags are unchanged.	O(size of state)

In the following table, X denotes a random distribution class returning objects of type T, u is a value of X, x is a (possibly const) value of X, and e is an lvalue of an arbitrary type that meets the requirements of a uniform random number generator, returning values of type U.

Random distribution requirements (in addition to number generator, CopyConstructible, and Assignable)
expression	return type	pre/post-condition	complexity
X::input_type	U	-	compile-time
u.reset()	void	subsequent uses of u do not depend on values produced by e prior to invoking reset.	constant
u(e)	T	the sequence of numbers returned by successive invocations with the same object e is randomly distributed with some probability density function p(x)	amortized constant number of invocations of e
os << x	std::ostream&	writes a textual representation for the parameters and additional internal data of the distribution x to os. post: The os.fmtflags and fill character are unchanged.	O(size of state)
is >> u	std::istream&	restores the parameters and additional internal data of the distribution u. pre: is provides a textual representation that was previously written by operator<< post: The is.fmtflags are unchanged.	O(size of state)

Additional requirements: The sequence of numbers produced by repeated invocations of x(e) does not change whether or not os << x is invoked between any of the invocations x(e). If a textual representation is written using os << x and that representation is restored into the same or a different object y of the same type using is >> y, repeated invocations of y(e) produce the same sequence of random numbers as would repeated invocations of x(e).

In the following subclauses, a template parameter named UniformRandomNumberGenerator shall denote a class that satisfies all the requirements of a uniform random number generator. Moreover, a template parameter named Distribution shall denote a type that satisfies all the requirements of a random distribution. Furthermore, a template parameter named RealType shall denote a type that holds an approximation to a real number. This type shall meet the requirements for a numeric type (26.1 [lib.numeric.requirements]), the binary operators +, -, *, / shall be applicable to it, a conversion from double shall exist, and function signatures corresponding to those for type double in subclause 26.5 [lib.c.math] shall be available by argument-dependent lookup (3.4.2 [basic.lookup.koenig]). [Note: The built-in floating-point types float and double meet these requirements.]

Header <random> synopsis

namespace std {
  template<class UniformRandomNumberGenerator, class Distribution>
  class variate_generator;

  template<class IntType, IntType a, IntType c, IntType m>
  class linear_congruential;

  template<class UIntType, int w, int n, int m, int r, UIntType a, int u,
  int s, UIntType b, int t, UIntType c, int l>
  class mersenne_twister;

  template<class IntType, IntType m, int s, int r>
  class subtract_with_carry;

  template<class RealType, int w, int s, int r>
  class subtract_with_carry_01;

  template<class UniformRandomNumberGenerator, int p, int r>
  class discard_block;

  template<class UniformRandomNumberGenerator1, int s1,
           class UniformRandomNumberGenerator2, int s2>
  class xor_combine;

  class random_device;

  template<class IntType = int>
  class uniform_int;

  template<class RealType = double>
  class bernoulli_distribution;

  template<class IntType = int, class RealType = double>
  class geometric_distribution;

  template<class IntType = int, class RealType = double>
  class poisson_distribution;

  template<class IntType = int, class RealType = double>
  class binomial_distribution;

  template<class RealType = double>
  class uniform_real;

  template<class RealType = double>
  class exponential_distribution;

  template<class RealType = double>
  class normal_distribution;

  template<class RealType = double>
  class gamma_distribution;

} // namespace std

Class template variate_generator

template<class Engine, class Distribution>
class variate_generator
{
public:
  typedef Engine engine_type;
  typedef /* implementation defined */ engine_value_type;
  typedef Distribution distribution_type;
  typedef typename Distribution::result_type result_type;

  variate_generator(engine_type eng, distribution_type d);

  result_type operator()();
  template<class T>  result_type operator()(T value);

  engine_value_type& engine();
  const engine_value_type& engine() const;

  distribution_type& distribution();
  const distribution_type& distribution() const;

  result_type min() const;
  result_type max() const;
};

Specializations of variate_generator satisfy the requirements of CopyConstructible. They also satisfy the requirements of Assignable unless the template parameter Engine is of the form U&.

The complexity of all functions specified in this section is constant. No function described in this section except the constructor throws an exception.

    variate_generator(engine_type eng, distribution_type d)

    result_type operator()()

    template<class T> result_type operator()(T value)

    engine_value_type& engine()

    const engine_value_type& engine() const

    distribution_type& distribution()

    const distribution_type& distribution() const

    result_type min() const

    result_type max() const

Random number engine class templates

The class templates specified in this section satisfy all the requirements of a pseudo-random number engine (given in tables in section x.x), except where specified otherwise. Descriptions are provided here only for operations on the engines that are not described in one of these tables or for operations where there is additional semantic information.

All members declared static const in any of the following class templates shall be defined in such a way that they are usable as integral constant expressions.

Class template linear_congruential

namespace std {
  template<class IntType, IntType a, IntType c, IntType m>
  class linear_congruential
  {
  public:
    // types
    typedef IntType result_type;

    // parameter values
    static const IntType multiplier = a;
    static const IntType increment = c;
    static const IntType modulus = m;

    //  constructors and member function
    explicit linear_congruential(IntType x0 = 1);
    template<class In> linear_congruential(In& first, In last);
    void seed(IntType x0 = 1);
    template<class In> void seed(In& first, In last);
    result_type min() const;
    result_type max() const;
    result_type operator()();
  };

  template<class IntType, IntType a, IntType c, IntType m>
  bool operator==(const linear_congruential<IntType, a, c, m>& x,
                  const linear_congruential<IntType, a, c, m>& y);
  template<class IntType, IntType a, IntType c, IntType m>
  bool operator!=(const linear_congruential<IntType, a, c, m>& x,
                  const linear_congruential<IntType, a, c, m>& y);

  template<class CharT, class traits,
           class IntType, IntType a, IntType c, IntType m>
  basic_ostream<CharT, traits>& operator<<(basic_ostream<CharT, traits>& os,
                                           const linear_congruential<IntType, a, c, m>& x);  
  template<class CharT, class traits,
           class IntType, IntType a, IntType c, IntType m>
  basic_istream<CharT, traits>& operator>>(basic_istream<CharT, traits>& is, 
                                           linear_congruential<IntType, a, c, m>& x);
}

The size of the state x(i) is 1.

    explicit linear_congruential(IntType x0 = 1)

    void seed(IntType x0 = 1)

    template<class In> linear_congruential(In& first, In last)

  template<class CharT, class traits,
           class IntType, IntType a, IntType c, IntType m>
  basic_ostream<CharT, traits>& operator<<(basic_ostream<CharT, traits>& os,
                                           const linear_congruential<IntType, a, c, m>& x);  

Class template mersenne_twister

• y(i) = (x(i-n) bitand um) | (x(i-(n-1)) bitand lm)
• If the lowest bit of the binary representation of y(i) is set, x(i) = x(i-(n-m)) xor (y(i) rshift 1) xor a; otherwise x(i) = x(i-(n-m)) xor (y(i) rshift 1).
• z1(i) = x(i) xor ( x(i) rshift u )
• z2(i) = z1(i) xor ( (z1(i) lshift s) bitand b )
• z3(i) = z2(i) xor ( (z2(i) lshift t) bitand c )
• o(x(i)) = z3(i) xor ( z3(i) rshift l )

namespace std {
  template<class UIntType, int w, int n, int m, int r, UIntType a, int u,
  int s, UIntType b, int t, UIntType c, int l>
  class mersenne_twister
  {
  public:
    // types
    typedef UIntType result_type;

    // parameter values
    static const int word_size = w;
    static const int state_size = n;
    static const int shift_size = m;
    static const int mask_bits = r;
    static const UIntType parameter_a = a;
    static const int output_u = u;
    static const int output_s = s;
    static const UIntType output_b = b;
    static const int output_t = t;
    static const UIntType output_c = c;
    static const int output_l = l;

    //  constructors and member function
    mersenne_twister();
    explicit mersenne_twister(UIntType value);
    template<class In> mersenne_twister(In& first, In last);
    void seed();
    void seed(UIntType value);
    template<class In> void seed(In& first, In last);
    result_type min() const;
    result_type max() const;
    result_type operator()();
  };

  template<class UIntType, int w, int n, int m, int r, UIntType a, int u,
           int s, UIntType b, int t, UIntType c, int l>
  bool operator==(const mersenne_twister<UIntType, w, n, m, r, a, u, s, b, t, c, l>& y,
                  const mersenne_twister<UIntType, w, n, m, r, a, u, s, b, t, c, l>& x);
  template<class UIntType, int w, int n, int m, int r, UIntType a, int u,
           int s, UIntType b, int t, UIntType c, int l>
  bool operator!=(const mersenne_twister<UIntType, w, n, m, r, a, u, s, b, t, c, l>& y,
                  const mersenne_twister<UIntType, w, n, m, r, a, u, s, b, t, c, l>& x);

  template<class CharT, class traits,
           class UIntType, int w, int n, int m, int r, UIntType a, int u,
           int s, UIntType b, int t, UIntType c, int l>
  basic_ostream<CharT, traits>& operator<<(basic_ostream<CharT, traits>& os,
                                           const mersenne_twister<UIntType, w, n, m, r, a, u, s, b, t, c, l>& x);
  template<class CharT, class traits,
           class UIntType, int w, int n, int m, int r, UIntType a, int u,
           int s, UIntType b, int t, UIntType c, int l>
  basic_istream<CharT, traits>& operator>>(basic_istream<CharT, traits>& is, 
                                           mersenne_twister<UIntType, w, n, m, r, a, u, s, b, t, c, l>& x);
}

The size of the state x(i) is n.

    mersenne_twister()

    explicit mersenne_twister(result_type value)

    template<class In> mersenne_twister(In& first, In last)

    void seed()

    void seed(result_type value)

    template<class In> void seed(In& first, In last)

    template<class UIntType, int w, int n, int m, int r, UIntType a, int u,
             int s, UIntType b, int t, UIntType c, int l>
    bool operator==(const mersenne_twister<UIntType, w, n, m, r, a, u, s, b, t, c, l>& y,
                    const mersenne_twister<UIntType, w, n, m, r, a, u, s, b, t, c, l>& x)

    template<class CharT, class traits,
             class UIntType, int w, int n, int m, int r, UIntType a, int u,
             int s, UIntType b, int t, UIntType c, int l>
    basic_ostream<CharT, traits>& operator<<(basic_ostream<CharT, traits>& os,
                                             const mersenne_twister<UIntType, w, n, m, r, a, u, s, b, t, c, l>& x)

Class template subtract_with_carry

namespace std {
  template<class IntType, IntType m, int s, int r>
  class subtract_with_carry
  {
  public:
    // types
    typedef IntType result_type;

    // parameter values
    static const IntType modulus = m;
    static const int long_lag = r;
    static const int short_lag = s;

    //  constructors and member function
    subtract_with_carry();
    explicit subtract_with_carry(IntType value);
    template<class In> subtract_with_carry(In& first, In last);
    void seed(IntType value = 19780503);
    template<class In> void seed(In& first, In last);
    result_type min() const;
    result_type max() const;
    result_type operator()();
  };
  template<class IntType, IntType m, int s, int r>
  bool operator==(const subtract_with_carry<IntType, m, s, r> & x,
                  const subtract_with_carry<IntType, m, s, r> & y);

  template<class IntType, IntType m, int s, int r>
  bool operator!=(const subtract_with_carry<IntType, m, s, r> & x,
                  const subtract_with_carry<IntType, m, s, r> & y);

  template<class CharT, class Traits,
           class IntType, IntType m, int s, int r>
  std::basic_ostream<CharT,Traits>& operator<<(std::basic_ostream<CharT,Traits>& os,
                                               const subtract_with_carry<IntType, m, s, r>& f);

  template<class CharT, class Traits,
          class IntType, IntType m, int s, int r>
  std::basic_istream<CharT,Traits>& operator>>(std::basic_istream<CharT,Traits>& is, 
                                               subtract_with_carry<IntType, m, s, r>& f);
}

The size of the state is r.

    subtract_with_carry()

    explicit subtract_with_carry(IntType value)

    template<class In> subtract_with_carry(In& first, In last)

    void seed(IntType value = 19780503)

    template<class In> void seed(In& first, In last)

    template<class IntType, IntType m, int s, int r>
    bool operator==(const subtract_with_carry<IntType, m, s, r> & x,
                    const subtract_with_carry<IntType, m, s, r> & y)

    template<class CharT, class Traits,
          class IntType, IntType m, int s, int r>
    std::basic_ostream<CharT,Traits>& operator<<(std::basic_ostream<CharT,Traits>& os,
                                                 const subtract_with_carry<IntType, m, s, r>& f)

Class template subtract_with_carry_01

namespace std {
  template<class RealType, int w, int s, int r>
  class subtract_with_carry_01
  {
  public:
    // types
    typedef RealType result_type;

    // parameter values
    static const int word_size = w;
    static const int long_lag = r;
    static const int short_lag = s;

    //  constructors and member function
    subtract_with_carry_01();
    explicit subtract_with_carry_01(unsigned int value);
    template<class In> subtract_with_carry_01(In& first, In last);
    void seed(unsigned int value = 19780503);
    template<class In> void seed(In& first, In last);
    result_type min() const;
    result_type max() const;
    result_type operator()();
  };
  template<class RealType, int w, int s, int r>
  bool operator==(const subtract_with_carry_01<RealType, w, s, r> x,
                  const subtract_with_carry_01<RealType, w, s, r> y);

  template<class RealType, int w, int s, int r>
  bool operator!=(const subtract_with_carry_01<RealType, w, s, r> x,
                  const subtract_with_carry_01<RealType, w, s, r> y);

  template<class CharT, class Traits,
           class RealType, int w, int s, int r>
  std::basic_ostream<CharT,Traits>& operator<<(std::basic_ostream<CharT,Traits>& os,
                                               const subtract_with_carry_01<RealType, w, s, r>& f);

  template<class CharT, class Traits,
           class RealType, int w, int s, int r>
  std::basic_istream<CharT,Traits>& operator>>(std::basic_istream<CharT,Traits>& is, 
                                               subtract_with_carry_01<RealType, w, s, r>& f);
}

The size of the state is r.

    subtract_with_carry_01()

    explicit subtract_with_carry_01(unsigned int value)

    template<class In> subtract_with_carry_01(In& first, In last)

    void seed(unsigned int value = 19780503)

    template<class In> void seed(In& first, In last)

    template<class RealType, int w, int s, int r>
    bool operator==(const subtract_with_carry<RealType, w, s, r> x,
                    const subtract_with_carry<RealType, w, s, r> y);

    template<class CharT, class Traits,
             class RealType, int w, int s, int r>
    std::basic_ostream<CharT,Traits>& operator<<(std::basic_ostream<CharT,Traits>& os,
                                                 const subtract_with_carry<RealType, w, s, r>& f);

Class template discard_block

namespace std {
  template<class UniformRandomNumberGenerator, int p, int r>
  class discard_block
  {
  public:
    // types
    typedef UniformRandomNumberGenerator base_type;
    typedef typename base_type::result_type result_type;
  
    // parameter values
    static const int block_size = p;
    static const int used_block = r;
  
    //  constructors and member function
    discard_block();
    explicit discard_block(const base_type & rng);
    template<class In> discard_block(In& first, In last);
    void seed();
    template<class In> void seed(In& first, In last);
    const base_type& base() const;
    result_type min() const;
    result_type max() const;
    result_type operator()();  
  private:
    // base_type b;                 exposition only
    // int n;                       exposition only
  };
  template<class UniformRandomNumberGenerator, int p, int r>
  bool operator==(const discard_block<UniformRandomNumberGenerator,p,r> & x,
                 (const discard_block<UniformRandomNumberGenerator,p,r> & y);
  template<class UniformRandomNumberGenerator, int p, int r,
    typename UniformRandomNumberGenerator::result_type val>
  bool operator!=(const discard_block<UniformRandomNumberGenerator,p,r> & x,
                 (const discard_block<UniformRandomNumberGenerator,p,r> & y);

  template<class CharT, class traits,
           class UniformRandomNumberGenerator, int p, int r>
  basic_ostream<CharT, traits>& operator<<(basic_ostream<CharT, traits>& os,
                                           const discard_block<UniformRandomNumberGenerator,p,r> & x);
  template<class CharT, class traits,
           class UniformRandomNumberGenerator, int p, int r>
  basic_istream<CharT, traits>& operator>>(basic_istream<CharT, traits>& is, 
                                           discard_block<UniformRandomNumberGenerator,p,r> & x);

}

    discard_block()

    explicit discard_block(const base_type & rng)

    template<class In> discard_block(In& first, In last)

    void seed()

    template<class In> void seed(In& first, In last)

    const base_type& base() const

    result_type operator()()

  template<class CharT, class traits,
           class UniformRandomNumberGenerator, int p, int r>
  basic_ostream<CharT, traits>& operator<<(basic_ostream<CharT, traits>& os,
                                           const discard_block<UniformRandomNumberGenerator,p,r> & x);

Class template xor_combine

namespace std {
  template<class UniformRandomNumberGenerator1, int s1,
           class UniformRandomNumberGenerator2, int s2>
  class xor_combine
  {
  public:
    // types
    typedef UniformRandomNumberGenerator1 base1_type;
    typedef UniformRandomNumberGenerator2 base2_type;
    typedef typename base_type::result_type result_type;
  
    // parameter values
    static const int shift1 = s1;
    static const int shift2 = s2;
  
    //  constructors and member function
    xor_combine();
    xor_combine(const base1_type & rng1, const base2_type & rng2);
    template<class In> xor_combine(In& first, In last);
    void seed();
    template<class In> void seed(In& first, In last);
    const base1_type& base1() const;
    const base2_type& base2() const;
    result_type min() const;
    result_type max() const;
    result_type operator()();  
  private:
    // base1_type b1;               exposition only
    // base2_type b2;               exposition only
  };
  template<class UniformRandomNumberGenerator1, int s1,
           class UniformRandomNumberGenerator2, int s2>
  bool operator==(const xor_combine<UniformRandomNumberGenerator1, s1, 
                                    UniformRandomNumberGenerator2, s2> & x,
                 (const xor_combine<UniformRandomNumberGenerator1, s1,
                                    UniformRandomNumberGenerator2, s2> & y);
  template<class UniformRandomNumberGenerator1, int s1,
           class UniformRandomNumberGenerator2, int s2>
  bool operator!=(const xor_combine<UniformRandomNumberGenerator1, s1,
                                    UniformRandomNumberGenerator2, s2> & x,
                 (const xor_combine<UniformRandomNumberGenerator1, s1,
                                    UniformRandomNumberGenerator2, s2> & y);

  template<class CharT, class traits,
           class UniformRandomNumberGenerator1, int s1,
           class UniformRandomNumberGenerator2, int s2>
  basic_ostream<CharT, traits>& operator<<(basic_ostream<CharT, traits>& os,
                                           const xor_combine<UniformRandomNumberGenerator1, s1,
                                                             UniformRandomNumberGenerator2, s2> & x);
  template<class CharT, class traits,
           class UniformRandomNumberGenerator1, int s1,
           class UniformRandomNumberGenerator2, int s2>
  basic_istream<CharT, traits>& operator>>(basic_istream<CharT, traits>& is, 
                                           xor_combine<UniformRandomNumberGenerator1, s1,
                                                       UniformRandomNumberGenerator2, s2> & x);

}

    xor_combine()

    xor_combine(const base1_type & rng1, const base2_type & rng2)

    template<class In> xor_combine(In& first, In last)

    void seed()

    template<class In> void seed(In& first, In last)

    const base1_type& base1() const

    const base2_type& base2() const

    result_type operator()()

  template<class CharT, class traits,
           class UniformRandomNumberGenerator1, int s1,
           class UniformRandomNumberGenerator2, int s2>
  basic_ostream<CharT, traits>& operator<<(basic_ostream<CharT, traits>& os,
                                           const xor_combine<UniformRandomNumberGenerator1, s1,
                                                             UniformRandomNumberGenerator2, s2> & x);

Engines with predefined parameters

namespace std {
  typedef linear_congruential</* implementation defined */, 16807, 0, 2147483647> minstd_rand0;
  typedef linear_congruential</* implementation defined */, 48271, 0, 2147483647> minstd_rand;

  typedef mersenne_twister</* implementation defined */,32,624,397,31,0x9908b0df,11,7,0x9d2c5680,15,0xefc60000,18> mt19937;

  typedef subtract_with_carry_01<float, 24, 10, 24> ranlux_base_01;
  typedef subtract_with_carry_01<double, 48, 10, 24> ranlux64_base_01;

  typedef discard_block<subtract_with_carry</* implementation defined */, (1<<24), 10, 24>, 223, 24> ranlux3;
  typedef discard_block<subtract_with_carry</* implementation defined */, (1<<24), 10, 24>, 389, 24> ranlux4;

  typedef discard_block<subtract_with_carry_01<float, 24, 10, 24>, 223, 24> ranlux3_01;
  typedef discard_block<subtract_with_carry_01<float, 24, 10, 24>, 389, 24> ranlux4_01;
}

For a default-constructed mt19937 object, x(10000) = 3346425566.

For a default-constructed ranlux3 object, x(10000) = 5957620. For a default-constructed ranlux4 object, x(10000) = 8587295. For a default-constructed ranlux3_01 object, x(10000) = 5957620 * 2^-24. For a default-constructed ranlux4_01 object, x(10000) = 8587295 * 2^-24.

Class random_device

If implementation limitations prevent generating non-deterministic random numbers, the implementation can employ a pseudo-random number engine.

namespace std {
  class random_device
  {
  public:
    // types
    typedef unsigned int result_type;

    // constructors, destructors and member functions
    explicit random_device(const std::string& token = /* implementation-defined */);
    result_type min() const;
    result_type max() const;
    double entropy() const;
    result_type operator()();
  
  private:
    random_device(const random_device& );
    void operator=(const random_device& );
  };
}

    explicit random_device(const std::string& token = /* implementation-defined */)

    result_type min() const

    result_type max() const

    double entropy() const

    result_type operator()()

Random distribution class templates

A template parameter named IntType shall denote a type that represents an integer number. This type shall meet the requirements for a numeric type (26.1 [lib.numeric.requirements]), the binary operators +, -, *, /, % shall be applicable to it, and a conversion from int shall exist. [Footnote: The built-in types int and long meet these requirements.]

Given an object whose type is specified in this subclause, if the lifetime of the uniform random number generator referred to in the constructor invocation for that object has ended, any use of that object is undefined.

No function described in this section throws an exception, unless an operation on values of IntType or RealType throws an exception. [Note: Then, the effects are undefined, see [lib.numeric.requirements]. ]

The algorithms for producing each of the specified distributions are implementation-defined.

Class template uniform_int

A uniform_int random distribution satisfies all the requirements of a uniform random number generator (given in tables in section x.x).

namespace std {
  template<class IntType = int>
  class uniform_int
  {
  public:
    // types
    typedef IntType input_type;
    typedef IntType result_type;

    //  constructors and member function
    explicit uniform_int(IntType min = 0, IntType max = 9);
    result_type min() const;
    result_type max() const;
    void reset();
    template<class UniformRandomNumberGenerator>
    result_type operator()(UniformRandomNumberGenerator& urng);
    template<class UniformRandomNumberGenerator>
    result_type operator()(UniformRandomNumberGenerator& urng, result_type n);
  };
}

    uniform_int(IntType min = 0, IntType max = 9)

    result_type min() const

    result_type max() const

    result_type operator()(UniformRandomNumberGenerator& urng, result_type n)

Class template bernoulli_distribution

namespace std {
  template<class RealType = double>
  class bernoulli_distribution
  {
  public:
    // types
    typedef int input_type;
    typedef bool result_type;

    //  constructors and member function
    explicit bernoulli_distribution(const RealType& p = RealType(0.5));
    RealType p() const;
    void reset();
    template<class UniformRandomNumberGenerator>
    result_type operator()(UniformRandomNumberGenerator& urng);
  };
}

    bernoulli_distribution(const RealType& p = RealType(0.5))

    RealType p() const

Class template geometric_distribution

namespace std {
  template<class IntType = int, class RealType = double>
  class geometric_distribution
  {
  public:
    // types
    typedef RealType input_type;
    typedef IntType result_type;

    //  constructors and member function
    explicit geometric_distribution(const RealType& p = RealType(0.5));
    RealType p() const;
    void reset();
    template<class UniformRandomNumberGenerator>
    result_type operator()(UniformRandomNumberGenerator& urng);
  };
}

    geometric_distribution(const RealType& p = RealType(0.5))

   RealType p() const

Class template poisson_distribution

namespace std {
  template<class IntType = int, class RealType = double>
  class poisson_distribution
  {
  public:
    // types
    typedef RealType input_type;
    typedef IntType result_type;

    //  constructors and member function
    explicit poisson_distribution(const RealType& mean = RealType(1));
    RealType mean() const;
    void reset();
    template<class UniformRandomNumberGenerator>
    result_type operator()(UniformRandomNumberGenerator& urng);
  };
}

    poisson_distribution(const RealType& mean = RealType(1))

   RealType mean() const

Class template binomial_distribution

namespace std {
  template<class IntType = int, class RealType = double>
  class binomial_distribution
  {
  public:
    // types
    typedef /* implementation-defined */ input_type;
    typedef IntType result_type;

    //  constructors and member function
    explicit binomial_distribution(IntType t = 1, const RealType& p = RealType(0.5));
    IntType t() const;
    RealType p() const;
    void reset();
    template<class UniformRandomNumberGenerator>
    result_type operator()(UniformRandomNumberGenerator& urng);
  };
}

    binomial_distribution(IntType t = 1, const RealType& p = RealType(0.5))

   IntType t() const

   RealType p() const

Class template uniform_real

A uniform_real random distribution satisfies all the requirements of a uniform random number generator (given in tables in section x.x).

namespace std {
  template<class RealType = double>
  class uniform_real
  {
  public:
    // types
    typedef RealType input_type;
    typedef RealType result_type;

    //  constructors and member function
    explicit uniform_real(RealType min = RealType(0), RealType max = RealType(1));
    result_type min() const;
    result_type max() const;
    void reset();
    template<class UniformRandomNumberGenerator>
    result_type operator()(UniformRandomNumberGenerator& urng);
  };
}

    uniform_real(RealType min = RealType(0), RealType max = RealType(1))

    result_type min() const

    result_type max() const

Class template exponential_distribution

namespace std {
  template<class RealType = double>
  class exponential_distribution
  {
  public:
    // types
    typedef RealType input_type;
    typedef RealType result_type;

    //  constructors and member function
    explicit exponential_distribution(const result_type& lambda = result_type(1));
    RealType lambda() const;
    void reset();
    template<class UniformRandomNumberGenerator>
    result_type operator()(UniformRandomNumberGenerator& urng);
  };
}

    exponential_distribution(const result_type& lambda = result_type(1))

    RealType lambda() const

Class template normal_distribution

namespace std {
  template<class RealType = double>
  class normal_distribution
  {
  public:
    // types
    typedef RealType input_type;
    typedef RealType result_type;

    //  constructors and member function
    explicit normal_distribution(const result_type& mean = 0,
                                 const result_type& sigma = 1);
    RealType mean() const;
    RealType sigma() const;
    void reset();
    template<class UniformRandomNumberGenerator>
    result_type operator()(UniformRandomNumberGenerator& urng);
  };
}

    explicit normal_distribution(const result_type& mean = 0,
                                 const result_type& sigma = 1);

    RealType mean() const

    RealType sigma() const

Class template gamma_distribution

namespace std {
  template<class RealType = double>
  class gamma_distribution
  {
  public:
    // types
    typedef RealType input_type;
    typedef RealType result_type;

    //  constructors and member function
    explicit gamma_distribution(const result_type& alpha = result_type(1));
    RealType alpha() const;
    void reset();
    template<class UniformRandomNumberGenerator>
    result_type operator()(UniformRandomNumberGenerator& urng);
  };
}

    explicit gamma_distribution(const result_type& alpha = result_type(1));

    RealType alpha() const

V. Acknowledgements

• Thanks to Walter Brown, Mark Fischler and Marc Paterno from Fermilab for input about the requirements of high-energy physics.
• Thanks to David Abrahams for additional comments on the design.
• Thanks to the Boost community for a platform for experimentation.

VI. References

• William H. Press, Saul A. Teukolsky, William A. Vetterling, Brian P. Flannery, "Numerical Recipes in C: The art of scientific computing", 2nd ed., 1992, pp. 274-328
• Bruce Schneier, "Applied Cryptography", 2nd ed., 1996, ch. 16-17. [I haven't read this myself. Yet.]
• D. H. Lehmer, "Mathematical methods in large-scale computing units", Proc. 2nd Symposium on Large-Scale Digital Calculating Machines, Harvard University Press, 1951, pp. 141-146
• P.A. Lewis, A.S. Goodman, J.M. Miller, "A pseudo-random number generator for the System/360", IBM Systems Journal, Vol. 8, No. 2, 1969, pp. 136-146
• Stephen K. Park and Keith W. Miller, "Random Number Generators: Good ones are hard to find", Communications of the ACM, Vol. 31, No. 10, October 1988, pp. 1192-1201
• Makoto Matsumoto and Takuji Nishimura, "Mersenne Twister: A 623-dimensionally equidistributed uniform pseudo-random number generator", ACM Transactions on Modeling and Computer Simulation: Special Issue on Uniform Random Number Generation, Vol. 8, No. 1, January 1998, pp. 3-30. http://www.math.keio.ac.jp/matumoto/emt.html
• Donald E. Knuth, "The Art of Computer Programming, Vol. 2", 3rd ed., 1997, pp. 1-193.
• Carter Bays and S.D. Durham, "Improving a poor random number generator", ACM Transactions on Mathematical Software, Vol. 2, 1979, pp. 59-64.
• Martin Lüscher, "A portable high-quality random number generator for lattice field theory simulations.", Computer Physics Communications, Vol. 79, 1994, pp. 100-110.
• William J. Hurd, "Efficient Generation of Statistically Good Pseudonoise by Linearly Interconnected Shift Registers", Technical Report 32-1526, Volume XI, The Deep Space Network Progress Report for July and August 1972, NASA Jet Propulsion Laboratory, 1972 and IEEE Transactions on Computers Vol. 23, 1974.
• Pierre L'Ecuyer, "Efficient and Portable Combined Random Number Generators", Communications of the ACM, Vol. 31, pp. 742-749+774, 1988.
• Pierre L'Ecuyer, "Maximally equidistributed combined Tausworthe generators", Mathematics of Computation Vol. 65, pp. 203-213, 1996.
• Pierre L'Ecuyer, "Good parameters and implementations for combined multple recursive random number generators", Operations Research Vol. 47, pp. 159-164, 1999.
• S. Kirkpatrick and E. Stoll, "A very fast shift-register sequence random number generator", Journal of Computational Physics, Vol. 40, pp. 517-526, 1981.
• R. C. Tausworthe, "Random numbers generated by iinear recurrence modulo two", Mathematics of Computation, Vol. 19, pp. 201-209, 1965.
• George Marsaglia and Arif Zaman, "A New Class of Random Number Generators", Annals of Applied Probability, Vol. 1, No. 3, 1991.