【C++】随机数生成原理

本文主要讨论 C++ 中最常用的几种随机数生成方法。本文不会去讲 random 库的使用，而是着力讲解其背后的原理。如果你只对 random 库的使用方法感兴趣，请移步 https://en.cppreference.com/w/cpp/numeric/random

我们在写 C++ 时，使用 rand() 之前总会 srand() 一个种子。种子相同，得到的随机数也一模一样。这是因为随机数是根据这个种子计算出来的，而不是真正的随机。通常，这样的随机数被称为“伪随机数”。

随机数分为真随机数和伪随机数。

真随机数利用某些自然因素（如熵）的随机性生成。Linux 中的 /dev/random 生成的就是真随机数。

伪随机数则利用一些生成算法来产生。通常来讲，C++ 中生成的随机数就是伪随机数。

LCG

rand() 利用的是线性同余法（LCG）生成伪随机数，即计算方法为
$$x_{n+1}=\left(a{x}_n+c\right)\mathrm{m}\mathrm{o}\mathrm{d}m$$
为了使得循环周期等于 $m$，需要满足

1. $\left(m,c\right)=1$
2. $a-1$ 可以被 $m$ 的所有素因子整除
3. 若 $4\mid m$，则 $4\mid a-1$

以上三条被称为 Hull–Dobell Theorem。

$a$ 和 $c$ 的值由编译器决定。例如在 gcc 编译器中，使用的是
$$x_{n+1}=\left(1103515245x_n+12345\right)\mathrm{m}\mathrm{o}\mathrm{d}{2}^{31}$$
各个编译器的值可以参考 https://en.wikipedia.org/wiki/Linear_congruential_generator#Parameters_in_common_use

LCG 的优势在于速度快、占用内存小。但是由于其不够随机，所以并不能把它用在蒙特卡洛之类的算法上。

RANLUX

ranlux 本质上还是 LCG，不过它使用的是带进位减法（Subtract-With-Carry）。

其算法为：
$${\Delta }_n=x_{n-s}-x_{n-r}-c_{n-1}\left(s<r\right)$$

