【硬技术】随机数和随机分布

本文出自商汤 SenseParrots 内部技术分享中“实用算法及工程经验”专题下的一次分享。

首先我们从  （伪）随机性 谈起，例举生活中对随机性的应用，以此引入随机数生成器。

然后我们将介绍两种简单的随机数生成器算法：线性同余法和线性反馈移位寄存器，并给出简单的 Python 代码示例。

拥有随机数生成器后，我们就可以用它生成各类经典的概率分布，这里不仅介绍了逆变法和拒绝采样法，而且还给出了八种经典的概率分布。

同时，我们聚焦在生成分布的方法上，列举出 Python 标准库 random 模块和 PyTorch 中对应分布的源码。

最后，为了大批量生成随机数，我们将介绍 4 种常见的并行策略，并基于线性同余法简单给出其中两种并行策略（块分法和跳跃法）的 Python 代码示例。

1 随机性和伪随机性

本节讲述（伪）随机性及其应用，以及随机数生成器。

1.1 基本概念

• 随机性 (Randomness)：可预测性缺失的特性。无法准确预测。
• 伪随机性 (Pseudorandomness)：看似可预测性缺失（看似具有随机性），实际没有。可以准确预测。

1.2 随机性应用

1. 仿真：  使用计算机模拟自然现象，引入随机数可以变得逼真。比如：

• 核物理中模拟粒子的随机碰撞；
• 运筹学中模拟在一段时间内机场的人流量；
• 3D 游戏、动画等各种逼真模拟。

2. 抽样：  随机抽样，用样本估计总体。

3. 数值分析与程序设计：蒙特卡洛模拟、随机数生成等。

4. 决策：  使用抛硬币或掷飞镖来做出选择。

5. 美学：  在图形设计和音乐创作中引入一些随机性，让作品更加生动。

6. 娱乐：  摇色子、转轮盘等娱乐活动。

1.3 随机数生成器

随机数生成器 (Random number generator)  是通过一些算法、物理讯号、环境噪音等来产生看起来似乎没有关联性的数列的方法或装置。丢硬币、掷色子、洗牌就是生活中常见的随机数产生方式。

这里我们主要关注基于计算机的随机数生成器，主要包括真随机数生成器（True Random Number Generator）和伪随机数生成器（Pseudo Random Number Generator）。

大部分计算机上的伪随机数，并不是真正的随机数，只是重复的周期比较大的数列，是按一定的算法和种子值生成的。

1.3.1 真随机数生成器

随机数生成器（True Random Number Generator，缩写：TRNG）有两种类型，一种是基于特定的硬件来生成真随机数，一种是捕获用户活动来生成真随机数。

• 基于硬件的真随机数生成：  硬件随机数生成器（hardware random number generator），比如通过 CPU 内置的放大电路的热噪声来产生随机数；
• 基于用户活动的真随机数生成：  比如 Linux 下的 /dev/random，利用系统的“熵池”收集随机性（比如：当前进程数、IO 错误次数、用户键入等等“客观”随机过程选择随机 bit）来产生随机数。

真随机数生成器通常每秒只能产生很有限的随机比特，这意味着它是相对较慢的。为了提高数据产生效率，它们都常被用来生成伪随机数生成器的“种子”。

1.3.2 伪随机数生成器

伪随机数生成器（pseudo random number generator，PRNG）是一个生成数字序列的算法。PRNG 生成的序列并不是真随机，它完全由一个初始值决定，这个初始值被称为 PRNG 的随机种子（seed，这个种子可能包含真随机数）。

注意，一般伪随机数生成器不能用于安全性或加密用途。像 Python 标准库中的 random 模块就有警告：

“The pseudo-random generators of this module should not be used for security purposes. For security or cryptographic uses, see the secrets module.”

1.3.3 随机数生成器对比

P.S.：另外还有密码学伪随机数生成器（CSPRNG，Cryptographically Secure Pseudo-Random Number Generator），其本质还是伪随机数生成器，只是因为借助了一些技术使其更难预测。对于安全性较高（比如生成银行账户的验证码）的需求可以采用该类生成器。

2 随机数生成算法

本节内容仅涉及伪随机数的生成算法。这里介绍两个简单的伪随机数生成器。

2.1 线性同余法

1. 基本形式：

线性同余法 (Linear Congruential Generator, LCG) 是 D.H.Lehmer 于 1949 年提出的，其具体形式如下：

