C++标准库中的随机数生成

一、伪随机与真随机

数字计算机的结果可以说是固定的、必然的。都是根据现有数据的状态得出接下来的状态。除非硬件损坏，计算机不会产生真正的随机和无法预料的事。在生活中随手抛一个硬币也受到出手动作、状态，以及风速等环境影响。一只蟑螂的运动路线也可能不是随机的。但是在隔绝的环境中一个分子的运动方向可以说是随机的。在量子力学中的测不准原理和波函数坍缩，都说明随机是一个概率，这个概率是定值可估计的，但是具体在哪儿却是随机的，人类并不能掌控计算它。所以人造出来的计算机也无法实现真正的随机，仅仅是模拟出来接近随机的效果，所以叫做伪随机。即在初始条件一定的情况下，产生的随机数结果是固定的。在玩许多可以生成随机地图的游戏时，很多游戏会提供一个种子值，固定的种子值会生成一模一样的地图。所以伪随机也不一定没有优点，至少程序员不期望bug的出现是随机的。在非硬件问题（比如太阳耀斑之比特翻转）的逻辑结果上，如果出现的结果与推导得出的结果不同，就说明推导过程错误。在这一点上，计算机模拟的世界就和人所在的世界大有不同，所以人工智能的什么就再等下辈子了，现有的AI只是建立在数学统计学的模型之上的，从理论上我就觉得不会有多么智能，除非换一种原理出发，比如从大脑怎么运作的，人类如何思考的。

所以，程序上的随机数一般都是通过生成一个种子值作为初始状态，再通过递归算法根据这个值计算出下一个值。只是按这个算法产生的数字序列分布看起来是随机的。比如高中的数学课本后面就有伪随机数表，几页的数字序列虽然本身是固定的。但我们可以随便选一个数字（比如掷骰子），然后取百位往后得到另一个数字。比如百位是8就往后得到往后8个的数字，按这种方法也就取出了新的序列，虽然实际是确定的，但只要起点和规则不同，这个新的数列也可以作为一般的随机用途了。

二、常用实现

C和C++的许多代码都如下生成随机数：

//以距1970年1月1号的0点的秒数作为种子值（注意long转到了uint）
srand(time(0));//返回0至RAND_MAX，RAND_MAX被定义为2147483647
int num = rand();//分布到 [min, max)
num = min + num % (max - min);

一直以来我都是使用的HGE引擎的算法，如下：

//省略了参数检查
UINT32 Math::g_seed = 0;
int DND::Math::GetRandInt(int min, int max)
{g_seed = 214013 * g_seed + 2531011;return min + (g_seed ^ g_seed >> 15) % (max - min + 1);
}float DND::Math::GetRandFloat(float min, float max)
{g_seed = 214013 * g_seed + 2531011;return min + (g_seed >> 16)*(1.0f / 65535.0f)*(max - min);
}void DND::Math::SetSeed(UINT32 seed)
{if (!seed) g_seed = (UINT32)time(0);else g_seed = seed;
}

三、标准库的使用

在实际的使用中，我严重怀疑了上述代码的可靠性。至于C语言的rand函数我就感觉更有局限性了。于是我决定花点时间用上新一点的、更强大的标准库实现。标准库提供了21个概率分布类，我只列举几个常用的：

1.离散均匀分布（整数）

//返回[min, max]的概率相等，返回整型//定义种子发生器
random_device rd;//定义随机数序列生成器rng，传入无符号整型作为种子，rd()会返回一个随机值（或许比time(0)好）
default_random_engine rng {rd()};//指定类型和范围 [1, 100]
uniform_int_distribution<int> dist {1, 100};//返回随机数 [1, 100]
dist(rng);

2.连续均匀分布（实数）

//返回[min, max)的概率相等，返回浮点数，注意前闭后开
//忽略了 序列生成器和种子值//指定类型和范围 [0, 50)
uniform_real_distribution<float> dist {0.0f, 50.0f};//返回随机数 [0, 50)
dist(rng);

