# The C++ Programming Language (2013)

### Part IV: The Standard Library

### 40. Numerics

*The purpose of computing is insight, not numbers.*

*– R. W. Hamming*

*... but for the student,numbers are often the best road to insight.*

*– A. Ralston*

• *Introduction*

• *Numerical Limits*

*Limit Macros*

• *Standard Mathematical Functions*

• *complex** Numbers*

• *A Numerical Array: **valarray*

*Constructors and Assignments*; *Subscripting*; *Operations*; *Slices*; ** slice_array**;

*Generalized Slices*

• *Generalized Numerical Algorithms*

** accumulate()**;

**;**

*inner_product()*

*partial_sum()**and*

**;**

*adjacent_difference()*

*iota()*• *Random Numbers*

*Engines*; *Random Device*; *Distributions*; *C-Style Random Numbers*

• *Advice*

**40.1. Introduction**

C++ was not designed primarily with numeric computation in mind. However, numeric computation typically occurs in the context of other work – such as database access, networking, instrument control, graphics, simulation, and financial analysis – so C++ becomes an attractive vehicle for computations that are part of a larger system. Furthermore, numeric methods have come a long way from being simple loops over vectors of floating-point numbers. Where more complex data structures are needed as part of a computation, C++’s strengths become relevant. The net effect is that C++ is widely used for scientific, engineering, financial, and other computation involving sophisticated numerics. Consequently, facilities and techniques supporting such computation have emerged. This chapter describes the parts of the standard library that support numerics. I make no attempt to teach numeric methods. Numeric computation is a fascinating topic in its own right. To understand it, you need a good course in numerical methods or at least a good textbook – not just a language manual and tutorial.

In addition to the standard-library facilities described here, *Chapter 29* is an extended example of numerical programming: an N-dimensional matrix.

**40.2. Numerical Limits**

To do anything interesting with numbers, we typically need to know something about the general properties of built-in numeric types. To allow the programmer to best take advantage of hardware, these properties are implementation-defined rather than fixed by the rules of the language itself (§*6.2.8*). For example, what is the largest **int**? What is the smallest positive **float**? Is a **double** rounded or truncated when assigned to a **float**? How many bits are there in a **char**?

Answers to such questions are provided by the specializations of the **numeric_limits** template presented in **<limits>**. For example:

**void f(double d, int i){ char classification[numeric_limits<unsigned char>::max()];**

**if (numeric_limits<unsigned char>::digits==numeric_limits<char>::digits) {**

**//**

*chars are unsigned*

**}**

**if (i<numeric_limits<short>::min() || numeric_limits<short>::max()<i) {**

**//**

*i cannot be stored in a short without loss of digits*

**}**

**if (0<d && d<numeric_limits<double>::epsilon()) d = 0;**

**if (numeric_limits<Quad>::is_specialized) {**

**//**

*limits information is available for type Quad*

**}**

}

}

Each specialization provides the relevant information for its argument type. Thus, the general **numeric_limits** template is simply a notational handle for a set of constants and **constexpr** functions:

**template<typename T>class numeric_limits {public: static const bool is_specialized = false; //**

*is information available for numeric_limits<T>?*

**//**

*... uninteresting defaults ...*

**};**

The real information is in the specializations. Each implementation of the standard library provides a specialization of **numeric_limits** for each fundamental numeric type (the character types, the integer types, the floating-point types, and **bool**) but not for any other plausible candidates such as **void**, enumerations, or library types (such as **complex<double>**).

For an integral type such as **char**, only a few pieces of information are of interest. Here is **numeric_limits<char>** for an implementation in which a **char** has 8 bits and is signed:

**template<>class numeric_limits<char> {public: static const bool is_specialized = true; //**

*yes, we have information*

**static const int digits = 7; //**

*number of bits ("binary digits") excluding sign*

**static const bool is_signed = true; //**

*this implementation has char signed*

**static const bool is_integer = true; //**

*char is an integral type*

**static constexpr char min() noexcept { return –128; } //**

*smallest value*

**static constexpr char max() noexcept { return 127; } //**

*largest value*

**//**

*lots of declarations not relevant to a char*

**};**

The functions are **constexpr**, so that they can be used where a constant expression is required and without run-time overhead.

Most members of **numeric_limits** are intended to describe floating-point numbers. For example, this describes one possible implementation of **float**:

**template<>class numeric_limits<float> {public: static const bool is_specialized = true;**

**static const int radix = 2; //**

*base of exponent (in this case, binary)*

**static const int digits = 24; //**

*number of radix digits in mantissa*

**static const int digits10 = 9; //**

*number of base 10 digits in mantissa*

**static const bool is_signed = true;**

static const bool is_integer = false;

static const bool is_exact = false;

static const bool is_integer = false;

static const bool is_exact = false;

**static constexpr float min() noexcept { return 1.17549435E–38F; } //**

*smallest positive*

**static constexpr float max() noexcept { return 3.40282347E+38F; } //**

*largest positive*

**static constexpr float lowest() noexcept { return –3.40282347E+38F; } //**

*smallest value*

**static constexpr float epsilon() noexcept { return 1.19209290E–07F; }**

static constexpr float round_error() noexcept { return 0.5F; } //

static constexpr float round_error() noexcept { return 0.5F; } //

*maximum rounding error*

**static constexpr float infinity() noexcept { return /***

*some value*

***/; }**

**static constexpr float quiet_NaN() noexcept { return /***

*some value*

***/; }**

static constexpr float signaling_NaN() noexcept { return /*

static constexpr float signaling_NaN() noexcept { return /*

*some value*

***/; }**

static constexpr float denorm_min() noexcept { return min(); }

static constexpr float denorm_min() noexcept { return min(); }

**static const int min_exponent = –125;**

static const int min_exponent10 = –37;

static const int max_exponent = +128;

static const int max_exponent10 = +38;

static const int min_exponent10 = –37;

static const int max_exponent = +128;

static const int max_exponent10 = +38;

**static const bool has_infinity = true;**

static const bool has_quiet_NaN = true;

static const bool has_signaling_NaN = true;

static const float_denorm_style has_denorm = denorm_absent;

static const bool has_denorm_loss = false;

static const bool has_quiet_NaN = true;

static const bool has_signaling_NaN = true;

static const float_denorm_style has_denorm = denorm_absent;

static const bool has_denorm_loss = false;

**static const bool is_iec559 = true; //**

*conforms to IEC-559*

**static const bool is_bounded = true;**

static const bool is_modulo = false;

static const bool traps = true;

static const bool tinyness_before = true;

static const bool is_modulo = false;

static const bool traps = true;

static const bool tinyness_before = true;

**static const float_round_style round_style = round_to_nearest;**

};

};