通过上述公式不断迭代所获得的随机数序列
$<X_n>$ 就称为线性同余序列。

2. 周期性：

在上述式子中，如线性同余法的参数选择不好，是很容易产生明显周期性的。例如：取
$m=10$，
$X_0= a = c = 7$，那么我们计算可获得如下序列：

这显然不是我们所期望的随机性。上述例子中的周期长度为 4，我们希望周期能够尽可能的大，最好能大于计算机能表示的最大整数。所以线性同余的参数选择需要格外小心。

一般形式为
$X_{n+1}=f(X_n)$ 的伪随机数生成器都存在周期性，合适的参数配置才可以获得较大的周期。

3. 最大周期：

为了达到最大的周期，线性同余公式中的四个参数需要满足以下条件（具体原因参见 Donald E.Knuth 《计算机程序设计艺术（第二卷）半数值算法（第3版）》3.21. 线性同余法 p.10-p.22）：

• 
  $c$ 与
  $m$ 互质；
• 
  $m$ 的所有质因数都能整除
  $a-1$；
• 若
  $m$ 是 4 的倍数，那么
  $a-1$ 也需要是；
• 
  $a$,
  $c$,
  $X_0$ 都要比
  $m$ 小；
• 
  $a$,
  $c$ 是正整数。

下表摘自维基百科，可以看出 LCG 受欢迎的程度，以及推荐的参数取值。

2.2 线性反馈移位寄存器

1. 基本概念：

线性反馈移位寄存器（Linear Feedback Shift Register，LFSR）是一个移位寄存器，其输入位是其先前状态的线性函数。具体拆析如下：

• 移位寄存器（Shift Register，SR）：

移位寄存器的每一位都存着一个二进制数，每次运算下移位寄存器的整体会向右移动一位并输出数据，空缺的位默认填 0。

比如下面 5 个 bits 的移位寄存器：

其依次移动 5 位后输出的数据分是：1、0、1、0、0。

• 反馈移位寄存器（Feedback Shift Register，FSR）：

在移位寄存器每次向右移位一位之后，最左边就会空出一位，默认情况该位置会填充 0，但如果这时引入一个反馈函数，以寄存器中已有的某些位作为反馈函数的输入，在函数中经过一定的运算后，就会将输出的结果填充到移位寄存器的最左端。

像这样拥有反馈函数的移位寄存器就被称为反馈移位寄存器（Feedback Shift Register, FSR）。在如下图所示的五位的反馈移位寄存器中就引入了反馈函数，该函数的输入是：

• 线性反馈移位寄存器（Linear Feedback Shift Register，LFSR）：

如果反馈移位寄存器中的反馈函数是线性函数，那么这种寄存器就被称为线性反馈移位寄存器（Linear Feedback Shift Register，LFSR）。

如下图所示，LFSR 中每一位的数据参与异或运算就记为：
$c_i = 1$，不参与异或运算就记为：
$c_i = 0$。我们把参与异或的位称为 抽头（Taps） 。LFSR 的反馈函数就是简单地对移位寄存器中的某些位进行异或运算，并将异或的结果填充到 LFSR 的最左端。注意：  这里的位序和二进制的高低位序相反。

如上图中所示，移位寄存器中的值表示为 b1,b2,…,bn ，则最左端第 1 位的新值 b1new 可以表示为：

其中

• 
  $b_i, i \in [1, n]$ 表示移位寄存器中的数据（0 或 1）；
• 
  $c_i, i \in [1, n]$ 表示第
  $i$ 位是否是 抽头，如果是，则
  $c_i = 1$，表示该位将参与运算；如果不是，则
  $c_i=0$，表示该位将不参与运算；
• 
  $\oplus$ 表示异或运算。

上述递推关系公式的对应特征多项式是：

• 这里的
  $c_i$ 与上个公式的
  $c_i$ 相同；
• 该公式十分像二进制转十进制的公式，但位序和二进制的高低位序相反；
• +1 是因为去除全 0 的情况，这意味着隐藏了
  $b_0 = 1$, 同时它不对应任何抽头。

一个
$n$ 位的 LFSR 最多能存储
$2^n-1$ 种状态。LFSR 中全为 0 的情况会导致反馈函数的返回值永远是 0，输出序列将一直是 0。这样的反馈移位寄存器没有太大意义，因此要减 1 去除这种情况。

