伪随机数算法
本文主要讨论 C++ 中最常用的几种随机数生成方法。本文不会去讲 random 库的使用，而是着力讲解其背后的原理。如果你只对 random 库的使用方法感兴趣，请移步 https://en.cppreference.com/w/cpp/numeric/random
我们在写 C++ 时，使用
rand()之前总会srand()一个种子。种子相同，得到的随机数也一模一样。这是因为随机数是根据这个种子计算出来的，而不是真正的随机。通常，这样的随机数被称为“伪随机数”。
随机数分为真随机数和伪随机数。
真随机数利用某些自然因素（如熵）的随机性生成。Linux 中的 /dev/random 生成的就是真随机数。
伪随机数则利用一些生成算法来产生。通常来讲，C++ 中生成的随机数就是伪随机数。
minstd
rand() 利用的是线性同余法（LCG）生成伪随机数，即计算方法为
为了使得循环周期等于 \(m\)，需要满足
- \(\left(m,c\right)=1\);
- \(a-1\) 可以被 \(m\) 的所有素因子整除;
- 若 \(4\mid m\)，则 \(4\mid a-1\)
以上三条被称为 Hull–Dobell Theorem。
\(a\) 和 \(c\) 的值由编译器决定。例如在 gcc 编译器中，使用的是
\[x_{n+1}=\left(1103515245 x_n+12345\right)\mod 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=\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.\]在以上操作前，我们要先对前 $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{matrix}0&I_{w-1}\\a_{w-1}&\left(a_{w-2},\cdots,a_0\right)\end{matrix}\right)\]其中，\(x_k^u\) 为 \(x_k\) 的高 \(w-r\) 位，\(x_{k+1}^l\) 为 \(x_{k+1}\) 的低 \(r\) 位，\(I\) 为单位矩阵。
实际计算中，我们使用等价的式子：
\[\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.\]其中，\(x_0\) 为 \(\boldsymbol{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，取
我们可以使用 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() 本身的算法已经没有太大关系了。
Comments