3.正态分布

也叫高斯分布，它的曲线由期望和期望差决定。期望为所有数据的平均值，也是最有概率出现的值。标准差为方差的算术平方根，代表着随差值范围出现概率的变化程度。比如期望为50，标准差为10。那么最可能出现的数字是50。小于40和大于60的概率相同，小于32%。随机到小于30和大于70的值的概率小于5%。在确定一颗果树掉落几颗苹果的时候，填入(10,2)。那么玩家最可能获得10个苹果，然后获得小于8个或大于12个的概率大概三分之一。当然也有二十分之一概率获得小于6个或者大于14个。（注意正态分布只支持实数，至于转换为整数或许四舍五入会更符合需求）

//指定期望和标准差 (10.0, 2.0)
normal_distribution<double> dist {10.0, 2.0};//返回10的概率最高
round(dist(rng));

4.对数分布

对数分布与正态分布类似，即随机变量的对数符合正态分布，即log x符合正态分布，x的概率分布就是对数分布。曲线特点是较快的到达高峰然后拖一条长尾巴。

//指定期望和标准差 (5.0, 0.5)
lognormal_distribution<> dist {5.0, 0.5};dist(rng);

5.抽样离散分布

假设我们解决这样一个问题：一个农场出现3种不同动物的概率百分比分别是30、30、40，然后一共有N只动物，需要随机生成每一种。一般想到的代码实现，只需要用离散均匀分布生成0到99的数字，然后判断范围即可知道生成哪一种，然后迭代N次即可。不过当新加一种动物的时候，或者改变某一种动物的出现概率时，就需要改动很多的代码，并且一大串的if-else看上去也很不美观。用抽样分布即可，所有项的概率之和不需要为1，也就是说用权值代替每一项的出现概率。

//指定权值
discrete_distribution<size_t> dist {10, 20, 15, 5, 99};//注意这个返回的范围是索引 [0, 4]
//所以返回4的概率最高，返回3的概率最低
dist(rng);

一般来说权值的值和个数应该是动态变化的，如下改变：

//使用新的权值（个数不限）
dist.param({2,2,2,3,4,7,9});

6.抽样分段常数分布

这个分布需要指定两个列表，一个说明了数轴如何分段，第二个指定了每一段的权值。由于n个值可以分成n-1段，所以第二个列表长度应该是第一个列表长度减一。其中每一段内的数字出现概率相同。

//定义分段，分为了 [0, 60)，[60, 90), [90, 100)
//连续分布均为前闭后开区间
vector<float> b {0, 60, 90, 100};//指定3段的权值
vector<float> w{7, 5, 9};//应用数据
piecewise_constant_distribution<> d {begin(b), end(b), begin(w)};

7.分段线性分布

这个和前一个类似，只不过权值列表多一个，和分段的边界一一对应，因为它表示的不是区间的权值了，而是区间边界的权值。意思就是0分出现的权值是7，而60分出现的权值是5，那么30分的权值可计算得出是6。而前一个分段常数分布的权值指的是[0,60)的权值为7。

//定义分段，分为了 [0, 60)，[60, 90), [90, 100)
//连续分布均为前闭后开区间
vector<float> b {0, 60, 90, 100};//指定4个边界的权值，区间的权值由边界权值线性计算
//例如可以假想 30分的权值为 6
vector<float> w{7, 5, 9, 2};//应用数据
piecewise_linear_distribution<> d {begin(b), end(b), begin(w)};

8.其它分布

其他分布的类名如下，每一种都有各自适用的模型，使用形式都差不多，想要了解更多的可以再查资料^_^。

泊松分布：poisson_distribution

几何分布：geometric_distribution

指数分布：exponetial_distribution

伽马分布：gamma_distribution

威尔分布：weibull_distribution

二项式分布：binomial_distribution

