C++ 中的随机数

C++ 中的随机数

伪随机数

一般用法

C++11 新增的 <random> 完善了伪随机数生成的功能，其设计理念和用法来源于 Boost.Random。与 Boost.Random 一样，<random> 严格区分了“引擎”和“分布”两个概念，

• 引擎（Engine）负责产生随机源，反复调用引擎对象能输出一系列无符号整型数，每个数字出现的可能性都是一样的。
• 分布（Distribution）将引擎输出的随机数转换成符合特定统计分布的整型数、浮点数或布尔变量。

std::default_random_engine eng{};
std::uniform_int_distribution dist{1, 6};

构造出引擎对象 eng 和分布对象 dist 后，反复执行 dist(eng) 就等得到一列 1 到 6 之间的随机整数。

一种常见的错误使用方式是绕过分布对象 dist ，直接用 1 + eng() % 6 生成一组数字。尽管这种方法能生成性质不错的随机数，而且性能更高，但在对随机性要求高的场合并不推荐这样做。查看 STL 源代码就能知道，uniform_int_distribution 为生成 1 到 6 之间的随机整数所进行的考量和计算远比 1 + eng() % 6 复杂得多。

标准库中的引擎

STL 根据复杂度和自由度由低到高提供了三种层次的随机数引擎类型，

• 首先是最简单的 default_random_engine，作为默认推荐使用的引擎类 。default_random_engine 的具体实现因标准库的版本而异，Visual Studio 提供的实现中 default_random_engine 是 mt19937 的别名。
• 其次是一系列预定义的引擎类型：
  • minstd_rand0/minstd_rand，属于 Linear congruential 引擎类
  • mt19937/mt19937_64，属于 Mersenne twister 引擎类
  • ranlux24_base/ranlux48_base，属于 Subtract with carry 引擎类
  • ranlux24/ranlux48，属于 Discard block 引擎类
  • knuth_b，属于 Shuffle order 引擎类
  • philox4x32/philox4x64，属于 Philox 引擎类（C++26 引入）
• 最后是最底层的模板类，用户可以调整模板参数实例化自定义的随机数引擎：
  • linear_congruential_engine
  • mersenne_twister_engine
  • subtract_with_carry_engine
  • discard_block_engine
  • shuffle_order_engine
  • philox_engine（C++26 引入）

引擎的初始化与状态

引擎对象的初始化有四种方式：

1. 隐式默认初始化，default_random_engine eng;
2. 显示默认初始化，default_random_engine eng{};
3. 用随机数种子显示初始化，default_random_engine eng{13607};，需要注意的是随机数种子的类型必须是 default_random_engine::result_type
4. 前三种初始化可以保证产生的随机数具有“可复现性”，如果需要打破可复现性可以借助特殊的类型 random_device 来初始化:

std::random_device rdev{};
std::default_random_engine eng{rdev()};
std::uniform_int_distribution dist{1, 6};

std::cout << dist(eng) << std::endl;

这样每次执行程序将得到不同的结果。事实上，ISO C++ 标准允许但不要求 random_device 是非确定性的，具体情况要视操作系统、硬件和编译器而定。但在 Visual Studio 的实现中，random_device 是非确定性的。

引擎的状态可以通过 seed 方法重置，

1. eng.seed(13607)，eng 重置为种子是 13607 的引擎；
2. eng.seed()，eng 重置为默认初始化状态；
3. eng.seed(rdev)，eng 重置为不可复现的状态；

标准库中的分布

目前 <random> 提供的统计分布如下：

• 均匀分布族
  • uniform_int_distribution，产生均匀分布的整数
  • uniform_real_distribution，产生均匀分布的实数
• 伯努利分布族
  • bernoulli_distribution，产生伯努利分布上的 bool 值
  • binomial_distribution， 产生二项分布的整数
  • negative_binomial_distribution，产生负二项分布的整数
  • geometric_distribution，产生几何分布的整数
