Ch3nyang's blog

home

home

person

about

hive

project

collections_bookmark

mindclip

rss_feed

rss

伪随机数算法

calendar_month 2022-04
archive 编程
tag random tag algorithm tag cpp

本文主要讨论 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() 本身的算法已经没有太大关系了。

Comments

Share This Post