例如，一个三位的 LFSR 如下图所示，它最多可以包含 001、010、011、100、101、110、111 等 7 种状态。

其反馈函数为：

f(b1,b2,b3)=b1⊕b3

对应特征多项式为：

f(x)=x3+x+1

因为
$b_2$ 没有对应的抽头，即
$c_2 = 0$，所以相应的反馈函数和特征多项式都没有
$c_2$ 这一项。

2. 最大周期：

最大长度 LFSR 的形式是由 Solomon W. Golomb 在其 1967 年的书《Shift register sequences》中开发的。特征多项式的不同，即抽头的不同，亦是该随机数生成器的参数不同。

参数的不同决定了随机数生成器的周期，这里给出了每种长度下的 LFSR 达到最大周期的特征多项式，移位寄存器的长度不超过 24。

2.3 代码示例

我们基于 Python 对以上两种方法进行简单实现示例：

import math
import matplotlib.pyplot as plt

def LCG(x0, a, c, m, num):
    x, y = [0], [x0/m]
    for i in range(1,num):
        x0 = (a*x0+c)%m
        y.append(x0/m)
        x.append(i)
    return x, y

def LFSR(num):
    x, y = [], []
    start_state = 1 << 15 | 1
    lfsr = start_state
    period = 0

    for i in range(num):
        #taps: 16 15 13 4; feedback polynomial: x^16 + x^15 + x^13 + x^4 + 1
        bit = (lfsr ^ (lfsr >> 1) ^ (lfsr >> 3) ^ (lfsr >> 12)) & 1
        lfsr = (lfsr >> 1) | (bit << 15)
        period += 1
        x.append(i)
        y.append(lfsr/(2**16-1))
    return x, y

w = 20
x, y = LCG(12, 34, 12, 345346, 2000)
plt.figure(figsize=(w, 6))
plt.plot(x,y)
plt.ylabel("random value"), plt.xlabel("number of iterations")
plt.show()

plt.hist(y)
plt.ylabel("frequency(number)"), plt.xlabel("frequency distribution")
plt.show()

x, y = LFSR(2000)
plt.figure(figsize=(w, 6))
plt.plot(x,y)
plt.ylabel("random value"), plt.xlabel("number of iterations")
plt.show()

plt.hist(y)
plt.ylabel("frequency(number)"), plt.xlabel("frequency distribution")
plt.show()

输出如下：

3 随机分布与生成策略

本节介绍代表性随机分布特性及其生成策略简介。

3.1 生成策略

3.1.1 逆变法（连续随机变量）

模拟连续型随机变量的一种一般方法称为逆变换方法，它是基于下列命题实现的。

命题：

设
$U$ 为
$(0, 1)$ 上均匀随机变量，
$F$ 为任意一个连续分布函数，如果定义随机变量

则
$Y$ 具有分布函数
$F$。[
$F^{-1}(x)$ 是方程
$F(y)=x$ 的解。]

证明：

因为
$F(x)$ 是一个单调函数，所以
$F^{-1}(U)\leq a$ 成立的充要条件是
$U\leq F(a)$，因此，由上式可得

由上面的命题可知，我们可以通过产生一个随机数
$U$ 并令
$X=F^{-1}(U)$ 来模拟具有连续分布函数
$F$ 的随机变量
$X$。

3.1.2 逆变法（离散随机变量）

例如，如果我们想要模拟一个随机变量
$X$，它具有如下概率分布列：

可利用下面的方法，它是逆变换方法的离散版本。

设
$U$ 为随机数，令

由于

我们可以看出，所产生的随机变量
$X$ 具有离散分布列
$\{P_j, j=1, 2, \dots\}$。

3.1.3 拒绝采样法

拒绝采样法是一种相对简单的方法。

如下图所示，它经过两次采样，分别记作 值和 值，当点 处于即蓝线以下时则保留，否则舍弃。输出时丢弃 值，只输出 值，这样就能保证输出是符合曲线所代表的分布的。但这种方法会浪费大量的采样点（红色点都被丢弃了）。

优化后的拒绝采样 ——Ziggurat 法：

Ziggurat 法是优化版的拒绝采样，其主要目的是减少无效运算。

如下图所示，首先沿着分布曲线生成面积相等的矩形，将采样区域限制在这些面积相等的矩形框内。采样点如果落在绿色区域则保留，落在红色区域则根据拒绝采样法进行二次判断，最后只保留蓝线以下部分，丢弃蓝线以上部分。