Note that **min()** is the smallest *positive* normalized number and that **epsilon** is the smallest positive floating-point number such that **1+epsilon–1** is larger than **0**.

When defining a scalar type along the lines of the built-in ones, it is a good idea also to provide a suitable specialization of **numeric_limits**. For example, if I write a quadruple-precision type **Quad**, a user could reasonably expect me to provide **numeric_limits<Quad>**. Conversely, if I use a nonnumeric type, **Dumb_ptr**, I would expect for **numeric_limits<Dumb_ptr<X>>** to be the primary template that has **is_specialized** set to **false**, indicating that no information is available.

We can imagine specializations of **numeric_limits** describing properties of user-defined types that have little to do with floating-point numbers. In such cases, it is usually better to use the general technique for describing properties of a type than to specialize **numeric_limits** with properties not considered in the standard.

**40.2.1. Limit Macros**

From C, C++ inherited macros that describe properties of integers. They are found in **<climits>**:

Analogously named macros for **signed char**s, **long long**, etc., are also provided.

Similarly, **<cfloat>** and **<float.h>** define macros describing properties of floating-point numbers:

Analogously named macros for **long double** are also provided.

**40.3. Standard Mathematical Functions**

In **<cmath>** we find what are commonly called the *standard mathematical functions*:

There are versions taking **float**, **double**, **long double**, and **complex** (§*40.4*) arguments. For each function, the return type is the same as the argument type.

Errors are reported by setting **errno** from **<cerrno>** to **EDOM** for a domain error and to **ERANGE** for a range error. For example:

**void f(){ errno = 0; //**

*clear old error state*

**sqrt(–1);**

if (errno==EDOM) cerr << "sqrt() not defined for negative argument";

if (errno==EDOM) cerr << "sqrt() not defined for negative argument";

**pow(numeric_limits<double>::max(),2);**

if (errno == ERANGE) cerr << "result of pow() too large to represent as a double";

}

if (errno == ERANGE) cerr << "result of pow() too large to represent as a double";

}

For historical reasons, a few mathematical functions are found in **<cstdlib>** rather than in **<cmath>**:

The **l*()** versions are there because C does not support overloading. The results of the **ldiv()** functions are **struct**s **div_t**, **ldiv_t**, and **lldiv_t**, respectively. These **struct**s have members **quot** (for quotient) and **rem** (for remainder) of implementation-defined types.

There is a separate ISO standard for *special mathematical functions* [C++Math,2010]. An implementation may add these functions to **<cmath>**:

If you don’t know these functions, you are unlikely to need them.

**40.4. complex Numbers**

The standard library provides complex number types **complex<float>**, **complex<double>**, and **complex<long double>**. A **complex<Scalar>** where **Scalar** is some other type supporting the usual arithmetic operations usually works but is not guaranteed to be portable.

**template<typename Scalar>class complex {**

**//**

*a complex is a pair of scalar values, basically a coordinate pair*

**Scalar re, im;**

public:

complex(const Scalar & r = Scalar{}, const Scalar & i = Scalar{}) :re(r), im(i) { }

public:

complex(const Scalar & r = Scalar{}, const Scalar & i = Scalar{}) :re(r), im(i) { }

**Scalar real() const { return re; } //**

*real part*

**void real(Scalar r) { re=r; }**

**Scalar imag() const { return im; } //**

*imaginary part*

**void imag(Scalar i) { im = i; }**

**template<typename X>**

complex(const complex<X>&);

complex(const complex<X>&);

**complex<T>& operator=(const T&);**

complex& operator=(const complex&);

template<typename X>

complex<T>& operator=(const complex<X>&);

complex& operator=(const complex&);

template<typename X>

complex<T>& operator=(const complex<X>&);

**complex<T>& operator+=(const T&);**

template<typename X>

complex<T>& operator+=(const complex<X>&);

template<typename X>

complex<T>& operator+=(const complex<X>&);

**//**

*similar for operators -=, *=, /=*

**};**

The standard-library **complex** does not protect against narrowing:

**complex<float> z1 = 1.33333333333333333; //** *narrows***complex<double> z2 = 1.33333333333333333;//** *narrows***z1=z2; //** *narrows*

To protect against accidental narrowing, use **{}** initialization:

**complex<float> z3 {1.33333333333333333}; //** *error: narrowing conversion*

In addition to the members of **complex**, **<complex>** offers a host of useful operations:

The standard mathematical functions (§*40.3*) are also available for complex numbers. Note that **complex** does not provide **<** or **%**. For more details, see §*18.3*.

**40.5. A Numerical Array: valarray**

Much numeric work relies on relatively simple single-dimensional vectors of floating-point values. In particular, such vectors are well supported by high-performance machine architectures, libraries relying on such vectors are in wide use, and very aggressive optimization of code using such vectors is considered essential in many fields. The **valarray** from **<valarray>** is a single-dimensional numerical array. It provides the usual numeric vector arithmetic operations for an array type plus support for slices and strides:

The fundamental idea of **valarray** was to provide Fortran-like facilities for dense multidimensional arrays with Fortran-like opportunities for optimization. This can only be achieved with the active support of compiler and optimization suppliers and the addition of more library support on top of the very basic facilities provided by **valarray**. So far, that has not happened for all implementations.

**40.5.1. Constructors and Assignments**

The **valarray** constructors allow us to initialize **valarray**s from the auxiliary numeric array types and from single values:

For example:

**valarray<double> v0; //** *placeholder, we can assign to v0 later***valarray<float> v1(1000); //** *1000 elements with value float()==0.0F***valarray<int> v2(–1,2000); //** *2000 elements with value –1***valarray<double> v3(100,9.8064); //** *bad mistake: floating-point valarray size***valarray<double> v4 = v3; //** *v4 has v3.size() elements***valarray<int> v5 {–1,2000}; //** *two elements*

In the two-argument constructors, the value comes before the number of elements. This differs from the convention for standard containers (§*31.3.2*).

The number of elements of an argument **valarray** to a copy constructor determines the size of the resulting **valarray**.

Most programs need data from tables or input. In addition to initializer lists, this is supported by a constructor that copies elements from a built-in array. For example:

