17 数值计算

17.1 引言

当初设计C++语言时，数值计算并非其关注的重点。然而，数值计算在很多领域所发挥的作用，却日益显著，比如科学计算、数据库访问、网络安全、设备控制、图形图像处理、工业系统仿真、音视频、金融分析，等等。因此，C++业已成为许多大型系统中执行计算任务的强有力的工具之一。数值计算远不止于在简单的循环中处理浮点数序列，其间往往涉及非常复杂的数据结构和算法逻辑。所需数据结构和算法逻辑的复杂度越高，越能发挥出C++语言在数值计算方面的威力。迄今为止，C++已被广泛应用于科学、工程、金融等许多涉及复杂数值的计算任务中，而支持这类计算的各种功能和技术也逐渐发展起来，并成为标准库的一部分，以方便用户使用。

17.2 数学函数

标准库在<cmath>头文件中提供了许多数学函数，如参数类型为float、double和long double的sqrt、log和sin函数等，如下表所示：

这些函数的返回值类型与参数类型相同，其复数版本的声明位于<complex>头文件中。这些函数一旦出错，会将错误号设置在全局变量errno中，该全局变量的外部声明位于<cerrno>头文件中。定义域错误的错误号为EDOM，值域错误的错误号为ERANGE。这些函数在不发生错误的情况下，并不会自动清除错误号，因此需要在调用它们前手动清除错误号，以避免误判。例如：

xxxxxxxxxx
9
1
errno = 0; // 清除错误号
2
double x = sqrt(-1);
3
if (errno == EDOM)
4
    cerr << "sqrt not defined for negative argument" << endl;
5

6
errno = 0; // 清除错误号
7
double y = pow(mumeric_limits<double>::max(), 2);
8
if (errno == ERANGE)
9
    cerr << "result of pow too large to represent as a double" << endl;

在<cmath>和<cstdlib>头文件中，还有更多数学函数。尤其是<cmath>头文件，里面不乏一些特殊数学函数，如beta（$\beta$$\beta$函数，第一类欧拉积分）、rieman_zeta（黎曼$\zeta$$\zeta$函数）、sph_bessel（第一类球贝塞尔函数），等等。

C++17：特殊数学函数

17.3 数值算法

标准库在<numeric>头文件中提供了一系列泛型数值算法函数，如accumulate等。下表列出了其中的一部分：

其中的b、e、b1、e1、b2、e2都是迭代器或指向数组元素的指针，v为序列中元素类型的值，f、f1、f2都是函数指针、匿名函数、函数对象等可调用对象。

这些算法可以适配各种类型的序列，因而可以用于概括常见操作（如累加、内积、部分和、邻差、递增等）的一般形态，甚至可以将对元素的操作参数化（如f、f1、f2等）。对于每个算法而言，除了带有操作参数的通用版本（如accumulate(b,e,i,f)）以外，还有针对常见操作的简化版本（如accumulate(b,e,i)）。

17.3.1 并行数值算法

在<numeric>头文件中，除了提供前述串行数值算法外，还提供了一套并行数值算法。并行数值算法允许以未指定的顺序操作序列中的元素，并可以接受表示执行策略的参数，如seq、unseq、par、par_unseq等。下表列出了其中的一部分：

这里省略了带有操作参数的版本。

表示执行策略的参数，如seq、unseq、par、par_unseq等，位于<execution>头文件的std::execution命名空间中。

在决定使用并行数值算法前，最好先行测量，以评估这样做是否真的值得。

C++17：并行算法

17.4 复数

标准库提供了名为complex的复数类型，同时为了支持单精度浮点数（float）和双精度浮点数（double）类型的实虚部，将complex定义成一个类模板。类似下面这个样子：

xxxxxxxxxx
6
1
template<typename Scalar>
2
class complex {
3
public:
4
    complex(const Scalar& re = {}, const Scalar& im = {});
5
    ...
6
};

complex支持常见的算术运算和数学函数。例如：

xxxxxxxxxx
14
1
complex<float> c1{ 1.2f, 3.4f };
2
cout << c1 << endl; // (1.2,3.4)
3