3.2 代表性分布

下面我们列举几个具有代表性的分布，以及在 Python 标准库和 PyTorch 中的实现方法，这些方法基本都是基于逆变法来实现的。

逆变法的关键在于求出累积分布函数的反函数，而下文中的分位函数（Quantile function）就是累积分布函数的反函数。

对于没有在图中给出分位函数的分布，我们进行了一些简单的推导示例。当然不是所有生成分布的方法都是逆变法，感兴趣的读者可以进一步查看源码，这里不再例举。

3.2.1 均匀分布

定义：

均匀分布的概率密度函数（Probability Density Function, PDF）、累积分布函数（Cumulative Distribution Function, CDF）及其图像如下图所示：

分位函数：

根据上图中的累积分布函数，我们可以轻松推导出它的反函数为：

代码：

Python 标准库 random 模块：

def uniform(self, a, b):
    "Get a random number in the range [a, b) or [a, b] depending on rounding."
    return a + (b-a) * self.random()

PyTorch 源码：

这里展示了用到逆变法的两个成员方法：

class Uniform(Distribution):
    def rsample(self, sample_shape=torch.Size()):
        shape = self._extended_shape(sample_shape)
        rand = torch.rand(shape, dtype=self.low.dtype, device=self.low.device)
        return self.low + rand * (self.high - self.low)

    def icdf(self, value):
        result = value * (self.high - self.low) + self.low
        return result

3.2.2 指数分布

定义：

指数分布的概率密度函数，累积分布函数及其图像如下图所示，其中，分位函数（Quantile function）是累积分布函数的反函数。

代码：

Python 标准库 random 模块：

def expovariate(self, lambd):
    """Exponential distribution. lambd is 1.0 divided by the desired mean. It should be nonzero. (The parameter would be called "lambda", but that is a reserved word in Python.) Returned values range from 0 to positive infinity if lambd is positive, and from negative infinity to 0 if lambd is negative. """
    # lambd: rate lambd = 1/mean
    # ('lambda' is a Python reserved word)

    # we use 1-random() instead of random() to preclude the
    # possibility of taking the log of zero.
    return -_log(1.0 - self.random())/lambd

PyTorch 源码：

这里展示了用到逆变法的两个成员方法：

class Exponential(ExponentialFamily):
    def rsample(self, sample_shape=torch.Size()):
        shape = self._extended_shape(sample_shape)
        if torch._C._get_tracing_state():
            # [JIT WORKAROUND] lack of support for ._exponential()
            u = torch.rand(shape, dtype=self.rate.dtype, device=self.rate.device)
            return -(-u).log1p() / self.rate
        return self.rate.new(shape).exponential_() / self.rate

    def icdf(self, value):
        return -torch.log(1 - value) / self.rate

3.2.3 三角分布

定义：

三角分布的概率密度函数，累积分布函数及其图像如下图所示。

其分位函数推导很简单，读者可自行尝试。

代码：

Python 标准库 random 模块：

def triangular(self, low=0.0, high=1.0, mode=None):
    """Triangular distribution. Continuous distribution bounded by given lower and upper limits, and having a given mode value in-between. http://en.wikipedia.org/wiki/Triangular_distribution """
    u = self.random()
    try:
        c = 0.5 if mode is None else (mode - low) / (high - low)
    except ZeroDivisionError:
        return low
    if u > c:
        u = 1.0 - u
        c = 1.0 - c
        low, high = high, low
    return low + (high - low) * (u * c) ** 0.5

未在 PyTorch 中找到三角分布的代码，这里没有对应展示。

3.2.4 正态分布

定义：

正态分布的概率密度函数，累积分布函数及其图像如下图所示，其中，分位函数（Quantile function）是累积分布函数的反函数。

生成方法：

（a）直接采用分位函数

上图中，分位函数中涉及到一个误差函数（error function），也叫高斯误差函数（Gauss error function）。它不是一个初等函数，定义如下：

$\operatorname{erf}(z) = \frac{2}{\sqrt\pi}\int_0^z e^{-t^2}\,dt$

采用该方法需要一些数值技巧，计算相对比较麻烦。

（b）Box-Muller 方法

