C++编程语言：标准库：第40章——数值(Numerics)(Bjarne Stroustrup)

原创于 2025-12-07 21:17:32 发布 · 631 阅读

20 ·

CC 4.0 BY-SA版权

文章标签：

#c++ #算法 #开发语言 #数值运算

c++ 同时被 2 个专栏收录

80 篇文章

订阅专栏

VC++

37 篇文章

订阅专栏

第 40 章数值(Numerics)

40.1 引言(Introduction)

40.2 数值极值(Numerical Limits)

40.2.1 极值宏(Limit Macros)

40.3 标准数学函数(Standard Mathematical Functions)

40.4 复数complex(Complex Numbers)

40.5 数值数组：valarray(A Numerical Array: valarray)

40.5.1 构造函数和赋值(Constructors and Assignments)

40.5.2 数组下标(Subscripting)

40.5.3 数组运算(Operations)

40.5.4 分片( Slices )

40.5.5 slice_array

40.5.6 扩展的分片(Generalized Slices)

40.6 扩展数值算法(Generalized Numerical Algorithms)

40.6.1 accumulate()

40.6.2 inner_product()

40.6.3 partial_sum()和adjacent_difference()

40.6.4 itoa

40.7 随机数(Random Numbers)

40.7.1 引擎

40.7.2 随机设备(Random Device)

40.7.3 分布(Distributions)

40.7.4 C风格随机数(C-Style Random Numbers)

40.8 建议(Advice)

40.1 引言(Introduction)

C++ 最初的设计并非以数值计算为主要目标。然而，数值计算通常出现在其他工作中——例如数据库访问、网络通信、仪器控制、图形处理、仿真和金融分析——因此，C++ 成为大型系统中进行计算的理想工具。此外，数值方法已经发展到远非简单的浮点数数组循环。当计算需要更复杂的数据结构时，C++ 的优势就得以充分发挥。最终结果是，C++ 被广泛用于科学、工程、金融以及其他涉及复杂数值计算的领域。因此，支持此类计算的工具和技术也应运而生。本章将介绍标准库中支持数值计算的部分。我并不打算教授数值方法。数值计算本身就是一个引人入胜的课题。要理解它，你需要一门优秀的数值方法课程，或者至少一本优秀的教科书——而不仅仅是一本语言手册和教程。

除了这里描述的标准库功能外，第 29 章是数值编程的一个扩展示例：一个N 维矩阵。

40.2 数值极值(Numerical Limits)

要对数进行任何有趣的操作，我们通常需要了解内置数值类型的一些基本属性。为了让程序员能够最大限度地利用硬件，这些属性是由实现定义的，而不是由语言本身的规则固定的(§6.2.8)。例如，最大的int是多少？最小的正float是多少？当一个double赋值给一个float时，它是四舍五入还是截断？一个char有多少二进制位？

此类问题的答案由 <limits> 中提供的 numeric_limits 模板的特化版本提供。例如：

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<shor t>::min() || numeric_limits<shor t>::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_specializ ed) {

// limits infor mation is available for type Quad

}

每个特化版本都为其参数类型提供相关信息。因此，通用的 numeric_limits 模板仅仅是一组常量的符号句柄和 constexpr 函数：

template<typename T>

class numeric_limits {

public:

static const bool is_specialized = false; // is information available for numer ic_limits<T>?

// ... uninteresting defaults ...

};

真正的关键在于特化。标准库的每一个实现都为每一个基本数值类型(字符类型、整数类型、浮点类型和布尔类型)提供了 numeric_limits 的特化，但没有为任何其他可能的候选类型(例如 void、枚举或库类型(例如 complex<double>))提供特化。

对于像 char 这样的整数类型，我们只需要关注少数几个信息。以下是 numeric_limits<char> 的一个实现示例，其中 char 类型为 8 位有符号数：

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

};

这些函数是 constexpr 函数，因此可以在需要常量表达式的地方使用，而不会产生运行时开销。

numeric_limits 中的大多数成员旨在描述浮点数。例如，以下描述了 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 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; } // maximum rounding error

static constexpr float infinity() noexcept { return /* some value */; }

1162 Numerics Chapter 40

static constexpr float quiet_NaN() noexcept { return /* some value */; }

static constexpr float signaling_NaN() noexcept { return /* some value */; }

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 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 is_iec559 = true; // confor ms 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 float_round_style round_style = round_to_nearest;

};

请注意，min() 是最小的正归一化数，epsilon 是最小的正浮点数，使得 1+epsilon−1 大于 0 。

当定义类似于内置类型的标量类型时，最好也提供 numeric_limits 的适当特化版本。例如，如果我定义了一个四精度类型 Quad，用户理所当然地会期望我提供 numeric_limits<Quad>。相反，如果我使用非数值类型 Dumb_ptr，我期望 numeric_limits<Dumb_ptr<X>> 作为主要模板，并将 is_specialized 设置为 false，表示没有可用信息。

我们可以设想，numeric_limits 函数可以有一些特化，用于描述用户自定义类型的属性，而这些属性与浮点数本身关系不大。在这种情况下，通常最好使用描述类型属性的通用方法，而不是针对标准中未考虑的属性来专门定义 numeric_limits 函数。

40.2.1 极值宏(Limit Macros)

C++ 从 C 语言继承了描述整数属性的宏。这些宏位于 <climits> 中：