4
complex<double> c2{ 5.6l, 7.8l };
5
cout << c2 << endl; // (5.6,7.8)
6

7
complex<long double> c3 = sqrt(c2);
8
cout << c3 << endl; // (2.757,1.41458)
9

10
auto c4 = pow(c3, 2);
11
cout << c4 << endl; // (5.6,7.8)
12

13
auto c5 = c4 + c3;
14
cout << c5 << endl; // (8.357,9.21458)

sqrt（平方根）、pow（幂）等数学函数的复数版本，声明于<complex>头文件中。

17.5 随机数

随机数在诸如测试、游戏、模拟、安全等领域中的应用颇多。为了满足多样化的应用需求，标准库在<random>头文件中提供了各种各样的随机数生成器。随机数生成器由两部分组成：

• 引擎：用于产生随机或伪随机的数值序列

• 分布：所产生的随机数服从哪种概率分布

  • 平均分布：uniform_int_distribution

  • 正态分布：normal_distribution

  • 指数分布：exponential_distribution

  • ······

在选择分布的同时，还可以指定随机数的生成范围。例如：

xxxxxxxxxx
9
1
default_random_engine eng; // 引擎
2
uniform_int_distribution dist{ 0, 5 }; // 分布
3

4
int dies[6]{};
5
for (int i = 0; i < 10000; ++i) // 掷一万次骰子
6
    ++dies[dist(eng)]; // 统计每个点数出现的次数
7

8
for (int i = 0; i < 6; ++i)
9
    cout << i + 1 << ": " << dies[i] << endl;

运行上面的程序，输出如下统计结果：

xxxxxxxxxx
6
1
1: 1680
2
2: 1613
3
3: 1648
4
4: 1712
5
5: 1662
6
6: 1685

标准库始终秉持通用性和性能毫不妥协的准则，因此它的随机数组件很难被视为“新手友好”的。适当使用类型别名和匿名函数，有望使代码的可读性有所提升。例如：

xxxxxxxxxx
13
1
using rengine = default_random_engine eng;
2
using uniform = uniform_int_distribution<>;
3

4
rengine eng; // 引擎
5
uniform dist{ 0, 5 }; // 分布
6
auto die = [&]() { return dist(eng); }; // 掷骰子
7

8
int dies[6]{};
9
for (int i = 0; i < 10000; ++i) // 掷一万次骰子
10
    ++dies[die()]; // 统计每个点数出现的次数
11

12
for (int i = 0; i < 6; ++i)
13
    cout << i + 1 << ": " << dies[i] << endl;

对于（任何背景的）新手而言，随机数库的通用接口有可能成为一个严重的障碍。定义一个简单的随机数工具类，不失为一种更好的选择。例如：

xxxxxxxxxx
17
1
class Rand {
2
public:
3
    Rand(int lower, int upper)
4
        : m_dist{ lower, upper } {}
5

6
    void seed(int seed) {
7
        m_eng.seed(seed);
8
    }
9

10
    int operator()() {
11
        return m_dist(m_eng);
12
    }
13

14
private:
15
    default_random_engine m_eng; // 引擎
16
    uniform_int_distribution<> m_dist; // 分布
17
};

而使用它生成随机数的代码会变得非常简单。例如：

xxxxxxxxxx
8
1
Rand rand{ 0, 5 };
2

3
int dies[6]{};
4
for (int i = 0; i < 10000; ++i) // 掷一万次骰子
5
    ++dies[rand()]; // 统计每个点数出现的次数
6

7
for (int i = 0; i < 6; ++i)
8
    cout << i + 1 << ": " << dies[i] << endl;

为了重复获得相同的随机数序列，或者为了确保每次产生的随机数序列不同，可以为引擎设置随机种子。种子相同序列相同，种子不同序列不同。例如：

xxxxxxxxxx
17
1
Rand r1{ 10, 99 };
2
r1.seed(1);
3
for (int i = 0; i < 10; ++i)
4
    cout << r1() << ' ';