Box-Muller 方法 (Box & Muller 1958b; Pike 1965) 由一对均匀分布随机数得到一对高斯随机数，是一种最早的精确转换方法，并且在很长时间内该方法都是生成正态分布随机数的“标准”算法。

该算法的特点是效率高，计算过程比较简单。该方法引入了概率论的一些技巧，但背后还是用了逆变法，其中用到的指数分布和均匀分布都是逆变法。

二维标准正态分布的概率密度函数为

在极坐标系下又可写成

由于
$p(r)$ 是关于极坐标对称的，我们只需要考虑半径
$r$，而角度
$\theta$ 可以用
$[0,2\pi]$ 区间上的均匀分布（
$\theta=2\pi u_2$）。由于
$p(r)$是对于变量
$r^2$ 的指数分布，其中参数
$a=1/2$，
$r$ 可以通过上述指数分布生成

其中
$u_1$ 和
$u_2$ 是满足均匀分布的随机数，进而生成两个满足标准正态分布且完全独立的随机数

代码：

Python 标准库 random 模块：

def gauss(self, mu, sigma):
    """Gaussian distribution. mu is the mean, and sigma is the standard deviation. This is slightly faster than the normalvariate() function. Not thread-safe without a lock around calls. """

    # When x and y are two variables from [0, 1), uniformly
    # distributed, then
    #
    # cos(2*pi*x)*sqrt(-2*log(1-y))
    # sin(2*pi*x)*sqrt(-2*log(1-y))
    #
    # are two *independent* variables with normal distribution
    # (mu = 0, sigma = 1).
    # (Lambert Meertens)
    # (corrected version; bug discovered by Mike Miller, fixed by LM)

    # Multithreading note: When two threads call this function
    # simultaneously, it is possible that they will receive the
    # same return value. The window is very small though. To
    # avoid this, you have to use a lock around all calls. (I
    # didn't want to slow this down in the serial case by using a
    # lock here.)

    random = self.random
    z = self.gauss_next
    self.gauss_next = None
    if z is None:
        x2pi = random() * TWOPI
        g2rad = _sqrt(-2.0 * _log(1.0 - random()))
        z = _cos(x2pi) * g2rad
        self.gauss_next = _sin(x2pi) * g2rad

    return mu + z*sigma

PyTorch 源码：

下面代码中，成员函数 icdf 直接采用了分位函数公式。

class Normal(ExponentialFamily):
    def rsample(self, sample_shape=torch.Size()):
        shape = self._extended_shape(sample_shape)
        eps = _standard_normal(shape, dtype=self.loc.dtype, device=self.loc.device)
        return self.loc + eps * self.scale

    def icdf(self, value):
        return self.loc + self.scale * torch.erfinv(2 * value - 1) * math.sqrt(2)

3.2.5 对数正态分布

定义：

对数正态分布的概率密度函数，累积分布函数及其图像如下图所示，其中，分位函数（Quantile function）是累积分布函数的反函数。

代码：

下面两段代码示例，都是在正态分布的基础上，套了一个指数函数来实现。

Python 标准库 random 模块：

def lognormvariate(self, mu, sigma):
    """Log normal distribution. If you take the natural logarithm of this distribution, you'll get a normal distribution with mean mu and standard deviation sigma. mu can have any value, and sigma must be greater than zero. """
    return _exp(self.normalvariate(mu, sigma))

PyTorch 源码：

class LogNormal(TransformedDistribution):
    def __init__(self, loc, scale, validate_args=None):
        base_dist = Normal(loc, scale, validate_args=validate_args)
        super(LogNormal, self).__init__(base_dist, ExpTransform(), validate_args=validate_args)

上述代码就是正态分布的结果套了一个指数函数 ExpTransform：

3.2.6 帕累托分布

定义：

帕累托分布的概率密度函数，累积分布函数及其图像如下图所示，其中，分位函数（Quantile function）是累积分布函数的反函数。

代码：

下面两段代码示例，都是直接采用帕累托分布的分位函数（Quantile）来实现。

Python 标准库 random 模块：

def paretovariate(self, alpha):
    """Pareto distribution. alpha is the shape parameter."""
    # Jain, pg. 495

    u = 1.0 - self.random()
    return 1.0 / u ** (1.0/alpha)

PyTorch 源码：