**void f(const int* p, int n){ const double vd[] = {0,1,2,3,4}; const int vi[] = {0,1,2,3,4};**

**valarray<double> v1{vd,4}; //**

*4 elements: 0,1,2,3*

**valarray<double> v2{vi,4}; //**

*type error: vi is not pointer to double*

**valarray<double> v3{vd,8}; //**

*undefined: too few elements in initializer*

**valarray<int> v4{p,n}; //**

*p had better point to at least n ints*

**}**

The **valarray** and its auxiliary facilities were designed for high-speed computing. This is reflected in a few constraints on users and by a few liberties granted to implementers. Basically, an implementer of **valarray** is allowed to use just about every optimization technique you can think of. The **valarray** operations are assumed to be free of side effects (except on their explicit arguments, of course), **valarray**s are assumed to be alias free, and the introduction of auxiliary types and the elimination of temporaries is allowed as long as the basic semantics are maintained. There is no range checking. The elements of a **valarray** must have the default copy semantics (§*8.2.6*).

Assignment can be with another **valarray**, a scalar, or a subset of a **valarray**:

A **valarray** can be assigned to another of the same size. As one would expect, **v1=v2** copies every element of **v2** into its corresponding position in **v1**. If **valarray**s have different sizes, the result of assignment is undefined.

In addition to this conventional assignment, it is possible to assign a scalar to a **valarray.** For example, **v=7** assigns **7** to every element of the **valarray v**. This may be surprising to some programmers and is best understood as an occasionally useful degenerate case of the operator assignment operations. For example:

**valarray<int> v {1,2,3,4,5,6,7,8};v *= 2; //**

*v=={2,4,6,10,12,14,16}*

**v = 7; //**

*v=={7,7,7,7,7,7,7,7}*

**40.5.2. Subscripting**

Subscripting can be used to select an element of a **valarray** or a subset of its elements:

Each **operator[]** returns a subset of the elements from a **valarray**. The return type (the type of the object representing the subset) depends on the argument type.

For **const** arguments, the result contains copies of elements. For non-**const** arguments, the result holds references to elements. Since C++ doesn’t directly support arrays of references (e.g., we can’t say **valarray<int&>**), the implementation will somehow simulate this. This can be done efficiently. An exhaustive list, with examples (based on §iso.26.6.2.5), is in order. In each case, the subscript describes the elements to be returned, and **v1** must be a **vallaray** with an appropriate length and element type:

• A **slice** of a **const valarray**:

**valarray<T> operator[](slice) const;//** *copy of elements*

**//** *...*

**const valarray<char> v0 {"abcdefghijklmnop",16}; valarray<char> v1 {v0[slice(2,5,3)]}; //**

*{"cfilo",5}*

• A **slice** of a non-**const valarray**:

**slice_array<T> operator[](slice); //** *references to elements*

**//** *...*

**valarray<char> v0 {"abcdefghijklmnop",16}; valarray<char> v1 {"ABCDE",5}; v0[slice(2,5,3)] = v1; //**

*v0=={"abAdeBghCjkDmnEp",16}*

• A **gslice** of a **const valarray**:

**valarray<T> operator[](const gslice&) const; //** *copies of elements*

**//** *...*

**const valarray<char> v0 {"abcdefghijklmnop",16}; const valarray<size_t> len {2,3}; const valarray<size_t> str {7,2}; valarray<char> v1 {v0[gslice(3,len,str)]}; //**

*v1=={"dfhkmo",6}*

• A **gslice** of a non-**const valarray**:

**gslice_array<T> operator[](const gslice&); //** *references to elements*

**//** *...*

**valarray<char> v0 {"abcdefghijklmnop",16}; valarray<char> v1 {"ABCDE",5}; const valarray<size_t> len {2,3}; const valarray<size_t> str {7,2}; v0[gslice(3,len,str)] = v1; //**

*v0=={"abcAeBgCijDlEnFp",16}*

• A **valarray<bool>** (a mask) of a **const valarray**:

**valarray<T> operator[](const valarray<bool>&) const; //** *copies of elements*

**//** *...*

**const valarray<char> v0 {"abcdefghijklmnop",16}; const bool vb[] {false, false, true, true, false, true}; valarray<char> v1 {v0[valarray<bool>(vb, 6)]}; //**

*v1=={"cdf",3}*

• A **valarray<bool>** (a mask) of a non-**const valarray**:

**mask_array<T> operator[](const valarray<bool>&); //** *references to elements*

**//** *...*

**valarray<char> v0 {"abcdefghijklmnop", 16}; valarray<char> v1 {"ABC",3}; const bool vb[] {false, false, true, true, false, true}; v0[valarray<bool>(vb,6)] = v1; //**

*v0=={"abABeCghijklmnop",16}*

• A **valarray<size_t>** (a set of indices) of a **const valarray**:

**valarray<T> operator[](const valarray<size_t>&) const; //** *references to elements*

**//** *...*

**const valarray<char> v0 {"abcdefghijklmnop",16}; const size_t vi[] {7, 5, 2, 3, 8}; valarray<char> v1 {v0[valarray<size_t>(vi,5)]}; //**

*v1=={"hfcdi",5}*

• A **valarray<size_t>** (a set of indices) of a non-**const valarray**:

**indirect_array<T> operator[](const valarray<size_t>&); //** *references to elements*

**//** *...*

**valarray<char> v0 {"abcdefghijklmnop",16}; valarray<char> v1 {"ABCDE",5}; const size_t vi[] {7, 5, 2, 3, 8}; v0[valarray<size_t>(vi,5)] {v1}; //**

*v0=={"abCDeBgAEjklmnop",16}*

Note that subscripting with a mask (a **valarray<bool>**) yields a **mask_array**, and subscripting with a set of indices (a **valarray<size_t>**) yields an **indirect_array**.

**40.5.3. Operations**

The purpose of **valarray** is to support computation, so a host of basic numerical operations are directly supported:

There is no range checking: the effect of using a function that tries to access an element of an empty **valarray** is undefined.

Note that **resize()** does not retain any old values.

The binary operations are defined for **valarray**s and for combinations of a **valarray** and its scalar type. A scalar type is treated as a **valarray** of the right size with every element having the scalar’s value. For example:

