第14章 数值计算

计算的目的在于洞悉事物，而非数字本身。

—— 理查德·卫斯里·汉明

……但对学生而言，

数据往往是开启洞察力的最佳途径。

—— Anthony Ralston

14.1 导言

C++在设计之初并未做数值计算方面着重考虑过。 但数值计算却常穿插于其它业务——比如数据库访问、网络系统、仪器控制、图形学、 仿真及金融分析等等——因此对于较大系统中的计算部分，C++就成了香饽饽。 此外，数值方法早已远非遍历浮点数向量这样简单的任务。 在参与计算的数据结构日益复杂之处，C++的威力变得举足轻重。 这导致C++广泛用于科学、工程、金融及其它涉及复杂计算的领域。 因此，此类计算的辅助构件和技术则应运而生。 本章讲述标准库中有关数值计算的部分。

14.2 数学函数

头文件<cmath>提供了 标准数学函数（standard mathematical functions） ， 例如针对参数类型float、double以及long double的sqrt()、log()和sin()等：

标准数学函数
abs(x)	绝对值
ceil(x)	>=x的最小整数
floor(x)	<=x的最大整数
sqrt(x)	平方根；x不能是负数
cos(x)	余弦函数
sin(x)	正弦函数
tan(x)	正切函数
acos(x)	反余弦函数；结果不为负
asin(x)	反正弦函数；返回最靠近0的结果
atan(x)	反正切函数
sinh(x)	双曲正弦函数
cosh(x)	双曲余弦函数
tanh(x)	双曲正切函数
exp(x)	e(自然常数)的x次幂
log(x)	自然对数，以e为底；x必须是正数
log10(x)	以10为底的对数

针对complex（§14.4）的版本在<complex>中。 以上函数的返回值类型与参数相同。

错误报告的方式是将errno设置为<cerrno>中的值， 定义域超出范围设为EDOM，值域超出范围设为ERANGE。例如：

void f()
{
    errno = 0; // 清除错误状态
    sqrt(-1);
    if (errno==EDOM)
        cerr << "sqrt() not defined for negative argument";

    errno = 0; // 清除错误状态
    pow(numeric_limits<double>::max(),2);
    if (errno == ERANGE)
        cerr << "result of pow() too large to represent as a double";
}

有些被称为 特殊数学函数（special mathematical functions） 的数学函数在<cstdlib>里， 还有几个在<cmath>比如 beta()、rieman_zeta()、sph_bessel()。

14.3 数值算法

<numeric>里有几个泛化过的数值算法，比如accumulate()。

数值算法
x=accumulate(b,e,i)	x是i与[b:e)间元素的和
x=accumulate(b,e,i,f)	调用accumulate时用f替换+
x=inner_product(b,e,b2,i)	x是[b:e)与 [b2:b2+(e-b))的内积， 即i与(*p1) * (*p2)的和， 其中p1是[b:e) 中的元素，且对应来自[b2:b2+(e-b)) 中的元素p2
x=inner_product(b,e,b2,i,f,f2)	调用inner_product时用f 和f2分别替换+和*
p=partial_sum(b,e,out)	[out:p)的第i个元素是 [b:b+i]间所有元素的和
p=partial_sum(b,e,out,f)	调用partial_sum时以f替换+
p=adjacent_difference(b,e,out)	i>0时，[out:p) 的第i个元素是 *(b+i)-* (b+i-1)；e-b>0时， *out就是* b
p=adjacent_difference(b,e,out,f)	调用adjacent_difference时以f替换-
iota(b,e,v)	把++v依次赋值给[b:e) 之间的元素，因此元素序列就变成v+1,v+2,……
x=gcd(n,m)	x是整数n和m的最大公约数
x=lcm(n,m)	x是整数n和m的最小公倍数

这些算法泛化了常见运算，比如求和运算被应用到所有类型的元素序列上了。 也把应用在元素序列上的操作参数化了。 对于每个算法，最常规的版本是将常规运算代入到通用版本得到的。例如：