(译注：在 C , C++ 等编程语言中，常常有宏指令(macro-instruction)，简称宏(macros)。它是一组程序指令，压缩成更简单的表示形式，在编程时以单条指令的形式出现。macro- 词缀的词义为 “长(long), 异常大（abnormally large）, 关于大尺寸的(on a large scale)”，而汉语中宏指的就是 “长，大” 之意。而现在程序中的宏还指用一个名称代替一个不好记的数或字符串。)

整数极值宏(_iso.diff.library_, abbreviated)
CHAR_BIT	一个char中的位数(通常为8位)
CHAR_MIN	最小的 char 值(可能为负值)
CHAR_MAX	最大的 char 值(若char为有符号则通常为127，若为无符号则为255 )
INT_MIN	最小int值
LONG_MAX	最大 long 值

同时还提供了名称类似的宏，用于表示有signed char ，long long 等类型。

类似地，<cfloat> 和 <float.h> 定义了描述浮点数属性的宏：

浮点数极值宏(_iso.diff.library_, abbreviated)
FLT_MIN	最小正float值(例如，1.175494351e−38F )
FLT_MAX	最大float值(例如，3.402823466e+38F )
FLT_DIG	浮点数的精度小数位数(例如，6 )
FLT_MAX_10_EXP	浮点数的最大十进制指数(例如，38)
DBL_MIN	最小 double 值
DBL_MAX	最大 double 值,例如，1.7976931348623158e+308
DBL_EPSILON	满足 1.0+DBL_EPSILON!=1.0 的最小double

long double 也有类似的宏。

40.3 标准数学函数(Standard Mathematical Functions)

在 <cmath> 中，我们可以找到通常所说的标准数学函数：

浮点数极值宏(_iso.diff.library_, abbreviated)
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 为底的指数
log(x)	自然对数，以 e 为底；x 必须为正数。
log10(x)	以 10 为底的对数

该函数有多种版本，分别接受 float、double、long double 和 complex(§40.4)类型的参数。每个函数的返回类型都与其参数类型相同。

错误报告方式为：将 errno 值从 <cerrno> 设置为 EDOM(域错误)或 ERANGE(范围错误)。例如：

void f()

{

errno = 0; // clear old error state

sqrt(−1);

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";

}

由于历史原因，一些数学函数位于 <cstdlib> 而不是 <cmath> 中：

更多数学函数(§iso.26.8)
n2=abs(n)	绝对值；n 为 int、long 或 long long 类型；n2 与 n 类型相同。
n2=labs(n)	“长整型绝对值”；n 和 n2 是long
n2=llabs(n)	“长长整型绝对值”；n 和 n2 是long long
p=div(n,d)	p=div(n,d)
p=ldiv(n,d)	n 除以 d；p 是 {商，余数}；n 和 d 是long型数。
p=lldiv(n,d)	n 除以 d；p 是 {商，余数}；n 和 d 是long long型数。

之所以存在 l∗() 版本，是因为 C 语言不支持重载。ldiv() 函数的结果分别是 div_t、ldiv_t 和 lldiv_t 结构体。这些结构体包含成员 quot(表示商)和 rem(表示余数)，其类型由实现定义。

针对特殊数学函数，有一个单独的 ISO 标准 [C++Math,2010]。实现过程中可以将这些函数添加到 <cmath> 中：

数学特殊函数(可选)
assoc_laguerre()	assoc_legendre()	beta()	comp_ellint_1()
comp_ellint_2()	comp_ellint_3()	cyl_bessel_i()	cyl_bessel_j()
cyl_bessel_k()	cyl_neumann()	ellint_1()	ellint_2()
ellint_3()	expint()	hermite()	laguerre()
legendre()	riemann_z eta()	sph_bessel()	sph_legendre()
sph_neumann()

如果你不知道这些功能，你可能也用不到它们。

40.4 复数complex(Complex Numbers)

标准库提供了复数类型 complex<float>、complex<double> 和 complex<long double>。complex<Scalar>(其中 Scalar 是支持常用算术运算的其他类型)通常可以工作，但不能保证其可移植性。

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) { }

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

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

Scalar imag() const { return im; } // imaginar y par t

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

template<typename X>

complex(const complex<X>&);

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

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>&);

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

};

标准库complex并不能防止窄化：

complex<float> z1 = 1.33333333333333333; // narrows

complex<double> z2 = 1.33333333333333333;// narrows

z1=z2; //narrows

为防止意外窄化，请使用 {} 初始化：

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

除了 complex 类的成员之外，<complex> 还提供了一系列有用的操作：

Complex运算符
z1+z2	加
z1−z2	减
z1∗z2	乘
z1/z2	除
z1==z2	相等
z1!=z2	不等
norm(z)	abs(z) 的平方
conj(z)	共轭：{z.re,−z.im}
polar(x,y)	给定极坐标(ρ，θ)，构造复数。
real(z)	实部
imag(z)	虚部
abs(z)	距离 (0,0) 的距离：sqrt(z.re∗z.re+z.im∗z.im)；也称为ρ
arg(z)	从正实轴出发的角度：atan2(z.im/z.re)；也称为θ 。
out<<z	复数输出
in>>z	复数输入

标准数学函数(§40.3)也适用于复数。请注意，comeplex不支持小于号 (<) 或百分比 (%)。更多详情请参见 §18.3。

40.5 数值数组：valarray(A Numerical Array: valarray)

许多数值运算都依赖于相对简单的一维浮点值数组(vectors)。特别是，高性能机器架构对这类数组的支持非常出色，基于这类向量的库也得到了广泛应用，并且在许多领域，对使用这类数组的代码进行高度优化被认为是至关重要的。valarray 是 <valarray> 类中的一维数值数组。它为数组类型提供了常用的数值数组算术运算，并支持分片(slices)和步长(strides)：