• 泊松分布族
  • poisson_distribution，产生泊松分布的整数
  • exponential_distribution，产生指数分布的实数
  • gamma_distribution，产生 $\Gamma$ 分布的实数
  • weibull_distribution，产生威布尔分布的实数
  • extreme_value_distribution，产生极值分布的实数
• 正态分布族
  • normal_distribution，产生标准正态分布的实数
  • lognormal_distribution，产生对数正态分布的实数
  • chi_squared_distribution，产生 $\chi^2$ 分布的实数
  • cauchy_distribution，产生柯西分布的实数
  • fisher_f_distribution，产生费舍尔 F 分布的实数
  • student_t_distribution，产生t 分布的实数
• 采样分布族
  • discrete_distribution，产生离散分布的整数
  • piecewise_constant_distribution，产生分段常数分布的实数（密度函数是分段常数函数）
  • piecewise_linear_distribution，产生分布线性分布的实数（密度函数是分段线性函数）

分布的初始化与参数

分布的初始化有两种方式，

1. 直接用若干参数初始化分布对象，std::normal_distribution<double> d{0.0, 1.0};
2. 用特定的 param_type 对象初始化分布对象，

std::normal_distribution<double>::param_type param{0.0, 0.2};
std::normal_distribution<double> dist{param};

除了初始化分布对象外，param_type 对象还能用来修改分布的参数。

• dist.param(param) 可以重置分布的参数；
• dist(eng, param) 使得分布对象临时使用参数 param 来生成随机数。

拟随机数

目前标准库的 <random> 仅包含了伪随机数引擎，要产生拟随机数需要依赖 Boost.Random。拟随机数引擎的用法和前面的伪随机数引擎完全一样，仅需要用维度参数初始化引擎对象，再和分布对象组合起来使用。

#include <boost/random/sobol.hpp>
#include <iostream>
#include <random>

int main() {
    boost::random::sobol eng{3};

    std::uniform_real_distribution<> dist{0.0, 1.0};

    double x, y, z;

    for (int i = 0; i < 10; ++i) {
        x = dist(eng);
        y = dist(eng);
        z = dist(eng);

        std::cout << x << ", " << y << ", " << z << std::endl;
    }

    return EXIT_SUCCESS;
}

0.5, 0.5, 0.5
0.75, 0.25, 0.25
0.25, 0.75, 0.75
0.375, 0.375, 0.625
0.875, 0.875, 0.125
0.625, 0.125, 0.875
0.125, 0.625, 0.375
0.1875, 0.3125, 0.9375
0.6875, 0.8125, 0.4375
0.9375, 0.0625, 0.6875

和 scipy 中的 Sobol 不同，boost::random::sobol 略去了起始的“零点”，其余结果是一致的。

然而事情并不像表面上看起来那么简单。对于拟随机数来说，顺序至关重要。以 3 维 Sobol 序列为例，连续三次调用 eng() 得到结果作为一个整体成为一个 3 维拟随机向量，如果跳过第二次 eng() 调用，将第一、三和四次调用的结果作为一个 3 维随机向量，将得到错误的结果。查看源代码可以知道无论是 STL 还是 Boost，分布类生成特定随机数的计算过程相当复杂，无法保证不会出现跳过某次 eng() 的情况。因此，最稳妥的方式是自定义均匀分布类，避免跳过 eng() 情况的发生，再用逆分布法得到特定分布的随机数。

和伪随机数的情况类似，Boost.Random 既提供了最底层的类模板，又提供了若干预定义引擎类：

• 预定义类
  • niederreiter_base2，生成 Niederreiter 序列，最高支持维度是 4720；
  • sobol，生成 Sobol 序列，最高支持维度是 3667；
  • faure，生成 Faure 序列，最高支持维度是 1117；
• 底层类模板
  • niederreiter_base2_engine
  • sobol_engine
  • faure_engine

拟随机数的最高支持维度由上述三个类模板的第三个模板参数决定，自定义这个模板参数就能达到扩展维度的目的。

扩展到更高维度