class Pareto(TransformedDistribution):
    def __init__(self, scale, alpha, validate_args=None):
        self.scale, self.alpha = broadcast_all(scale, alpha)
        base_dist = Exponential(self.alpha, validate_args=validate_args)
        transforms = [ExpTransform(), AffineTransform(loc=0, scale=self.scale)]
        super(Pareto, self).__init__(base_dist, transforms, validate_args=validate_args)

上面 PyTorch 的实现也是用的分位函数，只不过拆成了指数分布 + ExpTransform 和 AffineTransform 两个变换组合起来实现：

3.2.7 韦伯分布

定义：

韦伯分布的概率密度函数，累积分布函数及其图像如下图所示。

分位函数：

根据上图中的累积分布函数，我们可以轻松推导出它的反函数为：

代码：

Python 标准库 random 模块：

def weibullvariate(self, alpha, beta):
    """Weibull distribution. alpha is the scale parameter and beta is the shape parameter. """
    # Jain, pg. 499; bug fix courtesy Bill Arms

    u = 1.0 - self.random()
    return alpha * (-_log(u)) ** (1.0/beta)

PyTorch 源码：

class Weibull(TransformedDistribution):
    def __init__(self, scale, concentration, validate_args=None):
        self.scale, self.concentration = broadcast_all(scale, concentration)
        self.concentration_reciprocal = self.concentration.reciprocal()
        base_dist = Exponential(torch.ones_like(self.scale), validate_args=validate_args)
        transforms = [PowerTransform(exponent=self.concentration_reciprocal),
                      AffineTransform(loc=0, scale=self.scale)]
        super(Weibull, self).__init__(base_dist,
                                      transforms,
                                      validate_args=validate_args)

上面 PyTorch 的实现也是用的分位函数，将其拆成了指数分布和一系列的组合变换来实现：

3.2.8 几何分布

定义：

几何分布是一种离散分布，其概率质量函数（Probability Mass Function），累积分布函数及其图像如下图所示。

分位函数：

假设有一独立重复试验，每次成功的概率为
$p$，
$0<p<1$，试验进行到出现成功为止，记
$X$ 为试验的次数，则

$X=i$ 表示前
$i-1$ 次试验均失败，而第
$i$ 次成功。称随机变量
$X$ 是参数为
$p$ 的几何随机变量。因为

所以
$X$ 可以由下列方式产生：产生一个随机数
$U$，当

时，
$X$ 取为
$j$，上式与

是等价的，又由于
$U$ 与
$1-U$ 具有相同的分布，因此，我们可以定义
$X$ 为

代码：

Python 标准库 random 模块中没有涉及到几何分布，所以这里没有展示。

PyTorch 源码：

class Geometric(Distribution):
    def sample(self, sample_shape=torch.Size()):
        shape = self._extended_shape(sample_shape)
        tiny = torch.finfo(self.probs.dtype).tiny
        with torch.no_grad():
            if torch._C._get_tracing_state():
                # [JIT WORKAROUND] lack of support for .uniform_()
                u = torch.rand(shape, dtype=self.probs.dtype, device=self.probs.device)
                u = u.clamp(min=tiny)
            else:
                u = self.probs.new(shape).uniform_(tiny, 1)
            return (u.log() / (-self.probs).log1p()).floor()

4 随机数并行生成策略

在并行应用中，假设我们需要生成随机数流
$t_{j,i}$，该随机流由
$j=0,1,\cdots,p-1$ 编号，其中
$p$ 是总的进程数。并且每个随机流内以及流间的
$t_{j,i}$ 是统计独立的。

实际情况下，一般有四种不同的并行化技术。

4.1 随机播种（Random seeding）

所有进程使用相同的 PRNG，但使用不同的“随机”种子。

我们希望它们能生成原始 PRNG 的非重叠和不相关子序列。然而，这种方法没有理论基础。并且随机播种违反了 Donald Knuth 的建议——“不应随机选择随机数生成器”。

4.2 参数化（Parameterization）

所有进程都使用相同类型的PRNG，但每个生成器的参数不同。

例如：对于第
$j$ 个随机流 具有加法常数
$b_{j}$ 的线性同余生成器：
$t_{j,i} = a \cdot t_{j,i−1} + b_{j} \mod 2^e$ , 其中
$b_{j}$ 是第
$(j + 2)$ 个素数。另一种变体对不同的流使用不同的乘数
$a$。

这些方法的理论基础薄弱，经验测试表明流之间存在严重的相关性。在大规模并行系统上，可能需要数千个并行流，并且找到一种具有数千个“经过良好测试”的参数集的 PRNG 并非易事。