数值数组类(§iso.26.6.1)
valarray<T>	类型为 T 的数值数组
slice	类似 BLAS 的分片(起始位置、长度和步幅)；§40.5.4
slice_array<T>	由分片标识的子数组；§40.5.5
gslice	推广分片用于描述矩阵
gslice_array<T>	由推广分片标识的子矩阵；§40.5.6
mask_array<T>	由掩码标识的数组子集；§40.5.2
indirect_array<T>	由索引列表标识的数组子集；§40.5.2

valarray 的基本理念是为稠密多维数组提供类似 Fortran 的功能，并具备类似 Fortran 的优化能力。这只有在编译器和优化库供应商的积极支持下，并在 valarray 提供的基本功能基础上增加更多库支持的情况下才能实现。然而，目前并非所有实现都做到了这一点。

40.5.1 构造函数和赋值(Constructors and Assignments)

valarray 构造函数允许我们使用辅助数值数组类型和单个值来初始化 valarray：

valarray<T> 构造函数(§iso.26.6.2.2)
valarray va{};	不包含任何元素的 valarray
valarray va {n};	包含 n 个元素且值为 T{} 的 valarray；explicit
valarray va {t,n};	包含 n 个元素且值为 t 的 valarray。
valarray va {p,n};	包含 n 个元素的 valarray，其值从 [p:p+n) 复制
valarray va {v2};	移动和复制构造函数
valarray va {a};	使用数组 a 中的元素构造数组 va； a 可以是 slice_array、gslice_array、mask_array 或 indirect_array；元素个数等于数组 a 中的元素个数。
valarray va {args};	根据初始化列表 {args} 构建；元素个数等于 {args} 中的元素个数。
va.˜valarray()	析构函数

例如：

valarray<double> v0; // 点位, 我们可以在后面再对 v0 赋值

valarray<float> v1(1000); // 1000 个值为0.0F 的元素

valarray<int> v2(−1,2000); // 2000 个值为 −1的元素

valarray<double> v3(100,9.8064); // 严重错误:用浮点作为 valarray 的大小

valarray<double> v4 = v3; // v4 有 v3.size() 个元素

valarray<int> v5 {−1,2000}; // 2个元素

在双参数构造函数中，值位于元素数量之前。这与标准容器的约定不同(§31.3.2)。

复制构造函数的参数 valarray 的元素数量决定了生成的 valarray 的大小。

大多数程序都需要从表格或输入中获取数据。除了初始化列表之外，还可以通过构造函数从内置数组复制元素。例如：

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 个元素: 0,1,2,3

valarray<double> v2{vi,4}; // 类型错误: vi 不是一个指向 double 的指针

valarray<double> v3{vd,8}; // 未定义: 初始化器中元素太少

valarray<int> v4{p,n}; // p 最好至少指向 n 个整数。

}

valarray 及其辅助功能专为高速计算而设计。这体现在对用户的一些限制以及对实现者的一些自由度上。基本上，valarray 的实现者可以使用几乎所有你能想到的优化技术。valarray 操作被假定为无副作用(当然，显式参数除外)，valarray 假定为无别名，并且允许引入辅助类型和消除临时变量，只要保持基本语义即可。没有范围检查。valarray 的元素必须具有默认的复制语义(§8.2.6)。

赋值对象可以是另一个 valarray、一个标量或 valarray 的子集：

valarray<T> 赋值(§iso.26.6.2.2)
va2=va	复制赋值：va2.size() 变为 va.size()
va2=move(va)	移动分配：va 变为空
va=t	标量赋值：va 的每个元素都是 t 的副本。
va={args}	从初始化列表 {args} 中赋值； va 的元素个数变为 {args}.size()
va=a	从 a 赋值；a.size() 必须等于 va.size()； a 可以是 slice_array，gslice_array，mask_array 或 indirect_array。
va@=va2	对于 va 的每一个元素，v[i]@=va2[i]；@ 可以是 /，%，+，−，ˆ，&，\|，<< 或 >>
va@=t	对于 va 中的每一个元素，v[i]@=t；@ 可以是 /，%，+，−，ˆ，&，\|，<< 或 >>

可以将一个 valarray 赋值给另一个大小相同的 valarray。正如预期的那样，v1=v2 会将 v2 中的每一个元素复制到 v1 中对应的位置。如果 valarray 的大小不同，则赋值的结果未定义。

除了这种常规赋值之外，还可以将标量赋值给 valarray。例如，v=7 会将 7 赋值给 valarray v 的每一个元素。这可能会让一些程序员感到惊讶，最好将其理解为运算符赋值操作的一种偶尔有用的退化情况。例如：

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)

可以使用下标来选择 valarray 中的一个元素或其元素的子集：

valarray<T> Subscripting (§iso.26.6.2.4, §iso.26.6.2.5)
t=va[i]	下标：t 是指向 va 的第 i 个元素的引用；不进行范围检查
a2=va[x]	子集：x 是一个分片，一个 gslice，一个 valarray<bool> 或一个 valarray<size_t>

每一个运算符 [] 都返回一个valarray元素的子集。返回类型(表示该子集的对象类型)取决于参数类型。

对于 const 参数，结果包含元素的副本。对于非常量参数，结果包含对元素的引用。由于 C++ 不直接支持引用数组(例如，我们不能写成 valarray<int&> )，实现会以某种方式模拟这种行为。这可以高效地完成。有必要提供一个包含示例的详尽列表(基于 §iso.26.6.2.5)。在每种情况下，下标描述要返回的元素，并且 v1 必须是具有适当长度和元素类型的 vallaray ：