以性能最好也是最常用的 Sobol 序列为例，boost::random::sobol 的第三个模板参数是 boost::random::default_sobol_table，也就是 boost::random::detail::qrng_tables::sobol 的别名。如果要将维度扩展到 Joe 和 Kuo 得到 21201 维，需要参照 default_sobol_table 自定义方向数信息表 struct JoeKou6。首先，非常显然的有，

struct JoeKuoD5 {
    static const unsigned int max_dimension = 21201;
    static const unsigned int num_polynomials = max_dimension - 1;

    ...
};

从网站上下载 new-joe-kuo-6.21201 文件，下一步根据文件内容构建静态函数 polynomial 和 minit，需要按照 Boost 的规定表示本原多项式的系数和初始方向数。文件中 d 是维度，s 是本原多项式阶数，a 是表示本原多项式系数的整数（本原多项式的系数只有 0 和 1 两个数字，可以很方便的用一个二进制的整数表示），m 是初始方向数序列。

Joe 和 Kuo 规定的本原多项式具有如下形式：

\[x^s + a_1 x^{s-1} + a_2 x^{s-2} + \cdots + a_{s-1}x + 1\]

二进制整数 $a_1 a_2 \cdots a_{s-1}$ 就是文件中列 a 的数字。以第 10 行为例

10      5       7       1 1 7 11 19

对应的本原多项式是 5 阶的，二进制数 $a_1 a_2 a_3 a_4 = 0111$（也就是十进制数 7），因此本原多项式为，

\[x^5 + x^3 + x^2 + x + 1\]

Boost 规定的本原多项式具有如下形式：

\[a_0 + a_1 z + a_2 z^2 + \cdots + a_{31} z^{31}\]

因此二进制数 $a_5 a_4 a_3 a_2 a_1 a_0 = 101111$，也就是十进制数 47。

按照上述规律可以根据文件内容得到函数 polynomial，内部包含一个长度 21200 的整数数组，最后一个数字是 524263。$\log_2 524263$ 取整等于 18，而且 524263 适合用 std::uint32_t 表示，显然

struct JoeKuoD5 {
    static const unsigned int max_dimension = 21201;
    static const unsigned int num_polynomials = max_dimension - 1;

    // log2(polynomial(num_polynomials - 1)), i.e., integer log2 of the last polynomial in the table
    static const unsigned int max_degree = 18;

    using value_type = std::uint32_t;

    static value_type polynomial(std::size_t n) {...}
    ...
};

minit 的构造简单多了，将文件中 m 列的每个不等长序列用 0 补全为长度为 18 的序列，例如 1 1 7 11 19 补全为

1 1 7 11 19 0 0 0 0 0 0 0 0 0 0 0 0 0

按照上述规律可以根据文件得到函数 minit，以及函数中长度 18 x 21200 的数组 sobol_minit。

JoeKuo6 替换 default_sobol_table 就可以得到支持 21201 维的 Sobol 序列引擎。

并行化

金融工程和科学计算通常要求能够并行生成随机数。对于 Sobol 序列来说，并行化并不是难事，因为 Sobol 序列存在高效的跳转机制，即类 sobol 的成员函数 discard 计算成本很低。因此，可以采用“分块策略”实施并行化。

尽管标准库的伪随机数引擎也有成员函数 discard，但两方面的原因使得这些引擎难以借助 discard 实施分块策略的并行化：

1. discard 函数本身的计算成本过高；
2. 更严重的是，分布对象在生成随机数的时候可能会跳过某次 eng() 调用，这会导致分块策略彻底失败。

好消息是 C++26 计划引入的 Philox 引擎能够解决并行化的问题。

此外，第三方库 TRNG 是实现并行化的非常理想的选择，TRNG 提供了若干可并行的随机数生成器，并且能够适配多种并行计算架构，例如 MPI、OpenMP 甚至 CUDA。

扩展阅读

• Sobol 序列并行化的实践经验
• Random Number Generation in C++11