4.3 块分法（Block splitting）

设
$M$ 为每个进程调用 PRNG 的最大次数，并且
$p$ 为进程数。我们可以将顺序 PRNG 的序列
$r_{i}$ 拆分为长度为
$M$ 的连续块，使得

这种方法只有在我们预先知道
$M$ 或者至少可以靠谱地估计出
$M$ 的上限时才有效。要应用块分法，必须可以保证从第
$i$ 个随机数跳到第
$(i + M)$ 个数而不计算之间的所有数字，这对于许多 PRNG 来说都无法有效地完成（因为大多数 PRNG 是这样的：
$X_{n+1}=f(X_n)$）。

这种方法的一个潜在缺点是，通常在连续模拟中不会出现的长程相关性，反而可能会在子流之间的呈现短程相关性。

块拆分如下图所示：

4.4 跳跃法（Leapfrog）

跳跃方法通过抽取
$r_{i}$ 基本序列，将随机数序列
$r_{i}$ 分布在
$p$ 个进程上，使得：

跳跃法如下图所示。它是相对通用和强大的并行化方法，它不需要先验估计每个处理器将消耗多少随机数（避免了上述分块法需要
$M$ 的限制）。

要有效实现该方法，往往要求一个 PRNG 可以被修改为：只直接生成原始序列中每个第
$p$ 位置的元素。这也会使得一些比较流行的 PRNG 不支持该方式。

4.5 示例：并行线性同余法

线性同余法的公式如下：

现在对其进行展开：

可以看到公式右侧中的两项，第一项中的
$r_{i-P}$ 和左侧项中的
$r_i$ 之间间隔了
$P$ 次迭代，而右侧的第二个求和项不涉及
$r(i)$ 相关内容。

所以：

1. 这个公式可以采用跳跃法来实现并行。跳跃并行法的关键就是要求 PRNG 可以实现跳跃计算，被跳过的数值可以交给其他进程来进行计算。

2. 这个公式也可以采用块分法来实现并行。块分并行法的关键是要求 PRNG 可以实现跳跃计算，并给出线程中待计算的个数。跳过的初值在初始化中计算获得，剩下每个进程计算一个块。

更进一步观察到上述展开式子中的第二项是等比数列求和项：

下面我们简单实现这两种并行方法：

• 跳跃法：

from multiprocessing.dummy import Pool as ThreadPool

def LCG(x0, a, c, m, num, P):
    '''
    线性同余（LCG），并行技术示例：【跳跃法】
        * x0:  第一个开始值；
        * a:   乘数；
        * c:   增量；
        * m:   模数；
        * num: 输出个数；
        * P:   并行个数；
    '''
    items = [(0, x0)]            # 初始值：序号和数值
    out = [x0]*num               # 输出值

    # 【1】计算几个并行的初值 （跳跃点起始值，是连续的）
    for i in range(1,P):
        x0 = (a*x0+c)%m          # x(n+1)=a*x(n)+c mod m
        items.append((i, x0))
        out[i] = x0

    en = num//P                  # 线程中的元素个数
    a_p = a**P
    k = c*(1-a_p)//(1-a)         # 整数等比数列求和必是整数 
    # 【2】定义并行的函数
    def fun(v):
        j, x = v
        for i in range(1, en):
            x = (a_p*x + k)%m    # x(n+j)=a*x(n)+k mod m
            out[i*P+j] = x

    # 【3】启动多线程：
    pool = ThreadPool()
    pool.map(fun, items)
    pool.close()
    pool.join()
    return out

y = LCG(12, 34, 12, 345346, 20, 5)
print(y)

• 块分法：

from multiprocessing.dummy import Pool as ThreadPool

def LCG(x0, a, c, m, num, P):
    '''
    线性同余（LCG），并行技术示例：【块分法】
        * x0:  第一个开始值；
        * a:   乘数；
        * c:   增量；
        * m:   模数；
        * num: 输出个数；
        * P:   并行个数；
    '''
    items = [(0, x0)]            # 初始值：序号和数值
    out = [x0]*num               # 输出值

    # 【1】计算几个并行的初值 （块的起始值）
    en = num//P                  # 线程中的元素个数
    a_p = a**en
    k = c*(1-a_p)//(1-a)         # 整数等比数列求和必是整数
    for i in range(1, P):
        x0 = (a_p*x0 + k)%m      # x(n*i)=a*x(n)+k mod m
        out[i*en] = x0
        items.append((i*en, x0))

    # 【2】定义并行的函数
    def fun(v):
        j, x = v
        for i in range(1, en):
            x = (a*x+c)%m        # x(n+1)=a*x(n)+c mod m
            out[j+i] = x

    # 【3】启动多线程：
    pool = ThreadPool()
    pool.map(fun, items)
    pool.close()
    pool.join()
    return out