• 一个const valarray 的一个 slice:

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

// ...

const valarray<char> v0 {"abcdefghijklmnop",16};

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

• 一个非const valarray 的一个 slice:

slice_array<T> operator[](slice); //元素的引用

// ...

valarray<char> v0 {"abcdefghijklmnop",16};

valarray<char> v1 {"ABCDE",5};

v0[slice(2,5,3)] = v1; // v0=={"abAdeBghCjkDmnEp",16}

• 一个const valarray 的一个 gslice:

valarray<T> operator[](const gslice&) const; // 元素复制

// ...

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}

• 一个非const valarray 的一个 gslice:

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}

• 一个const valarray 的一个 valarray<bool>( 一个掩码 ):

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}

• 一个非const valarray 的一个 valarray<bool>( 一个掩码 ):

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}

• 一个const valarray 的一个 valarray<size_t>( 一个索引集 ):

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}

• 一个非const valarray 的一个 valarray<size_t>( 一个索引集 ):

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}

请注意，使用掩码(valarray<bool>)进行下标运算会产生 mask_array，而使用一组索引( valarray<size_t> )进行下标运算会产生 indirect_array 。

40.5.3 数组运算(Operations)

valarray 的目的是支持计算，因此它直接支持许多基本的数值运算：

valarray<T> 成员运算 (§iso.26.6.2.8)
va.swap(va2)	下标：t 是指向 va 的第 i 个元素的引用；不进行范围检查
n=va.size()	子集：x 是一个分片，一个 gslice，一个 valarray<bool> 或一个 valarray<size_t>
t=va.sum()	t 是使用 += 计算的 va 元素的总和
t=va.min()	t 是使用 < 找到的 va 的最小元素
t=va.max()	t 是使用 < 找到的 va 的最大元素
va2=va.shift(n)	元素线性移动，n > 0 左移，n < 0 右移
va2=va.cshift(n)	元素循环移动，n > 0 左移，n < 0 右移
va2=va.apply(f)	应用函数 f ：每个元素 va2[i] 的值为 va[f(i)]
va.resize(n,t)	创建一个包含 n 个元素的 valarray，其值为 t。
va.resize(n)	va.resiz e(n,T{})

没有范围检查：使用尝试访问空 valarray 元素的函数的效果是未定义的。

请注意，resize() 函数不会保留任何旧值。

valarray<T> 运算 (§iso.26.6.2.6，§iso.26.6.2.7) v 或 v2(但不能同时是两者)可以是标量；对于算术运算，结果是一个 valarray<T>。
swap(va,va2)	va.swap(va2)
va3=va@va2	对 va 和 va2 的元素执行 @ 操作，得到 va3； @ 可以是 + , −,∗ , / , % , & , \| , ˆ , << , >> , && , \|\|
vb=v@v2	对 v 和 v2 的元素执行 @ 操作，得到一个 valarray<bool>； @ 可以是 == , != , < , <= , > , >=
v2=@(v)	对 v 的元素执行 @() 操作，得到 v2； @ 可以是 abs , acos , asin , atan , cos , cosh , exp , log 或 log10。
v3=atan2(v,v2)	对 v 和 v2 的元素执行 atan2() 函数
v3=pow(v,v2)	对 v 和 v2 执行 pow() 操作
p=begin(v)	p 是一个随机访问迭代器，指向 v 的第一个元素。
p=end(v)	p 是一个随机访问迭代器，指向 v 的倒数第二个元素。

二元运算适用于 valarray 以及 valarray 与其标量类型的组合。标量类型视为大小合适的 valarray，其中每一个元素都具有该标量的值。例如：

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

}

这些数组(vector)运算都会对操作数的每一个元素进行运算，具体方式如 ∗ 和 cos() 示例所示。当然，只有当相应的运算已为标量类型定义时，才能使用该运算。否则，编译器在尝试特化运算符或函数时会报错。

如果结果是一个 valarray 类型，则其长度与其 valarray 操作数的长度相同。如果两个数组的长度不同，则对两个 valarray 类型执行二元运算符的结果未定义。

这些 valarray 操作会返回新的 valarray，而不是修改其操作数。这可能会增加性能开销，但如果应用积极的优化技术，则不必如此。

例如，如果 v 是一个 valarray 类型，它可以这样缩放：v∗=0.2，也可以这样缩放：v/=1.3。也就是说，将标量应用于数组意味着将该标量应用于向量的每一个元素。通常，∗= 比 ∗ 和 = 的组合更简洁(参见§18.3.1)，也更容易优化。

请注意，非赋值操作会创建一个新的值数组。例如：

double incr(double d) { return d+1; }

void f(valarray<double>& v)

{

valarray<double> v2 = v.apply(incr); // produce incremented valarray

// ...

}

这不会改变 v 的值。遗憾的是，apply() 不接受函数对象(§3.4.3，§11.4)作为参数。

逻辑移位函数 shift() 和循环移位函数 cshift() 会返回一个新的 valarray，其中元素已进行适当的移位，而原始 valarray 保持不变。例如，循环移位函数 v2=v.cshift(n) 会生成一个新的 valarray，使得 v2[i]==v[(i+n)%v.size()]。逻辑移位函数 v3=v.shift(n)会生成一个新的 valarray，使得如果 i+n 是 v 的有效索引，则 v3[i] 等于 v[i+n]。否则，结果为元素的默认值。这意味着，shift() 和 cshift() 在给定正参数时向左移位，在给定负参数时向右移位。例如：

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

}