**void f(valarray<double>& v, valarray<double>& v2, double d){ valarray<double> v3 = v*v2; //**

*v3[i] = v[i]*v2[i] for all i*

**valarray<double> v4 = v*d; //**

*v4[i] = v[i]*d for all i*

**valarray<double> v5 = d*v2; //**

*v5[i] = d*v2[i] for all i*

**valarray<double> v6 = cos(v); //**

*v6[i] = cos(v[i]) for all i*

**}**

These vector operations all apply their operations to each element of their operand(s) in the way indicated by the ***** and **cos()** examples. Naturally, an operation can be used only if the corresponding operation is defined for the scalar type. Otherwise, the compiler will issue an error when trying to specialize the operator or function.

Where the result is a **valarray**, its length is the same as its **valarray** operand. If the lengths of the two arrays are not the same, the result of a binary operator on two **valarray**s is undefined.

These **valarray** operations return new **valarray**s rather than modifying their operands. This can be expensive, but it doesn’t have to be when aggressive optimization techniques are applied.

For example, if **v** is a **valarray**, it can be scaled like this: **v*=0.2**, and this: **v/=1.3**. That is, applying a scalar to a vector means applying the scalar to each element of the vector. As usual, ***=** is more concise than a combination of ***** and **=** (§*18.3.1*) and easier to optimize.

Note that the non-assignment operations construct a new **valarray**. For example:

**double incr(double d) { return d+1; }****void f(valarray<double>& v){ valarray<double> v2 = v.apply(incr); //**

*produce incremented valarray*

**//**

*...*

**}**

This does not change the value of **v**. Unfortunately, **apply()** does not accept a function object (§*3.4.3*, §*11.4*) as an argument.

The logical and cyclic shift functions, **shift()** and **cshift()**, return a new **valarray** with the elements suitably shifted and leave the original one unchanged. For example, the cyclic shift **v2=v.cshift(n)** produces a **valarray** so that **v2[i]==v[(i+n)%v.size()]**. The logical shift **v3=v.shift(n)**produces a **valarray** so that **v3[i]** is **v[i+n]** if **i+n** is a valid index for **v**. Otherwise, the result is the default element value. This implies that both **shift()** and **cshift()** shift left when given a positive argument and right when given a negative argument. For example:

**void f(){ int alpha[] = { 1, 2, 3, 4, 5 ,6, 7, 8 }; valarray<int> v(alpha,8); //**

*1, 2, 3, 4, 5, 6, 7, 8*

**valarray<int> v2 = v.shift(2); //**

*3, 4, 5, 6, 7, 8, 0, 0*

**valarray<int> v3 = v<<2; //**

*4, 8, 12, 16, 20, 24, 28, 32*

**valarray<int> v4 = v.shift(–2); //**

*0, 0, 1, 2, 3, 4, 5, 6*

**valarray<int> v5 = v>>2; //**

*0, 0, 0, 1, 1, 1, 1, 2*

**valarray<int> v6 = v.cshift(2); //**

*3, 4, 5, 6, 7, 8, 1, 2*

**valarray<int> v7 = v.cshift(–2); //**

*7, 8, 1, 2, 3, 4, 5, 6*

**}**

For **valarray**s, **>>** and **<<** are bit shift operators, rather than element shift operators or I/O operators. Consequently, **<<=** and **>>=** can be used to shift bits within elements of an integral type. For example:

**void f(valarray<int> vi, valarray<double> vd){ vi <<= 2; //**

*vi[i]<<=2 for all elements of vi*

**vd <<= 2; //**

*error: shift is not defined for floating-point values*

**}**

All of the operators and mathematical functions on **valarray**s can also be applied to **slice_array**s (§*40.5.5*), **gslice_array**s (§*40.5.6*), **mask_array**s (§*40.5.2*), **indirect_array**s (§*40.5.2*), and combinations of these types. However, an implementation is allowed to convert an operand that is not a **valarray** to a **valarray** before performing a required operation.

**40.5.4. Slices**

A **slice** is an abstraction that allows us to manipulate a one-dimensional array (e.g., a built-in array, a **vector,** or a **valarray**) efficiently as a matrix of arbitrary dimension. It is the key notion of Fortran vectors and of the BLAS (Basic Linear Algebra Subprograms) library, which is the basis for much numeric computation. Basically, a **slice** is every **n**th element of some part of an array:

**class std::slice {**

**//** *starting index, a length, and a stride***public: slice(); //**

*slice{0,0,0}*

**slice(size_t start, size_t size, size_t stride);**

**size_t start() const; //**

*index of first element*

**size_t size() const; //**

*number of elements*

**size_t stride() const; //**

*element n is at start()+n*stride()*

**};**

A *stride* is the distance (in number of elements) between two elements of the **slice**. Thus, a **slice** describes a mapping of non-negative integers into indices. The number of elements (the **size()**) doesn’t affect the mapping (addressing) but allows us to find the end of a sequence. This mapping can be used to simulate two-dimensional arrays within a one-dimensional array (such as **valarray**) in an efficient, general, and reasonably convenient way. Consider a 3-by-4 matrix (three rows, each with four elements):

**valarray<int> v {**

**{00,01,02,03}, //** *row 0*

**{10,11,12,13}, //** *row 1*

**{20,21,22,23} //** *row 2***};**

or graphically:

Following the usual C/C++ conventions, the **valarray** is laid out in memory with row elements first (*row-major* order) and contiguous:

**for (int x : v) cout << x << ' ' ;**

This produces:

**0 1 2 3 10 11 12 13 20 21 22 23**

or graphically:

Row **x** is described by **slice(x*4,4,1)**. That is, the first element of row **x** is the **x*4**th element of the vector, the next element of the row is the **(x*4+1)**th, etc., and there are **4** elements in each row. For example, **slice{0,4,1}** describes the first row of **v** (row **0**): **00**, **01**, **02**, **03**, and **slice{1,4,1}**describes the second row (row **1**).

Column **y** is described by **slice(y,3,4)**. That is, the first element of column **y** is the **y**th element of the vector, the next element of the column is the **(y+4)**th, etc., and there are **3** elements in each column. For example, **slice{0,3,4}** describes the first column (column **0**): **00**, **10**, **20**, and**slice{1,3,4}** describes the second column (column **1**).

In addition to its use for simulating two-dimensional arrays, a **slice** can describe many other sequences. It is a fairly general way of specifying very simple sequences. This notion is explored further in §*40.5.6*.

One way of thinking of a slice is as an odd kind of iterator: a **slice** allows us to describe a sequence of indices for a **valarray**. We could build an STL-style iterator based on that:

**template<typename T>class Slice_iter { valarray<T>*v; slice s; size_t curr; //**

*index of current element*

**T& ref(size_t i) const { return (*v)[s.start()+i*s.stride()]; }**

public:

public:

**Slice_iter(valarray<T>*vv, slice ss, size_t pos =0)**

**:v{vv}, s{ss}, curr{0} { }**

**Slice_iter end() const { return {this,s,s.size()}; }**

**Slice_iter& operator++() { ++curr; return**

***this; }**

**Slice_iter operator++(int) { Slice_iter t = *this; ++curr; return t; }**

**T& operator[](size_t i) { return ref(i); } //**

*C-style subscript*

**T& operator()(size_t i) { return ref(i); } //**

*Fortran-style subscript*

**T& operator*() { return ref(curr); } //**

*current element*

**bool operator==(const Slice_iter& q) const**

{

return curr==q.curr && s.stride()==q.s.stride() && s.start()==q.s.start();

}

{

return curr==q.curr && s.stride()==q.s.stride() && s.start()==q.s.start();

}

**bool operator!=(const Slice_iter& q) const**

{

return !(*this==q);

}

{

return !(*this==q);

}

**bool operator<(const Slice_iter& q) const**

{

return curr<q.curr && s.stride()==q.s.stride() && s.start()==q.s.start();

};

}

{

return curr<q.curr && s.stride()==q.s.stride() && s.start()==q.s.start();

};

}