5
cout << endl; // 47 99 74 93 10 21 37 99 23 31
6

7
Rand r2{ 10, 99 };
8
r2.seed(1);
9
for (int i = 0; i < 10; ++i)
10
    cout << r2() << ' ';
11
cout << endl; // 47 99 74 93 10 21 37 99 23 31
12

13
Rand r3{ 10, 99 };
14
r3.seed(2);
15
for (int i = 0; i < 10; ++i)
16
    cout << r3() << ' ';
17
cout << endl; // 49 26 12 93 59 95 49 53 47 38

重复序列对于确定性调试非常重要，而不重复序列则使伪随机数更相似于真随机数。当然，如果想得到真正的真随机数，则需要借助于random_device的支持。

C++11：随机分布和随机引擎

17.6 向量算术

vector被设计成一种通用机制，它可以存放不同类型的对象且足够灵活，能够适应容器、迭代器和算法的体系结构。它虽然名为vector（向量），但并不支持数学意义上的向量算术。为vector添加这类运算并不难，但其对通用性和灵活性的追求限制了为数值计算所需的优化。因此，标准库在<valarray>头文件中提供了名为valarray的模板类。与vector相比，valarray的通用性不强，但针对数值计算做了充分优化。valarray的定义形式类似下面这个样子：

xxxxxxxxxx
4
1
template<typename T>
2
class valarray {
3
    ...
4
};

valarray支持常见的算术运算和大多数数学函数。例如：

xxxxxxxxxx
19
1
valarray v1{ 0.1, 0.2, 0.3 };
2
for (double x : v1) cout << x << ' ';
3
cout << endl; // 0.1 0.2 0.3
4

5
valarray v2 = -v1;
6
for (double x : v2) cout << x << ' ';
7
cout << endl; // -0.1 -0.2 -0.3
8

9
valarray v3 = v1 + v2 * 3;
10
for (double x : v3) cout << x << ' ';
11
cout << endl; // -0.2 -0.4 -0.6
12

13
valarray v4 = abs(v3);
14
for (double x : v4) cout << x << ' ';
15
cout << endl; // 0.2 0.4 0.6
16

17
valarray v5 = sqrt(v4 / v1);
18
for (double x : v5) cout << x << ' ';
19
cout << endl; // 1.41421 1.41421 1.41421

这些操作都是向量操作，即作用于向量中的每个元素。

此外，valarray还支持跨步访问，以实现多维计算。

17.7 类型属性

在<limits>头文件中，标准库提供了名为numeric_limits的类模板，用于获取各种内置类型的属性，比如是否带有符号位，或者最大值是多少，等等。例如：

xxxxxxxxxx
1
1
static_assert(numeric_limits<T>::is_signed, "这个类型不带符号位！");

或者

xxxxxxxxxx
1
1
static_assert(numeric_limits<T>::max() > 100000, "这个类型的最大值不够大！");

甚至可以为自定义类型定义numeric_limits。

17.8 类型别名

基本类型，如short、int、long等的字长，是实现定义的。因此，它们在不同的C++实现中可能会不一致。如果需要显式指定整数的字长，可以使用定义在<stdint>头文件中的类型别名，如int16_t、int32_t、int64_t等。甚至可以用uint_least64_t定义至少64位字长的无符号整数。

形如“_t”形式的类型后缀，是源自C语言时代的历史遗迹。当时的人们认为，在名称上显式地反映出“这是一个类型的别名（而非别的什么东西）”很重要。

其它常见的类型别名，如size_t（sizeof操作符的结果类型）、ptrdiff_t（指针相减的结果类型）等，被收录在<stddef>头文件中。

C++11：整数类型别名，如int16_t、uint32_t、int_fast64_t等

17.9 数学常数

