【计算机程序设计艺术 TAOCP】卷2·3.1 引言：程序怎样产生“随机数”

对应原书：第2卷《半数值算法》第3版，§3.1，书页1至8。第3章从本卷开始，主题是随机数。本文按原书的问题顺序展开，例子和 Python 程序为讲解重新设计。

写一个计算平方的函数，我们希望同一个输入总得到同一个结果。写一个随机数生成程序，我们又希望它的输出没有容易利用的规律。这两个要求放在一起，似乎有些冲突。

§3.1 先让我们面对这个问题：如果计算步骤完全确定，所谓“随机”究竟从哪里来？理解这一点，后面的周期、统计检验和分布转换才不会变成互不相干的技巧。

1. 应用需要的是一串数的性质

原书列出随机数的许多用途，包括模拟、抽样、数值分析、程序测试，以及某些带随机选择的算法。

以排队模拟为例，程序可能需要生成顾客到达的间隔、服务所用的时间，再据此估计等待长度。以抽样为例，我们可能希望从一万条记录中选出一百条，让每条记录有相同的入选机会。

这些任务需要的“随机”并不只是一张杂乱的数字表。

若每次到达间隔都符合目标分布，但相邻两次间隔总是一长一短，队列的行为仍然可能被改变。若抽样时每个位置单独看都差不多公平，但总把相邻记录一起选中，样本的联合分布也可能不对。

因此，至少需要分清两个问题：

• **单个输出的分布：**哪些值可以出现，各自有多大概率？
• **多个输出的关系：**前面的值是否影响后面的值，哪些组合会被过多或过少地产生？

本章首先研究0与1之间的一致分布，也常称均匀分布。理想模型是

$$U_0,U_1,U_2,\dots \text{相互独立，且服从 }U[0,1).$$

在这个模型下，长度为 $b-a$ 的区间 $[a,b)$ 得到的概率就是 $b-a$。独立性还要求，例如连续两项都落在 $[0,1/2)$ 中的概率为 $1/4$。

计算机实际给出的只是有限精度的数，因此这个连续模型是分析目标；具体程序还需要说明它怎样把整数状态映射到有限网格上。

2. 可复现是伪随机数的一项用途

设生成程序保存状态 $S_n$，每调用一次执行

$$S_{n+1}=T(S_n),X_n=g(S_n).$$

初始状态由种子确定。如果 $T$ 和 $g$ 都是确定函数，相同初始状态就会重放相同的输出。

在调试模拟程序时，这很有用。一次运行出现异常后，保存初始状态和调用过程，就能回到当时的数据条件。若每次重跑都换成另一串数，我们反而很难判断，是程序修改起了作用，还是输入恰好变了。

但“能够重放”和“适合模拟”是两个需要分别验证的性质。一个计数器

$$0,1,2,3,\dots ,9,0,1,2,3,\dots$$

当然可以重放，每个数字在完整周期中出现得也一样多；它却能被轻易预测，而且相邻数字的关系完全固定。

种子也不是让糟糕方法突然变好的开关。更换种子可能只是在同一个循环上换一个出发位置。后面研究生成器时，我们会把“状态怎样演化”和“初始状态怎样选取”分开讨论。

3. 有限状态为什么必然重复

假设程序只有 $M$ 个可能状态，并且运行过程中没有加入新的外部输入。观察

$$S_0,S_1,\dots ,S_M$$

这 $M+1$ 个状态，按抽屉原理，至少有两个相同。若 $S_i=S_j$ 且 $i<j$，确定性就给出

$$S_{i+t}=S_{j+t},t\ge 0.$$

从这时开始，后面的演化也会重复。轨道因此可以分成两段：进入循环以前的长度 $\mu$，以及循环本身的长度 $\lambda$。

初始状态 → …… → 循环入口 → …… → 循环入口 → ……
             前段长 μ       循环长 λ

这里判断的是完整状态。一个输出值出现两次，并不一定意味着程序状态已经重复。例如，状态中还保存了其他数，或者输出只取了状态的一部分。

反过来，状态重复一定会使后续输出重复。输出序列的最小周期也可能比状态周期更短，因为不同状态可能被 $g$ 映射到同一个值。

这个结论直接说明两件事：无限期不重复不是固定有限状态程序能做到的；周期虽然重要，也只是有限状态结构的一个指标。

