29 Numerics library [numerics]

29.5 Random number generation [rand]

29.5.4 Random number engine class templates [rand.eng]

29.5.4.1 General [rand.eng.general]

Each type instantiated from a class template specified in [rand.eng] meets the requirements of a random number engine type.

Except where specified otherwise, the complexity of each function specified in [rand.eng] is constant.

Except where specified otherwise, no function described in [rand.eng] throws an exception.

Every function described in [rand.eng] that has a function parameter q of type Sseq& for a template type parameter named Sseq that is different from type seed_seq throws what and when the invocation of q.generate throws.

Descriptions are provided in [rand.eng] only for engine operations that are not described in [rand.req.eng] or for operations where there is additional semantic information.

In particular, declarations for copy constructors, for copy assignment operators, for streaming operators, and for equality and inequality operators are not shown in the synopses.

Each template specified in [rand.eng] requires one or more relationships, involving the value(s) of its non-type template parameter(s), to hold.

A program instantiating any of these templates is ill-formed if any such required relationship fails to hold.

For every random number engine and for every random number engine adaptor X defined in [rand.eng] and in [rand.adapt]:

(7.1)
if the constructor template<class Sseq> explicit X(Sseq& q); is called with a type Sseq that does not qualify as a seed sequence, then this constructor shall not participate in overload resolution;
(7.2)
if the member function template<class Sseq> void seed(Sseq& q); is called with a type Sseq that does not qualify as a seed sequence, then this function shall not participate in overload resolution.

The extent to which an implementation determines that a type cannot be a seed sequence is unspecified, except that as a minimum a type shall not qualify as a seed sequence if it is implicitly convertible to X::result_type.

29.5.4.2 Class template linear_congruential_engine [rand.eng.lcong]

A linear_congruential_engine random number engine produces unsigned integer random numbers.

The state x

_{i}

of a linear_congruential_engine object x is of size 1 and consists of a single integer.

The transition algorithm is a modular linear function of the form

T A (x_{i}) = (a \cdot x_{i} + c) mod m

; the generation algorithm is

G A (x_{i}) = x_{i + 1}

🔗