负二项式分布：negative_binomial_distribution

极值分布：extreme_value_distribution

PS：概率密度函数简称PDF（Probability Density Function）

9.其他随机数序列的生成器

看了这么多分布后，我们应该知道了生成一个随机数大致分为了三步。第一步是设定种子值起点，第二步生成随机数序列，第三步应用分布。而STL提供的default_random_engine生成器是默认的随机数序列生成器，它生成随机的无符号整型序列。虽然是无符号整型，但也可以用到所有分布上。至于分布类是如何完成和保证这些结果正确的，确实有些像变魔术。

所以额外的标准库还提供了不同的随机数序列生成器。如下：

线性同余引擎：minstd_rand（32位无符号整型）、minstd_rand0、knuth_b

马特赛特旋转算法引擎：mt19937（32位无符号整型）、mt19937_64（64位无符号整型）

带进位减法引擎：ranlux24_base（24位整数）、anlux48_base（48位整数）、ranlux24、ranlux48

其中minstd_rand是minstd_rand0的改进版本，ranlux24和ranlux48是base版本产生的序列舍弃一大部分得来的。不过我们用default_random_engine默认生成器就够了。

四、原样封装一下

由于已经有很多地方使用了现有的随机数接口，所以暂时不可能大幅度的修改接口。不过只替换一下实现倒是没问题，等以后重写的时候再写得好看点（再说标准库实现得这么好，也没必要封装了，我也不考虑封装进dll了，开源就好了）。

首先我们需要定义变量，用于保存种子值、随机数序列生成器、分布对象（放到函数定义）：

//引入头文件和命名空间
#include <random>
using namespace std;class Math
{
private://种子值static unsigned int g_seed;//随机数序列生成器static default_random_engine g_random;
};//静态成员需要类外定义
unsigned int Math::g_seed = 0;
default_random_engine Math::g_random;

对于种子值，为无符号整型。传入0时让程序自动生成随机种子，有时候我们需要知道随机产生的种子值是多少，所以有必要缓存一下：

//设置种子 0代表随机
static void SetSeed(unsigned int s)
{if (s == 0){random_device rd;g_random.seed(g_seed = rd());}elseg_random.seed(g_seed = s);
}
//返回种子值
static unsigned int GetSeed() { return g_seed; }

分布的一般化写法可以写成模板，将分布对象放入函数作为静态变量：

//返回[min,max]区间的 整型随机值
template <typename T>
static T GetRandInteger(T min, T max)
{static uniform_int_distribution<T> dist_int;//设定区间dist_int.param(uniform_int_distribution<T>::param_type{ min, max });return dist_int(g_random);
}//返回[min,max)区间的 实数随机值
template <typename T>
static T GetRandReal(T min, T max)
{static uniform_real_distribution<T> dist_real;//设定区间dist_real.param(uniform_real_distribution<T>::param_type{ min, max });return dist_real(g_random);
}

这下旧的接口特例化模板即可：

//返回[min,max]区间的随机int
static int GetRandInt(int min, int max)
{return GetRandInteger<int>(min, max);
}//返回[min,max)区间的随机float
static float GetRandFloat(float min, float max)
{return GetRandReal<float>(min, max);
}

满足了旧的需求，我再简单用一下正态分布吧，也很简洁：

//返回 期望mu，标准差sigma 的正态分布随机值
template <typename T>
static T GetRandNormal(T mu, T sigma)
{static normal_distribution<T> dist_normal;//设定期望与标准差dist_normal.param(normal_distribution<T>::param_type{ mu, sigma });return dist_normal(g_random);
}

五、更新实现后的地图效果

从HGE的算法改为STL实现后，并没有出现bug。但是生成的岛屿风格还是有那么一点点的不同，其中必定有一些原因。不过算是如愿以偿的解决了随机数的需求，以后再用随机数的时候就不用担心暗藏bug了。

六、蒙特卡洛方法估计圆周率

