【C++】随机数生成原理_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 = { Δ n , c n = 0 if Δ n ≥ 0 Δ n + m , c n = 1 if Δ n < 0 x_n=\left\{

$$\begin{array}{llll} \Delta_n,&c_n=0&\text{if}&\Delta_n\geq0\\ \Delta_n+m,&c_n=1&\text{if}&\Delta_n<0 \end{array}$$ \right. xn​={Δn​,Δn​+m,​cn​=0cn​=1​ifif​Δn​≥0Δn​<0​

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

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

显然，等级越高，随机数越难以被预测。综合考虑性能等因素，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 = ( 0 I w − 1 a w − 1 ( a w − 2 , ⋯   , a 0 ) ) A=\left(

$$\begin{matrix}0&I_{w-1}\\a_{w-1}&\left(a_{w-2},\cdots,a_0\right)\end{matrix}$$ \right) A=(0aw−1​​Iw−1​(aw−2​,⋯,a0​)​)

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

实际计算中，我们使用等价的式子：
x A = { x ≫ 1 x 0 = 0 ( x ≫ 1 ) ⊕ a x 0 = 1 \boldsymbol{x}A=\left\{

$$\begin{array}{ll}\boldsymbol{x}\gg1&x_0=0\\\left(\boldsymbol{x}\gg1\right)\oplus \boldsymbol{a}&x_0=1\end{array}$$ \right. xA={x≫1(x≫1)⊕a​x0​=0x0​=1​
其中，$x_0$ 为 $\mathbf{x}$ 的最低位。

由于 $A$ 为有理范式，故需要级联一个 tempering transform 来补偿
y ≡ x ⊕ ( ( x ≫ u ) & d ) y ≡ y ⊕ ( ( y ≪ s ) & b ) y ≡ y ⊕ ( ( y ≪ t ) & c ) y ≡ y ⊕ ( y ≫ l )

$$\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≡x⊕((x≫u)&d)y≡y⊕((y≪s)&b)y≡y⊕((y≪t)&c)y≡y⊕(y≫l)​
最后的 $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，取
( w , n , m , r ) = ( 32 , 624 , 397 , 31 ) a = 9908B0DF 16 ( u , d ) = ( 11 , FFFFFFFF 16 ) ( s , b ) = ( 7 , 9D2C5680 16 ) ( t , c ) = ( 15 , EFC60000 16 ) l = 18 f = 1812433253

$$\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}$$ (w,n,m,r)a(u,d)(s,b)(t,c)lf​=======​(32,624,397,31)9908B0DF16​(11,FFFFFFFF16​)(7,9D2C568016​)(15,EFC6000016​)181812433253​
我们可以使用 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()
1
2
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
28
29
30
31
32
33
34

TRNG

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

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

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

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

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