进行数学计算时，常常会用到一些数学常数，如$e$$e$（自然常数：$2.718281828459045\cdots$$2.718281828459045\cdots$）、$\pi$$\pi$（圆周率：$3.141592653589793\cdots$$3.141592653589793\cdots$）、$\gamma$$\gamma$（欧拉—马斯切罗尼常数：$0.577215664901532\cdots$$0.577215664901532\cdots$）、$\phi$$\varphi$（黄金分割比：$\frac{\sqrt{5}-1}{2}\approx 0.618$$\frac{\sqrt5-1}2\approx0.618$）等。标准库在<numbers>头文件的std::numbers命名空间中，为这些数学常数提供了两种形式的定义：

• 模板形式：通过模板参数显式指定确切类型，如pi_v<float>、pi_v<double>等

• 简略形式：不带模板参数的简短写法，如pi（相当于pi_v<double>）等

例如：

xxxxxxxxxx
3
1
float r = 2.34f; // 圆的半径
2
float c = 2 * pi_v<float> * r; // 圆的周长
3
double s = pi * r * r; // 圆的面积

从编程角度讲，float类型的$\pi$$\pi$即pi_v<float>，和double类型的$\pi$$\pi$即pi_v<double>，之间的差别很小，必须精确到小数点后16位左右才能看到，但在实际的科学计算中，这种差异所导致的结果却可能有着天壤之别。在诸如图形处理、人工智能等对数值精度要求较高的领域，浮点数的表达能力将变得越来越重要。

在<numbers>头文件的std::numbers命名空间中，定义了很多常用的数学常数，例如：

• e：自然常数$e=2.718281828459045\cdots$$e=2.718281828459045\cdots$

• pi：圆周率$\pi =3.141592653589793\cdots$$\pi=3.141592653589793\cdots$

• inv_pi：$\frac{1}{\pi }$$\frac1{\pi}$

• inv_sqrtpi：$\frac{1}{\sqrt{\pi }}$$\frac1{\sqrt\pi}$

• egamma：欧拉—马斯切罗尼常数$\gamma =0.577215664901532\cdots$$\gamma=0.577215664901532\cdots$

• phi：黄金分割比$\phi =\frac{\sqrt{5}-1}{2}\approx 0.618$$\varphi=\frac{\sqrt5-1}2\approx0.618$

• log2e：${\log }_{2}e$$\log_2e$

• log10e：${\log }_{10}e$$\log_{10}e$

• ln2：${\log }_{e}2$$\log_e2$

• ln10：${\log }_{e}10$$\log_e10$

• sqrt2：$\sqrt{2}$$\sqrt2$

• sqrt3：$\sqrt{3}$$\sqrt3$

• inv_sqrt3：$\frac{1}{\sqrt{3}}$$\frac1{\sqrt3}$

• ······

基于这些基本的数学常数，还可以组合出更多针对特定领域数学常数。例如：

xxxxxxxxxx
4
1
template<typename T>
2
constexpr T tau_v = 2 * pi_v<T>;
3

4
constexpr double tau = tau_v<double>;

之后，就可以在代码中直接使用tau_v或者tau，表示类似$2\pi$$2\pi$这样的常量了。

C++20：pi、log10e等数学常数

17.10 建议

• 数值问题通常都非常微妙，如果对某个问题中的数学含义不能百分之百地确定，一定要征询相关方面专家的建议、实验验证，或二者兼而有之

• 进行严肃的数学计算，一定要使用库，而不要试图直接用编程语言实现

• 如果想用循环从序列中计算某个结果，优先考虑使用accumulate、inner_product、partial_sum或者adjacent_difference等标准库函数

• 对于大量的数据，尝试并行和向量化算法，以从专门的优化中获益

• 用std::complex类表示复数，并参与运算

• 将引擎绑定到分布上以获得一个随机数发生器

• 确保随机数足够随机

• 不要再使用C语言标准库的rand函数产生随机数，对于实际用途来说，它显然不够随机

• 如果运行效率比操作和元素类型的灵活性更重要的话，应该使用valarray而非vector，参与向量计算

• 借助numeric_limits类模板，可以访问数值类型的属性

• 借助numeric_limits类模板，可以检查数值类型能否满足特定的计算需求

• 如果希望显式指定整型数据的字长，可以使用整数类型别名