向一块正方形区域投掷飞镖，假设飞镖随落在正方形内任意位置的概率相等，那么落在内接圆的概率等于圆面积除以正方形面积。通过统计圆内飞镖数就可以估算出pi值，这种方法就叫统计模拟方法。

即 $\frac{SumCircle}{SumQuad}=\frac{\pi r^{2}}{4r^{2}}=\frac{\pi}{4}$ ，落在圆内的数量比上总数的比值应该为四分之pi。

通过编写代码，测出了以下数据：

默认序列 + 连续均匀分布
总量	圆内	pi值
10^1	5	2
10^2	69	2.76
10^3	784	3.136
10^4	7900	3.16
10^5	78593	3.14372
10^6	785144	3.14058
10^7	7855222	3.14209
10^8	78535828	3.14143
10^9	一分钟内未得出结果	-

代码如下，得出的数据算是符合了。

#include <iostream>
#include <random>using namespace std;template <typename T, typename T2>
void one_test(unsigned all_num, T& dist, T2& rng)
{double x, y;unsigned circle_num = 0;for (unsigned i = 0; i != all_num; ++i){x = dist(rng);y = dist(rng);if (x*x + y*y < 1.0){++circle_num;}}cout << "随机了 " << all_num << "个点，在圆内的有: " << circle_num << "个" << endl;cout << "估计PI值为: " << 4.0 * circle_num / all_num << endl;cout << endl;
}int main()
{//随机种子random_device rd;//序列生成器default_random_engine rng{ rd() };//分布到 [-1.0, 1.0)uniform_real_distribution<double> dist{ -1.0, 1.0 };for (unsigned i = 1; i != 10; ++i){one_test(pow(10, i), dist, rng);}return 0;
}

不过标准库还有一种分布叫标准均匀分布（generate_canonical）专门产生[0,1)的连续分布，区别是可以指定尾数比特个数。不过这是一个模板函数，所以如下使用：

template <typename T>
void one_test2(unsigned all_num, T& rng)
{double x, y;unsigned circle_num = 0;for (unsigned i = 0; i != all_num; ++i){x = generate_canonical<double, 32>(rng) * 2.0 - 1.0;y = generate_canonical<double, 32>(rng) * 2.0 - 1.0;if (x*x + y*y < 1.0){++circle_num;}}cout << "随机了 " << all_num << "个点，在圆内的有: " << circle_num << "个" << endl;cout << "估计PI值为: " << 4.0 * circle_num / all_num << endl;cout << endl;
}

默认序列 + 标准均匀分布(32位)
总量	圆内	pi值
10^1	7	2.8
10^2	69	2.76
10^3	789	3.156
10^4	7898	3.1592
10^5	78478	3.13912
10^6	785271	3.14108
10^7	7853267	3.14131
10^8	78539235	3.14157
10^9	一分钟内未得出结果	-

这个标准均匀分布的结果在总量相同的时候，pi的精度也大概高一位，所以还是比前者好。不过梅森旋转算法被称为最好的随机数算法，所以我们也该试一下：

//梅森旋转随机数序列生成器
mt19937_64 rng_2{ rd() };for (unsigned i = 1; i != 10; ++i)
{one_test2(pow(10, i), rng_2);
}

这个的结果我就不贴图了，效果也一般。不过还说有个特别适合蒙特卡洛模拟的生成器（带进位减法），如下：

//带进位减法随机数序列生成器
ranlux48 rng_3{ rd() };//连续均匀分布
for (unsigned i = 1; i != 10; ++i)
{one_test(pow(10, i), dist, rng_3);
}//标准均匀分布
for (unsigned i = 1; i != 10; ++i)
{one_test2(pow(10, i), rng_3);
}

经测试，无论哪一种分布，在10^6数量生成的时候就超过了十几秒，并且精度也没有显著提高。所以使用默认序列生成器应该能满足一般的需求。

2019年11月21日，略游