$$x_n=\{\begin{array}{llll}{\Delta }_n,& c_n=0& \text{if}& {\Delta }_n\ge 0\\ {\Delta }_n+m,& c_n=1& \text{if}& {\Delta }_n<0\end{array}$$

在以上操作前，我们要先对前 $r$ 位进行初始化。初始化的方法为
$$b_n=b_{n-13}+b_{n-31}$$
这里，每个 $b$ 是一个比特。$b$ 的前 $31$ 位由一个整数确定，这个整数就是我们的随机数种子。

为了进一步保证其随机性，算法还会丢弃掉一部分生成的随机数。丢弃的数量与奢侈等级（Luxury level）有关。在原论文中，奢侈等级的定义如下：

奢侈等级	采用的数量 r	丢弃的数量 p-r	生成的数量 p
0	24	0	24
1	24	24	48
2	24	73	97
3	24	199	223
4	24	365	389

显然，等级越高，随机数越难以被预测。综合考虑性能等因素，3 级在现实中使用较多。

MT

梅森旋转（Mersenne Twister - C++ 中的 mt1993）也是一种伪随机数生成方法，不过可以生成比 LCG 质量高得多的随机数。

MT 得名于其周期为梅森素数 $2^{nw-r}$。其利用的是 LFSR（更准确地讲，是 GFSR）。具体算法如下。
$$x_{k+n}=x_{k+m}\oplus \left(\left(x_k^u\mid x_{k+1}^l\right)A\right)$$

$$A=\left(\begin{array}{cc}0& I_{w-1}\\ a_{w-1}& \left(a_{w-2},\cdots ,a_0\right)\end{array}\right)$$

其中，$x_k^u$ 为 $x_k$ 的高 $w-r$ 位，$x_{k+1}^l$ 为 $x_{k+1}$ 的低 $r$ 位，$I$ 为单位矩阵。

实际计算中，我们使用等价的式子：
$$\mathbf{x}A=\{\begin{array}{ll}\mathbf{x}\gg 1& x_0=0\\ \left(\mathbf{x}\gg 1\right)\oplus \mathbf{a}& x_0=1\end{array}$$
其中，$x_0$ 为 $\mathbf{x}$ 的最低位。

由于 $A$ 为有理范式，故需要级联一个 tempering transform 来补偿
$$\begin{array}{l}y\equiv x\oplus \left(\left(x\gg u\right)\&d\right)\\ y\equiv y\oplus \left(\left(y\ll s\right)\&b\right)\\ y\equiv y\oplus \left(\left(y\ll t\right)\&c\right)\\ y\equiv y\oplus \left(y\gg l\right)\end{array}$$
最后的 $y$ 即为得到的随机数。

在实行以上操作之前，我们需要提前对移位寄存器的前 $n-1$ 位进行初始化，其计算方法为
$$x_i=f\times \left(x_{i-1}\oplus \left(x_{i-1}\gg \left(w-2\right)\right)\right)+2$$
同样的，上面的变量由编译器决定。例如，对于 C++11 中的 MT19937，取
$$\begin{array}{rcl}\left(w,n,m,r\right)& =& \left(32,624,397,31\right)\\ a& =& {\text{9908B0DF}}_{16}\\ \left(u,d\right)& =& \left(11,{\text{FFFFFFFF}}_{16}\right)\\ \left(s,b\right)& =& \left(7,{\text{9D2C5680}}_{16}\right)\\ \left(t,c\right)& =& \left(15,{\text{EFC60000}}_{16}\right)\\ l& =& 18\\ f& =& 1812433253\end{array}$$
我们可以使用 python 简单模拟一下：

def _int32(x): # 截取 32 位
    return int(0xFFFFFFFF & x)

class MT19937:
    def __init__(self, seed):
        self.mt = [0] * 624
        self.mt[0] = seed
        self.mti = 0
        for i in range(1, 624): # 初始化移位寄存器的前 623 位
            self.mt[i] = _int32(1812433253 * (self.mt[i - 1] ^ self.mt[i - 1] >> 30) + i)


    def extract_number(self):
        if self.mti == 0:
            self.twist()
        y = self.mt[self.mti]
        y = y ^ y >> 11
        y = y ^ y << 7 & 0x9D2C5680
        y = y ^ y << 15 & 0xEFC60000
        y = y ^ y >> 18
        self.mti = (self.mti + 1) % 624
        return _int32(y)


    def twist(self):
        for i in range(0, 624):
            # 高位和低位级联
            y = _int32((self.mt[i] & 0x80000000) + (self.mt[(i + 1) % 624] & 0x7fffffff))
            self.mt[i] = (y >> 1) ^ self.mt[(i + 397) % 624]

            if y % 2 != 0: # 如果最低为不为零
                self.mt[i] = self.mt[i] ^ 0x9908b0df

MT19937(seed).extract_number()

TRNG

C++ 中也提供了真随机数 random_device。它在 Windows 下调用 rand_s，在 Linux 下调用 /dev/urandom。

用过 Linux 的肯定经历过这么一个场景：要求你随意敲键盘，让你停你再停。这其实就是一个 TRNG 的过程。程序根据你键盘的输入产生了独特的熵值。不单单是敲键盘，鼠标位置、环境噪音、CPU 温度等都可以作为熵的产生方法。

当然，产生的熵并不能直接使用，因为它并不随机（比如键盘上有些键的敲击频率更高），所以会经过一些处理。这里面涉及到很多硬件知识，已经超出了我的能力范围。

真随机数的优点是足够随机，但它会消耗很多系统资源，在某些情况下是不可接受的。

如果只要生成一个随机数，我们也可以使用伪随机数算法取生成真随机数。比如常用的 srand(time(nullptr))，就是利用未初始化内存的随机性作为种子生成随机数。不过这样的生成方法其实和 rand() 本身的算法已经没有太大关系了。