4. 平方取中法：复杂的数字变化也可能很快塌缩

原书介绍了平方取中法：保存一个固定长度的十进制整数，平方后补齐位数，取中间几位作为新状态。

以四位状态为例：

$$6100^2=37210000,\text{中间四位为 }2100.$$

继续计算会得到

6100 → 2100 → 4100 → 8100 → 6100

这正是原书展示的一种短循环。另一个问题是0：它的平方仍为0，一旦进入就无法离开。

下面把“补齐八位”明确写进程序，避免把前导零丢掉。orbit 返回重复前出现过的不同状态、前段长度和循环长度。它只适合状态数较少的教学检查，空间占用与访问过的状态数成正比。

def middle_square4(x):
    if not isinstance(x, int) or not 0 <= x < 10000:
        raise ValueError("x must be a four-digit state")
    digits = f"{x * x:08d}"
    return int(digits[2:6])


def orbit(seed, step, max_states):
    if max_states <= 0:
        raise ValueError("max_states must be positive")
    seen, states = {}, []
    x = seed
    while x not in seen:
        if len(states) >= max_states:
            raise ValueError("state bound was too small")
        seen[x] = len(states)
        states.append(x)
        x = step(x)
    mu = seen[x]
    return states, mu, len(states) - mu


states, mu, period = orbit(6100, middle_square4, 10000)
print(states, mu, period)
print(orbit(1000, middle_square4, 10000))
assert orbit(6100, middle_square4, 10000) == (states, mu, period)
assert middle_square4(0) == 0

输出：

[6100, 2100, 4100, 8100] 0 4
([1000, 0], 1, 1)

第二行说明：种子1000只经过一步就进入0的固定点。种子看上去有四位，并不能保证得到很长的轨道。

§3.1 还讨论了一个拼接许多操作的复杂生成方案。它的价值在于反例：增加条件分支、交换数字或组合算术操作，并不会自动获得可证明的好性质。判断方法仍然要回到状态结构与输出行为。

5. “每个数字一样多”为什么不够

考虑长度为100的二进制串：

0101010101……

它包含50个0和50个1。如果只检查单个数字的频数，这个结果完美得近乎刻意。但观察相邻两位，00 和 11 从来不出现；知道当前位，就知道下一位。

再看真正独立的100次公平抛硬币。令1的个数为 $H$，则

$$\mathbb{E}[H]=50,\mathrm{Var}(H)=25.$$

50是均值，不是每次实验必须恰好达到的指标。把输出强行修正成50个0、50个1，会引入新的约束：前99位一旦确定，最后一位便可能被迫取某个值。

所以，检验中既需要寻找明显偏差，也要警惕不合理的“过度整齐”。不过，一个样本恰好平衡本身不构成错误证据；要判断它在指定模型和重复实验中是否异常。

同样，一段很长的连续0也不能单独证明生成器坏了。它的可疑程度取决于观察长度、检查规则和我们提出的问题。先看结果再挑一个最显眼的模式，会改变判断的概率背景。

6. 带着哪些问题进入下一节

接下来可以给每个生成方案记一份简短记录：

要记录什么	它能帮助回答什么
完整状态与更新公式	哪些信息决定后续输出
种子的合法范围	是否会进入零状态或特殊短循环
周期与证明条件	至少多久以后一定重复
输出映射方式	是否只使用了较差的位，是否产生额外偏差
针对联合分布的检验	频数之外还有哪些依赖
实际调用方法	应用取相邻值、隔项值，还是按某种条件取值

最后一行容易遗漏。生成器的数只有进入应用以后，才真正形成到达时间、坐标、排列或样本。应用如何消费输出，也会影响结果。

§3.2 将从最容易写下、也最值得认真分析的一类递推开始：线性同余法。我们要寻找的并非看起来最复杂的公式，而是参数、周期和结构都能够说明白的方法。

---

阅读定位：§3.1 开头的应用与随机性讨论见书页1至3；平方取中法和复杂生成方案的反例见3至6页；随后习题可用于练习周期与有限状态分析。上一篇为卷1·2.6，下一篇为卷2·3.2。

专栏：计算机程序设计艺术 TAOCP。底本：Donald E. Knuth 著，苏运霖译，《计算机程序设计艺术》第2卷《半数值算法》第3版，国防工业出版社，2002年版。