y = LCG(12, 34, 12, 345346, 20, 5)
print(y)

输出结果：

非并行：
[12, 420, 14292, 140594, 290710, 214464, 39522, 307722, 102180, 20672, 12168, 68378, 252788, 306500, 60632, 334770, 331120, 207020, 131772, 336108]
块分法：
[12, 420, 14292, 140594, 290710, 214464, 39522, 307722, 102180, 20672, 12168, 68378, 252788, 306500, 60632, 334770, 331120, 207020, 131772, 336108]
跳跃法：
[12, 420, 14292, 140594, 290710, 214464, 39522, 307722, 102180, 20672, 12168, 68378, 252788, 306500, 60632, 334770, 331120, 207020, 131772, 336108]

参考资料

1. Knuth D E. 计算机程序设计艺术（第 2 卷）: 半数值算法[J]. 2002.
2. （美）罗斯著. 概率论基础教程 原书第9版[M]. 北京：机械工业出版社, 2014.01.
3. Bauke H. Tina’s random number generator library[J]. 2011.
4. A. J. Kinderman and J. F. Monahan. 1977. Computer Generation of Random Variables Using the Ratio of Uniform Deviates. ACM Trans. Math. Softw. 3, 3 (Sept. 1977), 257–260.
5. Marsaglia G. Xorshift rngsJ. Journal of Statistical Software, 2003, 8: 1-6.
6. Random number generation – Wikipedia
7. Pseudorandom generator – Wikipedia
8. Continuous uniform distribution – Wikipedia
9. Normal distribution – Wikipedia
10. Exponential distribution – Wikipedia
11. Triangular distribution – Wikipedia
12. Log-normal distribution – Wikipedia
13. Pareto distribution – Wikipedia
14. Weibull distribution – Wikipedia
15. Geometric distribution – Wikipedia
16. Cryptographically-secure pseudorandom number generator – Wikipedia
17. random — Generate pseudo-random numbers — Python 3.10.5 documentation
18. cpython/random.py at main · python/cpython (github.com)
19. Linear congruential generator – Wikipedia
20. Linear-feedback shift register – Wikipedia
21. The Book of Shaders: Random
22. B OX , G. E. P. AND M ULLER , M. E. 1958b. A note on the generation of random normal deviates. Annals Math.Stat. 29, 2, 610–611.
23. MARSAGLIA , G. AND B RAY , T. A. 1964. A convenient method for generating normal variables. SIAM Rev. 6, 3,260–264.
24. （九）详解线性反馈移位寄存器（LFSR） – come back的文章 – 知乎
25. 真正意义的随机数生成器存在吗？ – acalephs的回答 – 知乎
26. 工程上存在那么多不确定情况，为什么计算机不能利用它们产生真随机数，而只能根据逻辑产生伪随机数？ – invalid s的回答 – 知乎
27. 智能优化方法——伪随机数生成算法 – 知乎
28. 谈谈梅森旋转：算法及其爆破 – 孟晨的文章 – 知乎
29. 为什么游戏里的都是伪随机，做不出真随机？希望来个简单易懂的解释？ – 腾讯游戏学堂的回答 – 知乎
30. 采样理论概述（逆变换采样、拒绝采样） – 摸鱼的文章 – 知乎
31. 如何产生正态分布的随机数？ – 孙晓博的回答 – 知乎
32. Reversed Fibonacci LFSR taps – Cryptography Stack Exchange
33. Schneier on Security: Applied Cryptography: Second Edition Errata
34. 谈谈梅森旋转：算法及其爆破 | 始终 (liam.page)

感谢阅读，欢迎在评论区留言讨论哦~

P.S. 如果喜欢本篇文章，请多多 点赞，让更多的人看见我们 :D

关注 公众号「SenseParrots」，获取人工智能框架前沿业界动态与技术思考。