对于 valarray 类型，>> 和 << 是位移运算符，而不是元素位移运算符或 I/O 运算符。因此，<<= 和 >>= 可以用于移动整数类型元素内的位。例如：

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

}

所有作用于 valarray 的运算符和数学函数也适用于 slice_arrays(§40.5.5)、gslice_arrays(§40.5.6)、mask_arrays(§40.5.2)、indirect_arrays(§40.5.2)以及这些类型的组合。但是，实现允许先将非 valarray 类型的操作数转换为 valarray 类型，然后再执行所需的操作。

40.5.4 分片( Slices )

slice是一种抽象概念，它允许我们高效地将一维数组(例如，内置数组、vector或 valarray)作为任意维数的矩阵进行操作。它是 Fortran 向量和 BLAS(基本线性代数子程序)库的关键概念，而 BLAS 库是许多数值计算的基础。简而言之，slice是数组某个部分的 n 个元素集合：

class std::slice {

// star ting index, a length, and a stride

public:

slice(); //slice{0,0,0}

slice(size_t star t, siz e_t siz e, siz e_t stride);

size_t star t() const; // index of first element

size_t siz e() const; // number of elements

size_t stride() const; // element n is at start()+n*str ide()

};

步长是指分片中两个元素之间的距离(以元素个数为单位)。因此，分片描述了非负整数到索引的映射。元素个数size() 函数不影响映射(寻址)，但允许我们找到序列的末尾。这种映射可以高效、通用且相当方便地在一维数组(例如valarray )中模拟二维数组。考虑一个 3×4 的矩阵(三行，每行四个元素)：

(译注：以下内部括号在代码中是不应该加的，加括号只是用于示意矩阵形式)

valarray<int> v {

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

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

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

};

或者图形化为：

按照通常的 C/C++ 约定，valarray 在内存中以行元素优先(行优先顺序)且连续的方式布局：

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

输出：

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

或图形化为：

第 x 行由 slice(x∗4,4,1) 描述。也就是说，第 x 行的第一个元素是向量的第 x∗4 个元素，该行的下一个元素是第 (x∗4+1) 个元素，依此类推，每行有 4 个元素。例如，slice{0,4,1} 描述向量 v 的第一行(数组第 0 行)：00，01，02，03，而 slice{1,4,1} 描述第二行(数组第1行)。

第 y 列由 slice(y,3,4) 描述。也就是说，第 y 列的第一个元素是向量的第 y 个元素，该列的下一个元素是第 (y+4) 个元素，依此类推，每列有 3 个元素。例如，slice{0,3,4} 描述第一列(第 0 列)：00 , 10 , 20，而 slice{1,3,4} 描述第二列(第 1 列)。

除了用于模拟二维数组之外，分片还可以描述许多其他序列。它是一种相当通用的指定简单序列的方法。这一概念将在 §40.5.6 中进一步探讨。

可以将分片理解为一种特殊的迭代器：分片允许我们描述一个值数组的索引序列。我们可以基于此构建一个 STL 风格的迭代器：

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:

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.star t();

}

bool operator!=(const Slice_iter& q ) const

{

return !(∗this==q);

}

bool operator<(const Slice_iter& q) const

{

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

}

};

由于slice具有大小，我们甚至可以进行范围检查。这里，我利用了 slice::size() 方法，提供了一个 end() 操作，用于返回slice末尾元素之后的下一个元素。

由于slice可以描述行或列，因此 Slice_iter 允许我们按行或按列遍历 valarray。

40.5.5 slice_array

我们可以利用 valarray 和 slice 构建一个看起来像 valarray 的东西，但实际上它只是引用 slice 所描述的数组子集的一种方式。

slice_array<T> (§iso.26.6.5)
slice_array sa {sa2};	复制构造函数：sa 引用与 sa2 相同的元素
sa=sa2	将 sa[i] 所引用的元素赋值给 sa2[i] 所引用的每一个对应元素。
sa=va	将 va[i] 分配给 sa[i] 引用的每一个对应元素。
sa=v	将 v 赋值给 sa 所引用的每一个元素。
sa@=va	sa[i]@=va[i] 对 sa 的每个元素； @ 可以是 /, %, +, −, ˆ, &, \|, <<, 或 >>

用户无法直接创建 slice_array。相反，用户需要使用 valarray 作为下标来为给定的分片创建 slice_array。slice_array 初始化后，所有对其的引用都会间接指向创建它的 valarray。例如，我们可以创建一个表示数组中每隔一个元素的 slice_array，如下所示：

void f(valarray<double>& d)

{

slice_array<double>& v_even = d[slice(0,d.size()/2+d.siz e()%2,2)];

slice_array<double>& v_odd = d[slice(1,d.siz e()/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

}

分片数组可以复制。例如：

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)

一个slice(§29.2.2，§40.5.4)可以描述一个n维数组的行或列。然而，有时我们需要提取的子数组并非行或列。例如，我们可能想要从一个4×3矩阵的左上角提取一个3×2矩阵：

遗憾的是，这些元素的分配方式无法用单个分片来描述：

gslice 是一个“扩展分片”，它包含了(几乎)n 个分片的信息：

class std::gslice {

//代替1个步长和1个像 slice 的大小 , gslice 存储 n 个步长和 n 个大小

public:

gslice();

gslice(size_t sz, const valarray<size_t>& lengths, const valarray<size_t>& strides);

size_t start()const; //第一个元素之索引

valarray<size_t> size() const; // 元素按维数的数量

valarray<size_t> stride() const; // index[0], index[1], ... 的步长

};

这些额外的参数允许 gslice 指定 n 个整数与用于寻址数组元素的索引之间的映射。例如，我们可以使用一对 (长度, 步长) 参数来描述 3×2 矩阵的布局：

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};// 第一维中的两个元素