Since a **slice** has a size, we could even provide range checking. Here, I have taken advantage of **slice::size()** to provide an **end()** operation to provide an iterator for the one-past-the-end element of the **slice**.

Since a **slice** can describe either a row or a column, the **Slice_iter** allows us to traverse a **valarray** by row or by column.

**40.5.5. slice_array**

From a **valarray** and a **slice**, we can build something that looks and feels like a **valarray** but is really simply a way of referring to the subset of the array described by the slice.

A user cannot directly create a **slice_array**. Instead, the user subscripts a **valarray** to create a **slice_array** for a given slice. Once the **slice_array** is initialized, all references to it indirectly go to the **valarray** for which it is created. For example, we can create something that represents every second element of an array like this:

**void f(valarray<double>& d){ slice_array<double>& v_even = d[slice(0,d.size()/2+d.size()%2,2)]; slice_array<double>& v_odd = d[slice(1,d.size()/2,2)];**

**v_even *= v_odd; //**

*multiply element pairs and store results in even elements*

**v_odd = 0; //**

*assign 0 to every odd element of d*

**}**

A **slice_array** can be copied. For example:

**slice_array<double> row(valarray<double>& d, int i){ slice_array<double> v = d[slice(0,2,d.size()/2)]; //**

*...*

**return d[slice(i%2,i,d.size()/2)];**

}

}

**40.5.6. Generalized Slices**

A **slice** (§*29.2.2*, §*40.5.4*) can describe a row or a column of an **n**-dimensional array. However, sometimes we need to extract a subarray that is not a row or a column. For example, we might want to extract the 3-by-2 matrix from the top-left corner of a 4-by-3 matrix:

Unfortunately, these elements are not allocated in a way that can be described by a single **slice**:

A **gslice** is a “generalized slice” that contains (almost) the information from **n slices**:

**class std::gslice {**

**//** *instead of 1 stride and one size like slice, gslice holds n strides and n sizes***public: gslice(); gslice(size_t sz, const valarray<size_t>& lengths, const valarray<size_t>& strides);**

**size_t start() const; //**

*index of first element*

**valarray<size_t> size() const; //**

*number of elements in dimension*

**valarray<size_t> stride() const; //**

*stride for index[0], index[1], ...*

**};**

The extra values allow a **gslice** to specify a mapping between **n** integers and an index to be used to address elements of an array. For example, we can describe the layout of the 3-by-2 matrix by a pair of (length,stride) pairs:

**size_t gslice_index(const gslice& s, size_t i, size_t j) //** *max (i,j) to their corresponding index***{ return s.start()+i*s.stride()[0]+j*s.stride()[1];}**

**valarray<size_t> lengths {2,3};//**

*2 elements in the first dimension*

**//**

*3 elements in the second dimension*

**valarray<size_t> strides {3,1}; //**

*3 is the stride for the first index*

**//**

*1 is the stride for the second index*

**void f()**

{

gslice s(0,lengths,strides);

{

gslice s(0,lengths,strides);

**for (int i=0; i<3; ++i) //**

*for each row*

**for (int j=0; j<2; ++j) //**

*for each element in row*

**cout << "(" << i << "," << j << ")–>" << gslice_index(s,i,j) << "; "; //**

*print mapping*

**}**

This prints:

**(0,0)–>0; (0,1)–>1; (1,0)–>3; (1,1)–>4; (2,0)–>6; (2,1)–>7**

In this way, a **gslice** with two (length,stride) pairs describes a subarray of a two-dimensional array, a **gslice** with three (length,stride) pairs describes a subarray of a three-dimensional array, etc. Using a **gslice** as the index of a **valarray** yields a **gslice_array** consisting of the elements described by the **gslice.** For example:

**void f(valarray<float>& v){ gslice m(0,lengths,strides); v[m] = 0; //**

*assign 0 to v[0],v[1],v[3],v[4],v[6],v[7]*

**}**

The **gslice_array** offers the same set of members as **slice_array** (§*40.5.5*). A **gslice_array** is the result of using a **gslice** as the subscript of a **valarray** (§*40.5.2*).

**40.6. Generalized Numerical Algorithms**

In **<numeric>**, the standard library provides a few generalized numeric algorithms in the style of the non-numeric algorithms from **<algorithm>** (*Chapter 32*). These algorithms provide general versions of common operations on sequences of numerical values:

These algorithms generalize common operations such as computing a sum by letting them apply to all kinds of sequences and by making the operation applied to elements of those sequences a parameter. For each algorithm, the general version is supplemented by a version applying the most common operator for that algorithm.

**40.6.1. accumulate()**

The simple version of **accumulate()** adds elements of a sequence using their **+** operator:

**template<typename In, typename T>****T accumulate(In first, In last, T init){ for (; first!=last; ++first) //**

*for all elements in [first:last)*

**init = init + *first; //**

*plus*

**return init;**

}

}

It can be used like this:

**void f(vector<int>& price, list<float>& incr){ int i = accumulate(price.begin(),price.end(),0); //**

*accumulate in int*

**double d = 0;**

d = accumulate(incr.begin(),incr.end(),d); //

d = accumulate(incr.begin(),incr.end(),d); //

*accumulate in double*

**int prod = accumulate(price.begin,price.end(),1,[](int a, int b) { return a*b; });**

//

//

*...*

**}**

The type of the initial value passed determines the return type.

We can provide an initial value and an operation for “combining elements” as arguments to **accumulate()**, so **accumulate()** is not all about addition.

Extracting a value from a data structure is common operation for **accumulate()**. For example:

**struct Record {**

**//** *...*

**int unit_price; int number_of_units;};**

**long price(long val, const Record& r)**

{

return val + r.unit_price * r.number_of_units;

}

{

return val + r.unit_price * r.number_of_units;

}

**void f(const vector<Record>& v)**

{

cout << "Total value: " << accumulate(v.begin(),v.end(),0,price) << '\n';

}

{

cout << "Total value: " << accumulate(v.begin(),v.end(),0,price) << '\n';

}

Operations similar to **accumulate** are called **reduce**, **reduction**, and **fold** in some communities.

**40.6.2. inner_product()**

Accumulating from a sequence is very common, and accumulating from a pair of sequences is not uncommon:

**template<typename In, typename In2, typename T>T inner_product(In first, In last, In2 first2, T init){ while (first != last) init = init + *first++ **first2++; return init;}**

**template<typename In, typename In2, typename T, typename BinOp, typename BinOp2>**

T inner_product(In first, In last, In2 first2, T init, BinOp op, BinOp2 op2)

{

while (first != last)

init = op(init,op2(*first++,*first2++));

return init;

}

T inner_product(In first, In last, In2 first2, T init, BinOp op, BinOp2 op2)

{

while (first != last)

init = op(init,op2(*first++,*first2++));

return init;

}

As usual, only the beginning of the second input sequence is passed as an argument. The second input sequence is assumed to be at least as long as the first.

The key operation in multiplying a **Matrix** by a **valarray** is an **inner_product**:

**valarray<double> operator*(const Matrix& m, valarray<double>& v){ valarray<double> res(m.dim2());**

**for (size_t i=0; i<m.dim2(); i++) {**

auto& ri = m.row(i);

res[i] = inner_product(ri,ri.end(),&v[0],double(0));

}

return res;

}

auto& ri = m.row(i);

res[i] = inner_product(ri,ri.end(),&v[0],double(0));

}

return res;

}

**valarray<double> operator*(valarray<double>& v, const Matrix& m)**

{

valarray<double> res(m.dim1());

{

valarray<double> res(m.dim1());

**for (size_t i=0; i<m.dim1(); i++) {**

auto& ci = m.column(i);

res[i] = inner_product(ci,ci.end(),&v[0],double(0));

}

return res;

}

auto& ci = m.column(i);

res[i] = inner_product(ci,ci.end(),&v[0],double(0));

}

return res;

}

Some forms of **inner_product** are referred to as “dot product.”

**40.6.3. partial_sum() and adjacent_difference()**

The **partial_sum()** and **adjacent_difference()** algorithms are inverses of each other and deal with the notion of incremental change.

Given a sequence **a**, **b**, **c**, **d**, etc., **adjacent_difference()** produces **a**, **b–a**, **c–b**, **d–c**, etc.

Consider a vector of temperature readings. We could transform it into a vector of temperature changes like this:

**vector<double> temps;****void f(){ adjacent_difference(temps.begin(),temps.end(),temps.begin());}**

For example, **17**, **19**, **20**, **20**, **17** turns into **17**, **2**, **1**, **0**, **–3**.

Conversely, **partial_sum()** allows us to compute the end result of a set of incremental changes:

**template<typename In, typename Out, typename BinOp>Out partial_sum(In first, In last, Out res, BinOp op){ if (first==last) return res;**

***res = *first;**

**T val = *first;**

while (++first != last) {

val = op(val,*first);

while (++first != last) {

val = op(val,*first);

***++res = val;**

}

return ++res;

}

}

return ++res;

}

**template<typename In, typename Out>**

**Out partial_sum(In first, In last, Out res)**

{

return partial_sum(first,last,res,plus); //

{

return partial_sum(first,last,res,plus); //

*use std::plus (§33.4)*

**}**

Given a sequence **a**, **b**, **c**, **d**, etc., **partial_sum()** produces **a**, **a+b**, **a+b+c**, **a+b+c+d**, etc. For example:

**void f(){ partial_sum(temps.begin(),temps.end(),temps.begin());}**

Note the way **partial_sum()** increments **res** before assigning a new value through it. This allows **res** to be the same sequence as the input; **adjacent_difference()** behaves similarly. Thus,

**partial_sum(v.begin(),v.end(),v.begin());**

turns the sequence **a**, **b**, **c**, **d** into **a**, **a+b**, **a+b+c**, **a+b+c+d**, and

**adjacent_difference(v.begin(),v.end(),v.begin());**

reproduces the original value. Thus, **partial_sum()** turns **17**, **2**, **1**, **0**, **–3** back into **17**, **19**, **20**, **20**, **17**.

For people who think of temperature differences as a boring detail of meteorology or science lab experiments, I note that analyzing changes in stock prices or sea levels involves exactly the same two operations. These operations are useful for analyzing any series of changes.

**40.6.4. iota()**

A call **iota(b,e,n)** assigns **n+i** to the **i**th element of [**b**:**e**). For example:

**vector<int> v(5);iota(v.begin(),v.end(),50);vector<int> v2 {50,51,52,53,54}**

**if (v!=v2)**

error("complain to your library vendor");

error("complain to your library vendor");

The name **iota** is the Latin spelling of the Greek letter , which was used for that function in APL.

Do not confuse **iota()** with the non-standard, but not uncommon, **itoa()** (**int**-to-alpha; §*12.2.4*).

**40.7. Random Numbers**

Random numbers are essential to many applications, such as simulations, games, sampling-based algorithms, cryptography, and testing. For example, we might want to choose the TCP/IP address for a router simulation, decide whether a monster should attack or scratch its head, or generate a set of values for testing a square root function. In **<random>**, the standard library defines facilities for generating (pseudo-)random numbers. These random numbers are sequences of values produced according to mathematical formulas, rather than unguessable (“truly random”) numbers that could be obtained from a physical process, such as radioactive decay or solar radiation. If the implementation has such a truly random device, it will be represented as a **random_device** (§*40.7.1*).

Four kinds of entities are provided:

• A *uniform random number generator* is a function object returning unsigned integer values such that each value in the range of possible results has (ideally) equal probability of being returned.

• A *random number engine* (an engine) is a uniform random number generator that can be created with a default state **E{}** or with a state determined by a **seed E{s}**.

• A *random number engine adaptor* (an adaptor) is a random number engine that takes values produced by some other random number engine and applies an algorithm to those values in order to deliver a sequence of values with different randomness properties.

• A *random number distribution* (a distribution) is a function object returning values that are distributed according to an associated mathematical probability density function *p*(*z*) or according to an associated discrete probability function *P*(*zi*).

For details see §iso.26.5.1.

In simpler terms, the users’ terms, a random number generator is an engine plus a distribution. The engine produces a uniformly distributed sequence of values, and the distribution bends those into the desired shape (distribution). That is, if you take lots of numbers from a random number generator and draw them, you should get a reasonably smooth graph of their distribution. For example, binding a **normal_distribution** to the **default_random_engine** gives me a random number generator that produces a normal distribution:

**auto gen = bind(normal_distribution<double>{15,4.0},default_random_engine{});****for (int i=0; i<500; ++i) cout << gen();**

The standard-library function **bind()** makes a function object that will invoke its first argument given its second argument (§*33.5.1*).

Using ASCII graphics (§*5.6.3*), I got:

**3 ******4 *****5 *********6 ********7 ********8 **********9 ****************10 *******************************11 *******************************12 ******************************************13 **************************************************************14 *******************************************************15 ***********************************************************16 ************************************17 ******************************************************18 *****************************************19 **************************************20 *******************21 ****************22 ****************23 ***********24 *********25 ********26 *****27 ***

Most of the time, most programmers just need a simple uniform distribution of integers or floating-point numbers in a given range. For example:

**void test(){**

**Rand_int ri {10,20}; //**

*uniform distribution of ints in [10:20)*

**Rand_double rd {0,0.5}; //**

*uniform distribution of doubles in [0:0.5)*

**for (int i=0; i<100; ++i)**

cout << ri() << ' ';

for (int i=0; i<100; ++i)

cout << rd() << ' ';

}

cout << ri() << ' ';

for (int i=0; i<100; ++i)

cout << rd() << ' ';

}