namespace std { template<class UIntType, UIntType a, UIntType c, UIntType m> class linear_congruential_engine { public: // types using result_type = UIntType; // engine characteristics static constexpr result_type multiplier = a; static constexpr result_type increment = c; static constexpr result_type modulus = m; static constexpr result_type min() { return c == 0u ? 1u: 0u; } static constexpr result_type max() { return m - 1u; } static constexpr result_type default_seed = 1u; // constructors and seeding functions linear_congruential_engine() : linear_congruential_engine(default_seed) {} explicit linear_congruential_engine(result_type s); template<class Sseq> explicit linear_congruential_engine(Sseq& q); void seed(result_type s = default_seed); template<class Sseq> void seed(Sseq& q); // equality operators friend bool operator==(const linear_congruential_engine& x, const linear_congruential_engine& y); // generating functions result_type operator()(); void discard(unsigned long long z); // inserters and extractors template<class charT, class traits> friend basic_ostream<charT, traits>& operator<<(basic_ostream<charT, traits>& os, const linear_congruential_engine& x); template<class charT, class traits> friend basic_istream<charT, traits>& operator>>(basic_istream<charT, traits>& is, linear_congruential_engine& x); }; }

If the template parameter m is 0, the modulus m used throughout [rand.eng.lcong] is numeric_limits<result_type>::max() plus 1.

[Note 1:

m need not be representable as a value of type result_type.

— end note]

If the template parameter m is not 0, the following relations shall hold: a < m and c < m.

The textual representation consists of the value of x

_{i}

🔗

explicit linear_congruential_engine(result_type s);

Effects: If

c mod m

is 0 and

s mod m

is 0, sets the engine's state to 1, otherwise sets the engine's state to

s mod m

🔗

template<class Sseq> explicit linear_congruential_engine(Sseq& q);

Effects: With

k = ⌈ \frac{log_{2} m}{32} ⌉

and a an array (or equivalent) of length

k + 3

, invokes q.generate(

a + 0

a + k + 3

) and then computes

S = (\sum_{j = 0}^{k - 1} a_{j + 3} \cdot 2^{32 j}) mod m

c mod m

is 0 and S is 0, sets the engine's state to 1, else sets the engine's state to S.

29.5.4.3 Class template mersenne_twister_engine [rand.eng.mers]

A mersenne_twister_engine random number engine248 produces unsigned integer random numbers in the closed interval

[0, 2^{w} - 1]

The state x

_{i}

of a mersenne_twister_engine object x is of size n and consists of a sequence X of n values of the type delivered by x; all subscripts applied to X are to be taken modulo n.

The transition algorithm employs a twisted generalized feedback shift register defined by shift values n and m, a twist value r, and a conditional xor-mask a.

To improve the uniformity of the result, the bits of the raw shift register are additionally tempered (i.e., scrambled) according to a bit-scrambling matrix defined by values u, d, s, b, t, c, and ℓ.

The state transition is performed as follows:

(2.1)
Concatenate the upper $w - r$ bits of $X_{i - n}$ with the lower r bits of $X_{i + 1 - n}$ to obtain an unsigned integer value Y.
(2.2)
With $α = a \cdot (Y b i t a n d 1)$ , set $X_{i}$ to $X_{i + m - n} x o r (Y r s h i f t 1) x o r α$ .

The sequence X is initialized with the help of an initialization multiplier f.

The generation algorithm determines the unsigned integer values

z_{1}, z_{2}, z_{3}, z_{4}

as follows, then delivers

z_{4}

as its result:

(3.1)
Let $z_{1} = X_{i} x o r ((X_{i} r s h i f t u) b i t a n d d)$ .
(3.2)
Let $z_{2} = z_{1} x o r ((z_{1} {l s h i f t}_{w} s) b i t a n d b)$ .
(3.3)
Let $z_{3} = z_{2} x o r ((z_{2} {l s h i f t}_{w} t) b i t a n d c)$ .
(3.4)
Let $z_{4} = z_{3} x o r (z_{3} r s h i f t ℓ)$ .

🔗

namespace std { template<class UIntType, size_t w, size_t n, size_t m, size_t r, UIntType a, size_t u, UIntType d, size_t s, UIntType b, size_t t, UIntType c, size_t l, UIntType f> class mersenne_twister_engine { public: // types using result_type = UIntType; // engine characteristics static constexpr size_t word_size = w; static constexpr size_t state_size = n; static constexpr size_t shift_size = m; static constexpr size_t mask_bits = r; static constexpr UIntType xor_mask = a; static constexpr size_t tempering_u = u; static constexpr UIntType tempering_d = d; static constexpr size_t tempering_s = s; static constexpr UIntType tempering_b = b; static constexpr size_t tempering_t = t; static constexpr UIntType tempering_c = c; static constexpr size_t tempering_l = l; static constexpr UIntType initialization_multiplier = f; static constexpr result_type min() { return 0; } static constexpr result_type max() { return

2^{w} - 1

; } static constexpr result_type default_seed = 5489u; // constructors and seeding functions mersenne_twister_engine() : mersenne_twister_engine(default_seed) {} explicit mersenne_twister_engine(result_type value); template<class Sseq> explicit mersenne_twister_engine(Sseq& q); void seed(result_type value = default_seed); template<class Sseq> void seed(Sseq& q); // equality operators friend bool operator==(const mersenne_twister_engine& x, const mersenne_twister_engine& y); // generating functions result_type operator()(); void discard(unsigned long long z); // inserters and extractors template<class charT, class traits> friend basic_ostream<charT, traits>& operator<<(basic_ostream<charT, traits>& os, const mersenne_twister_engine& x); template<class charT, class traits> friend basic_istream<charT, traits>& operator>>(basic_istream<charT, traits>& is, mersenne_twister_engine& x); }; }

The following relations shall hold: 0 < m, m <= n, 2u < w, r <= w, u <= w, s <= w, t <= w, l <= w, w <= numeric_limits<UIntType>::digits, a <= (1u<<w) - 1u, b <= (1u<<w) - 1u, c <= (1u<<w) - 1u, d <= (1u<<w) - 1u, and f <= (1u<<w) - 1u.

The textual representation of x

_{i}

consists of the values of

X_{i - n}, \dots, X_{i - 1}

, in that order.

🔗

explicit mersenne_twister_engine(result_type value);

Effects: Sets

X_{- n}

value mod 2^{w}

Then, iteratively for

i = 1 - n, \dots, - 1

, sets

X_{i}

[f \cdot (X_{i - 1} x o r (X_{i - 1} r s h i f t (w - 2))) + i mod n] mod 2^{w} .

Complexity:

O (n)

🔗

template<class Sseq> explicit mersenne_twister_engine(Sseq& q);

Effects: With

k = ⌈ w / 32 ⌉

and a an array (or equivalent) of length

n \cdot k

, invokes q.generate(

a + 0

a + n \cdot k

) and then, iteratively for

i = - n, \dots, - 1

, sets

X_{i}

(\sum_{j = 0}^{k - 1} a_{k (i + n) + j} \cdot 2^{32 j}) mod 2^{w}

Finally, if the most significant

w - r

bits of

X_{- n}

are zero, and if each of the other resulting

X_{i}

is 0, changes

X_{- n}

2^{w - 1}

248)248)

The name of this engine refers, in part, to a property of its period: For properly-selected values of the parameters, the period is closely related to a large Mersenne prime number.

29.5.4.4 Class template subtract_with_carry_engine [rand.eng.sub]

A subtract_with_carry_engine random number engine produces unsigned integer random numbers.

The state x

_{i}

of a subtract_with_carry_engine object x is of size

O (r)

, and consists of a sequence X of r integer values

0 \leq X_{i} < m = 2^{w}

; all subscripts applied to X are to be taken modulo r.

The state x

_{i}

additionally consists of an integer c (known as the carry) whose value is either 0 or 1.

The state transition is performed as follows:

(3.1)
Let $Y = X_{i - s} - X_{i - r} - c$ .
(3.2)
Set $X_{i}$ to $y = Y mod m$ .

Set c to 1 if $Y < 0$ , otherwise set c to 0.

[Note 1:

This algorithm corresponds to a modular linear function of the form

T A (x_{i}) = (a \cdot x_{i}) mod b

, where b is of the form

m^{r} - m^{s} + 1

and

a = b - (b - 1) / m

— end note]

The generation algorithm is given by

G A (x_{i}) = y

, where y is the value produced as a result of advancing the engine's state as described above.

🔗

namespace std { template<class UIntType, size_t w, size_t s, size_t r> class subtract_with_carry_engine { public: // types using result_type = UIntType; // engine characteristics static constexpr size_t word_size = w; static constexpr size_t short_lag = s; static constexpr size_t long_lag = r; static constexpr result_type min() { return 0; } static constexpr result_type max() { return

m - 1

; } static constexpr uint_least32_t default_seed = 19780503u; // constructors and seeding functions subtract_with_carry_engine() : subtract_with_carry_engine(0u) {} explicit subtract_with_carry_engine(result_type value); template<class Sseq> explicit subtract_with_carry_engine(Sseq& q); void seed(result_type value = 0u); template<class Sseq> void seed(Sseq& q); // equality operators friend bool operator==(const subtract_with_carry_engine& x, const subtract_with_carry_engine& y); // generating functions result_type operator()(); void discard(unsigned long long z); // inserters and extractors template<class charT, class traits> friend basic_ostream<charT, traits>& operator<<(basic_ostream<charT, traits>& os, const subtract_with_carry_engine& x); template<class charT, class traits> friend basic_istream<charT, traits>& operator>>(basic_istream<charT, traits>& is, subtract_with_carry_engine& x); }; }

The following relations shall hold: 0u < s, s < r, 0 < w, and w <= numeric_limits<UIntType>::digits.

The textual representation consists of the values of

X_{i - r}, \dots, X_{i - 1}

, in that order, followed by c.

🔗

explicit subtract_with_carry_engine(result_type value);

Effects: Sets the values of

X_{- r}, \dots, X_{- 1}

, in that order, as specified below.

X_{- 1}

is then 0, sets c to 1; otherwise sets c to 0.

To set the values

X_{k}

, first construct e, a linear_congruential_engine object, as if by the following definition: linear_congruential_engine<uint_least32_t, 40014u,0u,2147483563u> e(value == 0u ? default_seed : value);

Then, to set each

X_{k}

, obtain new values

z_{0}, \dots, z_{n - 1}

from

n = ⌈ w / 32 ⌉

successive invocations of e.

Set

X_{k}

(\sum_{j = 0}^{n - 1} z_{j} \cdot 2^{32 j}) mod m

Complexity: Exactly

n \cdot r

invocations of e.

🔗

template<class Sseq> explicit subtract_with_carry_engine(Sseq& q);

Effects: With

k = ⌈ w / 32 ⌉

and a an array (or equivalent) of length

r \cdot k

, invokes q.generate(

a + 0

a + r \cdot k

) and then, iteratively for

i = - r, \dots, - 1

, sets

X_{i}

(\sum_{j = 0}^{k - 1} a_{k (i + r) + j} \cdot 2^{32 j}) mod m

X_{- 1}

is then 0, sets c to 1; otherwise sets c to 0.

29.5.4.5 Class template philox_engine [rand.eng.philox]

A philox_engine random number engine produces unsigned integer random numbers in the closed interval [0, m], where

m = 2^{w} - 1

and the template parameter w defines the range of the produced numbers.

The state of a philox_engine object consists of a sequence X of n unsigned integer values of width w, a sequence K of

n / 2

values of result_type, a sequence Y of n values of result_type, and a scalar i, where

(1.1)
X is the interpretation of the unsigned integer counter value $Z : = \sum_{j = 0}^{n - 1} X_{j} \cdot 2^{w j}$ of $n \cdot w$ bits,
(1.2)
K are keys, which are generated once from the seed (see constructors below) and remain constant unless the seed function ([rand.req.eng]) is invoked,
(1.3)
Y stores a batch of output values, and
(1.4)
i is an index for an element of the sequence Y.

The generation algorithm returns

Y_{i}

, the value stored in the

i^{t h}

element of Y after applying the transition algorithm.

The state transition is performed as if by the following algorithm: i = i + 1 if (i == n) { Y = Philox(K, X) // see below Z = Z + 1 i = 0 }

The Philox function maps the length-

n / 2

sequence K and the length-n sequence X into a length-n output sequence Y.

Philox applies an r-round substitution-permutation network to the values in X.

A single round of the generation algorithm performs the following steps:

(4.1)
The output sequence $X^{'}$ of the previous round (X in case of the first round) is permuted to obtain the intermediate state V: $V_{j} = X_{f_{n} (j)}^{'}$ where $j = 0, \dots, n - 1$ and $f_{n} (j)$ is defined in Table 124.
Table 124: Values for the word permutation $f_{n} (j)$ [tab:rand.eng.philox.f]
🔗
$f_{n} (j)$
j
🔗
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
🔗
n
2
0
1
🔗
4
0
3
2
1
🔗
8
2
1
4
7
6
5
0
3
🔗
16
0
9
2
13
6
11
4
15
10
7
12
3
14
5
8
1
🔗

[Note 1:
For $n = 2$ the sequence is not permuted.
— end note]
(4.2)
The following computations are applied to the elements of the V sequence: $X_{2 k + 0} = mullo (V_{2 k + 1}, M_{k}, w)$ $X_{2 k + 1} = mulhi (V_{2 k + 1}, M_{k}, w) x o r {k e y}_{k}^{q} x o r V_{2 k}$ where:
- (4.2.1)
  mullo(a, b, w) is the low half of the modular multiplication of a and b: $(a \cdot b) mod 2^{w}$ ,
- (4.2.2)
  mulhi(a, b, w) is the high half of the modular multiplication of a and b: $(⌊ (a \cdot b) / 2^{w} ⌋)$ ,
- (4.2.3)
  $k = 0, \dots, n / 2 - 1$ is the index in the sequences,
- (4.2.4)
  $q = 0, \dots, r - 1$ is the index of the round,
- (4.2.5)
  ${k e y}_{k}^{q}$ is the $k^{th}$ round key for round q, ${k e y}_{k}^{q} : = (K_{k} + q \cdot C_{k}) mod 2^{w}$ ,
- (4.2.6)
  $K_{k}$ are the elements of the key sequence K,
- (4.2.7)
  $M_{k}$ is multipliers[k], and
- (4.2.8)
  $C_{k}$ is round_consts[k].

After r applications of the single-round function, Philox returns the sequence

Y = X^{'}

🔗

namespace std { template<class UIntType, size_t w, size_t n, size_t r, UIntType... consts> class philox_engine { static constexpr size_t array-size = n / 2; // exposition only public: // types using result_type = UIntType; // engine characteristics static constexpr size_t word_size = w; static constexpr size_t word_count = n; static constexpr size_t round_count = r; static constexpr array<result_type, array-size> multipliers; static constexpr array<result_type, array-size> round_consts; static constexpr result_type min() { return 0; } static constexpr result_type max() { return m - 1; } static constexpr result_type default_seed = 20111115u; // constructors and seeding functions philox_engine() : philox_engine(default_seed) {} explicit philox_engine(result_type value); template<class Sseq> explicit philox_engine(Sseq& q); void seed(result_type value = default_seed); template<class Sseq> void seed(Sseq& q); void set_counter(const array<result_type, n>& counter); // equality operators friend bool operator==(const philox_engine& x, const philox_engine& y); // generating functions result_type operator()(); void discard(unsigned long long z); // inserters and extractors template<class charT, class traits> friend basic_ostream<charT, traits>& operator<<(basic_ostream<charT, traits>& os, const philox_engine& x); template<class charT, class traits> friend basic_istream<charT, traits>& operator>>(basic_istream<charT, traits>& is, philox_engine& x); }; }

Mandates:

(6.1)
sizeof...(consts) == n is true, and
(6.2)
n == 2 || n == 4 || n == 8 || n == 16 is true, and
(6.3)
0 < r is true, and
(6.4)
0 < w && w <= numeric_limits<UIntType>::digits is true.

The template parameter pack consts represents the

M_{k}

and

C_{k}

constants which are grouped as follows:

[M_{0}, C_{0}, M_{1}, C_{1}, M_{2}, C_{2}, \dots, M_{n / 2 - 1}, C_{n / 2 - 1}]

The textual representation consists of the values of

K_{0}, \dots, K_{n / 2 - 1}, X_{0}, \dots, X_{n - 1}, i

, in that order.

[Note 2:

The stream extraction operator can reconstruct Y from K and X, as needed.

— end note]

🔗

explicit philox_engine(result_type value);

Effects: Sets the

K_{0}

element of sequence K to

value mod 2^{w}

All elements of sequences X and K (except

K_{0}

) are set to 0.

The value of i is set to

n - 1

🔗

template<class Sseq> explicit philox_engine(Sseq& q);

Effects: With

p = ⌈ w / 32 ⌉

and an array (or equivalent) a of length

(n / 2) \cdot p

, invokes q.generate(a + 0, a + n / 2 * p) and then iteratively for

k = 0, \dots, n / 2 - 1

, sets

K_{k}

(\sum_{j = 0}^{p - 1} a_{k p + j} \cdot 2^{32 j}) mod 2^{w}

All elements of sequence X are set to 0.

The value of i is set to

n - 1

🔗

void set_counter(const array<result_type, n>& c);

Effects: For

j = 0, \dots, n - 1

sets

X_{j}

C_{n - 1 - j} mod 2^{w}

The value of i is set to

n - 1

[Note 3:

The counter is the value Z introduced at the beginning of this subclause.

— end note]

🔗	$f_{n} (j)$		j
🔗			0	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15
🔗	n	2	0	1
🔗		4	0	3	2	1
🔗		8	2	1	4	7	6	5	0	3
🔗		16	0	9	2	13	6	11	4	15	10	7	12	3	14	5	8	1
🔗