2022-05-01-伪随机数算法.md

layout	title	date	categories	tags	comments	mathjax	copyrights
post	伪随机数算法	2022-05-01 00:00:00 +0800	编程	random algorithm cpp	true	true	原创

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

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

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

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

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

minstd

rand() 利用的是线性同余法（LCG）生成伪随机数，即计算方法为

$$ x_{n+1}=\left(a x_n+c\right)\mod 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(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，取

$$ \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() 本身的算法已经没有太大关系了。