Unfortunately, **Ran_int** and **Rand_double** are not standard classes, but they are easy to build:

**class Rand_int {**

**Rand_int(int lo, int hi) : p{lo,hi} { } //** *store the parameters*

**int operator()() const { return r(); }private: uniform_int_distribution<>::param_type p; auto r = bind(uniform_int_distribution<>{p},default_random_engine{});};**

I store the parameters using the distribution’s standard **param_type** alias (§*40.7.3*) so that I can use **auto** to avoid having to name the result of the **bind()**.

Just for variation, I use a different technique for **Rand_double**:

**class Rand_double {public:**

**Rand_double(double low, double high)**

**:r(bind(uniform_real_distribution<>(low,high),default_random_engine())) { }**

double operator()() { return r(); }

private:

function<double()> r;

};

double operator()() { return r(); }

private:

function<double()> r;

};

One important use of random numbers is for sampling algorithms. In such algorithms we need to choose a *sample* of some size from a much larger *population*. Here is algorithm R (the simplest algorithm) from a famous old paper [Vitter,1985]:

**template<typename Iter, typename Size, typename Out, typename Gen>Out random_sample(Iter first, Iter last, Out result, Size n, Gen&& gen){ using Dist = uniform_int_distribution<Size>; using Param = typename Dist::param_type;**

**//**

*Fill the reservoir and advance first:*

**copy(first,n,result);**

advance(first,n);

advance(first,n);

**//**

*Sample the remaining values in [first+n:last) by selecting a random*

**//**

*number r in the range [0:k], and, if r<n, replace it.*

**//**

*k increases with each iteration, making the probability smaller.*

**//**

*For random access iterators, k = i-first (assuming we increment i and not first).*

**Dist dist;**

for (Size k = n; first!=last; ++first,++k) {

for (Size k = n; first!=last; ++first,++k) {

**Size r = dist(gen,Param{0,k});**

if(r < n)

if(r < n)

***(result + r) = *first;**

}

return result;

}

}

return result;

}

**40.7.1. Engines**

A uniform random number generator is a function object that produces an approximately uniformly distributed sequence of values of its **result_type**:

A random number engine is a uniform random number generator with additional properties to make it widely useful:

A seed is a value in the range [0:2^{32}) that can be used to initialize a particular engine. A seed sequence, **g**, is an object that provides a function **g.generate(b,e)** that when called fills [**b**:**e**) with newly generated seeds (§iso.26.5.1.2).

The **UI** parameter for a standard random number engine must be an unsigned integer type. For **linear_congruential_engine<UI,a,c,m>**, if the modulus **m** is **0**, the value **numeric_limits<result_type>::max()+1** is used. For example, this writes out the index of the first repetition of a number:

**map<int,int> m;linear_congruential_engine<unsigned int,17,5,0> linc_eng;for (int i=0; i<1000000; ++i) if (1<++m[linc_eng()]) cout << i << '\n';**

I was lucky; the parameters were not too bad and I got no duplicate values. Try **<unsigned int,16,5,0>** instead and see the difference. Use the **default_random_engine** unless you have a real need and know what you are doing.

A *random number engine adaptor* takes a random number engine as an argument and produces a new random number engine with different randomness properties.

For example:

**independent_bits_engine<default_random_engine,4,unsigned int> ibe;for (int i=0; i<100; ++i) cout << '0'+ibe() << ' ';**

This will produce 100 numbers in the range [48:63] ([’0’:’0’+2^{4}-1)).

A few aliases are defined for useful engines:

**using minstd_rand0 = linear_congruential_engine<uint_fast32_t, 16807, 0, 2147483647>;using minstd_rand = linear_congruential_engine<uint_fast32_t, 48271, 0, 2147483647>;using mt19937 = mersenne_twister_engine<uint_fast32_t, 32,624,397,**

**31,0x9908b0df,**

11,0xffffffff,

7,0x9d2c5680,

15,0xefc60000,

18,1812433253>

11,0xffffffff,

7,0x9d2c5680,

15,0xefc60000,

18,1812433253>

**using mt19937_64 = mersenne_twister_engine<uint_fast64_t, 64,312,156,**