list<double> lst {1, 2, 3, 4, 5, 9999.99999};
auto s = accumulate(lst.begin(),lst.end(),0.0); // 求和得到：10014.9999

这些算法适用于标准库中的所有元素序列，可以将操作以参数的形式传入（§14.3）。

14.3.1 并行算法

<numeric>中，数值算法具有略带差异的并行版本（§12.9）：

并行数值算法
x=reduce(b,e,v)	无序执行的x=accumulate(b,e,v)
x=reduce(b,e)	x=reduce(b,e,V{})，其中V 是b的值类型
x=reduce(pol,b,e,v)	采用执行策略pol的x=reduce(b,e,v)
x=reduce(pol,b,e)	x=reduce(pol,b,e,V{})，其中V 是b的值类型
p=exclusive_scan(pol,b,e,out)	按照pol策略执行p=partial_sum(b,e,out) 计算第i个和的时候，第i个输入元素不参与计算
p=inclusive_scan(pol,b,e,out)	按照pol策略执行p=partial_sum(b,e,out) 计算第i个和的时候，第i个输入元素参与计算

并行数值算法（续表）
p=transform_reduce(pol,b,e,f,v)	对[b:e)中的每个x执行f(x)， 而后执行reduce
p=transform_exclusive_scan(pol,b,e,out,f,v)	对[b:e)中的每个x执行f(x)， 而后执行exclusive_scan
p=transform_inclusive_scan(pol,b,e,f,v)	对[b:e)中的每个x执行f(x)， 而后执行inclusive_scan

为简化叙述，此处没有提及那些采用仿函数参数替代+和=算法版本。 除reduce()意外，采用默认（顺序）执行策略和缺省值的版本也未提及。

此处的算法和<algorithm>里的并行算法一样，可以指定执行策略：

vector<double> v {1, 2, 3, 4, 5, 9999.99999};
auto s = reduce(v.begin(),v.end());     // 以double的默认值为初值累加求和

vector<double> large;
// ... 以大量的值填充large ...
auto s2 = reduce(par_unseq,large.begin(),large.end());  // 求和，并行策略可用则用，不可用就顺序执行

（reduce()之类的）并行算法区别于顺序版本（即accumulate()）之处在于： 并行算法中针对元素的操作执行顺序不定。

14.4 复数

标准库提供了一系列的复数类型，它们符合§4.2.1中complex的描述。 为了让其中的标量能支持单精度浮点数（float）、双精度浮点数（double）等类型， 标准库的complex是个模板：

template<typename Scalar>
class complex {
public:
    complex(const Scalar& re ={}, const Scalar& im ={});    // 函数参数缺省值；参见 §3.6.1
    // ...
};

复数支持常见的算术操作和大多数的数学函数。例如：

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

sqrt()和pow()（幂运算）属于<complex>中定义的常见数学函数。

14.5 随机数

许多领域需要随机数，比如测试、游戏、仿真以及安全系统。 标准库在<random>中提供了种类繁多的随机数发生器，它们反映了应用领域的多样性。 随机数发生器由两部分组成：

• [1] 引擎（engine） ，负责生成随机值或伪随机值的序列

• [2] 分布器（distribution） ，负责将这些值映射到特定的数学分布

分布器的例子有：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};                   // 映射到整数 1..6 的分布器
auto die = [](){ return one_to_six(re); }           // 创建一个随机数生成器

int x = die();                                      // 掷骰子：x的值在闭区间[1:6]内

出于对标准库中随机数组件通用性和性能的持续关注， 一位专家称其为“每个随机数程序库成长的榜样”。 但是要论“新手之友”的称号，它可就愧不敢当了。 前述代码示例借助using语句和lambda表达式，稍稍提升了一点代码可读性。

对于（任何背景的）新手而言，随机数程序库那个完整的通用接口绝对是个大坑。 一个简洁统一的随机数生成器往往就足以起步了。例如：