// 第二维中的三个元素

valarray<size_t> strides {3,1}; // 3是第一个索引的步长

// 1是第二个索引的步长

void f()

{

gslice s(0,lengths,strides);

for (int i=0; i<3; ++i) // 对每一行

for (int j=0; j<2; ++j) // 对行中的每一个元素

cout << "(" << i << "," << j << ")−>" << gslice_index(s,i,j) << "; "; //打印映射

}

打印：

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

按这种方式，一个包含两个 (长度, 步长) 对的 gslice 描述了一个二维数组的子数组，一个包含三个 (长度, 步长) 对的 gslice 描述了一个三维数组的子数组，依此类推。将 gslice 用作 valarray 的索引，会生成一个 gslice_array，其中包含由该 gslice 描述的元素。例如：

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]

}

gslice_array 提供的成员集与 slice_array 相同(§40.5.5)。gslice_array 是将 gslice 用作 valarray 的下标的结果(§40.5.2)。

40.6 扩展数值算法(Generalized Numerical Algorithms)

在 <numeric> 中，标准库提供了一些通用的数值算法，其风格与 <algorithm>(第 32 章)中的非数值算法类似。这些算法提供了对数值序列进行常见运算的通用版本：

数值算法 (§iso.26.7)(续) 这些算法取输入迭代器
x=accumulate(b,e ,i)	x 是 i 与 [b:e) 中元素的和
x=accumulate(b,e ,i,f)	使用 f 而不是 + 进行累加
x=inner_product(b,e ,b2,i)	x 是 [b:e) 和 [b2:b2+(e−b)) 的内积，即对于 [b:e) 中的每一个 p1 和 [b2:b2+(e−b)) 中对应的 p2，i 与 (∗p1)∗(∗p2) 之和。
x=inner_product(b,e ,b2,i,f,f2)	使用 f 和 f2 代替 + 和 ∗ 进行内积运算
p=partial_sum(b,e,out)	[out:p) 中的元素 i 是元素 [b:b+i] 之和
p=partial_sum(b,e,out,f)	使用 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)	使用 f 而不是 − 来计算相邻差
iota(b,e,v)	对于 [b:e] 中的每个元素，赋值 ++v；因此序列变为 v+1, v+2, ...

这些算法通过使其适用于各种序列，并将应用于这些序列元素的运算作为参数，从而概括了诸如计算总和之类的常见运算。对于每一个算法，除了通用版本外，还有一个版本专门用于该算法中最常用的运算符。

40.6.1 accumulate()

accumulate() 函数的简单版本使用 + 运算符将序列中的元素相加：

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;

}

它可以这样使用：

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); // accumulate in double

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

// ...

}

传入的初始值的类型决定了返回值的类型。

我们可以将初始值和“合并元素”的操作作为参数传递给 accumulate()，因此 accumulate() 不仅仅是加法。

从数据结构中提取值是 accumulate() 函数的常见操作。例如：

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;

}

void f(const vector<Record>& v)

{

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

}

在某些社群中，与accumulate类似的操作称为 reduce，reduction和fold 。

40.6.2 inner_product()

从一个序列中累加非常常见，从一对序列中累积也并非罕见：

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;

}

与往常一样，只有第二个输入序列的开头部分作为参数传递。假定第二个输入序列的长度至少与第一个输入序列一样长。

将Matrix乘以valarray的关键操作是inner_product：

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

{

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

for (siz e_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;

}

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

{

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

for (siz e_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;

}

某些形式的内积称为“点积”。

40.6.3 partial_sum()和adjacent_difference()

partial_sum() 和 adjacent_difference() 算法互为逆算法，处理的是增量变化的概念。

给定一个序列 a , b , c , d 等，agnancy_difference() 生成 a , b−a , c−b , d−c 等。

考虑一个温度读数数组。我们可以将其转换为温度变化数组，如下所示：

vector<double> temps;

void f()

{

adjacent_difference(temps.begin(),temps.end(),temps.begin());

}

例如，17 , 19 , 20 , 20 , 17 变为 17 , 2 , 1 , 0 , -3。

相反，partial_sum() 函数允许我们计算一组增量变化的最终结果：

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);

∗++res = val;

}

return ++res;

}

template<typename In, typename Out>

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

{

return partial_sum(first,last,res,plus); // use std::plus (§33.4)

}

给定一个序列 a, b , c , d 等，partial_sum() 函数会生成 a , a+b , a+b+c, a+b+c+d 等结果。例如：

void f()

{

partial_sum(temps.begin(),temps.end(),temps.begin());

}

注意 partial_sum() 函数在给 res 赋新值之前会先递增它。这使得 res 可以与输入序列相同；agnancy_difference() 函数的行为类似。因此，

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

将序列 a，b，c，d 转换为 a，a+b，a+b+c，a+b+c+d 和

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

还原原始值。因此，partial_sum() 将 17，2，1，0，-3 还原为 17，19，20，20，17 。

对于那些认为温度差异只是气象学或科学实验室实验中枯燥细节的人来说，我想指出，分析股票价格或海平面变化也涉及完全相同的两个操作。这些操作对于分析任何一系列变化都非常有用。