**31,0xb5026f5aa96619e9,**

29, 0x5555555555555555,

17, 0x71d67fffeda60000,

37, 0xfff7eee000000000,

43, 6364136223846793005>;

using ranlux24_base = subtract_with_carry_engine<uint_fast32_t, 24, 10, 24>;

using ranlux48_base = subtract_with_carry_engine<uint_fast64_t, 48, 5, 12>;

using ranlux24 = discard_block_engine<ranlux24_base, 223, 23>;

using ranlux48 = discard_block_engine<ranlux48_base, 389, 11>;

using knuth_b = shuffle_order_engine<minstd_rand0,256>;

29, 0x5555555555555555,

17, 0x71d67fffeda60000,

37, 0xfff7eee000000000,

43, 6364136223846793005>;

using ranlux24_base = subtract_with_carry_engine<uint_fast32_t, 24, 10, 24>;

using ranlux48_base = subtract_with_carry_engine<uint_fast64_t, 48, 5, 12>;

using ranlux24 = discard_block_engine<ranlux24_base, 223, 23>;

using ranlux48 = discard_block_engine<ranlux48_base, 389, 11>;

using knuth_b = shuffle_order_engine<minstd_rand0,256>;

**40.7.2. Random Device**

If an implementation is able to offer a truly random number generator, that source of random numbers is presented as a uniform random number generator called **random_device**:

Think of **s** as the name of a random number source, such as a Geiger counter, a Web service, or a file/device containing the record of a truly random source. The **entropy()** is defined as

for a device with *n* states whose respective probabilities are *P*_{0},..., *P*_{n}_{–1}. The entropy is an estimate of the randomness, the degree of unpredictability, of the generated numbers. In contrast to thermodynamics, high entropy is good for random numbers because that means that it is hard to guess subsequent numbers. The formula reflects the result of repeatedly throwing a perfect *n*-sided dice.

The **random_device** is intended to be useful for cryptograpic applications, but it would be against all rules of that kind of application to trust an implementation of **random_device** without first studying it closely.

**40.7.3. Distributions**

A random number distribution is a function object that, when called with a random number generator argument, produces a sequence of values of its **result_type**:

In the following tables, a template argument **R** means a real is required in the mathematical formula and **double** is the default. An **I** means that an integer is required and **int** is the default.

The *precondition* field specifies requirements on the distribution arguments. For example:

**uniform_int_distribution<int> uid1 {1,100}; //** *OK***uniform_int_distribution<int> uid2 {100,1}; //** *error: a>b*

The *default* field specifies default arguments. For example:

**uniform_real_distribution<double> urd1 {}; //** *use a==0.0 and b==1.0***uniform_real_distribution<double> urd2 {10,20}; //** *use a==10.0 and b==20.0***uniform_real_distribution<> urd3 {}; //** *use double and a==0.0 and b==1.0*

The *result* field specifies the range of the results. For example:

**uniform_int_distribution<> uid3 {0,5};default_random_engine e;for (int i=0; i<20; ++i) cout << uid3(e) << ' ';**

The range for **uniform_int_distribution** is closed, and we see the six possible values:

**2 0 2 5 4 1 5 5 0 1 1 5 0 0 5 0 3 4 1 4**

For **uniform_real_distribution**, as for all other distributions with floating-point results, the range is half-open.

Bernoulli distributions reflect sequences of tosses of coins with varying degrees of loading:

Poisson distributions express the probability of a given number of events occurring in a fixed interval of time and/or space:

Normal distributions map real values into real values. The simplest is the famous “bell curve” that distributes values symmetrically around a peak (mean) with the distance of elements from the mean being controlled by a standard deviation parameter:

To get a feeling for these distributions, look at graphical representations for a variety of parameters. Such representations are easily generated and even more easily found on the Web.

Sampling distributions map integers into a specific range according to their probability density function **P**:

**40.7.4. C-Style Random Numbers**

In **<cstdlib>** and **<stdlib.h>**, the standard library provides a simple basis for the generation of random numbers:

**#define RAND_MAX implementation_defined /*** *large positive integer* ***/****int rand(); //** *pseudo-random number between 0 and RAND_MAX***void srand(unsigned int i); //** *seed random number generator by i*

Producing a good random number generator isn’t easy, and unfortunately not all systems deliver a good **rand()**. In particular, the low-order bits of a random number are often suspect, so **rand()%n** is not a good portable way of generating a random number between **0** and **n–1**. Often,**int((double(rand())/RAND_MAX)*n)** gives acceptable results. However, for serious applications, generators based on **uniform_int_distribution** (§*40.7.3*) will give more reliable results.

A call **srand(s)** starts a new sequence of random numbers from the *seed*, **s**, given as argument. For debugging, it is often important that a sequence of random numbers from a given seed be repeatable. However, we often want to start each real run with a new seed. In fact, to make games unpredictable, it is often useful to pick a seed from the environment of a program. For such programs, some bits from a real-time clock often make a good seed.

**40.8. Advice**

[1] Numerical problems are often subtle. If you are not 100% certain about the mathematical aspects of a numerical problem, either take expert advice, experiment, or do both; §*29.1*.

[2] Use variants of numeric types that are appropriate for their use; §*40.2*.

[3] Use **numeric_limits** to check that the numeric types are adequate for their use; §*40.2*.

[4] Specialize **numeric_limits** for a user-defined numeric type; §*40.2*.

[5] Prefer **numeric_limits** over limit macros; §*40.2.1*.

[6] Use **std::complex** for complex arithmetic; §*40.4*.

[7] Use **{}** initialization to protect against narrowing; §*40.4*.

[8] Use **valarray** for numeric computation when run-time efficiency is more important than flexibility with respect to operations and element types; §*40.5*.

[9] Express operations on part of an array in terms of slices rather than loops; §*40.5.5*.

[10] Slices is a generally useful abstraction for access of compact data; §*40.5.4*, §*40.5.6*.

[11] Consider **accumulate()**, **inner_product()**, **partial_sum()**, and **adjacent_difference()** before you write a loop to compute a value from a sequence; §*40.6*.

[12] Bind an engine to a distribution to get a random number generator; §*40.7*.

[13] Be careful that your random numbers are sufficiently random; §*40.7.1*.

[14] If you need genuinely random numbers (not just a pseudo-random sequence), use **random_device**; §*40.7.2*.

[15] Prefer a random number class for a particular distribution over direct use of **rand()**; §*40.7.4*.