Rand_int rnd {1,10};    // 创建一个[1:10]之间的随机数生成器
int x = rnd();          // x是闭区间[1:10]内的一个值

但到哪儿去找这个东西呢？我们得弄个跟die()差不多的东西， 把引擎和分布器撮合起来，装进一个Rand_int类：

class Rand_int {
public:
    Rand_int(int low, int high) :dist{low,high} { }
    int operator()() { return dist(re); }       // 抽一个 int
    void seed(int s) { re.seed(s); }            // 设置新的随机数引擎种子
private:
    default_random_engine re;
    uniform_int_distribution<> dist;
};

这个定义仍然是“专家级”的，但是Rand_int()的 使用 ， 学习C++第一周的新手就可以轻松掌握了。例如：

int main()
{
    constexpr int max = 9;
    Rand_int rnd {0,max};           // 创建一个统一的随机数生成器

    vector<int> histogram(max+1);   // 创建一个容量合适的vector
    for (int i=0; i!=200; ++i)
        ++histogram[rnd()];         // 以[0:max]之间的数字作为频率填充直方图

    for (int i = 0; i!=histogram.size(); ++i) { // 绘制柱状图
        cout << i << '\t';
        for (int j=0; j!=histogram[i]; ++j) cout << '*';
        cout << endl;
    }
}

输出是个（索然无味的）均匀分布（具有合理的统计波动）。

0   *********************
1   ****************
2   *******************
3   ********************
4   ****************
5   ***********************
6   **************************
7   ***********
8   **********************
9   *************************

C++没有标准的图形库，所以这里用了“ASCII图形”。 毫无疑问，C++有众多开源以及商业的图形和GUI库， 但我在本书中限定自己仅使用ISO标准的构件。

14.6 矢量算术

§11.2叙述的vector设计目的是作为一个承载值的通用机制， 要灵活，并且能融入容器、迭代器、算法这套架构。 可惜它不支持数学矢量（vector）的运算。 给vector添加这些运算没什么难度， 但它的通用性和灵活性跟繁重的数值作业所需的优化格格不入。 因此，标准库（在<valarray>中）提供了一个类似vector的模板， 被称为valarray，它的通用性不济，但在数值计算所需的优化方面却精益求精：

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

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

void f(valarray<double>& a1, valarray<double>& a2)
{
    valarray<double> a = a1*3.14+a2/a1;         // 数值数组运算符 *、+、/、和=
    a2 += a1*3.14;
    a = abs(a);
    double d = a2[7];
    // ...
}

除算术运算之外，valarray还支持跨步访问，以辅助多维计算的实现。

14.7 数值范围

标准库在<limits>中提供了一些类，用于描述内建类型的属性—— 比如float指数部分的最大值，或者每个int的字节数。 比如说，可以断言char是有符号的：

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

请留意，第二个断言能够运行的原因，来（且仅来）自于numeric_limits<int>::max()是个constexpr函数（§1.6）这一事实。

14.8 忠告

• [1] 数值计算问题通常很微妙，如果对这类问题的某一方面没有100%的确信， 要么尝试咨询专家建议，要么实践检验，或者干脆双管齐下；§14.1。

• [2] 别把繁重的数值计算建立在编程语言的基础构件上，请采用程序库；§14.1。

• [3] 如果要从序列里计算出一个值，尝试写循环之前，请先考虑accumulate()、inner_product()、partial_sum()、adjacent_difference() ；§14.3。

• [4] 为复数运算采用std::complex；§14.4。

• [5] 把随机数引擎和一个分布器组合起来创建随机数生成器；§14.5。

• [6] 确保你的随机数足够随机；§14.5。

• [7] 别用C标准库的rand()；对于实际应用而言，它的随机度不够；§14.5。

• [8] 如果运行时的效率压倒操作和元素类型方面的灵活性， 请为数值计算采用valarray；§14.6。

• [9] 数值类型的属性可以通过numeric_limits获取；§14.7。

• [10] 请使用numeric_limits查询数值类型的属性，确保它们够用；§14.7。