40.6.4 itoa

调用 iota(b,e,n) 将 n+i 赋值给集合 [b:e) 的第 i 个元素。例如：

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");

iota 这个名称是希腊字母 ι 的拉丁语拼写，在 APL 中用于表示该函数。

不要将 iota() 与非标准但并不罕见的 itoa()(int-to-alpha；§12.2.4)混淆。

40.7 随机数(Random Numbers)

随机数对于许多应用至关重要，例如模拟、游戏、基于采样的算法、密码学和测试。例如，我们可能需要为路由器模拟选择 TCP/IP 地址，决定怪物应该攻击还是挠头，或者生成一组值来测试平方根函数。在标准库 <random> 中，定义了用于生成(伪)随机数的工具。这些随机数是根据数学公式生成的数值序列，而不是像放射性衰变或太阳辐射等物理过程产生的无法预测的“真正随机”数。如果实现中包含这样的真正随机设备，它将表示为 random_device(§40.7.1)。

有四种实体：

• 均匀随机数生成器是一个函数对象，它返回无符号整数值，使得可能结果范围内的每个值(理想情况下)都有相等的概率返回。

• 随机数引擎(引擎)是一个均匀随机数生成器，可以用默认状态 E{} 或由种子 E{s} 确定的状态来创建。

• 随机数引擎适配器(适配器)是一种随机数引擎，它接收由其他随机数引擎产生的值，并将算法应用于这些值，以便提供具有不同随机性属性的值序列。

• 随机数分布(分布)是一个函数对象，它返回的值根据相关的数学概率密度函数 p(z) 或相关的离散概率函数 P(zi) 分布。

详情请参见§iso.26.5.1。

简单来说，随机数生成器就是一个引擎加上一个分布。引擎生成一个均匀分布的值序列，而分布则将这些值调整成所需的形状(分布)。也就是说，如果你从随机数生成器中取出大量数并绘制它们，你应该能得到一个相当平滑的分布图。例如，将 normal_distribution 绑定到 default_random_engine，就能得到一个生成正态分布的随机数生成器：

auto gen = bind(normal_distribution<double>{15,4.0},default_random_engine{});

for (int i=0; i<500; ++i) cout << gen();

标准库函数 bind() 创建一个函数对象，该对象将根据其第二个参数调用其第一个参数(§33.5.1)。

使用 ASCII 图形(§5.6.3)，我得到：

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 ∗

在大多数情况下，大多数程序员只需一个在给定范围内均匀分布的整数和浮点数。例如：

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() << ' ';

}

遗憾的是，Ran_int 和 Rand_double 不是标准类，但它们很容易构建：

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{});

};

我使用发行版的标准 param_type 别名(§40.7.3)存储参数，这样我就可以使用 auto 来避免为 bind() 的结果命名。

为了增加一些变化，我使用了不同的 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;

};

随机数的一个重要用途是用于抽样(采样)算法。这类算法需要从远大于总体大小的样本中抽取一定大小的样本。以下是来自一篇著名论文 [Vitter,1985] 的算法 R(最简单的算法)：

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<Siz e>;

using Param = typename Dist::param_type;

// Fill the reservoir and advance first:

copy(first,n,result);

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 (Siz e k = n; first!=last; ++first,++k) {

Size r = dist(g en,Param{0,k});

if(r < n)

∗(result + r) = ∗first;

}

return result;

}

40.7.1 引擎

均匀随机数生成器是一个函数对象，它生成一个近似均匀分布的 result_type 值序列：

均匀随机数生成器G<T> (§iso.26.5.1.3)
G::result_type	这个序列一个元素的类型
x=g()	应用运算符：x 是序列的下一个元素
x=G::min()	x 是g() 可以返回的最小元素
x=G::max()	x 是g() 可以返回的最大元素

随机数引擎是一种均匀随机数生成器，它具有一些附加特性，使其用途更加广泛：

随机数引擎E<T> (§iso.26.5.1.4)
E e {};	默认构造函数
E e {e2};	复制构造函数
E e {s};	e 将处于由种子(seed) s 确定的状态(译注：seed指的是调用伪随机数生成器的初始值)
E e {g};	e 的状态将由对种子序列 g 调用 generate() 函数决定。
e.seed()	e 将处于默认状态
e.seed(s)	e 将处于由种子 s 确定的状态
e.seed(g)	e 的状态将由对种子序列 g 调用 generate() 函数确定
e.discard(n)	跳过序列中的接下来的 n 个元素
e==e2	e 和 e2 会产生完全相同的序列吗？
e!=e2	!(e==e2)
os<<e	输出 e的一个表示到 os 中
is>>e	从 is读取一个先前由 << 输出的一个引擎表示到 e中

种子是介于 [0 : $2^{32}$ ) 之间的值，可用于初始化特定的引擎。种子序列 g 是一个对象，它提供了一个函数 g.generate(b,e)，调用该函数会将 [b:e) 填充为新生成的种子(§iso.26.5.1.2)。

标准随机数引擎的 UI 参数必须是无符号整数类型。对于 linear_congruential_engine<UI,a,c,m>，如果模数 m 为 0，则使用 numeric_limits<result_type>::max()+1 的值。例如，以下代码会输出一个数首次重复的索引：

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';

我运气不错，参数设置还算合理，没有出现重复值。试试 <unsigned int,16,5,0>，看看有什么区别。除非你真的有需要并且清楚自己在做什么，否则请使用 default_random_engine。

随机数引擎适配器接受一个随机数引擎作为参数，并生成一个具有不同随机属性的新随机数引擎。

