Lesson 17: Numerics

Lesson 0017 — A Tour of C++, Chapter 17 (§17.1–§17.10)

C++ was not designed primarily with numeric computation in mind — but numeric computation typically occurs in the context of other work (scientific computation, database access, networking, instrument control, graphics, simulation, financial analysis), so C++ becomes an attractive vehicle for computations that are part of a larger system. Where more complex data structures are needed as part of a computation, C++'s strengths become relevant. This chapter describes the parts of the standard library that support numerics: mathematical functions, numerical algorithms, complex numbers, random numbers, vector arithmetic, numeric limits, type aliases, and mathematical constants. The thread running through it: the library already knows how to do numerics correctly — subtle, and easy to get wrong by hand.

17.1 Introduction

C++ is widely used for scientific, engineering, financial, and other computation involving sophisticated numerics; consequently, facilities and techniques supporting such computation have emerged. 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.

17.2 Mathematical Functions

In <cmath>, we find the standard mathematical functions, such as sqrt(), log(), and sin(), for arguments of type float, double, and long double:

The versions of these functions for complex numbers (§17.4) are found in <complex>. 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:

errno = 0;                              // clear old error state
double d = sqrt(-1);
if (errno==EDOM) cerr << "sqrt() not defined for negative argument\n";
errno = 0;                              // clear old error state
double dd = pow(numeric_limits<double>::max(),2);
if (errno == ERANGE) cerr << "result of pow() too large to represent as a double\n";

More mathematical functions are found in <cmath> and <cstdlib>. The so-called special mathematical functions, such as beta(), riemann_zeta(), and sph_bessel(), are also in <cmath>.

17.3 Numerical Algorithms

In <numeric>, we find a small set of generalized numerical algorithms, such as accumulate():

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

list<double> lst {1, 2, 3, 4, 5, 9999.99999};
auto s = accumulate(lst.begin(),lst.end(),0.0);   // calculate the sum: 10014.9999

These algorithms work for every standard-library sequence and can have operations supplied as arguments.

17.3.1 Parallel Numerical Algorithms

The numerical algorithms have parallel versions that differ slightly from the sequential ones — in particular, the parallel versions allow operations on elements in unspecified order. The parallel numerical algorithms can take an execution policy argument (§13.6): seq, unseq, par, and par_unseq, hidden in namespace std::execution in <execution>.

Just as for the parallel algorithms in <algorithm> (§13.6), we can specify an execution policy:

vector<double> v {1, 2, 3, 4, 5, 9999.99999};
auto s = reduce(v.begin(),v.end());                 // calculate the sum using a double as the accumulator
vector<double> large;                               // ... fill large with lots of values ...
auto s2 = reduce(par_unseq,large.begin(),large.end());   // calculate the sum using available parallelism

Measure to verify that using a parallel or vectorized algorithm is worthwhile.

17.4 Complex Numbers

The standard library supports a family of complex number types along the lines of the complex class described in §5.2.1. To support complex numbers where the scalars are single-precision (floats), double-precision (doubles), etc., the standard-library complex is a template:

template<typename Scalar> class complex {
public:
    complex(const Scalar& re ={}, const Scalar& im ={});   // default function arguments; see §3.4.1
    // ...
};

The usual arithmetic operations and the most common mathematical functions are supported for complex numbers:

void f(complex<float> fl, complex<double> db) {
    complex<long double> ld {fl+sqrt(db)};
    db += fl * 3;
    fl = pow(1/fl,2);
    // ...
}

The sqrt() and pow() (exponentiation) functions are among the usual mathematical functions defined in <complex> (§17.2).

17.5 Random Numbers

Random numbers are useful in many contexts, such as testing, games, simulation, and security. The diversity of application areas is reflected in the wide selection of random number generators provided by the standard library in <random>. A random number generator consists of two parts:

• An engine that produces a sequence of random or pseudo-random values;
• A distribution that maps those values into a mathematical distribution in a range.

Examples of distributions are uniform_int_distribution (where all integers produced are equally likely), normal_distribution ("the bell curve"), and exponential_distribution (exponential growth); each for some specified range:

using my_engine = default_random_engine;             // type of engine
using my_distribution = uniform_int_distribution<>;  // type of distribution
my_engine eng {};               // the default version of the engine
my_distribution dist {1,6};     // distribution that maps to the ints 1..6
auto die = [&](){ return dist(eng); };   // make a generator
int x = die();                  // roll the die: x becomes a value in [1:6]

Thanks to its uncompromising attention to generality and performance, one expert has deemed the standard-library random number component "what every random number library wants to be when it grows up" — however, it can hardly be deemed "novice friendly." The using statements and the lambda make what is being done a bit more obvious. For novices, a simple uniform random number generator is often sufficient to get started:

Rand_int rnd {1,10};   // make a random number generator for [1:10]
int x = rnd();          // x is a number in [1:10]

To get that, we combine an engine with a distribution inside a class:

class Rand_int {
public:
    Rand_int(int low, int high) :dist{low,high} { }
    int operator()() { return dist(re); }    // draw an int
    void seed(int s) { re.seed(s); }         // choose new random engine seed
private:
    default_random_engine re;
    uniform_int_distribution<> dist;
};

That definition is still "expert level," but the use of Rand_int() is manageable in the first week of a C++ course for novices — e.g., a histogram of 200 draws from Rand_int rnd {0,max} produces a reassuringly boring uniform distribution.

To get a repeated or different sequence of values, we seed the engine — that is, give its internal state a new value:

Rand_int rnd {10,20};
for (int i = 0; i<10; ++i) cout << rnd() << ' ';   // 16 13 20 19 14 17 10 16 15 14
rnd.seed(999);
for (int i = 0; i<10; ++i) cout << rnd() << ' ';   // 11 17 14 19 20 13 20 14 16 19
rnd.seed(999);                                    // same seed, same sequence
for (int i = 0; i<10; ++i) cout << rnd() << ' ';   // 11 17 14 19 20 13 20 14 16 19

Repeated sequences are important for deterministic debugging; seeding with different values is important when we don't want repetition. If you need genuine random numbers, rather than a generated pseudo-random sequence, look to see how random_device is implemented on your machine.

17.6 Vector Arithmetic

The vector described in §12.2 was designed to be a general mechanism for holding values, to be flexible, and to fit into the architecture of containers, iterators, and algorithms — it does not support mathematical vector operations. Adding such operations to vector would be easy, but its generality and flexibility preclude optimizations that are often considered essential for serious numerical work. Consequently, the standard library provides (in <valarray>) a vector-like template, called valarray, that is less general and more amenable to optimization for numerical computation:

template<typename T> class valarray { /* ... */ };

void f(valarray<double>& a1, valarray<double>& a2) {
    valarray<double> a = a1 * 3.14+a2/a1;   // numeric array operators *, +, /, and =
    a2 += a1 * 3.14;
    a = abs(a);
    double d = a2[7];
    // ...
}

The operations are vector operations — applied to each element of the vectors involved. In addition to arithmetic operations, valarray offers stride access to help implement multidimensional computations.

17.7 Numeric Limits

In <limits>, the standard library provides classes that describe the properties of built-in types — such as the maximum exponent of a float or the number of bytes in an int:

static_assert(numeric_limits<char>::is_signed,"unsigned characters!");
static_assert(100000<numeric_limits<int>::max(),"small ints!");

The second assert works because numeric_limits<int>::max() is a constexpr function (§1.6). We can define numeric_limits for our own user-defined types.

17.8 Type Aliases

The size of fundamental types, such as int and long long, is implementation-defined — they may be different on different implementations of C++. If we need to be specific about the size of our integers, we can use aliases defined in <stdint>, such as int32_t and uint_least64_t (an unsigned integer with at least 64 bits). The curious _t suffix is a relic from the days of C when it was deemed important to have a name reflect that it named an alias. Other common aliases, such as size_t (the type returned by the sizeof operator) and ptrdiff_t (the type of the result of subtracting one pointer from another), can be found in <stddef>.

17.9 Mathematical Constants

When doing mathematical computations, we need common mathematical constants such as e, pi, and log2e. The standard library offers those and more, in <numbers>. They come in two forms: a template that allows us to specify the exact type (e.g., pi_v<T>) and a short name for the most common use (e.g., pi meaning pi_v<double>):

void area(float r) {
    using namespace std::numbers;        // this is where the mathematical constants are kept
    double d = pi * r * r;
    float f = pi_v<float> * r * r;
    // ...
}

Here the difference is small (we would have to print with precision 16 or so to see it), but in real physics calculations such differences quickly become significant. Other areas where the precision of constants matters are graphics and AI, where smaller representations of values are increasingly important. In <numbers> we find e (Euler's number), log2e, log10e, pi, inv_pi (1/pi), inv_sqrtpi (1/sqrt(pi)), ln2, ln10, sqrt2, sqrt3, inv_sqrt3, egamma (the Euler-Mascheroni constant), and phi (the golden ratio). Naturally, we would like more constants and constants for different domains — that's easily done because such constants are variable templates with specializations for double (or whatever type is most useful for a domain):

template<typename T> constexpr T tau_v = 2 * pi_v<T>;
constexpr double tau = tau_v<double>;

17.10 Advice

Here is a summary of the guidance from this chapter. All 12 items, with the section where each is introduced. The C++ Core Guidelines link each item to its recommended practice.

Retrieval Quiz

Notes

Primary source: Stroustrup, B. (2022). A Tour of C++, 3rd ed., Chapter 17: "Numerics." Addison-Wesley.
Reference: Chapter 1–12 Quick Reference & Glossary — keep it beside you while you study.
Recommended supplement: cppreference on <cmath>, <numeric>, and random number generation.

Questions? Ask your agent — your teacher — about anything unclear: why reduce can run out of order, when to seed a random engine, or which numeric type is big enough for your data. Follow-ups are expected, not optional.

Chapters

1. Lesson 1: The Basics (Ch. 1)
2. Lesson 2: User-Defined Types (Ch. 2)
3. Lesson 3: Separate Compilation (Ch. 3)
4. Lesson 4: Error Handling (Ch. 4)
5. Lesson 5: Classes (Ch. 5)
6. Lesson 6: Essential Operations (Ch. 6)
7. Lesson 7: Templates (Ch. 7)
8. Lesson 8: Concepts and Generic Programming (Ch. 8)
9. Lesson 9: Library Overview (Ch. 9)
10. Lesson 10: Strings and Regular Expressions (Ch. 10)
11. Lesson 11: Input and Output (Ch. 11)
12. Lesson 12: Containers (Ch. 12)
13. Lesson 13: Algorithms (Ch. 13)
14. Lesson 14: Ranges (Ch. 14)
15. Lesson 15: Pointers and Containers (Ch. 15)
16. Lesson 16: Utilities (Ch. 16)
17. Lesson 17: Numerics (Ch. 17 — this lesson)
18. Lesson 18: Concurrency (Ch. 18)
19. Lesson 19: History and Compatibility (Ch. 19)