Chapter 14: Numerics
14.1
C++ was not designed primarily with numeric computation in mind. (C was also not designed for numeric computation. Fortran was.) However, numeric computation typically occurs in the context of other work, so C++ becomes an attractive vehicle for computations that are part of a largeer system.
14.2
In <cmath>, we find the standard mathematical functions. They are reimlpementation of <math.h> from C library. After C++2x, these functions will becomes constexpr. The versions for comlpex are found in <complex>. Errors are reported by setting errno from <cerrno> to EDOM for a domain error and to ERANGE for a range error.
void f()
{
errno = 0; // clear old error state
sqrt(-1);
if (errno == EDOM){
cerr << "sqrt() not defined for negative argument":
}
errno = 0; // clear old error state
pow(numeric_limits<double>::max(), 2);
if (errno == ERANGE){
cerr << "result of power() too large to represents as a double";
}
}
A few more mathematical functions are found in <cstdlib> and the so-called special mathematical functions, such as beta(), rieman_zeta(), sph_bessel() are also in <cmath>.
14.3
In <numeric>, we find a small set of generalized numerical algorithsm, such as accumulate()
Numeric Algorithms
| algorithm | explanation |
| :--------------------------------- | :-------------------------------------------------------------------------------------------------------- |
| x=accumulate(b,e,i) | x is the sum of i and the elements of [b:e) |
| x=accumulate(b,e,i,f) | accumulate using f instead of + |
| x=inner_product(b,e,b2,i) | x is the inner product of [b:e) and [b2: b2+(e-b)), then add i |
| x=inner_product(b,e,b2,i,f,f2) | inner_product using f instead of + to aggregrate i, and use f2 instead of * for inner product |
| p=partial_sum(b,e,out) | element i of [out:p) is the sume of elements [b: b+i] |
| p=partial_sum(b,e,out,f) | partial_sum using f instead of + |
| p=adjacent_difference(b,e,out) | element i of [out:p) is *(b+i) - *(b+i-1) for i>0; if e-b>0, then *out is *b |
| p=adjacent_difference(b,e,out,f) | adjacent_difference using function f instead of - |
| iota(b,e,v) | for each element in [b:e) assign ++v; thus the sequence becomse v+1, v+2 |
|x=gcd(n,m)|xis the greatest common denominator of integersnandm|
|x=lcm(n,m)|xis the least common multiple of integersnandm` |
14.3.1
In <numeric>, the numerical algorithms have parallel versions that are slightly different
parallel numeric algorithms
| algorithm | explanation |
|---|---|
x=reduce(b,e,v) |
equivalent to x=accumulate(b,e,v), except out of order |
x=reduce(b,e) |
x=reduce(b,e, V{}), where V is b's value type |
x=reduce(pol,e,b,v) |
x=reduce(b,e,v) with execution policy pol |
x=reduce(pol,e,b) |
x=reduce(pol,b,e,V{}), where V is b's value type |
p=exclusive_scan(pol,b,e,out) |
p=partial_sum(b,e,out) according to the policy pol, excludes the ith input element from the ith sum |
p=inclusive_scan(pol,b,e,out) |
p=partial_sum(b,e,out) according to the policy pol, includes the ith input element from the ith sum |
p=transform_reduce(pol, b, e, f, v) |
f(x) for each x in [b:e), then reduce |
p=transform_exclusive_scan(pol,b,e,out,f,v) |
f(x) for each x in [b:e), then exclusive_scan |
p=transform_inclusive_scan(pol,b,e,out,f,v) |
f(x) for each x in [b:e), then inclusive_scan |
and there are more, such as those functions with functors
14.4
the standard library supports a family of complex number types along the lines of the comlex class. To support complex numbers where the scalars are single-precision floating-point numbers (floats), double-precision floating-point numbers (doubles), etc, the standard library complex is a template
template<typename Scalar>
class complex {
public:
comlpex(const Scalar& re={}, const Scalar& im={});
};
The usual arithmetic options are supported. sqrt() and power() functions are among the usual mathematical functions defined in <complex>.
14.5
In <random> library, 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 include uniform_int_distribution, normal_distribution, exponential_distribution.
using my_engine = default_random_engine;
using my_distribution = uniform_int_distribution<>;
my_engine re{};
my_distribution one_to_six{1, 6};
auto die = [](){ return one_to_six(re); } // make a generator
int x = die();
This solution is likely too general for novice, here is an example for a uniform random number generator in class syntax
class Rand_int{
public:
Rand_int(int low, int high): dist{low, high} {}
int operator()() {return dist(re); }
void seed(int s) { re.seed(s); }
private:
default_random_engine re;
uniform_int_distribution<> dist;
};
// use case
Rand_int rnd {1, 10}; // construct a random number generator from instantiating a class.
int x = rnd(); // generate a random number, x in [1:10]
The definition is still "expert level"; but the use is now manageble.
14.6
standard vector doesn't support mathematical 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. The usual arithmetic operations and the most common mathematical functions are supported for valarray.
template<typename T> // standard definitino
class valarray{
//...
};
void f(valarray<double>* a1, valarray<double>& a2){
valarray<double> a = a1 * 3.14 + a2 / a1; // element-wise operations
a2 += a2 * 3.14;
a = abs(a);
}
14.7
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(10000 < numeric_limits<int>::max(), "int type is too small");
17.8 Type Aliases from the 3rd edition
The size of fundamental types, such as int and long long are implementation defined. If we need to be specific about the size of our intergers, we can use aliases defined in <stdint>, such as int32_t and uint_least64_t (at least 64 bits).
The curious _t suffix is a relict from the days of C when it was deemed important to have a name reflect that it named an alias.
17.9 Mathematical Constants from the 3rd edition
When doing mathematical computations, we need common mathematical constants such as e, pi, and log2e. The standard library offers those and more. they come in two forms: a template that allows us to specify exact type (e.g. pi_v<T>) and a short name for the most common use (e.g. pi meaning pi_v<double>)
Naturally we can define type specific constants.