标准随机数引擎适配器(§iso.26.5.4)
discard_block_engine<E,p,r>	E是引擎；§iso.26.5.4.2
independent_bits_engine<E,w,UI>	在类型UI中生成 w 个位; §iso.26.5.4.3
shuffle_order_engine<E,k>	§iso.26.5.4.4

例如：

independent_bits_engine<default_random_engine,4,unsigned int> ibe;

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

cout << '0'+ibe() << ' ';

这将会生成范围 [48:63]([‘0’: ‘0’+ $2^{4}-1$ ))中的 100 个数。

为一些常用引擎定义了别名：

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>;

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>;

40.7.2 随机设备(Random Device)

如果某个实现能够提供真正的随机数生成器，则该随机数源将以名为 random_device 的均匀随机数生成器的形式呈现：

random_device (§iso.26.5.6)
random_device rd {s};	string s 标识随机数源；实现定义；explicit
d=rd.entropy()	d 是 double 型，对一个伪随机数生成器，d=0.0

将 s 视为随机数源的名称，例如Geiger计数器、Web 服务或包含真正随机源记录的文件/设备。对于一个具有其分别的概率为 $P_{0} ,..., P_{n-1}$ 的 n 个状态的设备，其 entropy() 定义为：

$S( P_{0} ,..., P_{n-1} ) = -\sum_{i=0}^{i=n-1}P_{i} \log{ P_{i}}$ 。

熵(entropy)是对生成数的随机性(即不可预测程度)的估计。与热力学相反，高熵(high entropy)对于随机数来说是好事，因为这意味着很难预测后续的数。该公式反映了反复投掷一个完美的n面骰子的结果。

random_device 的设计初衷是用于密码学应用，但未经仔细研究就信任 random_device 的实现，这违反了此类应用的所有规则。

40.7.3 分布(Distributions)

随机数分布是一个函数对象，当使用随机数生成器参数调用它时，会生成一个 result_type 类型的值序列：

随机数生成D (§iso.26.5.1.6)
D::result_type	D的一个元素之类型
D::param_type	构造 D 所需的参数集类型
D d {};	默认构造函数
D d {p};	从 param_type p 构造
d.reset()	重置为默认状态
p=d.param()	p 是 d 的 param_type 参数
d.param(p)	重置为由 param_type p 确定的状态
x=d(g)	x 是由生成器 g 生成的，由生成器 d 产生的值
x=d(g,p)	x 是由 d 在给定生成器 g 和参数 p 的情况下生成的值。
x=d.min()	x 是 d 可以返回的最小值
x=d.max()	x 是 d可以返回的最大值
d==d2	d 和 d2 会产生相同的元素序列吗？
d!=d2	!(d==d2)
os<<d	将 d 的状态输出到os，以便可以通过 >> 读取回来。
is>>d	将一个前面使用 << 从 is 输出的状态读入 d 中

在下表中，模板参数 R 表示数学公式中必须使用实数，默认值为 double。I 表示必须使用整数，默认值为 int 。

均匀分布(§iso.26.5.8.2)
分布	预置条件	默认	结果
uniform_int_distribution<I>(a,b)	a ≤b P(i\|a, b)=1/(b-a+1)	(0,max)	[a:b]
uniform_real_distribution<R>(a,b)	a ≤b P(x\|a, b)=1/(b-a)	(0.0,1.0)	[a:b)

预置条件字段指定了分布参数的要求。例如：

uniform_int_distribution<int> uid1 {1,100}; // OK

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

默认值字段指定缺省参数。例如：

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

结果字段指定结果的范围。例如：

uniform_int_distribution<> uid3 {0,5};

default_random_engine e;

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

cout << uid3(e) << ' ';

uniform_int_distribution 的取值范围是封闭的，我们可以看到六个可能的值：

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

对于 uniform_real_distribution，以及所有其他具有浮点结果的分布，其范围是半开的。

Bernoulli分布反映了一系列具有不同程度负载的硬币抛掷结果：

Poisson 分布表示在固定的时间间隔和/或空间内发生给定数量事件的概率：

正态分布将实数值映射到实数值。最简单的正态分布是著名的“钟形曲线”，它将数值对称地分布在一个峰值(均值)周围，元素与均值的距离由标准差参数控制：

为了更好地理解这些分布情况，可以查看各种参数的图形表示。这类图形很容易生成，而且在网上也更容易找到。

抽样分布根据概率密度函数 P 将整数映射到特定范围内：

40.7.4 C风格随机数(C-Style Random Numbers)

在 <cstdlib> 和 <stdlib.h> 中，标准库提供了一个简单的随机数生成基础：

#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

产生一个好的随机数生成器并非易事，而且并非所有系统都能提供可靠的 rand() 函数。特别是，随机数的低位通常不太可靠，因此 rand()%n 并不是生成 0 到 n-1 之间随机数的理想通用方法。通常，int((double(rand())/RAND_MAX)∗n) 可以得到可接受的结果。然而，对于更严肃的应用，基于均匀整数分布(§40.7.3)的生成器会提供更可靠的结果。

调用 srand(s) 函数会从给定的种子(seed) s 开始生成一个新的随机数序列，种子 s 由参数 en 提供。为了便于调试，通常需要保证给定种子生成的随机数序列具有可重复性。然而，我们往往希望每次实际运行都使用新的种子。事实上，为了使游戏结果不可预测，从程序环境中选择种子通常更实用。对于这类程序，实时时钟的某些比特位通常可以作为不错的种子。