数值计算相关论文&书籍中文翻译

7.0 引言

使用计算机——人类思维所构想出的最精确、最确定的机器——来生成“随机”数，乍看之下似乎有些荒谬。甚至可以说，这在概念上是不可能的。毕竟，任何程序所产生的输出都是完全可预测的，因此并非真正“随机”。

然而，实用的计算机“随机数生成器”已被广泛使用。我们将把这一悖论的深层解释留给计算机时代的哲学家们去解决（参见，例如 Knuth [1] §3.5 的讨论和参考文献）。有时人们会将计算机生成的序列称为伪随机（pseudo-random），而将“随机”一词保留给本质上随机的物理过程所产生的输出，例如将盖革计数器置于某种放射性元素样本旁，两次“咔哒”声之间的时间间隔。我们不会试图做如此精细的区分。

在计算机生成序列的语境下，我们可以给出一个操作性的随机性定义：生成随机序列的确定性程序应当与使用其输出的计算机程序不同，并且在所有可测量的方面都统计上不相关。换句话说，任意两个不同的随机数生成器在与你的特定应用程序结合使用时，应当产生统计上相同的结果。如果结果不同，那么至少其中一个（从你的角度来看）不是一个好的生成器。

上述定义看似循环，因为它是在将一个生成器与另一个进行比较。然而，确实存在大量随机数生成器，它们在极其广泛的应用程序类别中相互满足这一定义。经验上也发现，由物理过程产生的随机数会得到统计上相同的结果。因此，既然我们知道这样的生成器是存在的，那么如何定义它们的问题就可以留给哲学家们去处理了。

因此，从实用主义的角度来看，随机性取决于观察者（或程序员）的眼光。对某个应用来说足够随机的序列，对另一个应用来说可能就不够随机。尽管如此，我们并非完全迷失在不可通约的应用程序海洋中：存在一个被广泛接受的统计检验列表——其中一些检验是合理的，另一些则只是因历史原因而被沿用下来——总体而言，这些检验能非常有效地检测出应用程序（在本例中，就是你的程序）可能发现的任何非随机性。好的随机数生成器应当通过所有这些检验，

或者至少用户应当清楚知道哪些检验未能通过，从而能够判断这些失败是否与当前问题相关。

关于这一主题的参考文献，首推 Knuth [1]。对于 1995 年之前的资料要谨慎使用，因为该领域在随后的十年中取得了巨大进展。

参考文献

Knuth, D.E. 1997, Seminumerical Algorithms, 3rd ed., vol. 2 of The Art of Computer Programming (Reading, MA: Addison-Wesley), Chapter 3, especially §3.5.[1]

Gentle, J.E. 2003, Random Number Generation and Monte Carlo Methods, 2nd ed. (New York: Springer).

7.1 均匀分布随机数（Uniform Deviates）

均匀分布随机数（uniform deviates）就是落在指定范围内的随机数，对于浮点数通常为 0.0 到 1.0，对于整数则通常为 0 到 $2^{32}-1$ 或 $2^{64}-1$。在该范围内，每个数出现的概率都相同。换句话说，这就是你通常所认为的“随机数”。然而，我们要将均匀分布随机数与其他类型的随机数区分开来，例如从具有指定均值和标准差的正态（高斯）分布中抽取的数。我们将在后续章节看到，这些其他类型的随机数几乎总是通过对一个或多个均匀分布随机数进行适当运算而生成的。因此，本节所讨论的可靠均匀分布随机数源，是任何随机建模或蒙特卡洛计算机工作的基本构建模块。

在过去十年中，生成均匀分布随机数的技术取得了显著进步，如今已开始呈现出一个成熟领域的面貌。现在可以合理地期望，每个随机数最多只需十余次算术或逻辑运算即可生成“完美”的随机数，而快速、“足够好”的随机数所需的运算则更少。这一领域的进步得益于三个因素：第一，新的数学算法；第二，对实际陷阱的更好理解；第三，编程语言（尤其是整数运算）的标准化——特别是 C 和 C++ 中无符号 64 位算术的普遍可用性。最后一个因素看似微不足道，却如此重要，这或许有些讽刺。但正如我们将看到的，它确实至关重要。

如今用户面临的最大隐患是，许多过时且低劣的方法仍在普遍使用。以下是一些需要警惕的陷阱：

• 切勿使用主要基于线性同余生成器（LCG）或乘法线性同余生成器（MLCG）的生成器。我们将在下文进一步说明这一点。

• 切勿使用周期小于 $\sim 2^{64}\approx 2\times 10^{19}$ 的生成器，或任何未公开周期的生成器。

• 切勿使用警告不要将其低位比特视为完全随机的生成器。这种建议曾经是合理的，但现在表明该算法已过时（通常是 LCG）。

• 切勿使用 C 和 C++ 语言中的内置生成器，尤其是 rand 和 rand。这些函数没有标准实现，且常常存在严重缺陷。

如果所有因上述陷阱之一或多个而导致结果存疑的科学论文都从图书馆书架上消失，那么每个书架上留下的空隙大概有你的拳头那么大。

你可能还想留意生成器是否过度设计，从而浪费资源：

• 避免使用生成一个 64 位整数或双精度浮点结果需要超过（比如说）二十多次算术或逻辑运算的生成器。

• 避免使用为严肃的密码学用途而（过度）设计的生成器。

• 避免使用周期 $>10^{100}$ 的生成器。你真的永远用不到这么长的周期，而且在超过某个最小界限后，生成器的周期与其质量关系不大。

既然我们已经告诉了你过去应避免什么，那么接下来应立即介绍当前公认的智慧：

一个可接受的随机数生成器必须至少结合两种（理想情况下互不相关）的方法。所结合的方法应独立演化且不共享状态。结合方式应采用简单运算，且不会产生比其操作数更不随机的结果。

如果你不想阅读本节的其余内容，那么请使用以下代码生成你所需的所有均匀分布随机数。这是我们推荐的“背带加腰带、全身盔甲、绝无任何疑虑”的生成器；* 同时，它也符合上述避免浪费和过度设计方法的指导原则。（我们下面推荐的最快生成器，即使将其代码内联到应用程序中，也仅比它快约 2.5 倍。）

struct Ran {
    // 实现推荐的最高质量生成器。构造函数以一个整数种子调用，并创建生成器实例。
    // 成员函数 int64、doub 和 int32 返回随机序列中的下一个值，其变量类型由函数名指示。
    // 该生成器的周期约为 3.138 × 10^57。

    Ullong u, v, w;
    Ran(Ullong j) : v(4101842887655102017LL), w(1) {
        // 构造函数。以任意整数种子调用（但不要使用上面 v 的值）。
        u = j ^ v; int64();
        v = u; int64();
        w = v; int64();
    }
    inline Ullong int64() {
        // 返回 64 位随机整数。方法解释见正文。
        u = u * 2862933555777941757LL + 7046029254386353087LL;
        v ^= v >> 17; v ^= v << 31; v ^= v >> 8;
        w = 4294957665U * (w & 0xffffffff) + (w >> 32);
        Ullong x = u ^ (u << 21); x ^= x >> 35; x ^= x << 4;
        return (x + v) ^ w;
    }
    inline Doub doub() { return 5.42101086242752217E-20 * int64(); }
    // 返回范围在 0 到 1 之间的随机双精度浮点值。
    inline Unt int32() { return (Unt)int64(); }
    // 返回 32 位随机整数。
};

这里的基本前提是：由于随机数生成器在多次调用之间需要维护内部状态，因此它应当是一个对象（或结构体）。你可以声明多个该对象的实例（尽管很难想到这样做的理由），而且不同的实例之间完全不会相互影响。

构造函数 Ran() 接受一个整数参数，该参数将作为生成序列的种子。不同的种子会产生（就所有实际用途而言）完全不同的序列。一旦构造完成，Ran 的实例就提供了几种不同的随机输出格式。具体来说，假设你通过如下声明创建了一个实例：

Ran myran(17);

其中 myran 是该实例的名称，而 17 是其种子。那么，函数 myran.int64() 返回一个随机的 64 位无符号整数；函数 myran.int32() 返回一个 32 位无符号整数；函数 myran.doub() 返回一个双精度浮点数，范围在 0.0 到 1.0 之间。你可以随意混合调用这些函数，并将返回的随机位用于任何目的。例如，如果你需要一个介于 1 到 $n$（含）之间的随机整数，表达式 $1+\text{myran.int64()}\%(n-1)$ 是完全可行的（尽管使用 % 并不是最快的方式）。

在本节的其余部分，我们将简要回顾一些历史（线性同余生成器 LCG 的兴衰），然后详细介绍一些构成优质生成器的算法方法，以及如何组合这些方法。最后，除了上面提到的 Ran 之外，我们还会推荐一些其他生成器。

7.1.1 一些历史

事后看来，整个随机数生成领域似乎被如下简单的递推公式迷住了太长时间：

$$I_{i+1}=a{I}_i+c(\mathrm{mod}m)$$

其中 $m$ 称为模数，$a$ 是一个正整数，称为乘数，$c$（可以为零）是一个非负整数，称为增量。当 $c\ne 0$ 时，公式 (7.1.1) 被称为线性同余生成器（LCG）。当 $c=0$ 时，有时称为乘法 LCG（MLCG）。

递推式 (7.1.1) 最终必然重复，其周期显然不超过 $m$。如果 $m$、$a$ 和 $c$ 选择得当，则周期将达到最大长度，即长度为 $m$。在这种情况下，0 到 $m-1$ 之间的所有整数都会在某个时刻出现，因此任意初始种子 $I_0$ 的选择都和其他选择一样好：

序列从此点开始，后续的值 $I_j$ 即为返回的“随机”值。

LCG 的思想可追溯到计算的黎明时期，并在 1950 年代及之后被广泛使用。然而，早在 1960 年代中期，人们就开始注意到其中的问题（例如 [1]）：如果每次取 $k$ 个随机数用于在 $k$ 维空间中绘制点（每个坐标在 0 到 1 之间），这些点并不会“填满”整个 $k$ 维空间，而是会落在 $(k-1)$ 维的“平面上”。这样的平面最多大约有 $m^{1/k}$ 个。如果常数 $m$ 和 $a$ 没有非常仔细地选择，实际的平面数量会远少于这个上限。通常 $m$ 接近机器所能表示的最大整数，通常是 $\sim 2^{32}$。例如，在三维空间中，三点所构成的平面数量最多约为 $2^{32}$ 的立方根，即约 1600。如果你关注的物理过程只发生在总体积的一小部分内，那么这种平面离散性可能会非常明显。

更糟糕的是，许多早期的生成器恰好为 $m$ 和 $a$ 做了特别糟糕的选择。其中一个臭名昭著的例子是 RANDU，其 $a=65539$，$m=2^{31}$，在 IBM 大型机上流行多年，并被广泛复制到其他系统中。我们中的一位还记得自己读研究生时曾画出一个只有 11 个平面的“随机”图，结果被计算机中心的编程顾问告知他误用了随机数生成器：“我们保证每个数单独来看是随机的，但我们不保证多个数一起也是随机的。” 这至少耽误了我们一年的研究生学业！

LCG 和 MLCG 还有其他弱点：当 $m$ 选为 2 的幂（例如 RANDU）时，生成的低位几乎完全不随机。具体来说，最低有效位的周期最多为 2，第二位最多为 4，第三位最多为 8，依此类推。然而，如果你不将 $m$ 选为 2 的幂（实际上，通常选择 $m$ 为素数更好），那么你通常需要使用双倍长度的寄存器来执行公式 (7.1.1) 中的乘法和取模运算。而当时的计算机（甚至现在通常也）往往不具备这种能力。

随后，人们投入了大量精力来“修复”这些弱点。一种优雅的数论测试——谱测试（spectral test）被开发出来，用于刻画任意维空间中平面的密度。（参见 [2] 的最新综述，其中包含了一些历史上使用过的极其糟糕的生成器的图形展示，另见 [3]。）Schrage 方法 [4] 被发明出来，使得在 $m$ 高达 $2^{32}-1$ 的情况下，仅使用 32 位算术即可完成乘法 $a{I}_j$，但不幸的是，这种方法仅适用于某些特定的 $a$ 值，而这些值并不总是最优的。Park 和 Miller 的综述 [5] 很好地描绘了 LCG 在其鼎盛时期的状况。

回顾过去，很明显，该领域长期专注于 LCG 是有些误入歧途的。从技术上讲，过去十年中出现的更好的非 LCG 生成器本可以在几十年前就被发现；同样，人们本可以更早放弃对优雅“单一算法”生成器的不切实际的幻想（转而采用更务实的组合生成器拼接方法）。正如我们下面将要解释的，LCG 和 MLCG 仍然可以有用，但仅限于经过仔细控制的场景，并且必须充分注意其明显的弱点。

7.1.2 组合生成器中推荐使用的方法

如今，至少有十几种合理的算法值得认真考虑用于随机数生成器。我们对其中几种的选择既出于数学考量，也出于审美偏好。我们喜欢操作少而快、初始化简单可靠、状态足够小以便存放在寄存器或一级缓存中的算法（前提是编译器和硬件支持）。这意味着我们倾向于避免那些状态为某个长度数组的算法，尽管这类算法能相对简单地实现极其巨大的周期。关于更广泛方法的综述，参见 [6] 和 [7]。

要被推荐用于组合生成器，一种方法必须在理论上得到一定程度的理解，并通过一套相当广泛的实证测试（或者，即使失败，其弱点也必须被充分刻画）。我们的最低理论标准是：周期、返回值集合以及有效初始化集合都必须被完全理解。作为最低实证标准，我们使用了 Marsaglia 于 2003 年发布的第二版名为 Diehard 的统计测试套件 [8]。也可以使用替代测试套件 NIST-STS [9]，或与 Diehard 一起使用。

仅仅要求组合生成器通过 Diehard 或 NIST-STS 并不是一个足够严格的标准。这些测试套件仅对生成器调用约 $10^7$ 次，而用户程序可能会调用 $10^{12}$ 次或更多。更有意义的要求是：组合生成器中的每种方法都必须单独通过所选测试套件。这样，组合生成器（如果构造正确）的性能应远优于任何单一组件。在下表中，我们使用符号“$\ast$”表示该方法自身能通过 Diehard 测试。（对于 64 位量，该声明指高 32 位和低 32 位各自都能通过。）相应地，下文中的“可用作随机”并不意味着完美随机，而仅表示在不需要更好组合生成器的快速简易应用中达到最低水平。

我们现在转向具体方法，首先介绍使用 64 位无符号算术的方法（我们称之为 Ullong，即 Linux/Unix 中的 unsigned long long，或微软平台上的 unsigned __int64）。

(A) 64 位 XORshift 方法。该生成器由 Marsaglia [10] 发现并刻画。仅需三次 XOR 和三次移位（通常是快速操作），即可在 64 位上产生 $2^{64}-1$ 的完整周期。（缺失的值是零，它会自我延续，必须避免。）高、低位均能通过 Diehard 测试。生成器可以使用下面以 >> 开头的三行更新规则，或以 << 开头的规则。（这两种更新规则产生不同的序列，彼此通过位反转相关。）

$$\begin{array}{rl}& \text{state:}x(\text{unsigned 64-bit})\\ & \text{initialize:}x\ne 0\\ & \text{update:}x\leftarrow x\oplus (x\gg a_1),\\ & x\leftarrow x\oplus (x\ll a_2),\\ & x\leftarrow x\oplus (x\gg a_3);\\ & \text{or}x\leftarrow x\oplus (x\ll a_1),\\ & x\leftarrow x\oplus (x\gg a_2),\end{array}$$

可用作随机：$x\leftarrow x\oplus (x\ll a_3)$；

可用于位混合：$x$（所有位）$\ast$

可用于位混合：$x$（所有位）

可通过以下方式改进：输出 64 位 MLCG 的后继值

周期：$2^{64}-1$

以下是这些生成器背后理论的简要概述：将整数的 64 位视为长度为 64 的向量在模 2 线性空间中的分量。注意到 XOR（$\oplus$）等同于加法，更新规则中的每一行都可以写成一个 $64\times 64$ 矩阵对向量的作用，其中矩阵除对角线和恰好一条上对角线或下对角线（对应 << 或 >>）外全为零。将此矩阵记为 ${\mathbf{S}}_k$，其中 $k$ 是移位参数（例如，左移为正，右移为负）。那么，一次完整的更新步骤（上述更新规则的三行）对应于矩阵乘法 $\mathbf{T}\equiv {\mathbf{S}}_{k_1}{\mathbf{S}}_{k_2}{\mathbf{S}}_{k_3}$。

接下来需要找到能给出完整周期 $M\equiv 2^{64}-1$ 的整数三元组 $(k_1,k_2,k_3)$，例如 (21, -35, 4)。充要条件是 ${\mathbf{T}}^M=1$（单位矩阵），且对于以下七个 $N$ 值，${\mathbf{T}}^N\ne 1$：$M/6700417$、$M/65537$、$M/641$、$M/257$、$M/17$、$M/5$ 和 $M/3$，即 $M$ 除以其七个不同素因子的结果。所需的 $\mathbf{T}$ 的大幂次可通过连续平方快速计算，仅需约 $64^4$ 次操作。利用这一机制，可以通过穷举搜索以合理代价找到完整周期的三元组 $(k_1,k_2,k_3)$。

Brent [11] 指出，64 位 xorshift 方法在每一位上产生的比特序列，与某个特定的 64 位线性反馈移位寄存器（LFSR）所产生的序列完全相同。（我们将在 §7.5 中进一步了解 LFSR。）因此，xorshift 方法可能具有与 LFSR 类似的某些优点和缺点。然而，一个典型的 xorshift 生成器所对应的本原多项式具有大量非零项，这使其统计性质优于基于本原三项式（primitive trinomials）的 LFSR 生成器，从而在一定程度上缓解了上述问题。实际上，xorshift 生成器是一种仅使用六次快速的 64 位操作，同时推进 64 个非平凡的一比特 LFSR 寄存器的方法。虽然还有其他对 LFSR 进行快速推进的方法，以及组合多个此类生成器输出的方法 [12,13]，但都不如 xorshift 方法简单。

尽管 xorshift 生成器中每一位都遵循相同的递推关系，因此具有相同的周期为 $2^{64}-1$ 的序列，但该方法保证了各序列之间存在偏移，使得在一个完整周期内，所有非零的 64 位字都会在各位上出现（正如我们刚刚所见）。

文献 [10] 中列出了若干个具有完整周期的三元组。但其中只有很小一部分能通过 Diehard 测试。此外，某个三元组可能在其 <<-first 版本中通过测试，但在 >>-first 版本中失败，反之亦然。由于这两个版本仅产生比特顺序相反的序列，显然任一版本的失败都应被视为两者共同的失败（同时也暴露了 Diehard 测试的某种弱点）。以下推荐的参数集对 << 和 >> 规则均能通过 Diehard 测试。列表顶部附近的参数集可能略优于底部附近的参数集。ID 列为每个推荐的生成器分配了一个标识字符串，我们将在后文引用。

ID	$a_1$	$a_2$	$a_3$
A1	21	35	4
A2	20	41	5
A3	17	31	8
A4	11	29	14
A5	14	29	11
A6	30	35	13
A7	21	37	4
A8	21	43	4
A9	23	41	18

如果单独使用 xorshift 生成器，很容易设计出使其失败的测试。因为在第 $i+1$ 步中，每一位最多只依赖于第 $i$ 步中的 8 个比特，因此对两个时间步进行某些简单的逻辑组合（并配合适当的掩码）即可立即揭示其非随机性。此外，当状态经过汉明权重（即 1 的个数）很小的值时（这种状态最终必然会出现，称为低汉明权重状态），其恢复速度会比预期更慢。尽管如此，当与其他方法组合使用时，xorshift 生成器是一种异常强大且实用的方法。如果在 1949 年发现的是 xorshift 而不是 LCG，许多麻烦本可避免！

(B) 带进位乘法（Multiply with Carry, MWC），基数 $b=2^{32}$。该方法同样由 Marsaglia 发现，MWC 生成器的基数 $b$ 最方便地选为可用字长的一半（即对于 64 位字，$b=2^{32}$）。MWC 由其乘数 $a$ 定义。

状态：$x$（无符号 64 位整数）

初始化：$1\le x\le 2^{32}-1$

更新：$x\leftarrow a\cdot (x\&(2^{32}-1))+(x\gg 32)$

可用作随机数：$x$ 的低 32 位 $\ast$

可用于比特混合：$x$ 的全部 64 位

可改进方式：输出 64 位 xorshift 后继值作为 64 位 $x$

周期：$(2^{32}a-2)/2$（一个素数）

理论上，参数为 $b$ 和 $a$ 的 MWC 生成器与模数 $m=ab-1$、乘数为 $a$ 的 LCG 相关 [14]，但并不完全相同。很容易找到使 $m$ 为素数的 $a$ 值，从而仅使用 2 的幂模运算即可获得素数模数的好处。虽然无法选择 $a$ 使得周期达到最大值 $m$，但如果选择 $a$ 使得 $m$ 和 $(m-1)/2$ 均为素数，则 MWC 的周期为 $(m-1)/2$，几乎同样优秀。这样选出的部分 $a$ 值能通过标准统计测试套件；谱测试 [14] 是一个有前景的发展方向，但本文未予采用。

尽管只有状态 $x$ 的低 $b$ 位可视为算法意义上的随机，但代表乘积 $ab$ 的所有 $x$ 的比特中仍包含大量随机性。这在组合生成器中非常方便，允许将整个状态 $x$ 作为组件使用。事实上，下面推荐的前两个 $a$ 值使得 $ab$ 非常接近 $2^{64}$（误差约 2 ppm），以至于 $x$ 的高比特位实际上能通过 Diehard 测试。（这是一个很好的例子，说明任何测试套件都可能无法发现少量高度非随机的行为，例如此处高 32 位中可能缺失多达 8000 个值。）

除上述考虑外，以下数值推荐使用，无特定排序。

ID	$a$
B1	4294957665
B2	4294963023
B3	4162943475
B4	3947008974
B5	3874257210
B6	2936881968
B7	2811536238
B8	2654432763
B9	1640531364

(C) 模 $2^{64}$ 的 LCG。既然上面已对其大加贬斥，为何此处还要包含这种生成器？对于下面给出的参数（它们能很好地通过谱测试），其高 32 位几乎（但并未完全）通过 Diehard 测试，而低 32 位则完全失败。然而，正如我们在讨论组合生成器构造时将看到的，它仍有一席之地。以下推荐的乘数 $a$ 具有良好的谱特性 [15]。

状态：$x$（无符号 64 位整数）

初始化：任意值

更新：$x\leftarrow ax+c(\mathrm{mod}{2}^{64})$

可用作随机数：$x$ 的高 32 位（需谨慎）

可用于比特混合：$x$ 的高 32 位

可改进方式：输出 64 位 xorshift 后继值

周期：$2^{64}$

ID	a	c（任意奇数值均可）
C1	3935559000370003845	2691343689449507681
C2	3202034522624059733	4354685564936845319
C3	2862933555777941757	7046029254386353087

(D) 模 $2^{64}$ 的 MLCG（乘法同余生成器）。与前述生成器类似，该生成器的有用角色也极为有限。其低比特位高度非随机。推荐的乘数具有良好的谱特性（部分来自 [15]）。

状态：$x$（无符号 64 位整数）

初始化：$x\ne 0$

更新：$x\leftarrow ax(\mathrm{mod}{2}^{64})$

可用作随机数：$x$（高32位，需谨慎使用）

可用于位混合：$x$（高32位）

可改进方式：输出64位 xorshift 后继值

周期：$2^{62}$

ID	a
D1	2685821657736338717
D2	7664345821815920749
D3	4768777513237032717
D4	1181783497276652981
D5	702098784532940405

(E) 模数 $m\gg 2^{32}$ 且 $m$ 为素数的 MLCG。当支持64位无符号整数运算时，使用素数模数和具有良好谱特性的大乘数的 MLCG 是不错的32位生成器。其主要缺点是：为了仅获得32位（左右）的结果，却需要执行代价高昂的64位乘法和64位取余运算。

$$\begin{array}{rlr}& \mathrm{state}:& x(\mathrm{unsigned}64\text{-}\mathrm{bit})\\ & \mathrm{initialize}:& 1\le x\le m-1\\ & \mathrm{update}:& x\leftarrow ax(\mathrm{mod}m)\\ & \mathrm{can}\mathrm{use}\mathrm{as}\mathrm{random}:& x(1\le x\le m-1)\mathrm{or}\mathrm{low}32\mathrm{bits}\ast \\ & \mathrm{can}\mathrm{use}\mathrm{in}\mathrm{bit}\mathrm{mix}:& (\mathrm{same})\\ & \mathrm{period}:& m-1\end{array}$$

下表中的参数值由 P. L'Ecuyer 为我们精心计算得出。对于接近2的幂次的素数模数，这些乘数几乎是所能获得的最佳选择。尽管推荐仅使用低32位（这些位均通过了 Diehard 测试），但可以看出（取决于模数），最多可获得43位质量尚可的随机位，代价是执行64位乘法和取余运算。

ID	$m$	$a$
E1	$2^{39}-7=549755813881$	10014146 30508823 25708129
E2 E3
E4	$2^{41}-21=2199023255531$	5183781 1070739 6639568
E5 E6
E7	$2^{42}-11=4398046511093$	1781978 2114307 1542852
E8 E9
E10	$2^{43}-57=8796093022151$	2096259 2052163 2006881

(F) 模数 $m\gg 2^{32}$、$m$ 为素数，且满足 $a(m-1)\approx 2^{64}$ 的 MLCG。一种变体用于组合生成器：选择 $m$ 和 $a$，使得 $a(m-1)$ 尽可能接近 $2^{64}$，同时仍要求 $m$ 为素数，且 $a$ 通过谱测试。这样做的目的是使 $ax$ 成为一个64位值，其高位具有良好随机性，适用于组合生成器。然而，乘法和取余运算的高昂开销仍是主要缺点。$x$ 的低32位的随机性并不显著劣于前述 MLCG 生成器 E1–E12。

状态：$x$（无符号64位）

初始化：$1\le x\le m-1$

更新：$x\leftarrow ax(\mathrm{mod}m)$

可用作随机数：$x(1\le x\le m-1)$ 或低32位 $\ast$

可用于位混合：$ax$（但不要同时使用 $ax$ 和 $x$）$\ast$

可改进方式：输出 $ax$ 的64位 xorshift 后继值

周期：$m-1$

ID	$m$	$a$
F1	$1148\times 2^{32}+11=4930622455819$	3741260
F2	$1264\times 2^{32}+9=5428838662153$	3397916
F3	$2039\times 2^{32}+3=8757438316547$	2106408

7.1.3 如何构造组合生成器

虽然构造组合生成器是一门艺术，但仍需以数学原理为指导。关于组合生成器的严格定理通常仅在所组合的生成器在算法上相关时才可能成立；但根据“不要把所有鸡蛋放在一个篮子里”的一般原则，这种做法本身通常是不可取的。因此，我们只能依赖一些指导原则和经验法则。

被组合的方法应相互独立。它们不得共享状态（尽管其初始化可源自某个方便的共同种子）。它们应具有不同且不可通约的周期。理想情况下，它们在算法上应尽可能彼此不同。后一标准必然涉及一定的艺术性。

组合生成器的输出绝不应干扰各独立方法的演化过程，组合操作本身也不应产生任何副作用。

组合方法应通过二元运算进行，这些运算需满足：当一个输入固定时，输出的随机性不低于另一个输入。在32位或64位无符号整数运算中，实践中这意味着只能使用 + 和 ∧ 运算符。例如，乘法是禁止使用的：若一个操作数是2的幂，则无论另一个操作数多么随机，其乘积末尾都会出现若干个零。

组合输出中的每一位都应依赖于至少两种方法的高质量位，也可额外依赖其他方法的低质量位。在上表中，标记为“可用作随机数”的位被视为高质量；标记为“可用于位混合”的位被视为低质量，除非它们也通过了如 Diehard 等统计测试套件。

我们还有一种技巧可用，即把某个方法用作后继关系，而非独立的生成器。上述每种方法都是从某个64位状态 $x_i$ 到唯一后继状态 $x_{i+1}$ 的映射。若一个方法能通过良好的统计测试套件，则其状态与其后继之间应无任何可检测的相关性。此外，若该方法的周期为 $2^{64}$ 或 $2^{64}-1$，则所有值（可能除零外）都会恰好作为后继状态出现一次。

假设我们取一个周期为 $2^{64}$ 的生成器（例如上面的 C1）的输出，并将其作为后继关系送入周期为 $2^{64}-1$ 的生成器 A6。这可方便地记作 “A6(C1)”，我们称之为复合生成器。注意，复合输出绝不会反馈回 C1 的状态，C1 的演化不受任何干扰。复合生成器 A6(C1) 的周期等于 C1 的周期，不幸的是，并非两个周期的乘积。但 A6 对 C1 输出值的随机映射有效地修正了 C1 低有效位短周期的问题（若使用先左移形式的 A6，则效果更佳）。同时，A6(C1) 也弥补了 A6 的弱点——即某一位仅依赖于前一状态的少数几位。因此，我们将精心构造的复合生成器视为一种组合生成器，其地位与通过 + 或 ∧ 直接组合相当。

复合方式相比直接组合的劣势在于：其计算开销几乎相同，却未增加状态大小或周期长度。但其优势在于能更有效地混合差异较大的位位置。在前述例子中，我们不会接受 A6 + C1 作为组合生成器，因为 C1 的低有效位质量太差，对组合几乎无贡献；但 A6(C1) 无此缺陷，且优点显著。在前述各方法的汇总表中，我们已在“可改进方式”一栏中指明了推荐用于复合生成器的组合。

现在，我们可以用如下伪等式完整描述前面提到的 Ran 生成器：

$$\mathrm{Ran}=[{A1}_l(C3)+{A3}_r]\wedge B1$$

也就是说，这是四种不同生成器的组合和/或复合。对于方法 A1 和 A3，下标 $l$ 和 $r$ 表示首先执行的是左移还是右移操作。Ran 的周期是 C3、A3 和 B1 各自周期的最小公倍数。

我们能够轻松推荐的最简单且最快的生成器是

$$Ranq1\equiv D1(\mathrm{A}{1}_{\mathrm{r}})$$

其实现如下：

struct Ranq1 {
    // 推荐用于日常使用的生成器。周期约为 $1.8 \times 10^{19}$。
    // 调用约定与上面的 Ran 相同。
    Ullong v;
    Ranq1(Ullong j) : v(4101842887655102017LL) {
        v ^= j;
        v = int64();
    }
    inline Ullong int64() {
        v ^= v >> 21; v ^= v << 35; v ^= v >> 4;
        return v * 2685821657736338717LL;
    }
    inline Doub doub() { return 5.42101086242752217E-20 * int64(); }
    inline Unt int32() { return (Uint)int64(); }
};

Ranq1 在 3 次移位、3 次异或和 1 次乘法操作后生成一个 64 位随机整数，或者再额外进行一次乘法操作生成一个双精度浮点值。其实现足够简洁，可以轻松内联到应用程序中。它的周期“仅有” $1.8\times 10^{19}$，因此不应被调用次数超过 $\sim 10^{12}$ 的应用程序使用。在此限制下，我们认为 Ranq1 足以满足 99.99% 的用户应用程序需求，而 Ran 可保留给剩下的 0.01% 使用。

如果你对 Ranq1 的“较短”周期感到不安（其实大可不必），可以改用

$$Ranq2\equiv {A3}_r\wedge B1$$

其周期为 $8.5\times 10^{37}$。

// 如果 Ranq1 周期太短而 Ran 又太慢时的备用生成器。
// 周期约为 $8.5 \times 10^{37}$。调用约定与上面的 Ran 相同。
Ullong v, w;
Ranq2(Ullong j) : v(4101842887655102017LL), w(1) {
    v ^= j;
    w = int64();
    v = int64();
}
inline Ullong int64() {
    v ^= v >> 17; v ^= v << 31; v ^= v >> 8;
    w = 4294957665U * (w & 0xffffffff) + (w >> 32);
    return v ^ w;
}
inline Doub doub() { return 5.42101086242752217E-20 * int64(); }
inline Unt int32() { return (Uint)int64(); }

7.1.4 随机哈希与随机字节

偶尔，你会需要一个随机序列 $H_i$，其值可以按任意顺序访问或重新访问。也就是说，你需要对整数 $i$ 进行随机哈希，即使对于高度有序的 $i$ 序列，该哈希也能通过严格的随机性测试。用我们已有的术语来说，你需要一个完全无状态的生成器，它完全基于后继关系构建，并以值 $i$ 作为起点。

一个轻松通过 Diehard 测试的例子是：

$$\mathrm{Ranhash}\equiv A{2}_l(D3(A{7}_r(C1(i))))$$

注意其中交替使用了基于 64 位乘法的后继关系和基于移位与异或的后继关系。

struct Ranhash {
    // 将整数高质量地哈希为多种数值类型。
    inline Ullong int64(Ullong u) {
        // 返回 u 的 64 位整数哈希值。
        Ullong v = u * 3935559000370003845LL + 2691343689449507681LL;
        v ^= v >> 21; v ^= v << 37; v ^= v >> 4;
        v *= 4768777513237032717LL;
        v ^= v << 20; v ^= v >> 41; v ^= v << 5;
        return v;
    }
    inline Unt int32(Ullong u) {
        // 返回 u 的 32 位整数哈希值。
        return (Uint)(int64(u) & 0xffffffff);
    }
    inline Doub doub(Ullong u) {
        // 返回 u 的双精度浮点哈希值（范围在 0 到 1 之间）。
        return 5.42101086242752217E-20 * int64(u);
    }
};

由于 Ranhash 无状态，因此没有构造函数。你只需在需要时直接调用其 int64(i) 函数（或其他函数），传入你的 $i$ 值即可。

随机字节。在另一种情况下，你可能希望逐字节地生成随机整数。当然，你可以从上述任何推荐的组合生成器中提取字节，因为它们在所有位上都具有同等质量。以下代码可添加到上述任意生成器中，为其增加一个 int8() 方法。（务必在构造函数中将 bc 初始化为零。）

Ullong breg;
Int bc;
inline unsigned char int8() {
    if (bc >= 0) return (unsigned char)(breg >> (8 * bc--));
    breg = int64();
    bc = 7;
    return (unsigned char)breg;
}

如果你想要一个更面向字节（尽管不一定更快）的算法，一个有趣的选择是 Rivest 的 RC4，它被广泛用于许多互联网应用中。RC4 最初是 RSA 公司的专有算法，但仅作为商业秘密保护，并未申请专利或版权。结果是，当该秘密于 1994 年被匿名发布到互联网上后，RC4 在几乎所有方面都成为了公共财产。RC4 这个名称仍受保护，是 RSA 的商标。因此，为了严谨起见，我们将以下实现命名为 Ranbyte。

struct Ranbyte {
    // 使用通常称为 RC4 的算法生成随机字节。
    Int s[256], i, j, ss;
    Uint v;
    Ranbyte(Int u) {
        // 构造函数。使用任意整数种子调用。
        v = 2244614371U ^ u;
        for (i = 0; i < 256; i++) { s[i] = i; }
        for (j = 0, i = 0; i < 256; i++) {
            ss = s[i];
            j = (j + ss + (v >> 24)) & 0xff;
            s[i] = s[j]; s[j] = ss;
            v = (v << 24) | (v >> 8);
        }
        i = j = 0;
        for (Int k = 0; k < 256; k++) int8();
    }
    inline unsigned char int8() {
        // 返回序列中的下一个随机字节。
        i = (i + 1) & 0xff;
        ss = s[i];
        j = (j + ss) & 0xff;
        s[i] = s[j]; s[j] = ss;
        return (unsigned char)(s[(s[i] + s[j]) & 0xff]);
    }
    Uint int32() {
        // 由 4 个随机字节构造一个随机 32 位整数。较慢。
        v = 0;
        for (int k = 0; k < 4; k++) {
            i = (i + 1) & 0xff;
            ss = s[i];
            j = (j + ss) & 0xff;
            s[i] = s[j]; s[j] = ss;
            v = (v << 8) | s[(s[i] + s[j]) & 0xff];
        }
        return v;
    }
    Doub doub() {
        // 返回一个 0 到 1 之间的随机双精度浮点值。很慢！
        return 2.32830643653869629E-10 * (int32() +
               2.32830643653869629E-10 * int32());
    }
};

注意，Ranbyte 实例的初始化开销很大，因此不应在频繁执行的循环内部创建其实例。返回 32 位整数或双精度浮点值的方法相比上述其他生成器较慢，但提供这些方法是为了让你在需要时可以用 Ranbyte 作为其他（可能有问题的）生成器的测试替代品。

如果你在 Ranbyte 中发现任何非随机性，请不要告诉我们。但有几个国家密码机构可能会（也可能不会）想和你谈谈！

7.1.5 更快的浮点值生成

上述将 64 位整数转换为双精度浮点值的步骤涉及非平凡的类型转换和 64 位浮点乘法，这会成为性能瓶颈。另一种方法是直接通过联合体（union）、掩码和一些 64 位逻辑操作将随机位放入双精度浮点数的正确位置；但根据我们的经验，这种方法并没有显著更快。

如果对浮点值生成速度有绝对要求，我们需要稍微放宽一些设计规则。以下是“Knuth 的减法生成器”的一个变体，它是一种在包含 55 个值的循环列表上、滞后量为 24 和 55 的所谓延迟斐波那契生成器。其有趣之处在于，新值直接以浮点形式生成，通过对两个先前值进行浮点减法得到。

struct Ranfib {
    // 仅使用浮点运算实现 Knuth 的减法生成器。注意事项见正文。
    Doub dtab[55], dd;
    Int inext, inextp;
    Ranfib(Ullong j) : inext(0), inextp(31) {
        // 构造函数。使用任意整数种子调用。使用 Ranq1 进行初始化。
        Ranq1 init(j);
        for (int k = 0; k < 55; k++) dtab[k] = init.doub();
    }
    Doub doub() {
        // 返回 0 到 1 之间的随机双精度浮点值。
        if (++inext == 55) inext = 0;
        if (++inextp == 55) inextp = 0;
        dd = dtab[inext] - dtab[inextp];
        if (dd < 0) dd += 1.0;
        return (dtab[inext] = dd);
    }
    inline unsigned long int32() {
        // 返回随机 32 位整数。仅推荐用于测试目的。
        return (unsigned long)(doub() * 4294967295.0);
    }
};

int32 方法仅用于测试或偶然用途。另请注意，我们使用 Ranq1 来初始化 Ranfib 的 55 个随机值表。参见 Knuth 或《Numerical Recipes》早期版本中（略显笨拙的）纯内部初始化方法。

Ranfib 未能通过 Diehard 的“生日测试”，该测试能够识别滞后 0、24 和 55 处三个值之间的简单关系。除此之外，它是一个不错的（但并非卓越的）生成器，其主要优势在于速度。

7.1.6 计时结果

计时结果高度依赖于特定的硬件和编译器细节，因此很难判断单组测试是否具有任何实际意义。对于组合生成器尤其如此，因为优秀的编译器或具备复杂指令前瞻能力的 CPU 可以交错并流水线化各个方法的操作，直至最终的组合操作。此外，目前台式计算机正处于从 32 位向 64 位过渡的阶段，这将影响 64 位操作的计时结果。因此，你应该熟悉 C 语言中的 clock_t clock(void) 功能，并自行进行实验。

尽管如此，下表给出了本节中各例程在 2004 年款 3.4 GHz Pentium CPU 上的典型结果（已归一化）。单位为每秒返回值数量（$10^6$）。数值越大越好。

生成器	int64()	doub()	int8()
Ran	19	10	51
Ranq1	39	13	59
Ranq2	32	12	58
Ranfib		24
Ranbyte

表中 Ran、Ranq1 和 Ranq2 的 int8() 计时结果指的是如上所述经过增强的版本。

7.1.7 当你仅有 32 位算术运算时

我们最好的建议是：换一个更好的编译器！但如果你确实必须生活在仅有无符号 32 位算术运算的世界中，那么以下是一些备选方案。这些方案单独使用时均无法通过 Diehard 测试。

(G) 32 位 XORshift 随机数生成器

状态：$x$（无符号 32 位）

初始化：$x\ne 0$

更新： $x\leftarrow x\oplus (x\gg b_1)$,
$x\leftarrow x\oplus (x\ll b_2)$,
$x\leftarrow x\oplus (x\gg b_3)$；

或
$x\leftarrow x\oplus (x\ll b_1)$,
$x\leftarrow x\oplus (x\gg b_2)$,
$x\leftarrow x\oplus (x\ll b_3)$；

可用作随机数：$x$（32 位，需谨慎使用）

可用于位混合：$x$（32 位）

可改进方法：输出 32 位 MLCG 的后继值

周期：$2^{32}-1$

ID	$b_1$	$b_2$	$b_3$
G1	13	17	5
G2	7	13	3
G3	9	17	6
G4	6	13	5
G5	9	21	2
G6	17	15	5
G7	3	13	7
G8	5	13	6
G9	12	21	5

(H) 基数 $b=2^{16}$ 的 MWC

状态：$x$, $y$（无符号 32 位）

初始化：$1\le x$, $y\le 2^{16}-1$

更新： $x\leftarrow a\cdot (x\&(2^{16}-1))+(x\gg 16)$
$y\leftarrow b\cdot (y\&(2^{16}-1))+(y\gg 16)$

可用作随机数：$(x<<16)+y$

可在位混合中使用：相同，或（谨慎使用）$x$ 或 $y$

可通过以下方式改进：输出 32 位 xorshift 后继值

周期：$(2^{16}a-2)(2^{16}b-2)/4$（两个素数的乘积）

ID	a	b
H1	62904	41874
H2	64545	34653
H3	34653	64545
H4	57780	55809
H5	48393	57225
H6	63273	33378

(I) 模 $2^{32}$ 的线性同余生成器（LCG）

状态：$x$（无符号 32 位整数）
初始化：任意值
更新：$x\leftarrow ax+c(\mathrm{mod}{2}^{32})$
是否可用于随机数生成：不推荐
是否可用于位混合：不推荐
可通过以下方式改进：输出 32 位 xorshift 后继值
周期：$2^{32}$

ID	$a$	$c$（任意奇数均可）
I1	1372383749	1289706101
I2	2891336453	1640531513
I3	2024337845	797082193
I4	32310901	626627237
I5	29943829	1013904223

(J) 模 $2^{32}$ 的乘法线性同余生成器（MLCG）

状态：$x$（无符号 32 位整数）
初始化：$x\ne 0$
更新：$x\leftarrow ax(\mathrm{mod}{2}^{32})$
是否可用于随机数生成：不推荐
是否可用于位混合：不推荐
可通过以下方式改进：输出 32 位 xorshift 后继值
周期：$2^{30}$

ID	$a$
J1	1597334677
J2	741103597
J3	1914874293
J4	990303917
J5	747796405

一个高质量（尽管稍慢）的组合生成器为：

$$Ranlim32\equiv [{G3}_l(I2)+{G1}_r]\wedge [{G6}_l({H6}_h)+{H5}_h]$$

实现如下：

Uint u, v, w1, w2;
Ranlim32(Uint j) : v(2244614371U), w1(521288629U), w2(362436069U) {
    u = j ^ v; int32();
    v = u; int32();
}

inline Uint int32() {
    u = u * 2891336453U + 1640531513U;
    v ^= v >> 13; v ^= v << 17; v ^= v >> 5;
    w1 = 33378 * (w1 & 0xffff) + (w1 >> 16);
    w2 = 57225 * (w2 & 0xffff) + (w2 >> 16);
    Uint x = u ^ (u << 9); x ^= x >> 17; x ^= x << 6;
    Uint y = w1 ^ (w1 << 17); y ^= y >> 15; y ^= y << 5;
    return (x + v) ^ (y + w2);
}

inline Doub doub() { return 2.32830643653869629E-10 * int32(); }
inline Doub truedoub() {
    return 2.32830643653869629E-10 * ( int32() + 2.32830643653869629E-10 * int32() );
}

注意：doub() 方法返回的浮点数仅有 32 位精度。若需完整精度，请使用较慢的 truedoub() 方法。

参考文献

Gentle, J.E. 2003, Random Number Generation and Monte Carlo Methods, 2nd ed. (New York: Springer), Chapter 1.

Marsaglia, G 1968, "Random Numbers Fall Mainly in the Planes", Proceedings of the National Academy of Sciences, vol. 61, pp. 25–28.[1]

Entacher, K. 1997, "A Collection of Selected Pseudorandom Number Generators with Linear Structures", Technical Report No. 97, Austrian Center for Parallel Computation, University of Vienna. [Available on the Web at multiple sites.][2]

Knuth, D.E. 1997, Seminumerical Algorithms, 3rd ed., vol. 2 of The Art of Computer Programming (Reading, MA: Addison-Wesley).[3]

Schrage, L. 1979, "A More Portable Fortran Random Number Generator," ACM Transactions on Mathematical Software, vol. 5, pp. 132-138.[4]

Park, S.K., and Miller, K.W. 1988, "Random Number Generators: Good Ones Are Hard to Find," Communications of the ACM, vol. 31, pp. 1192-1201.[5]

L'Ecuyer, P. 1997 "Uniform Random Number Generators: A Review," Proceedings of the 1997 Winter Simulation Conference, Andradottir, S. et al., eds. (Piscataway, NJ: IEEE).[6]

Marsaglia, G. 1999, "Random Numbers for C: End, at Last?", posted 1999 January 20 to sci.stat.math.[7]

Marsaglia, G. 2003, "Diehard Battery of Tests of Randomness v0.2 beta," 2007+ at http://www.cs.hku.hk/~diehard./[8]

Rukhin, A. et al. 2001, "A Statistical Test Suite for Random and Pseudorandom Number Generators", NIST Special Publication 800-22 (revised to May 15, 2001).[9]

Marsaglia, G. 2003, "Xorshift RNGs", Journal of Statistical Software, vol. 8, no. 14, pp. 1-6.[10]

Brent, R.P. 2004, "Note on Marsaglia's Xorshift Random Number Generators", Journal of Statistical Software, vol. 11, no. 5, pp. 1-5.[11]

L'Ecuyer, P. 1996, "Maximally Equidistributed Combined Tausworththe Generators," Mathematics of Computation, vol. 65, pp. 203-213.[12]

L'Ecuyer, P. 1999, "Tables of Maximally Equidistributed Combined LSFR Generators," Mathematics of Computation, vol. 68, pp. 261-269.[13]

Couture, R. and L'Ecuyer, P. 1997, "Distribution Properties of Multiply-with-Carry Random Number Generators," Mathematics of Computation, vol. 66, pp. 591-607.[14]

L'Ecuyer, P. 1999, "Tables of Linear Congruential Generators of Different Sizes and Good Lattice Structure", Mathematics of Computation, vol. 68, pp. 249-260.[15]

7.2 对大型数组进行完全哈希

我们在 §7.1.4 中介绍了随机哈希或哈希函数的概念。有时我们可能需要一个不仅作用于单个字，而是作用于长度为 $M$ 的整个数组的哈希函数。作为完美主义者，我们希望哈希输出数组中的每一位都依赖于输入数组中的每一位。实现这一目标的一种方法是借鉴看似毫不相关的算法（如数据加密标准 DES 和快速傅里叶变换 FFT）中的结构思想！

DES 与其前身密码系统 LUCIFER 一样，是一种所谓的“分组乘积密码”[1]。它对 64 位输入进行操作，通过迭代（实际上是 16 次）应用一种高度非线性的位混合函数。图 7.2.1 展示了 DES 在此混合过程中信息的流动。将 32 位输入映射为 32 位输出的函数 $g$ 被称为“密码函数”。Meyer 和 Matyas [1] 讨论了密码函数必须是非线性的以及其他设计准则的重要性。

DES 通过一系列复杂的位置换和短连续位序列的查表操作来构造其密码函数 $g$。就我们的目的而言，更倾向于使用一种能在高级计算机语言中快速计算的不同函数 $g$。这样的函数可能会在密码学意义上削弱算法。然而，我们的目的并非密码学：我们希望找到最快的 $g$，以及图 7.2.1 中混合过程所需的最少迭代次数，使得输出的随机序列能够通过通常用于随机数生成器的测试。由此得到的算法不是 DES，而是一种更适合当前用途的“伪 DES”。

图 7.2.1. 数据加密标准（DES）以图示方式（引自 Meyer 和 Matyas [1]）对两个 32 位字迭代应用非线性函数 $g$。

根据上述准则，即 $g$ 应为非线性函数，我们必须在 $g$ 中赋予整数乘法操作重要地位。我们将自己限制在将 16 位操作数相乘得到 32 位结果的情形。$g$ 的基本思想是：计算输入高 16 位和低 16 位半字的三个不同的 32 位乘积，然后通过快速操作（例如加法或异或）将这些乘积（可能还包括一些固定常数）组合成一个 32 位结果。

实现这一通用方案的方式数量有限，因此可以系统地探索各种替代方案。通过对输出随机性的实验和测试，我们得到了图 7.2.2 所示的操作序列。图中一些新元素需要解释：值 $C_1$ 和 $C_2$ 是固定常数，随机选择，但要求它们恰好包含 16 个 1 位和 16 个 0 位；通过异或组合这些常数可确保整体函数 $g$ 对 0 位或 1 位没有偏向性。图 7.2.2 中的“反转半字”操作被证明是必不可少的；否则，最低位和最高位无法通过三次乘法得到充分混合。

接下来需要确定我们所能接受的最小迭代次数 $N_{it}$。在本节中，我们推荐 $N_{it}=2$。我们尚未发现由此方案生成的多达 $10^9$ 个随机数序列存在任何统计上的非随机性偏差。不过，我们仍为 $N_{it}\le 4$ 的情形提供了 $C_1$ 和 $C_2$ 常数。

void psdes(Uint &lword, Uint &rword) {
    // 对 64 位字 (lword, rword) 进行伪 DES 哈希。
    // 两个 32 位参数的所有位均被哈希后返回。
    const int NITER = 2;
    static const Uint c1[4] = {
        0xbaa96887L, 0x1e17d32cL, 0x03bcdc3cL, 0x0f33d1b2L
    };
    static const Uint c2[4] = {
        0x4b0f3b58L, 0xe874f0c3L, 0x6955c5a6L, 0x55a7ca46L
    };
    Uint i, ia, ib, iswap, itmp1, itmph;
    for (i = 0; i < NITER; i++) {
        // 执行 NITER 次 DES 逻辑迭代，使用一个更简单（非密码学）的非线性函数替代 DES 的函数。
        ia = (iswap = rword) ^ c1[i];
        itmp1 = ia & 0xffff;
        itmph = ia >> 16;
        // 位丰富的常数 c1 和（下方的）c2 保证了充分的非线性混合。
        ib = itmp1 * itmp1 + ~(itmph * itmph);
        rword = lword ^ ((ia = (ib >> 16) | ((ib & 0xffff) << 16)) ^ c2[i]) + itmp1 * itmph;
        lword = iswap;
    }
}

到目前为止，这似乎与对大型数组进行完全哈希关系不大。然而，psdes 为我们提供了一个基本构件：一个用于对两个任意 32 位整数进行相互哈希的例程。现在我们转向 FFT 中的“蝶形”概念，将哈希扩展到整个数组。

蝶形是一种特定的算法结构，适用于长度为 $N$（2 的幂）的数组。它通过大约 $N{\log }_{2}N$ 次操作使每个元素与其他所有元素相互“通信”。一个有用的比喻是：假设数组中的一个元素患有一种疾病，它会感染所有与之接触的其他元素。那么蝶形在此具有两个有趣的性质：(i) 经过 ${\log }_{2}N$ 个阶段后，所有人都被感染；(ii) 经过 $j$ 个阶段后，有 $2^j$ 个元素被感染；通信路径从未出现“瓶颈”或“缩颈”现象。

图 7.2.2. 例程 psdes 所使用的非线性函数 $g$

蝶形的描述非常简单：在第一阶段，数组前半部分的每个元素与其在后半部分对应的元素相互通信。然后对每个半部分递归地执行相同的操作，依此类推。通过归纳法可以证明，现在每个元素都与其他所有元素建立了通信路径：显然当 $N=2$ 时成立；若对 $N$ 成立，则对 $2N$ 也成立，因为第一步使每个元素都能与自身所在半区及另一半区通信，之后根据假设，它就能到达任何位置。

我们需要对蝶形稍作修改，使得数组大小 $M$ 不必是 2 的幂。令 $N$ 为大于等于 $M$ 的最小 2 的幂。我们在（虚拟的）大小 $N$ 上执行蝶形操作，忽略与大于 $M$ 的不存在元素之间的通信。仅这样做还不够，因为前 $N/2$ 个元素中的靠后部分无法“感染”后 $N/2$ 个元素（在后续递归层级中也存在类似问题）。然而，如果在最后额外进行一次前 $N/2$ 与后 $N/2$ 元素之间的通信，那么所有缺失的通信路径都可以通过前 $N/2$ 个元素恢复。

以下代码的第三行是一个惯用表达式，将 n 设置为大于等于 m 的最小 2 的幂，这是 S.E. Anderson [2] 的一个微型杰作。如果你仔细观察，会发现它本身也是一种蝶形操作，只不过是在比特位上进行的！

void hashall(Vec<Uint> &arr) {
    // 将数组 arr 替换为一个同样大小的哈希数组，其所有位都依赖于 arr 中的所有位。
    // 使用 psdes 对两个 32 位字进行相互哈希。
    Int m = arr.size(), n = m - 1;
    n |= n >> 1; n |= n >> 2; n |= n >> 4; n |= n >> 8; n |= n >> 16; n++;
    // 令人难以置信的是，n 现在是大于等于 m 的最小 2 的幂。
    Int nb = n, nb2 = n >> 1, j, jb;
    if (n < 2) throw("size must be > 1");
    while (nb > 1) {
        for (jb = 0; jb < n - nb + 1; jb += nb)
            for (j = 0; j < nb2; j++)
                if (jb + j + nb2 < m) psdes(arr[jb + j], arr[jb + j + nb2]);
        nb = nb2;
        nb2 >>= 1;
    }
    nb2 = n >> 1;
    if (m != n) for (j = nb2; j < m; j++) psdes(arr[j], arr[j - nb2]);
    // 仅当 m 不是 2 的幂时才需要最后的混合。
}

参考文献

Meyer, C.H. 与 Matyas, S.M. 1982，《密码学：计算机数据安全的新维度》（纽约：Wiley）。[1]

Zonst, A.E. 2000，《理解FFT》，第2版修订版（佛罗里达州泰特斯维尔：Citrus Press）。

Anderson, S.E. 2005，“位操作技巧（Bit Twiddling Hacks）”，2007年及以后，见 http://graphics.stanford.edu/~seander/bit hacks.html。[2]

《数据加密标准》，1977年1月15日，联邦信息处理标准出版物第46号（华盛顿：美国商务部国家标准局）。

《NBS数据加密标准的实施与使用指南》，1981年4月1日，联邦信息处理标准出版物第74号（华盛顿：美国商务部国家标准局）。

7.3 从其他分布生成随机数

在§7.1中，我们学习了如何生成在0到1之间服从均匀分布的随机数，记作 $U(0,1)$。生成一个介于 $x$ 与 $x+dx$ 之间的数的概率为

$$p(x)dx=\{\begin{array}{ll}dx& 0\le x<1\\ 0& \text{其他情况}\end{array}$$

我们写作

$$x\sim \mathrm{U}(0,1)$$

如§6.14所述，符号 $\sim$ 可读作“从分布中抽取”。

在本节中，我们将学习如何生成从其他概率分布中抽取的随机数，包括§6.14中讨论的所有分布。对具体分布的讨论将与所用通用方法的讨论穿插进行。

7.3.1 指数分布随机数

假设我们生成一个均匀分布的随机数 $x$，然后对其应用某个给定的函数 $y(x)$。$y$ 的概率分布记作 $p(y)dy$，由概率的基本变换法则决定，即

$$|p(y)dy|=|p(x)dx|$$

或

$$p(y)=p(x)\left|\frac{dx}{dy}\right|$$

举个例子，令

$$\gamma (x)=-\ln (x)$$

其中 $x\sim \mathrm{U}(0,1)$。那么

$$p(y)dy=\left|\frac{dx}{dy}\right|dy=e^{-y}dy$$

这就是单位均值的指数分布 $\text{Exponential}(1)$，如§6.14.5中所述。这种分布在现实生活中经常出现，通常用于描述独立泊松随机事件之间的等待时间，例如原子核的放射性衰变。你也可以很容易地看出（由式7.3.6），量 $y/\beta$ 的概率分布为 $\beta e^{-\beta y}$，因此

$$y/\beta \sim \mathrm{Exponential}\left(\beta \right)$$

因此，我们每次调用只需消耗大约一个均匀分布随机数和一次对数运算，即可生成一个指数分布的随机数。

deviates.h

struct Expondev : Ran {
    // 用于生成指数分布随机数的结构体。
    Doub beta;
    Expondev(Doub bbbeta, Ullong i) : Ran(i), beta(bbbeta) {}
    // 构造函数参数为 β 和随机数序列的种子。
    Doub dev() {
        // 返回一个指数分布的随机数。
        Doub u;
        do u = doub(); while (u == 0.);
        return -log(u)/beta;
    }
};

本节及后续小节中，我们约定：每种分布对应的类都从均匀随机数生成器类 Ran 派生而来。我们使用构造函数设置分布的参数并初始化生成器的种子，然后提供一个 dev() 方法，用于返回该分布的一个随机数。

7.3.2 一般情况下的变换法

让我们看看如何使用上述变换法来生成任意指定的 $y$ 的分布，例如某个满足 $p(y)=f(y)$ 的分布，其中 $f$ 是一个积分为1的正函数。根据式（7.3.4），我们需要求解如下微分方程：

图7.3.1：从已知概率分布 $p(y)$ 生成随机数 $y$ 的变换法。必须已知 $p(y)$ 的不定积分且该积分可逆。在 $0$ 到 $1$ 之间选取一个均匀分布的随机数 $x$，其在定积分曲线上对应的 $y$ 值即为所需的随机数。

$$\frac{dx}{dy}=f(y)$$

但该方程的解就是 $x=F(y)$，其中 $F(y)$ 是 $f(y)$ 的不定积分。因此，将均匀分布的随机数转换为服从分布 $f(y)$ 的随机数所需的变换为

$$\gamma (x)=F^{-1}(x)$$

其中 $F^{-1}$ 是 $F$ 的反函数。公式 (7.3.9) 是否可行，取决于 $f(y)$ 的积分的反函数是否能够解析或数值计算。有时可以，有时则不行。

顺便提一下，(7.3.9) 有一个直接的几何解释：由于 $F(y)$ 表示概率曲线在 $y$ 左侧的面积，(7.3.9) 的做法就是：选取一个均匀随机数 $x$，然后找到使得其左侧概率面积恰好为 $x$ 的值 $y$，并返回该 $y$ 值。（见图 7.3.1。）

7.3.3 Logistic 分布随机数

如 §6.14.4 中所述，Logistic 分布的随机数可通过变换法利用公式 (6.14.15) 轻松生成。每次生成一个 Logistic 随机数的代价主要是一次均匀随机数生成和一次对数运算。

struct Logisticdev : Ran {
    // 用于生成 Logistic 分布随机数的结构体。
    Doub mu, sig;
    Logisticdev(Doub mmu, Doub ssig, Ullong i) : Ran(i), mu(mmu), sig(ssig) {}
    // 构造函数参数为 $\mu$、$\sigma$ 和随机数序列种子。
    Doub dev() {
        // 返回一个 Logistic 分布的随机数。
        Doub u;
        do u = doub(); while (u*(1.-u) == 0);
        return mu + 0.551328895421792050*sig*log(u/(1.-u));
    }
};

7.3.4 通过变换法生成正态分布随机数（Box-Muller 方法）

变换法可以推广到多维情形。若 $x_1,x_2,\dots$ 是具有联合概率分布 $p(x_1,x_2,\dots )d{x}_{1}d{x}_2\dots$ 的随机变量，且 $y_1,y_2,\dots$ 均为所有 $x$ 的函数（$y$ 的个数与 $x$ 相同），则 $y$ 的联合概率分布为

$$p(y_1,y_2,\dots )d{y}_{1}d{y}_2\dots =p(x_1,x_2,\dots )\left|\frac{\partial (x_1,x_2,\dots )}{\partial (y_1,y_2,\dots )}\right|d{y}_{1}d{y}_2\dots$$

其中 $\left|\partial (\mathbf{x})/\partial (\mathbf{y})\right|$ 是 $x$ 关于 $y$ 的雅可比行列式（或 $y$ 关于 $x$ 的雅可比行列式的倒数）。

历史上一个重要的例子是利用 (7.3.10) 的 Box-Muller 方法，用于生成服从正态（高斯）分布的随机数（见 §6.14.1）：

$$p(y)dy=\frac{1}{\sqrt{2\pi }}{e}^{-y^2/2}dy$$

考虑在区间 (0,1) 上的两个均匀随机数 $x_1$、$x_2$ 与两个量 $y_1$、$y_2$ 之间的变换：

$$\begin{array}{r}{y}_1=\sqrt{-2\ln x_1}\cos 2\pi x_2\\ y_2=\sqrt{-2\ln x_1}\sin 2\pi x_2\end{array}$$

等价地，也可以写成：

$$\begin{array}{r}{x}_1=\exp \left[-\frac{1}{2}(y_1^2+y_2^2)\right]\\ x_2=\frac{1}{2\pi }\arctan \frac{y_2}{y_1}\end{array}$$

现在可以轻松计算雅可比行列式（试试看！）：

$$\frac{\partial (x_1,x_2)}{\partial (y_1,y_2)}=\left|\begin{array}{ll}\frac{\partial x_1}{\partial y_1}& \frac{\partial x_1}{\partial y_2}\\ \frac{\partial x_2}{\partial y_1}& \frac{\partial x_2}{\partial y_2}\end{array}\right|=-\left[\frac{1}{\sqrt{2\pi }}{e}^{-y_1^2/2}\right]\left[\frac{1}{\sqrt{2\pi }}{e}^{-y_2^2/2}\right]$$

由于该结果是仅关于 $y_1$ 的函数与仅关于 $y_2$ 的函数的乘积，因此每个 $y$ 均独立地服从正态分布 (7.3.11)。

在应用 (7.3.12) 时还有一个有用的技巧：假设我们不在单位正方形内选取均匀随机数 $x_1$ 和 $x_2$，而是在以原点为中心的单位圆内随机选取一点，其纵坐标和横坐标分别为 $v_1$ 和 $v_2$。那么它们的平方和 $R^2\equiv v_1^2+v_2^2$ 就是一个均匀分布的随机数，可用作 $x_1$；而向量 $(v_1,v_2)$ 与 $v_1$ 轴之间的夹角可作为随机角度 $2\pi x_2$。这样做的优势在于，(7.3.12) 中的余弦和正弦函数现在可以写成 $v_1/\sqrt{R^2}$ 和 $v_2/\sqrt{R^2}$，从而避免了三角函数调用！（下一节我们将大幅推广这一技巧。）

以下是使用 Box-Muller 方法生成正态分布随机数的代码。请注意，此代码仅用于教学目的，因为 §7.3.9 将介绍一种显著更快的正态分布随机数生成方法。

struct Normaldev_BM : Ran {
    // 用于生成正态分布随机数的结构体。
    Doub mu, sig;
    Doub storedval;
    Normaldev_BM(Doub mmu, Doub ssig, Ullong i)
    : Ran(i), mu(mmu), sig(ssig), storedval(0.) {}
    // 构造函数参数为 $\mu$、$\sigma$ 和随机数序列种子。
    Doub dev() {
        // 返回一个正态分布的随机数。
        Doub v1, v2, rsq, fac;
        if (storedval == 0.) {
            do {
                v1 = 2.0 * doub() - 1.0;
                v2 = 2.0 * doub() - 1.0;
                rsq = v1 * v1 + v2 * v2;
            } while (rsq >= 1.0 || rsq == 0.0);
            fac = sqrt(-2.0 * log(rsq) / rsq);
            storedval = v1 * fac;
            return mu + sig * v2 * fac;
        } else {
            fac = storedval;
            storedval = 0;
            return mu + sig * fac;
        }
    }
};

7.3.5 瑞利（Rayleigh）分布随机数

瑞利分布对正数 $z$ 定义为

$$p(z)dz=z\exp \left(-\frac{1}{2}{z}^2\right)dz(z>0)$$

由于其不定积分可以解析求出，且结果易于求反函数，因此可直接通过均匀分布随机数 $x$ 进行简单变换得到：

$$z=\sqrt{-2\ln x}.x\sim \mathrm{U}(0,1)$$

瑞利分布随机数 $z$ 也可以通过两个标准正态分布随机数 $y_1$ 和 $y_2$ 生成：

$$z=\sqrt{y_1^2+y_2^2},y_1,y_2\sim \mathrm{N}(0,1)$$

事实上，如果我们对 Box-Muller 方法（公式 (7.3.12)）中生成的两个正态偏差 $y_1$ 和 $y_2$ 进行平方并求和，公式 (7.3.16) 与 (7.3.17) 之间的关系就显而易见了。

7.3.6 拒绝法

拒绝法是一种强大而通用的技术，用于生成服从已知且可计算的概率密度函数 $p(x)dx$（即变量取值在 $x$ 到 $x+dx$ 之间的概率）的随机偏差。拒绝法不要求累积分布函数（即 $p(x)$ 的不定积分）易于计算，更不要求该函数的反函数——而上一节介绍的变换法恰恰需要这一条件。

拒绝法基于一个简单的几何论证（见图 7.3.2）：

图 7.3.2. 从已知概率分布 $p(x)$（处处小于另一函数 $f(x)$）中生成随机偏差 $x$ 的拒绝法。首先使用变换法生成服从分布 $f$ 的随机偏差 $x$（参见图 7.3.1）。再利用一个均匀分布的随机数决定是否接受该 $x$。若被拒绝，则重新生成一个新的 $f$ 的偏差，依此类推。被接受点与被拒绝点的比例等于 $p$ 曲线下面积与 $p$ 和 $f$ 之间区域面积之比。

绘制你希望生成的概率分布 $p(x)$ 的图像，使得曲线下任意 $x$ 区间内的面积对应于在该区间内生成 $x$ 的期望概率。如果我们有一种方法可以在该曲线下二维区域内均匀地选择一个随机点，那么该随机点的 $x$ 坐标就具有所需的分布。

现在，在同一张图上绘制另一条曲线 $f(x)$，要求其总面积有限（非无穷），并且处处位于原始概率分布 $p(x)$ 之上。（这总是可行的，因为根据概率的定义，原始曲线下的总面积为 1。）我们将此 $f(x)$ 称为比较函数。假设我们有一种方法可以在比较函数曲线下二维区域内均匀地选择一个随机点。每当该点落在原始概率分布曲线下区域之外时，我们就拒绝它并重新选择一个随机点；每当它落在原始概率分布曲线下区域之内时，我们就接受它。

显然，被接受的点在被接受区域内是均匀分布的，因此它们的 $x$ 值具有所需的分布。同样显然的是，被拒绝点的比例仅取决于比较函数面积与概率分布函数面积之比，而与任一函数的具体形状无关。例如，即使某个比较函数在某些 $x$ 值处对概率函数的近似效果很差（例如在 $p(x)$ 为零的区域仍保持有限值），只要其总面积小于 2，被拒绝的点就少于一半。

接下来只需说明如何在比较函数 $f(x)$ 下二维区域内均匀地选择一个随机点。变换法（§7.3）的一个变体可以很好地完成这一任务：确保所选的比较函数的不定积分在解析上是已知的，并且其反函数也能解析地给出 $x$ 作为“$x$ 左侧比较函数曲线下面积”的函数。现在，在 0 到 $A$ 之间生成一个均匀随机数（其中 $A$ 是 $f(x)$ 下的总面积），并用它得到对应的 $x$ 值。然后，在 0 到 $f(x)$ 之间生成另一个均匀随机数作为该二维点的 $y$ 坐标。最后，根据该点是否分别小于或大于 $p(x)$ 来决定接受或拒绝。

综上所述，对于给定的 $p(x)$，拒绝法要求我们一次性找到一个相当不错的比较函数 $f(x)$。此后，

每次生成一个偏差需要两个均匀随机数，一次 $f$ 的求值（用于得到坐标 $y$）和一次 $p$ 的求值（用于决定是否接受点 $(x,y)$）。图 7.3.1 展示了整个过程。当然，平均而言，该过程可能需要重复 $A$ 次才能最终获得一个偏差。

7.3.7 柯西偏差

在 Box-Muller 方法中，公式 (7.3.14) 后描述的“进一步技巧”现在被看作是一种用于获取均匀随机角度的三角函数的拒绝法。如果我们将此方法与柯西分布累积分布函数（cdf）反函数的显式公式（公式 (6.14.6)，参见 §6.14.2）结合起来，就可以非常高效地生成柯西偏差。

struct Cauchydev : Ran {
    // 用于生成柯西偏差的结构体。
    Doub mu, sig;
    Cauchydev(Doub mmu, Doub ssig, Ullong i) : Ran(i), mu(mmu), sig(ssig) {}
    // 构造函数参数为 $\mu$、$\sigma$ 和随机序列种子。
    Doub dev() {
        // 返回一个柯西偏差。
        Doub v1, v2;
        do {
            v1 = 2.0 * doub() - 1.0;
            v2 = doub();
        } while (SQR(v1) + SQR(v2) >= 1. || v2 == 0.);
        return mu + sig * v1 / v2;
    }
};

7.3.8 均匀比方法

在生成柯西偏差时，我们取了两个均匀偏差的比值，这两个偏差被限制在单位圆内。如果我们推广到单位圆以外的其他形状，并结合拒绝法的原理，就会得到一种强大的变体。Kinderman 和 Monahan [1] 指出，几乎任何概率分布 $p(x)$ 的偏差都可以通过以下令人惊讶的方法生成：

• 在 $(u,v)$ 平面上构造由 $0\le u\le [p(v/u){]}^{1/2}$ 所界定的区域。
• 在该区域内均匀地选择两个偏差 $u$ 和 $v$。
• 返回 $v/u$ 作为偏差。

证明：我们可以用 $(x,p)$ 平面上的方程表示普通的拒绝法：

$$p(x)dx={\int }_{p^{\prime }=0}^{p^{\prime }=p(x)}d{p}^{\prime }dx$$

由于被积函数为 1，只要 $p^{\prime }$ 在积分限内（即 $0<p^{\prime }<p(x)$），我们就有理由在 $(x,p^{\prime })$ 上均匀采样。现在进行变量替换：

$$\begin{array}{r}\frac{v}{u}=x\\ u^2=p\end{array}$$

图 7.3.3. 均匀比方法。该泪滴形状的内部是正态分布的接受区域：如果在该区域内随机选择一个点，则比值 $v/u$ 将是一个正态偏差。

于是公式 (7.3.18) 变为：

$$p(x)dx={\int }_{p^{\prime }=0}^{p^{\prime }=p(x)}d{p}^{\prime }dx={\int }_{u=0}^{u=\sqrt{p(x)}}\frac{\partial (p,x)}{\partial (u,v)}dudv=2{\int }_{u=0}^{u=\sqrt{p(v/u)}}dudv$$

因为（如你所计算的）雅可比行列式为常数 2。由于新的被积函数为常数，在 $(u,v)$ 上按所示 $u$ 的限值进行均匀采样等价于在 $(x,p)$ 上使用拒绝法。

上述对 $u$ 的限制通常定义了一个“泪滴”形状的区域。要理解这一点，请注意 $x=v/u$ 为常数的轨迹是径向直线。沿每条径向线，接受区域从原点延伸到 $u^2=p(x)$ 的点。由于大多数概率分布在 $x$ 很大或很小时都趋于零，因此接受区域沿径向线向原点收缩，形成泪滴形状。当然，重要的是这个泪滴的确切形状。图 7.3.3 显示了正态分布情况下接受区域的形状。

通常，当目标区域可以被矩形、平行四边形或其他易于均匀采样的形状紧密包围时，就会使用这种均匀比方法。然后，我们通过拒绝目标区域外的点，从对简单形状的采样过渡到对目标区域的采样。

均匀比方法的一个重要补充是“挤压”（squeeze）的概念。挤压是指任何易于计算的形状，它从内部或外部紧密地包围拒绝法的接受区域。最理想的情况是同时拥有内外两个挤压。这样，你可以立即拒绝位于外挤压之外的点，并立即接受位于内挤压之内的点。只有当你不幸抽到位于两个挤压之间的点时，才需要进行更耗时的实际拒绝边界比较计算。挤压在普通拒绝法和均匀比方法中都很有用。

7.3.9 使用均匀比方法生成正态偏差

Leva [2] 提出了一种使用均匀比（ratio-of-uniforms）方法生成正态偏差的算法，并取得了极大成功。他使用二次曲线在 $(u,v)$ 平面上同时提供内侧和外侧的“挤压”（squeeze），紧密贴合目标区域（见图 7.3.3）。仅有约 1% 的情况下需要计算精确边界（此时需计算对数）。

所得代码看起来如此简单且“非超越性”，以至于很难相信它能生成精确的正态偏差。但事实确实如此！

struct Normaldev : Ran {
    // 用于生成正态偏差的结构体。
    Doub mu, sig;
    Normaldev(Doub mmu, Doub ssig, Ullong i)
        : Ran(i), mu(mmu), sig(ssig) {}
    // 构造函数参数为 μ、σ 和随机序列种子。
    Doub dev() {
        // 返回一个正态偏差。
        Doub u, v, x, y, q;
        do {
            u = doub();
            v = 1.7156 * (doub() - 0.5);
            x = u - 0.449871;
            y = abs(v) + 0.386595;
            q = SQR(x) + y * (0.19600 * y - 0.25472 * x);
        } while (q > 0.27597
                 && (q > 0.27846 || SQR(v) > -4. * log(u) * SQR(u)));
        return mu + sig * v / u;
    }
};

注意，while 子句利用了 C（及 C++）语言对逻辑表达式求值的短路特性：若第一个操作数足以确定结果，则不会对第二个操作数求值。根据这一规则，仅当 q 介于 0.27597 与 0.27846 之间时，才会计算对数。

平均而言，每个正态偏差需使用 2.74 个均匀偏差。顺便指出，尽管各种常数仅给出六位有效数字，但该方法是精确的（达到双精度浮点数的全部精度）。边界曲线的微小扰动不会产生任何影响，因为精度隐含在（极少发生的）精确边界计算中。

7.3.10 伽马偏差

伽马分布 $\mathrm{Gamma}(\alpha ,\beta )$ 在 §6.14.9 中已有描述。参数 $\beta$ 仅作为尺度因子出现：

$$\mathrm{Gamma}(\alpha ,\beta )\cong \frac{1}{\beta }\mathrm{Gamma}(\alpha ,1)$$

（说明：要生成一个 $\mathrm{Gamma}(\alpha ,\beta )$ 偏差，只需先生成一个 $\mathrm{Gamma}(\alpha ,1)$ 偏差，再除以 $\beta$ 即可。）

若 $\alpha$ 是一个小的正整数，则生成 $x\sim \mathrm{Gamma}(\alpha ,1)$ 的一种快速方法是利用如下事实：它等价于单位均值泊松随机过程中第 $\alpha$ 次事件的等待时间。由于连续两次事件之间的间隔服从指数分布 $\mathrm{Exponential}(1)$，因此只需将 $\alpha$ 个指数分布的等待时间相加，即对 $\alpha$ 个均匀偏差取对数后求和。更进一步，由于对数之和等于乘积的对数，实际上只需计算 $\alpha$ 个均匀偏差的乘积，再取一次对数即可。然而，由于这是一个非常特殊的情形，我们并未将其包含在下面的代码中。

当 $\alpha <1$ 时，伽马分布的概率密度函数无界，这带来不便。然而，文献 [4] 指出，若

$$\gamma \sim \mathrm{Gamma}(\alpha +1,1),u\sim \mathrm{Uniform}(0,1)$$

则

$$\gamma u^{1/\alpha }\sim \mathrm{Gamma}(\alpha ,1)$$

我们将在下面的代码中使用这一性质。

对于 $\alpha >1$ 的情形，Marsaglia 和 Tsang [5] 提出了一种优雅的拒绝采样方法，该方法基于对伽马分布的简单变换并结合“挤压”技术。变换后，伽马分布可被一个高斯曲线所界定，且该高斯曲线的面积至多比伽马曲线大 5%。因此，生成一个伽马偏差的代价仅略高于用于采样比较函数的正态偏差。以下代码给出了精确的实现；完整解释请参阅原始论文。

// deviates.h
struct Gammadev : Normaldev {
    // 用于生成伽马偏差的结构体。
    Doub alpha, alph, bet;
    Doub a1, a2;
    Gammadev(Doub aalph, Doub bbet, Ullong i)
        : Normaldev(0., 1., i), alpha(aalph), alph(aalph), bet(bbet) {
        // 构造函数参数为 α、β 和随机序列种子。
        if (alpha <= 0.) throw("bad alpha in Gammadev");
        if (alpha < 1.) alpha += 1.;
        a1 = alpha - 1. / 3.;
        a2 = 1. / sqrt(9. * a1);
    }
    Doub dev() {
        // 使用 Marsaglia 和 Tsang 的方法返回一个伽马偏差。
        Doub u, v, x;
        do {
            do {
                x = Normaldev::dev();
                v = 1. + a2 * x;
            } while (v <= 0);
            v = v * v * v;
            u = doub();
        } while (u > 1. - 0.331 * SQR(SQR(x)) && 
                 log(u) > 0.5 * SQR(x) + a1 * (1. - v + log(v))); // 很少执行
        if (alpha == alph) return a1 * v / bet;
        else {
            do u = doub(); while (u == 0);
            return pow(u, 1. / alph) * a1 * v / bet;
        }
    }
};

伽马偏差满足一个求和法则。若有一组相互独立的偏差 $y_i$，其 ${\alpha }_i$ 可能不同，但共享相同的 $\beta$ 值，

$${\gamma }_i\sim \mathrm{Gamma}({\alpha }_i,\beta )$$

则它们的和也是一个伽马偏差：

$$y\equiv \sum _i{y}_i\sim \mathrm{Gamma}({\alpha }_T,\beta ),{\alpha }_T=\sum _i{\alpha }_i$$

若各 ${\alpha }_i$ 为整数，则此性质与前述泊松等待时间的讨论密切相关。

7.3.11 可由其他偏差轻松生成的分布

利用正态偏差、伽马偏差和均匀偏差，我们可以免费获得许多其他分布。重要提示：当你需要组合多个结果时，请确保所有不同的 Normaldev、Gammadev 和 Ran 实例均使用不同的随机种子！（Ran 及其派生类足够健壮，使用种子 $i,i+1,\dots$ 是可以接受的。）

卡方偏差（参见 §6.14.8）

这很简单：

$$\mathrm{Chisquare}(\nu )\cong \mathrm{Gamma}\left(\frac{\nu }{2},\frac{1}{2}\right)\cong 2\mathrm{Gamma}\left(\frac{\nu }{2},1\right)$$

学生 t 偏差（参见 §6.14.3）

偏离 Student-t 分布的随机数可以通过一种与 Box-Muller 方法非常相似的方法生成。其对应于公式 (7.3.12) 的表达式为：

$$y=\sqrt{\nu (u_1^{-2/\nu }-1)\cos 2\pi u_2}$$

若 $u_1$ 和 $u_2$ 是相互独立的均匀分布随机变量，即 $U(0,1)$，则有：

$$\gamma \sim \mathrm{Student}(\nu ,0,1)$$

或

$$\mu +\sigma \gamma \sim \mathrm{Student}(\nu ,\mu ,\sigma )$$

不幸的是，你无法像 Box-Muller 方法那样一次生成两个偏离值，因为类似于公式 (7.3.14) 的雅可比行列式无法分解。不过，你仍可能希望使用极坐标方法，仅用于计算 $\cos 2\pi u_2$，但此时其优势已不再那么显著。

另一种方法利用正态分布与伽马分布偏离值的比值。如果我们有：

$$x\sim \mathrm{N}(0,1),y\sim \mathrm{Gamma}(\frac{\nu }{2},\frac{1}{2})$$

那么：

$$x\sqrt{\nu /y}\sim \mathrm{Student}(\nu ,0,1)$$

贝塔分布偏离值（参见 §6.14.11）

若：

$$x\sim \mathrm{Gamma}(\alpha ,1),y\sim \mathrm{Gamma}(\beta ,1)$$

则：

$$\frac{x}{x+y}\sim \mathrm{Beta}(\alpha ,\beta )$$

F 分布偏离值（参见 §6.14.10）

若：

$$x\sim \mathrm{Beta}(\frac{1}{2}{\nu }_1,\frac{1}{2}{\nu }_2)$$

（见公式 7.3.33），则：

$$\frac{{\nu }_{2}x}{{\nu }_1(1-x)}\sim \mathrm{F}({\nu }_1,{\nu }_2)$$

图 7.3.4. 拒绝法应用于整数值分布的情形。该方法作用于图中虚线所示的阶梯函数，产生一个实数值偏离量。随后将该偏离量向下取整为下一个较小的整数，并作为输出。

7.3.12 泊松分布偏离值

泊松分布 Poisson($\lambda$)，如前所述（§6.14.13），是一种离散分布，因此其偏离值为整数 $k$。为了使用前述方法，可通过以下技巧将泊松分布转换为连续分布：将有限概率 $p(k)$ 均匀地铺展到区间 $[k,k+1)$ 上。这定义了一个连续分布 $q_{\lambda }(k)dk$，其形式为：

$$q_{\lambda }(k)dk=\frac{{\lambda }^{\lfloor k\rfloor }{e}^{-\lambda }}{\lfloor k\rfloor !}dk$$

其中 $\lfloor k\rfloor$ 表示不大于 $k$ 的最大整数。若我们使用拒绝法或其他任何方法，从 (7.3.36) 生成一个（非整数）偏离值，再取其整数部分，则结果等价于从离散泊松分布中抽取的样本。（参见图 7.3.4。）该技巧适用于任何整数值概率分布。除了使用“向下取整”（floor）操作符外，也可以同样使用“向上取整”（ceiling）或“四舍五入”（nearest）——只要能将概率均匀铺展到单位区间即可。

当 $\lambda$ 足够大时，分布 (7.3.36) 在定性上呈钟形（尽管是由许多小的矩形台阶构成的钟形）。此时，均匀比方法（ratio-of-uniforms method）效果良好。在 $(u,v)$ 平面上，不难找到形如 $v^2=Q(u)$ 的简单内、外挤压函数，其中 $Q(u)$ 是关于 $u$ 的简单多项式。唯一的技巧在于，要确保挤压函数之间留有足够大的间隙，以包含所有 $\lambda$ 值下真实、锯齿状的边界。（可提前参考图 7.3.5 中的类似示例。）

对于中等大小的 $\lambda$，锯齿效应过于显著，使得挤压函数变得不切实际，但未经修饰的均匀比方法仍然效果不错。

对于较小的 $\lambda$，我们可以采用类似于前述伽马分布在整数形状参数 $a$ 情况下的思路：当独立指数分布偏离值的累加和首次超过 $\lambda$ 时，其累加项数减 1 即为一个泊松偏离值 $k$。此外，正如伽马分布所述，我们可以将 $U(0,1)$ 的均匀偏离值相乘，而不必累加指数分布 Exponential(1) 的偏离值。

这些思想产生了以下例程。

struct Poissondev : Ran {
    // 用于生成泊松分布随机数的结构体。
    Doub lambda, sqlam, loglam, lamexp, lambold;
    VecDoub logfact;
    Int swch;
    Poissondev(Doub lambda, Ullong i) : Ran(i), lambda(lambda),
        logfact(1024, -1.), lambold(-1.) {}
    // 构造函数参数为 $\lambda$ 和一个随机序列种子。
    Int dev() {
        // 使用最近设置的 $\lambda$ 值返回一个泊松分布随机数。
        Doub u, u2, v, v2, p, t, lfac;
        Int k;
        if (lambda < 5.) {
            if (lambda != lambold) {
                lambold = lambda;
                lamexp = exp(-lambda);
            }
            k = -1;
            t = 1;
            do {
                ++k;
                t *= doub();
            } while (t > lamexp);
        } else {
            if (lambda != lambold) {
                lambold = lambda;
                sqlam = sqrt(lambda);
                loglam = log(lambda);
            }
            for (;;) {
                u = 0.64 * doub();
                v = -0.68 + 1.28 * doub();
                if (lambda > 13.5) {
                    v2 = SQR(v);
                    if (v >= 0.) {
                        if (v2 > 6.5 * u * (0.64 - u) * (u + 0.2)) continue;
                    } else {
                        if (v2 > 9.6 * u * (0.66 - u) * (u + 0.07)) continue;
                    }
                }
                k = Int(floor(sqlam * (v / u) + lambda + 0.5));
                if (k < 0) continue;
                u2 = SQR(u);
                if (lambda > 13.5) {
                    if (v >= 0.) {
                        if (v2 < 15.2 * u2 * (0.61 - u) * (0.8 - u)) break;
                    } else {
                        if (v2 < 6.76 * u2 * (0.62 - u) * (1.4 - u)) break;
                    }
                }
                if (k < 1024) {
                    if (logfact[k] < 0.) logfact[k] = gammln(k + 1);
                    lfac = logfact[k];
                } else lfac = gammln(k + 1);
                p = sqlam * exp(-lambda + k * loglam - lfac);
                if (u2 < p) break;
            }
        }
        return k;
    }
    Int dev(Doub lambda) {
        // 重置 $\lambda$，然后返回一个泊松分布随机数。
        lambold = -1.; // 强制重新计算
        this->lambda = lambda;
        return dev();
    }
};

图 7.3.5. 应用于二项分布随机数生成的均匀比方法。点在 $(\mathbf{u},\mathbf{v})$ 平面上随机选择。平滑曲线是内层和外层挤压边界。锯齿状曲线对应于各种二项分布，其中 $\mathbf{n}>64$ 且 $\mathbf{n}\mathbf{p}>30$。只有当随机点落在平滑曲线之间时，才需要计算二项概率。

在 $\lambda >13.5$ 的情形下，上述代码每生成一个泊松分布随机数大约使用 3.3 个均匀分布随机数，并进行 0.4 次精确概率的计算（包括一次指数运算，以及对较大的 $k$ 调用 gammln）。

如果你使用无参的 dev 函数多次生成相同 $\lambda$ 值的随机数，Poissondev 会稍快一些；而如果你每次调用都改变 $\lambda$（使用单参数重载的 dev 函数，该函数正是为此目的提供的），则会稍慢。差别仅在于额外的一次指数运算（$\lambda <5$ 时）或平方根与对数运算（$\lambda \ge 5$ 时）。另外请注意对象中存储的已计算对数阶乘表。如果你的 $\lambda$ 值达到 $\sim 10^3$ 量级，你可能需要增大该表的大小。

7.3.13 二项分布随机数

生成二项分布随机数 $k\sim \text{Binomial}(n,p)$ 涉及许多与泊松分布随机数相同的思想。该分布同样是整数值的，因此我们使用相同的技巧将其转换为阶梯状连续分布。我们总可以将注意力限制在 $p\le 0.5$ 的情形，因为分布的对称性使我们可以轻松地从 $p>0.5$ 的情形恢复结果。

当 $n>64$ 且 $np>30$ 时，我们使用均匀比方法，其挤压边界如图 7.3.5 所示。每次生成一个二项分布随机数的成本约为 3.2 个均匀分布随机数，外加 0.4 次精确概率的计算。

对于 $n>64$ 且 $np<30$ 的情形，若花费大量精力去优化是愚蠢的，因为只需简单地制表累积分布函数（cdf），例如对 $0\le k<64$，然后循环遍历 $k$ 直到找到正确的值即可。（下面实现的二分查找甚至更好。）使用长度为 64 的 cdf 表时，表尾被忽略的概率永远不会超过 $\sim 10^{-20}$。（以每秒 $10^9$ 个随机数的速度，你需要运行 3000 年才会丢失一个随机数。）

剩下的就是有趣的情形 $n<64$，我们将详细探讨，因为它展示了位并行随机比较这一重要概念。

类似于小整数 $a$ 的伽马分布随机数和小 $\lambda$ 的泊松分布随机数的方法，二项分布随机数也有如下直接方法：生成 $n$ 个 $U(0,1)$ 均匀分布随机数，统计其中小于 $p$ 的个数，并将该计数作为 $k\sim \mathrm{Binomial}(n,p)$ 返回。这实际上就是二项过程的基本定义！

直接方法的问题在于，即使 $k$ 的期望值很小，它似乎也需要 $n$ 个均匀分布随机数。如果我们告诉你，对于 $n\le 64$，你平均最多只需七个 64 位均匀分布随机数就能达到相同目标，你会感到惊讶吗？方法如下。

将 $p<1$ 展开为其前 5 位二进制位加上一个余项： $$p=b_1{2}^{-1}+b_2{2}^{-2}+\cdots +b_5{2}^{-5}+p_r{2}^{-5}$$ 其中每个 $b_i$ 为 0 或 1，且 $0\le p_r\le 1$。

现在假设你已生成并存储了 64 个均匀分布 $U(0,1)$ 随机数，且 64 位字 $P$ 仅显示这 64 个数的首位。将 $P$ 的每一位与 $b_1$ 比较。如果位相同，则我们尚不知道该均匀随机数是小于还是大于 $p$；但如果位不同，则我们知道该随机数小于 $p$（当 $b_1=1$ 时）或大于 $p$（当 $b_1=0$ 时）。如果我们保留一个“已知”与“未知”情形的掩码，就可以通过位逻辑运算以位并行方式执行这些比较（参见下面的代码了解具体实现）。接着以同样方式处理第二位 $b_2$。在每个阶段，我们将剩余未知情形的一半变为已知。经过五个阶段（对于 $n=64$），平均剩下两个未知情形，我们通过生成一个新的均匀随机数并与 $p_r$ 比较来完成处理。（这需要遍历 64 位；但由于 C++ 没有位计数（"popcount"）操作，我们无论如何都必须进行这样的循环。如果你能进行位计数，你可能更愿意进行更多阶段，直到未知掩码为零。）

关键在于，这五个阶段使用的位实际上并非 64 个生成器的前五位，而只是五个独立的 64 位随机整数。选择数字五是因为它最小化了 $64\times 2^{-j}+j$，即所需随机数的期望数量。

因此，二项分布随机数的代码最终包含三种独立的方法：位并行直接法、cdf 查表法（通过二分查找）和挤压均匀比法。

struct Binomialdev : Ran {
    // 用于生成二项分布随机数的结构体。
    Doub pp, p, pb, expnp, np, glnp, plog, pclog, sq;
    Int n, swch;
    Ullong uz, uo, unfin, diff, rltp;
    Int pbits[5];
    Doub cdf[64];
    Doub logfact[1024];
    Binomialdev(Int nn, Doub ppp, Ullong i) : Ran(i), pp(ppp), n(nn) {
        // 构造函数参数为 $n$、$p$ 和一个随机序列种子。
        Int j;
        pb = p = (pp <= 0.5 ? pp : 1.0 - pp);
        if (n <= 64) {
            uz = 0;
            uo = 0xffffffffffffffffLL;
            rltp = 0;
            for (j = 0; j < 5; j++) pbits[j] = 1 & ((Int)(pb * (1 << (j + 1))));
            for (j = 0; j < 5; j++) pb -= pbits[j] * (1 << -(j + 1));
            swch = 0;
        } else if (n * p < 30.) {
            cdf[0] = exp(n * log(1 - p));
            for (j = 1; j < 64; j++) 
                cdf[j] = cdf[j - 1] + exp(gammln(n + 1) - gammln(j + 1.) - gammln(n - j + 1.) + j * log(p) + (n - j) * log(1. - p));
            swch = 1;
        } else {
            np = n * p;
            glnp = gammln(n + 1);
            plog = log(p);
            pclog = log(1. - p);
            sq = sqrt(np * (1. - p));
            if (n < 1024) 
                for (j = 0; j <= n; j++) logfact[j] = gammln(j + 1);
            swch = 2;
        }
    }
    Int dev() {
        // 返回一个二项分布随机数。
        Int j, k, kl = 0, km;
        Doub y, u, v, u2, v2, b;
        if (swch == 0) {
            unfin = uo;
            rltp = 0;
            for (j = 0; j < 5; j++) {
                diff = unfin & (pbits[j] ? uo : uz);
                if (pbits[j]) rltp |= diff;
                else rltp &= ~diff;
                unfin &= ~diff;
            }
            k = 0;
            for (j = 0; j < n; j++) {
                if (unfin & 1) {
                    if (doub() < pb) k++;
                } else {
                    if (rltp & 1) k++;
                }
                unfin >>= 1;
                rltp >>= 1;
            }
        } else if (swch == 1) {
            y = doub();
            k = 64;
            kl = 0;
            while (k - kl > 1) {
                km = (kl + k) / 2;
                if (y < cdf[km]) k = km;
                else kl = km;
            }
        } else {
            for (;;) {
                u = 0.645 * doub();
                v = -0.63 + 1.25 * doub();
                v2 = SQR(v);
                // 尝试挤压以快速拒绝：
                if (v >= 0.) {
                    if (v2 > 6.5 * u * (0.645 - u) * (u + 0.2)) continue;
                } else {
                    if (v2 > 8.4 * u * (0.645 - u) * (u + 0.1)) continue;
                }
                k = Int(floor(sq * (v / u) + np + 0.5));
                if (k < 0 || k > n) continue;
                u2 = SQR(u);
                // 尝试挤压以快速接受：
                if (v >= 0.) {
                    if (v2 < 12.25 * u2 * (0.615 - u) * (0.92 - u)) break;
                } else {
                    if (v2 < 7.84 * u2 * (0.615 - u) * (1.2 - u)) break;
                }
                b = sq * exp(glnp + k * plog + (n - k) * pclog
                    - (n < 1024 ? logfact[k] + logfact[n - k]
                       : gammln(k + 1.) + gammln(n - k + 1.)));
                if (u2 < b) break;
            }
        }
        if (p != pp) k = n - k;
        return k;
    }
};

如果你处于这样一种情形：对许多不同的 $n$ 和/或 $p$ 值，每次只生成一个或少数几个随机数，那么你需要重构代码，使得 $n$ 和 $p$ 可以在不创建新对象实例且不重新初始化底层 Ran 生成器的情况下进行更改。

7.3.14 当你需要更高速度时

在特定情况下，你可以通过一些简化来获得更高的速度。以下是一些建议：

• 本节中的所有算法都可以通过使用 §7.1 中的 Ranq1 而非 Ran 来显著加速。我们不知道有任何理由不这样做。如果你将 Ranq1 的算法内联编码，从而消除函数调用，还可以进一步提高速度。

• 如果你使用 Poissondev 或 Binomialdev 且 $\lambda$ 或 $n$ 的值很大，上述代码会回退到调用较慢的 gammln。你可以改为增加存储表的长度。

• 对于 $\lambda <20$ 的泊松分布随机数，你可能希望使用存储的 cdf 表结合二分查找来确定 $k$ 的值。Binomialdev 中的代码展示了如何实现这一点。

• 如果你需要的是小 $n$ 的二项分布随机数，你可以轻松修改 Binomialdev 中的代码，使得每次执行位并行代码时生成多个随机数（实际上约为 $64/n$ 个）。

• 你需要精确的随机数，还是近似值就足够了？如果你感兴趣的分布可以用正态分布近似，考虑使用上面的 Normaldev，特别是如果你也将均匀随机数生成内联编码的话。

• 如果你将恰好 12 个均匀分布随机数 $U(0,1)$ 相加然后减去 6，就能得到一个相当好的标准正态分布随机数 $N(0,1)$ 的近似。在通用 CPU 上，这肯定比 Normaldev 慢（更不用说精度更低了）。然而，据报道，某些专用信号处理芯片可以并行地使用整数算术完成所有这些操作。

参见 Gentle [3]、Ripley [4]、Devroye [6]、Bratley [7] 和 Knuth [8] 以了解更多的算法。

参考文献

Kinderman, A.J. 和 Monahan, J.F. 1977, “利用均匀分布随机变量之比生成随机变量的计算机方法”，《ACM 数学软件汇刊》，第 3 卷，第 257–260 页。[1]

Leva, J.L. 1992, “一种快速的正态随机数生成器”，《ACM 数学软件汇刊》，第 18 卷，第 4 期，第 449–453 页。[2]

Gentle, J.E. 2003, 《随机数生成与蒙特卡罗方法》，第 2 版（纽约：Springer），第 4–5 章。[3]

Ripley, B.D. 1987, 《随机模拟》（纽约：Wiley）。[4]

Marsaglia, G. 和 Tsang W-W. 2000, “一种生成 Gamma 分布随机变量的简单方法”，《ACM 数学软件汇刊》，第 26 卷，第 3 期，第 363–372 页。[5]

Devroye, L. 1986, 《非均匀随机变量生成》（纽约：Springer）。[6]

Bratley, P.、Fox, B.L. 和 Schrage, E.L. 1983, 《模拟指南》，第 2 版（纽约：Springer）。[7]

Knuth, D.E. 1997, 《半数值算法》，第 3 版，《计算机程序设计艺术》第 2 卷（马萨诸塞州雷丁：Addison-Wesley），第 125 页及以后。[8]

7.4 多元正态随机变量

一个维度为 $M$ 的多元随机变量是 $M$ 维空间中的一个点。其坐标是一个向量，该向量的 $M$ 个分量均为随机变量——但通常既非相互独立，也非同分布。多元正态随机变量这一特殊情况由多维高斯密度函数定义：

$$N(\mathbf{x}\mid \mu ,\Sigma )=\frac{1}{(2\pi {)}^{M/2}\det (\Sigma {)}^{1/2}}\exp \left[-\frac{1}{2}(\mathbf{x}-\mu )\cdot {\Sigma }^{-1}\cdot (\mathbf{x}-\mu )\right]$$

其中参数 $\mathbf{\mu }$ 是一个向量，表示分布的均值；参数 $\mathbf{\Sigma }$ 是一个对称正定矩阵，表示分布的协方差矩阵。

有一种通用的方法可以构造具有指定协方差 $\mathbf{\Sigma }$ 和均值 $\mathbf{\mu }$ 的向量随机变量 $\mathbf{x}$：首先从一个各分量相互独立、均值为零、方差为 1 的随机向量 $\mathbf{y}$ 出发；然后利用 Cholesky 分解（§2.9）将 $\mathbf{\Sigma }$ 分解为一个下三角矩阵 $\mathbf{L}$ 与其转置的乘积：

$$\mathbf{\Sigma }=\mathbf{L}{\mathbf{L}}^T$$

由于 $\mathbf{\Sigma }$ 是正定矩阵，这种分解总是可行的，并且对于每个不同的 $\mathbf{\Sigma }$ 只需执行一次。接下来，每当需要一个新的随机向量 $\mathbf{x}$ 时，用独立的标准正态随机变量填充 $\mathbf{y}$，然后构造：

$$\mathbf{x}=\mathbf{Ly}+\mathbf{\mu }$$

证明是直接的，其中尖括号表示期望值：由于 $y_i$ 的分量相互独立且方差为 1，我们有

$$⟨\mathbf{y}\otimes \mathbf{y}⟩=1$$

其中 $1$ 是单位矩阵。于是，

$$\begin{array}{rl}& ⟨(\mathbf{x}-\mathbf{\mu })\otimes (\mathbf{x}-\mathbf{\mu })⟩=⟨(\mathbf{Ly})\otimes (\mathbf{Ly})⟩\\ & =⟨\mathbf{L}(\mathbf{y}\otimes \mathbf{y}){\mathbf{L}}^T⟩=\mathbf{L}⟨\mathbf{y}\otimes \mathbf{y}⟩{\mathbf{L}}^T\\ & ={\mathbf{LL}}^T=\Sigma \end{array}$$

尽管该方法具有一般性，但它通常仅适用于多元正态随机变量。原因在于，虽然 $\mathbf{x}$ 的分量确实具有正确的均值和协方差结构，但其具体的分布形式通常并不“良好”。$x_i$ 是 $y_i$ 的线性组合，而一般而言，随机变量的线性组合的分布是其各自分布的复杂卷积。

然而，对于高斯分布，我们确实拥有“良好”的性质：所有正态随机变量的线性组合本身仍服从正态分布，并且完全由其均值和协方差结构决定。因此，只要我们将 $\mathbf{y}$ 的每个分量都填充为标准正态随机变量：

$$y_i\sim \mathrm{N}(0,1)$$

那么由式 (7.4.3) 生成的随机向量 $\mathbf{x}$ 将服从式 (7.4.1) 所定义的分布。

实现起来非常直接，因为 Cholesky 分解不仅完成了矩阵分解，还提供了一种高效进行矩阵乘法的方法，充分利用了 $\mathbf{L}$ 的三角结构。正态随机变量的生成以内联方式实现以提高效率，与 §7.3 中的 Normaldev 完全相同。

struct Multinormaldev : Ran {
    // 用于生成多元正态随机变量的结构体。
    Int mm;
    VecDoub mean;
    MatDoub var;
    Cholesky chol;
    VecDoub spt, pt;

    Multinormaldev(Ullong j, VecDoub &mmean, MatDoub &vvar) 
        : Ran(j), mm(mmean.size()), mean(mmean), var(vvar), chol(var), spt(mm), pt(mm) {
        // 构造函数。参数为随机数生成器种子、（向量）均值和（矩阵）协方差。
        // 在此处对协方差矩阵进行 Cholesky 分解。
        if (var.ncols() != mm || var.nrows() != mm) throw("bad sizes");
    }

    VecDoub &dev() {
        // 返回一个多元正态随机向量。
        Int i;
        Doub u, v, x, y, q;
        for (i = 0; i < mm; i++) {
            do {
                u = doub();
                v = 1.7156 * (doub() - 0.5);
                x = u - 0.449871;
                y = abs(v) + 0.386595;
                q = SQR(x) + y * (0.19600 * y - 0.25472 * x);
            } while (q > 0.27597 && (q > 0.27846 || SQR(v) > -4. * log(u) * SQR(u)));
            spt[i] = v / u;
        }
        chol.emult(spt, pt); // 应用公式 (7.4.3)。
        for (i = 0; i < mm; i++) pt[i] += mean[i];
        return pt;
    }
};

7.4.1 多个随机变量的去相关

虽然与随机数的生成没有直接关系，但这里恰好可以指出如何以相反的方式使用 Cholesky 分解，即寻找相关随机变量的线性组合，使其彼此不相关。在此应用中，我们给定一个向量 $\mathbf{x}$，其分量具有已知的协方差矩阵 $\Sigma$ 和均值 $\mu$。将 $\Sigma$ 按照公式 (7.4.2) 进行分解，我们断言

$$\mathbf{y}={\mathbf{L}}^{-1}(\mathbf{x}-\mu )$$

的各个分量互不相关，且每个分量的方差均为 1。证明如下：

$$\begin{array}{rl}& ⟨\mathbf{y}\otimes \mathbf{y}⟩=⟨({\mathbf{L}}^{-1}[\mathbf{x}-\mathbf{\mu }])\otimes ({\mathbf{L}}^{-1}[\mathbf{x}-\mathbf{\mu }])⟩\\ & ={\mathbf{L}}^{-1}⟨(\mathbf{x}-\mathbf{\mu })\otimes (\mathbf{x}-\mathbf{\mu })⟩{\mathbf{L}}^{-1T}\\ & ={\mathbf{L}}^{-1}\Sigma {\mathbf{L}}^{-1T}={\mathbf{L}}^{-1}\mathbf{L}{\mathbf{L}}^T{\mathbf{L}}^{-1T}=1\end{array}$$

需要注意的是，这种线性组合并非唯一。事实上，一旦你得到了一个分量互不相关的向量 $\mathbf{y}$，你可以对其执行任意旋转操作，结果仍然保持分量互不相关。特别地，若 $\mathbf{K}$ 是一个正交矩阵，即满足

$${\mathbf{K}}^T\mathbf{K}=\mathbf{K}{\mathbf{K}}^T=1$$

则有

$$⟨(\mathbf{Ky})\otimes (\mathbf{Ky})⟩=\mathbf{K}⟨\mathbf{y}\otimes \mathbf{y}⟩{\mathbf{K}}^T=\mathbf{K}{\mathbf{K}}^T=1$$

Cholesky 分解的一个常见（但较慢）的替代方法是使用 Jacobi 变换（§11.1）将 $\Sigma$ 分解为

$$\mathbf{\Sigma }=\mathbf{V}\mathrm{diag}({\sigma }_i^2){\mathbf{V}}^T$$

其中 $\mathbf{V}$ 是由特征向量组成的正交矩阵，${\sigma }_i$ 是（新的）不相关变量的标准差。此时，$\mathbf{V}\mathrm{diag}({\sigma }_i)$ 在上述证明中扮演了 $\mathbf{L}$ 的角色。

§16.1.1 节进一步讨论了 Cholesky 分解在多元随机变量中的一些其他应用。

7.5 线性反馈移位寄存器

线性反馈移位寄存器（LFSR）由一个状态向量和一种特定的更新规则组成。状态向量通常是 32 位或 64 位字中的位集合，但有时也可以是一个数组中的多个字。要成为 LFSR，其更新规则必须对当前状态中的位（或字）进行线性组合，然后将结果移入状态向量的一端，而状态向量另一端的最旧值则被移出并丢弃。LFSR 的输出是由这些新移入的位（或字）组成的序列。

对于单个比特，“线性”意味着模 2 运算，这等价于使用逻辑异或（XOR）作为加法运算（$+$），逻辑与（AND）作为乘法运算（$\times$）。然而，使用算术符号书写方程更为方便。因此，对于长度为 $n$ 的 LFSR，上文所述内容可转化为如下形式：

$$\begin{array}{rl}& a_1^{\prime }=\left({\sum }_{j=1}^{n-1}{c}_j{a}_j\right)+a_n\\ & a_i^{\prime }=a_{i-1},i=2,\dots ,n\end{array}$$

图 7.5.1. 从移位寄存器和模 2 的本原多项式生成随机比特的两种相关方法。(a) 选定抽头位的内容通过 XOR（模 2 加法）组合，结果从右侧移入。该方法在硬件中实现最为简便。(b) 选定比特通过与最左侧比特进行 XOR 修改，然后将最左侧比特从右侧移入。该方法在软件中实现最为简便。

此处 ${\mathbf{a}}^{\prime }$ 是由更新规则从 $\mathbf{a}$ 推导出的新状态向量。上述第一行中单独列出 $a_n$ 的原因是其系数 $c_n$ 必须 $\equiv 1$。否则，该 LFSR 的长度将不是 $n$，而仅能达到 $c_i$ 中最后一个非零系数的位置。

此外，将比特编号从 1 开始（而非更习惯的 0）也有其原因（此后我们仅考虑比特向量的情形，而非字向量）。方程 (7.5.1) 的数学性质源于整数模 2 上的多项式性质。与 (7.5.1) 相关联的多项式为

$$P(x)=x^n+c_{n-1}{x}^{n-1}+\cdots +c_2{x}^2+c_{1}x+1$$

其中每个 $c_i$ 的取值为 0 或 1。因此，$c_0$ 与 $c_n$ 一样，虽然存在但隐含地 $\equiv 1$。描述此类特定多项式（如 (7.5.2)）有几种记法。一种是简单列出所有 $c_i$ 非零的 $i$ 值（按惯例包括 $c_n$ 和 $c_0$）。例如，多项式

$$x^{18}+x^5+x^2+x+1$$

可简写为

$$(18,5,2,1,0)$$

另一种记法（假设已知 $n$ 值（此处为 18），且 $c_n=c_0=1$）是从二进制字 $c_{n-1}{c}_{n-2}\cdots c_2{c}_1$（按惯例此时排除 $c_n$ 和 $c_0$）构造一个“序列号”。对于 (7.5.3)，该序列号为 19，即 $2^4+2^1+2^0$。非零的 $c_i$ 通常被称为 LFSR 的“抽头”。

图 7.5.1(a) 展示了多项式 (7.5.3) 和 (7.5.4) 如何作为 18 位寄存器上的更新过程。比特 0 是一个临时位置，用于计算将成为新比特 1 的值。

一个 $n$ 位 LFSR 在输出开始重复之前的最长周期为 $2^n-1$。这是因为其最多有 $2^n$ 个不同的状态，但全零的特殊状态会以周期 1 自我重复。如果你随机选择一个多项式 $P(x)$，那么所构造的生成器通常不会达到最大周期。在模 2 整数上的多项式中，有一部分是不可约的，即无法被因式分解；而在不可约多项式中，又有一部分是本原的，即能生成最大周期的 LFSR。例如，多项式 $x^2+1=(x+1)(x+1)$ 不是不可约的，因此也不是本原的。（记住系数运算需模 2。）多项式 $x^4+x^3+x^2+x+1$ 是不可约的，但并非本原。而多项式 $x^4+x+1$ 既是不可约的，也是本原的。

最大周期 LFSR 常被用作硬件设备中的随机比特源，因为图 7.5.1(a) 所示的逻辑仅需少量门电路，且可实现极高的运行速度。然而，LFSR 在软件应用中并无太大优势，因为在代码中实现公式 (7.5.1) 时，每个非零 $c_i$ 至少需要两次完整的字逻辑运算，而所有这些工作仅产生一个比特的输出。我们称这种方法为“方法 I”。一种更优的软件方法——“方法 II”——乍看之下似乎并非 LFSR，但实际上在数学上等价于 LFSR。如图 7.5.1(b) 所示，该方法在代码中可按如下方式从本原多项式实现：

设 maskp 和 maskn 为两个位掩码，

$$\begin{array}{rl}\mathrm{maskp}\equiv & (0\cdots 0c_{n-1}{c}_{n-2}\cdots c_2{c}_1)\\ \mathrm{maskn}\equiv & (0\cdots 100\cdots 00)\end{array}$$

那么，一个字 $\mathbf{a}$ 按如下方式更新：

$$\begin{array}{r}\mathrm{if}(\mathrm{a}\&\mathrm{maskn})\mathrm{a}=((\mathrm{a}\&\mathrm{maskp})<<1)|1;\\ \mathrm{else}\mathrm{a}<<=1;\end{array}$$

你应该仔细推导上述规则，以验证它与图中所示完全一致。此更新的输出（仍仅为一位）可取为 (a & maskn)，或者实际上也可以取 $\mathbf{a}$ 中任意固定的一位。

LFSR（无论是方法 I 还是方法 II）有时通过拼接 $m$ 次连续更新的输出位来生成随机的 $m$ 位字（或者对于方法 I，等价地，在每 $m$ 次更新后取状态的低 $m$ 位）。这通常是个糟糕的主意，因为生成的字通常无法通过一些标准的随机性统计检验。尤其当 $m$ 与 $2^n-1$ 不互质时，情况更糟，此时该方法甚至不能均匀地生成所有 $m$ 位字。

接下来，我们将发展一些理论，以揭示方法 I 与方法 II 之间的关系，并由此导出一个例程，用于检验任意给定的多项式（以系数 $c_i$ 的位串形式表示）是否为本原多项式。但如果你现在只需要一张本原多项式表以便立即使用，下一页就提供了一张。

由于更新规则（7.5.1）是线性的，它可以写成一个矩阵 $\mathbf{M}$，该矩阵从左侧乘以一个位列向量 $\mathbf{a}$，从而产生一个更新后的状态 ${\mathbf{a}}^{\prime }$。（注意，$\mathbf{a}$ 的低位从列向量的顶部开始。）我们可以直接写出：

$$\mathbf{M}=\left[\begin{array}{llllll}{c}_1& c_2& \cdots & c_{n-2}& c_{n-1}& 1\\ 1& 0& \cdots & 0& 0& 0\\ 0& 1& \cdots & 0& 0& 0\\ ⋮& ⋮& & ⋮& ⋮& ⋮\\ 0& 0& \cdots & 1& 0& 0\\ 0& 0& \cdots & 0& 1& 0\end{array}\right]$$

一些模 2 的本原多项式（引自 Watson [1]）
(1, 0)				(51, 6, 3, 1, 0)
(2, 1, 0)				(52, 3, 0)
(3, 1, 0)				(53, 6, 2, 1, 0)
(4, 1, 0)				(54, 6, 5, 4, 3, 2, 0)
(5, 2, 0)				(55, 6, 2, 1, 0)
(6, 1, 0)				(56, 7, 4, 2, 0)
(7, 1, 0)				(57, 5, 3, 2, 0)
(8, 4, 3, 2, 0)				(58, 6, 5, 1, 0)
(9, 4, 0)				(59, 6, 5, 4, 3, 1, 0)
(10, 3, 0)				(60, 1, 0)
(11, 2, 0)				(61, 5, 2, 1, 0)
(12, 6, 4, 1, 0)				(62, 6, 5, 3, 0)
(13, 4, 3, 1, 0)				(63, 1, 0)
(14, 5, 3, 1, 0)				(64, 4, 3, 1, 0)
(15, 1, 0)				(65, 4, 3, 1, 0)
(16, 5, 3, 2, 0)				(66, 8, 6, 5, 3, 2, 0)
(17, 3, 0)				(67, 5, 2, 1, 0)
(18, 5, 2, 1, 0)				(68, 7, 5, 1, 0)
(19, 5, 2, 1, 0)				(69, 6, 5, 2, 0)
(20, 3, 0)				(70, 5, 3, 1, 0)
(21, 2, 0)				(71, 5, 3, 1, 0)
(22, 1, 0)				(72, 6, 4, 3, 2, 1, 0)
(23, 5, 0)				(73, 4, 3, 2, 0)
(24, 4, 3, 1, 0)				(74, 7, 4, 3, 0)
(25, 3, 0)				(75, 6, 3, 1, 0)
(26, 6, 2, 1, 0)				(76, 5, 4, 2, 0)
(27, 5, 2, 1, 0)				(77, 6, 5, 2, 0)
(28, 3, 0)				(78, 7, 2, 1, 0)
(29, 2, 0)				(79, 4, 3, 2, 0)
(30, 6, 4, 1, 0)				(80, 7, 5, 3, 2, 1, 0)
(31, 3, 0)				(81, 4, 0)
(32, 7, 5, 3, 2, 1, 0)			(82, 8, 7, 6, 4, 1, 0)
(33, 6, 4, 1, 0)			(83, 7, 4, 2, 0)
(34, 7, 6, 5, 2, 1, 0)			(84, 8, 7, 5, 3, 1, 0)
(35, 2, 0)				(85, 8, 2, 1, 0)
(36, 6, 5, 4, 2, 1, 0)			(86, 6, 5, 2, 0)
(37, 5, 4, 3, 2, 1, 0)			(87, 7, 5, 1, 0)
(38, 6, 5, 1, 0)			(88, 8, 5, 4, 3, 2, 1, 0)
(39, 4, 0)			(89, 6, 5, 3, 0)
(40, 5, 4, 3, 0)			(90, 5, 3, 2, 0)
(41, 3, 0)			(91, 7, 6, 5, 3, 2, 0)
(42, 5, 4, 3, 2, 1, 0)			(92, 6, 5, 2, 0)
(43, 6, 4, 3, 0)			(93, 2, 0)
(44, 6, 5, 2, 0)			(94, 6, 5, 1, 0)
(45, 4, 3, 1, 0)			(95, 6, 5, 4, 2, 1, 0)
(46, 8, 5, 3, 2, 1, 0)			(96, 7, 6, 4, 3, 2, 0)
(47, 5, 0)			(97, 6, 0)
(48, 7, 5, 4, 2, 1, 0)			(98, 7, 4, 3, 2, 1, 0)
(49, 6, 5, 4, 0)			(99, 7, 5, 4, 0)
(50, 4, 3, 2, 0)			(100, 8, 7, 2, 0)

矩阵 $\mathbf{M}$ 需满足什么条件才能构成一个全周期生成器，从而证明系数为 $c_i$ 的多项式是本原的？显然我们必须有

$${\mathbf{M}}^{(2^n-1)}=1$$

其中 $1$ 是单位矩阵。这表明周期（或其某个倍数）为 $2^n-1$。但唯一可能的倍数是 $2^n-1$ 的因数。为了排除这些因数并确保全周期，我们只需验证

$${\mathbf{M}}^{q_k}\ne 1,q_k\equiv (2^n-1)/f_k$$

其中 $f_k$ 是 $2^n-1$ 的每个素因子。（这正是我们在 §7.1.2 中描述但未证明的矩阵 $\mathbf{T}$ 测试背后的逻辑。）

乍看之下，计算方程 (7.5.8) 和 (7.5.9) 中 $\mathbf{M}$ 的巨大幂次似乎令人望而生畏。但通过 $\mathbf{M}$ 的重复平方法，每个这样的幂次仅需大约 $n$（例如 32 或 64）次矩阵乘法。而且，由于所有运算都是模 2 进行的，因此绝不会发生溢出！条件 (7.5.8) 和 (7.5.9) 实际上是一种高效检验多项式是否本原的方法。以下代码实现了该检验。注意，你必须根据所选的 $n$（代码中称为 $\mathbb{N}$）自定义构造函数中的常量，特别是 $2^n-1$ 的素因子。这里展示了 $n=32$ 的情况。除了这一自定义外，所给代码对 $n\le 64$ 均有效。检验的输入是待测多项式的“序列号”，如方程 (7.5.4) 后所定义。声明一个 Primpolytest 结构的实例后，你可以反复调用其 test() 方法来检验多个多项式。为了使 Primpolytest 完全自包含，矩阵被实现为线性数组，且该结构从头构建了所需的少量矩阵运算。这虽然不够优雅，但很有效。

struct Primpolytest {
    // 在模 2 整数上检验多项式是否本原。
    Int N, nfactors;
    VecUllong factors;
    VecInt t, a, p;

    Primpolytest() : N(32), nfactors(5), factors(nfactors), t(N*N),
        a(N*N), p(N*N) {
        // 构造函数。常量针对 32 位 LFSR。
        Ullong factordata[5] = {3,5,17,257,65537};
        for (Int i=0; i<nfactors; i++) factors[i] = factordata[i];
    }

    Int ispident() {
        Int i,j;
        for (i=0; i<N; i++) for (j=0; j<N; j++) {
            if (i == j) { if (p[i*N+j] != 1) return 0; }
            else { if (p[i*N+j] != 0) return 0; }
        }
        return 1;
    }

    void mattimesseq(VecInt &a, VecInt &b) {
        // 矩阵乘法工具：a *= b。
        Int i,j,k, sum;
        VecInt tmp(N*N);
        for (i=0; i<N; i++) for (j=0; j<N; j++) {
            sum = 0;
            for (k=0; k<N; k++) sum += a[i*N+k] * b[k*N+j];
            tmp[i*N+j] = sum & 1;
        }
        for (k=0; k<N*N; k++) a[k] = tmp[k];
    }

    void matpow(Ullong n) {
        // 通过连续平方计算矩阵 p = a^n。
        Int k;
        for (k=0; k<N*N; k++) p[k] = 0;
        for (k=0; k<N; k++) p[k*N+k] = 1;
        while (1) {
            if (n & 1) mattimesseq(p,a);
            n >>= 1;
            if (n == 0) break;
            mattimesseq(a,a);
        }
    }

    Int test(Ullong n) {
        // 主检验例程。若序列号为 n 的多项式（见正文）是本原的则返回 1，否则返回 0。
        Int i,k,j;
        Ullong pow, tnm1, nn = n;
        tnm1 = ((Ullong)1 << N) - 1;
        if (n > (tnm1 >> 1)) throw("not a polynomial of degree N");
        for (k=0; k<N*N; k++) t[k] = 0;
        // 在 t 中构造更新矩阵。
        for (i=1; i<N; i++) t[i*N+(i-1)] = 1;
        j=0;
        while (nn) {
            if (nn & 1) t[j] = 1;
            nn >>= 1;
            j++;
        }
        t[N-1] = 1;
        for (k=0; k<N*N; k++) a[k] = t[k];
        matpow(tnm1);
        if (ispident() != 1) return 0;
        for (i=0; i<nfactors; i++) {
            pow = tnm1 / factors[i];
            for (k=0; k<N*N; k++) a[k] = t[k];
            matpow(pow);
            if (ispident() == 1) return 0;
        }
        return 1;
    }
};

将此方法推广到 $n>64$ 或模数 $p$ 为 2 以外的素数是直接的。如果 $p^n>2^{64}$，你需要用多字二进制表示整数 $p^n-1$ 及其与素因子的商，以便 matpow 仍能通过连续平方求幂。注意，计算时间大致按 $O(n^4)$ 增长，因此 $n=64$ 很快，而 $n=1024$ 则会相当耗时。

一些 $n=32$ 位的随机本原多项式（以十进制序列表示）为：2046052277、1186898897、221421833、55334070、1225518245、216563424、1532859853、1735381519、2049267032、1363072601 和 130420448。一些 $n=64$ 位的随机本原多项式为：926773948609480634、3195735403700392248、4407129700254524327、256457582706860311、5017679982664373343 和 1723461400905116882。

给定一个满足方程 (7.5.8) 和 (7.5.9) 的矩阵 $\mathbf{M}$，存在一些相关的矩阵也满足这些关系。一个例子是 $\mathbf{M}$ 的逆矩阵，你可以轻松验证其为

$${\mathbf{M}}^{-1}=\left[\begin{array}{llllll}0& 1& 0& \dots & 0& 0\\ 0& 0& 1& \dots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& 0& \dots & 0& 1\\ 1& c_1& c_2& \dots & c_{n-2}& c_{n-1}\end{array}\right]$$

这是将状态 ${\mathbf{a}}^{\prime }$ 回退到其前驱状态 $\mathbf{a}$ 的更新规则。你可以轻松地将 (7.5.10) 转换为类似于方程 (7.5.1) 或图 7.5.1(a) 的规则。

另一个满足保证全周期关系的矩阵是 $\mathbf{M}$ 的逆矩阵的转置（或转置矩阵的逆）。

$${\left({\mathbf{M}}^{-1}\right)}^T=\left[\begin{array}{llllll}0& 0& \dots & 0& 0& 1\\ 1& 0& \dots & 0& 0& c_1\\ 0& 1& \dots & 0& 0& c_2\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& \dots & 1& 0& c_{n-2}\\ 0& 0& \dots & 0& 1& c_{n-1}\end{array}\right]$$

惊喜！这正是方法 II，如图 7.5.1(b) 所示。（请自行验证。）

更具体地说，基于本原多项式 $P(x)$ 的方法 II LFSR 输出的比特序列，与使用互反多项式 $x^{n}P(1/x)$ 的方法 I LFSR 输出的序列完全相同。该结论的证明略超出了我们的范围，但其本质原因在于：矩阵 $\mathbf{M}$ 及其转置都是特征多项式（即公式 (7.5.2)）的根，而逆矩阵 ${\mathbf{M}}^{-1}$ 及其转置则都是互反多项式的根。根据定义，互反多项式只是将非零系数的位置首尾对调。例如，公式 (7.5.3) 的互反多项式为 (18, 17, 16, 13, 1)。若一个多项式是本原的，则其互反多项式也是本原的。

不妨做个实验：先运行一个方法 II 的生成器一段时间。然后取其输出中连续的 $n$ 个比特（例如从最高位开始），并将它们作为初始值填入一个方法 I 的移位寄存器中（最低位对应最近的比特）。接着，让这两个方法同步运行，其中方法 I 使用互反多项式。你会发现，两个生成器的输出完全一致。

参考文献

Knuth, D.E. 1997, Seminumerical Algorithms, 3rd ed., vol. 2 of The Art of Computer Programming (Reading, MA: Addison-Wesley), pp. 30ff.

Horowitz, P., and Hill, W. 1989, The Art of Electronics, 2nd ed. (Cambridge, UK: Cambridge University Press), §9.32 – §9.37.

Tausworthe, R.C. 1965, "Random Numbers Generated by Linear Recurrence Modulo Two," Mathematics of Computation, vol. 19, pp. 201–209.

Watson, E.J. 1962, “Primitive Polynomials (Mod 2)," Mathematics of Computation, vol. 16, pp. 368-369.[1]

7.6 哈希表与哈希存储

这是一个奇怪的梦。你身处一间类似邮件收发室的房间，四壁排列着编号的信箱格子。一位名叫哈彻（Mr. Hacher）的先生坐在一张桌子旁，而你站着。墙上挂着一个收件篮。你的任务是从收件篮中取出信件，并将它们分拣到对应的信箱格子里。

但如何分拣呢？信件是按收件人姓名投递的，而信箱却只有编号。这时哈彻先生就派上用场了。你把每封信给他看，他立刻告诉你该放入哪个编号的信箱。对于同一个名字，他总是给出相同的编号；而不同的名字则总是对应不同的编号（从而对应唯一的信箱）。

随着时间推移，收件人数量不断增加，空信箱越来越少，直到最终一个不剩。只要新来的信件都是寄给已有信箱持有者的，这并不会造成问题。但有一天，你发现一封信上是一个从未见过的名字。你忐忑不安地把它放到哈彻先生面前……然后你就醒了！

哈彻先生和他的桌子就是一个哈希表。哈希表的行为就好像它维护着一份所有曾见过的哈希键（即收件人姓名）的动态账本，为每个键分配一个唯一的编号，并且在每次收到新查询时，都能遍历所有名字：若该键已存在，则返回之前分配的编号；若为新键，则分配一个新编号。通常还支持删除某个键的操作。

实现哈希表的目标是让所有这些操作（插入、查询、删除等）都只需执行很少的计算机指令，甚至不需要 $O(\log N)$ 的时间复杂度。如果你仔细想想，这其实相当巧妙。即使你设法维护一个有序或按字母排序的键列表，通过二分查找定位某个键的位置仍需 $O(\log N)$ 次操作。哈希表的核心思想是利用随机数技术（§7.1）将哈希键映射为一个介于 0 到 $N-1$ 之间的伪随机整数，其中 $N$ 是信箱的总数。这里我们明确需要的是伪随机而非真正的随机整数，因为同一个键每次必须产生相同的整数。

在理想情况下（大多数时候），这个初始的伪随机整数——即哈希函数的输出，或简称为该键的“哈希值”——就是哈希表返回的结果，也就是哈彻先生给出的编号。然而，有可能（事实上随着不同键的数量接近 $N$，这种可能性会越来越大；一旦超过 $N$，根据鸽巢原理，就必然发生）两个不同的键具有相同的哈希值。因此，哈希表的实现必须包含一种冲突解决策略，以确保即使不同键的哈希值相同，也能返回唯一的整数。

许多厂商对 C++ 标准模板库（STL）的实现提供了名为 hash_map 的哈希表类。不幸的是，截至本文撰写时，hash_map 并非 STL 标准的正式组成部分，且各厂商实现的质量也参差不齐。因此，我们在此自行实现一个哈希表；这样既能更深入理解其原理，又能加入一些本书后续章节（例如 §21.8 和 §21.6）所需的具体功能。

7.6.1 哈希函数对象

所谓哈希函数对象，是指将哈希算法（如 §7.1 所述）与构建哈希表所需的“粘合”逻辑封装在一起的结构。该对象应能将任意类型的键 keyT（其本身可能是一个包含多个数据成员的结构体）映射为我们实现中所用的 64 位伪随机整数。哈希函数对象真正需要知道的关于 keyT 的信息仅仅是其字节长度，即 sizeof(keyT)，因为它并不关心这些字节的具体含义，只关心它们是待哈希键的一部分。因此，我们为哈希函数对象提供一个构造函数，用于指定要哈希的字节数；并通过一个指向键地址的 void* 指针来访问键的内容，从而可以按任意方式处理这些字节。

作为哈希函数对象的第一个例子，我们只需对 §7.1.4 中的哈希函数算法进行简单封装。当 sizeof(keyT) 为 4 或 8 时，这种方法效率很高。

struct Hashfn1 {
    // 封装哈希函数的对象示例，供 Hashmap 类使用。
    Ranhash hashed;
    Int n;
    Hashfn1(Int nn) : n(nn) {}
    Ullong fn(const void *key) {
        Uint *k;
        Ullong *kk;
        switch (n) {
            case 4:
                k = (Uint *)key;
                return hashed.int64(*k);
            case 8:
                kk = (Ullong *)key;
                return hashed.int64(*kk);
            default:
                throw("Hashfn1 is for 4 or 8 byte keys only");
        }
    }
};

（由于 n 在对象生命周期内是常量，每次调用时都进行判断略显低效；若你事先已知 n 的值，应删除不必要的代码。）

更一般地，哈希函数对象可通过逐字节处理的方式支持任意长度的键。这里存在速度与随机性程度之间的权衡。历史上，哈希函数更倾向于速度，采用简单的组合规则，例如：

$$\begin{array}{rl}& h_0=\mathrm{某个固定常数}\\ & h_i=(m\cdot h_{i-1}\mathrm{op}{k}_i)\mathrm{mod}{2}^{32}(i=1\dots K)\end{array}$$

其中 $k_i$ 是键的第 $i$ 个字节（$1\le i\le K$），$m$ 是一个乘数，常用值包括 33、63689 和 $2^{16}+2^6-1$（在第一和第三种情况下，可通过移位和加法实现乘法），而 “op” 可以是加法或按位异或（XOR）。当你使用 32 位无符号整数运算时，模运算会自动完成。然而，由于现代机器上 64 位运算速度很快，我们认为使用小乘数或每次仅改变少数几位的多次操作的时代已经过去。我们更倾向于能通过良好随机性测试的哈希函数。（当你对键的特性有深入了解时，可以设计出比随机更好的哈希函数，但这超出了本文范围。）

哈希函数对象在创建时也可以进行一些初始化（例如构建查找表等）。但与随机数生成器不同，它不能在调用之间保存任何依赖历史的状态，因为对于相同的键，它必须每次都返回相同的哈希值。以下是一个适用于任意长度键的自包含哈希函数对象示例。这也是我们后续将使用的哈希函数对象。

struct Hashfn2 {
    // 另一个封装哈希函数的对象示例，支持任意固定长度的键或以空字符结尾的变长字符串。
    // 哈希算法是自包含的。
    static Ullong hashfntab[256];
    Ullong h;
    Int n;
    Hashfn2(Int nn) : n(nn) {
        if (n == 1) n = 0;
        h = 0x544B2FBACAAF1684LL;
        for (Int j = 0; j < 256; j++) {
            for (Int i = 0; i < 31; i++) {
                h = (h >> 7) ^ h;
                h = (h << 11) ^ h;
                h = (h >> 10) ^ h;
            }
            hashfntab[j] = h;
        }
    }

    Ullong fn(const void *key) {
        Int j;
        char *k = (char *)key;
        h = 0xBB40E64DA205B064LL;
        j = 0;
        while (n ? j++ < n : *k) {
            h = (h * 7664345821815920749LL) ^ hashfntab[(unsigned char)(*k)];
            k++;
        }
        return h;
    }
};

Ullong Hashfn2::hashfntab[256]; // 为查找表定义存储空间。

该方法本质上基于公式 (7.6.1)，但做了两点改进：(i) 使用了一个已知对模 $2^{64}$ 的线性同余随机数生成器而言效果良好的大常数乘数；(ii) 更重要的是，通过查表将每个字节值（0 到 255）替换为一个固定的（但看似随机的）64 位值。此外，还做了一个小调整，使得 Hashfn2 既可用于固定长度的键类型（构造时传入 n > 1），也可用于以空字符结尾的变长字节数组（构造时传入 n = 0 或 1）。

7.6.2 哈希表

所谓哈希表，是指具有梦境中哈彻先生（及其表格）功能的对象，即将任意键转换为指定范围内的唯一整数。让我们直接进入主题。哈希表对象的大致结构如下：

template<class keyT, class hfnT> struct Hashtable {
    // 实例化一个哈希表，用于维护任意键与指定范围内唯一整数之间的一一对应关系。
    Int nhash, nmax, nn, ng;
    VecInt htable, next, garbg;
    VecUllong thehash;
    hfnT hash;
    Hashtable(Int hh, Int nv);
    // 构造函数。参数为哈希表大小和可存储元素（键）的最大数量。
    Int ibget(const keyT &key);
    // 返回先前已设置键对应的整数。
    Int iset(const keyT &key);
    // 为新键返回唯一整数。
    Int ieraase(const keyT &key);
    // 删除一个键。
    Int irerserve();
    // 预留一个整数（无对应键）。
    Int ireliniquish(Int k);
    // 释放一个预留的整数。
};

template<class keyT, class hfnT>
Hashtable<keyT,hfnT>::Hashtable(Int hh, Int nv):
    // 构造函数。设置哈希表大小 nhash 和可容纳元素（键）的最大数量 nmax。适当分配数组。
    hash(sizeof(keyT)), nhash(nh), nmax(nv), nn(0), ng(0),
    htable(nh), next(nv), garbg(nv), thehash(nv) {
    for (Int j=0; j<nh; j++) { htable[j] = -1; }
    // 表示空。
}

Hashtable 对象由两个类名模板化：键的类（可以简单如 int，也可以复杂如多重继承类）和哈希函数对象的类（例如上面的 Hashfn1 或 Hashfn2）。注意哈希函数对象如何利用 keyT 的大小自动创建，因此用户无需知道该值。如果你打算使用以空字符结尾的变长字节数组作为键，则 keyT 的类型应为 char，而非 char*；参见 §7.6.5 的示例。

哈希表对象由两个整数参数创建。其中最重要的是 nm，即最多可存储的对象数量——在梦境中，这相当于房间里的鸽笼数量。目前，假设第二个参数 nh 与 nm 取相同值。

整体方案是通过将哈希函数的输出对 nh 取模，将任意键转换为范围 $0\dots nh-1$ 内的整数，用作数组 htable 的索引。该数组的索引元素要么为 $-1$（表示“空”），要么为范围 $0\dots nm-1$ 内的索引，指向数组 thehash 和 next。（从计算机科学角度看，也可以用指针链接的链表元素实现，但本着数值计算的精神，我们使用数组；两种方式效率大致相当。）

哈希表中的每个元素包含之前分配给该索引的键的 64 位哈希值。我们将两个哈希值的相等视为其键完全相同的充分证据。当然，这并非绝对正确。两个键因偶然产生相同哈希值的概率为 $2^{-64}\sim 5\times 10^{-20}$。为确保无差错运行，哈希表实际上必须存储原始键，而不仅仅是哈希值；但就我们的目的而言，可以接受极小概率下两个元素被混淆的情况。（如果你通常在单个哈希表中存储超过十亿个元素，请不要使用这些例程。不过你早就知道这一点了！）

这种 $10^{-20}$ 的巧合并非哈希冲突（hash collision）所指。哈希冲突是指两个哈希值对 $\mathsf{nh}$ 取模后结果相同，从而指向 $\mathsf{htable}$ 中的同一元素。这种情况并不罕见，我们必须提供处理机制。数组 next 中的元素包含回指到 $\mathsf{thehash}$ 和 next 的索引，即形成一个链表。因此，当两个或多个键映射到同一个值 $i$（$0\le i<\mathsf{nh}$）时，若要检索其中某个特定键，它要么位于 $\mathsf{thehash}[i]$，要么位于从该位置开始、由 $\mathsf{next}[i]$、$\mathsf{next}[\mathsf{next}[i]]$ 等链接而成的（希望很短的）链表中。

现在可以进一步说明参数 $\mathsf{nh}$ 的初始设定值。对于一个已填满所有 nm 个值的完整表，htable 每个元素所附带的链表长度服从泊松分布，其均值 $\lambda \equiv \text{nm/nh}$。因此，较大的 $\lambda$（即 $\mathsf{nh}$ 过小）意味着大量链表遍历，而较小的 $\lambda$（即 $\mathsf{nh}$ 过大）则意味着 htable 中空间浪费。传统经验建议选择 $\lambda \sim 0.75$，此时（假设哈希函数良好）htable 中 47% 为空，67% 的非空元素链表长度为 1（即首次尝试即命中正确键），平均间接寻址次数（遍历 next 指针的步数）为 0.42。若取 $\lambda =1$（即 $\mathsf{nh}=\mathsf{nm}$），则对应值为：37% 表项为空，58% 首次尝试命中，平均间接寻址次数为 0.58。因此，在此范围内，任何选择基本都可接受。通用公式如下：

$$\begin{array}{rl}& \mathrm{空表项比例}=P_{\lambda }(0)=e^{-\lambda }\\ & \mathrm{首次尝试命中率}=P_{\lambda }(1)/[1-P_{\lambda }(0)]=\frac{\lambda e^{-\lambda }}{1-e^{-\lambda }}\\ & \mathrm{平均间接寻址次数}={\sum }_{j=2}^{\infty }\frac{(j-1)P_{\lambda }(j)}{1-P_{\lambda }(0)}=\frac{e^{-\lambda }-1+\lambda }{1-e^{-\lambda }}\end{array}$$

其中 $P_{\lambda }(j)$ 为泊松概率函数。

接下来讨论 Hashtable 内部的具体实现。最容易理解的是“获取”（get）函数，它仅在键先前已被“设置”时返回索引值，否则按约定返回 $-1$。我们的数据结构旨在使该操作尽可能快速。

template<class keyT, class hfnT>
Int Hashtable<keyT,hfnT>::iget(const keyT &key) {
    // 返回键对应的 0..nmax-1 范围内的整数，若该键先前未存储则返回 -1。
    Int j,k;
    Ullong pp = hash(&key);
    j = (Int)(pp % nhash);
    for (k = htable[j]; k != -1; k = next[k]) {
        if (thehash[k] == pp) {
            return k;
        }
    }
    return -1;
}

此处需注意一个语言细节：iget 以常量引用方式接收 key，然后将其地址 &key 传递给哈希函数对象。C++ 允许这样做，因为哈希函数对象的 void 指针参数本身声明为 const。

“设置”（set）键的例程稍复杂些。若键先前已被设置，则返回首次设置时的相同值；若未设置，则初始化必要的链接以备将来使用。

template<class keyT, class hfnT>
Int Hashtable<keyT,hfnT>::iset(const keyT &key) {
    // 返回 0..nmax-1 范围内的整数，此后该整数将对应此键。
    // 若键先前已设置，则返回之前的相同整数。

    Int j,k,kprev;
    Ullong pp = hash(&key);
    j = (Int)(pp % nhash);

    if (htable[j] == -1) {
        k = ng ? garbg[--ng] : nn++;
        htable[j] = k;
    } else {
        for (k = htable[j]; k != -1; k = next[k]) {
            if (thehash[k] == pp) {
                return k;
            }
            kprev = k;
        }
        k = ng ? garbg[--ng] : nn++;
        next[kprev] = k;
    }

    if (k >= nmax) throw("storing too many values");
    thehash[k] = pp;
    next[k] = -1;
    return k;
}

此处简要说明垃圾回收机制。当一个键被删除（通过下面的例程）时，我们希望将其整数释放以供未来“设置”使用，从而确保始终能存储 nmax 个键。若分配一个垃圾数组（garbg）并将其用作可用整数的后进先出栈，则该机制极易实现。上述设置例程在需要新整数时总会检查此栈。（顺便提一下，若我们设计 Hashtable 时使用指针链接的链表元素而非数组，则高效垃圾回收的实现会更困难；参见 Stroustrup [1]。）

template<class keyT, class hfnT>
Int Hashtable<keyT,hfnT>::ierase(const keyT &key) {
    // 删除一个键，返回被删除的整数（范围在 0..nmax-1），若该键之前未设置则返回 -1。
    Int j,k,kprev;
    Ullong pp = hash_fn(&key);
    j = (Int)(pp % nhash);
    if (htable[j] == -1) return -1;       // 键之前未设置。
    kprev = -1;
    for (k = htable[j]; k != -1; k = next[k]) {
        if (thehash[k] == pp) {
            if (kprev == -1) htable[j] = next[k];
            else next[kprev] = next[k];
            garbg[ng++] = k;
            return k;
        }
        kprev = k;
    }
    return -1;
}

最后，Hashtable 提供了保留（reserve）和释放（relinquish）整数的例程，这些整数的范围是 0 到 nmax。当一个整数被保留后，可以保证哈希表不会使用它。下面我们将利用这一特性，方便地构建一种哈希内存，使其能在单个键下存储多个元素。

template<class keyT, class hfnT>
Int Hashtable<keyT,hfnT>::ireserve() {
    // 保留一个 0..max-1 范围内的整数，使其不会被 set() 使用，并返回该整数值。
    Int k = ng ? garbg[--ng] : nn++;
    if (k >= max) throw("reserving too many values");
    next[k] = -2;
    return k;
}

template<class keyT, class hfnT>
Int Hashtable<keyT,hfnT>::irelinquish(Int k) {
    // 将之前通过 reserve() 保留的索引归还到池中，并返回该索引；若该索引之前未被保留，则返回 -1。
    if (next[k] != -2) { return -1; }
    garbg[ng++] = k;
    return k;
}

7.6.3 哈希内存（Hash Memory）

上述 Hashtable 类实现了哈彻先生的任务。在此基础上，我们接下来实现你在梦中的任务，即通过任意键对任意对象进行实际的存储和检索。这被称为哈希内存（hash memory）。

当你向普通计算机内存中存储数据时，之前存储在该位置的值会被覆盖。如果你希望你的哈希内存具有相同的行为，那么从 Hashtable 派生出的哈希内存类 Hash 几乎是微不足道的。该类由三种结构类型模板化：keyT 表示键的类型；elT 表示存储在哈希内存中的元素类型；hfnT 与之前一样，表示封装了你所选哈希函数的对象。

template<class keyT, class elT, class hfnT>
struct Hash : Hashtable<keyT, hfnT> {
    // 扩展 Hashtable 类，增加类型为 elT 的元素存储，并提供存储、检索和删除元素的方法。
    // 所有方法中，key 均以地址形式传递。
    using Hashtable<keyT,hfnT>::iget;
    using Hashtable<keyT,hfnT>::iset;
    using Hashtable<keyT,hfnT>::ierase;
    vector<elT> els;

    Hash(Int nh, Int nm) : Hashtable<keyT, hfnT>(nh, nm), els(nm) {}
    // 构造函数语法与 Hashtable 相同。

    void set(const keyT &key, const elT &el) {
        // 存储一个元素 el。
        els[iset(key)] = el;
    }

    Int get(const keyT &key, elT &el) {
        // 将元素检索到 el 中。若键下未存储元素则返回 0，成功则返回 1。
        Int l1 = iget(key);
        if (l1 < 0) return 0;
        el = els[l1];
        return 1;
    }

    elT& operator[] (const keyT &key) {
        // 使用下标符号存储或检索元素。返回一个可作为左值使用的引用。
        Int l1 = iget(key);
        if (l1 < 0) {
            l1 = iset(key);
            els[l1] = elT();
        }
        return els[l1];
    }

    Int count(const keyT &key) {
        // 返回键下存储的元素数量，即 0 或 1。
        Int l1 = iget(key);
        return (l1 < 0 ? 0 : 1);
    }

    Int erase(const keyT &key) {
        // 删除一个元素。成功返回 1，若键下无元素则返回 0。
        return (ierase(key) < 0 ? 0 : 1);
    }
};

上面的 operator[] 方法旨在用于两种不同的场景。首先，它实现了直观的语法来存储和检索元素，例如：

$$\mathrm{myhash}[some\_key]=\mathrm{rhs}$$

用于存储，以及

$$lhs=myhash[some\_key]$$

用于检索。但请注意，这里引入了一个微小的低效：当首次设置某个元素时，会额外调用一次 get。其次，该方法返回一个非常量引用，不仅可以作为左值使用，还可以取地址，例如：

some_pointer = &myhash[some_key];

现在可以通过该指针多次引用存储的元素，而无需额外的键查找开销。在某些应用中，这可以显著提升效率。当然，你也可以直接使用 set 和 get 方法。

7.6.4 哈希多重映射内存（Hash Multimap Memory）

接下来考虑一种情况：你希望能在同一个键下存储多个元素。如果普通计算机内存具有这种行为，你就可以将一个变量设置为一系列值，并让它记住所有这些值！显然，这比 Hash 对 Hashtable 的扩展要复杂一些。我们将其称为 Mhash，其中 M 代表“多值”（multivalued）或“多重映射”（multimap）。一个需求是提供一种便捷的语法，以便逐个检索单个键对应的多个值。我们通过 getinit 和 getnext 函数实现这一点。此外，在下面的 Mhash 中，nmax 现在表示可存储的最大值数量，而不是键的数量（后者通常更小）。

带注释的代码应该无需太多额外解释即可理解。我们利用 Hashtable 的 reserve 和 relinquish 特性，为所有存储的元素（包括某个键的首次实例（Hashtable 必须知道）和后续实例（对 Hashtable 不可见，但由 Mhash 通过链表 nextsis 管理））建立统一的编号系统。

template<class keyT, class elT, class hfnT>
struct Mhash : Hashtable<keyT,hfnT> {
    // 扩展 Hashtable 类，增加类型为 elT 的元素存储，允许在单个键下存储多个元素。
    using Hashtable<keyT,hfnT>::iget;
    using Hashtable<keyT,hfnT>::iset;
    using Hashtable<keyT,hfnT>::ierase;
    using Hashtable<keyT,hfnT>::ireserve;
    using Hashtable<keyT,hfnT>::irelinquish;
    vector<elT> els;
    VecInt nextsis;
    Int nextget;

    Mhash(Int nh, Int nm);
    Int store(const keyT &key, const elT &el); // 在键下存储一个元素。
    Int erase(const keyT &key, const elT &el); // 删除键下指定的元素。
    Int count(const keyT &key);
    Int getinit(const keyT &key);
    Int getnext(elT &el);
};

template<class keyT, class elT, class hfnT>
Mhash<keyT, elT, hfnT>::Mhash(Int nh, Int nm)
: Hashtable<keyT, hfnT>(nh, nm), nextget(-1), els(nm), nextsis(nm) {
    for (Int j=0; j<nm; j++) { nextsis[j] = -2; } // 初始化为“空”。
}

template<class keyT, class elT, class hfnT>
Int Mhash<keyT, elT, hfnT>::store(const keyT &key, const elT &el) {
    // 在键下存储一个元素 el。返回 0..nmax-1 范围内的索引，表示所使用的存储位置。
    Int j,k;
    j = iset(key);
    if (nextsis[j] == -2) {
        // 首次为该键存储元素
        els[j] = el;
        nextsis[j] = -1; // 标记为链表头
        return j;
    } else {
        // 该键已有元素，追加到链表
        k = ireserve(); // 保留一个新索引
        els[k] = el;
        nextsis[k] = nextsis[j]; // 插入到链表头部
        nextsis[j] = k;
        return k;
    }
}

els[j] = el;
nextsis[j] = -1;
return j;
} else {
    while (nextsis[j] != -1) { j = nextsis[j]; } // 获取一个新索引并将其链接到链表中。
}
}

template<class keyT, class elT, class hfnT>
Int Mhash<keyT, elT, hfnT>::erase(const keyT &key, const elT &el) {
    // 删除先前以 key 存储的元素 el。成功则返回 1，若未找到匹配元素则返回 0。
    // 注意：类型 elT 必须定义 == 操作符。
    Int j = -1, kpp = -1;
    Int k = iget(key);
    while (k >= 0) {
        if (j < 0 && el == els[k]) j = k;
        kpp = kp;
        kp = k;
        k = nextsis[k];
    }
    if (j < 0) return 0;
    if (kpp < 0) {
        ierase(key);
        nextsis[j] = -2;
    } else {
        if (j != kp) els[j] = els[kp];
        nextsis[kpp] = nextsis[kp];
        ireliquish(kp);
        nextsis[kp] = -2;
    }
    return 1;
}

template<class keyT, class elT, class hfnT>
Int Mhash<keyT, elT, hfnT>::count(const keyT &key) {
    // 返回以 key 存储的元素数量，若无则返回 0。
    Int next, n = 1;
    if ((next = iget(key)) < 0) return 0;
    while ((next = nextsis[next]) >= 0) { n++; }
    return n;
}

template<class keyT, class elT, class hfnT>
Int Mhash<keyT, elT, hfnT>::getinit(const keyT &key) {
    // 初始化 nextget，使其指向以 key 存储的第一个元素。成功则返回 1，若无则返回 0。
    nextget = iget(key);
    return ((nextget < 0) ? 0 : 1);
}

template<class keyT, class elT, class hfnT>
Int Mhash<keyT, elT, hfnT>::getnext(elT &el) {
    // 若 nextget 有效，则将其元素复制到 el 中，将 nextget 更新为同一 key 的下一个元素，并返回 1。
    // 否则，不修改 el，并返回 0。
    if (nextget < 0) { return 0; }
    el = els[nextget];
    nextget = nextsis[nextget];
    return 1;
}

方法 getinit 和 getnext 的设计用于如下代码中，其中 myhash 是类型为 Mhash 的变量：

// 检索以单个 key 存储的所有元素 el，并对它们进行处理。
if (myhash.getinit(key)) {
    while (myhash.getnext(el)) {
        // 此处使用返回的元素 el。
    }
}

7.6.5 使用示例

尽管我们详细揭示了 Hash 和 Mhash 类的内部工作原理，但这并不意味着它们难以使用。恰恰相反。以下代码声明了一个用于整数的哈希表，并存储了一些人物的出生年份：

Hash<string, Int, Hashfn2> year(1000, 1000);

year[string("Marie Antoinette")] = 1755;
year[string("Ludwig van Beethoven")] = 1770;
year[string("Charles Babbage")] = 1791;

如上声明，year 最多可容纳 1000 个条目。我们使用 C++ 的 string 类作为键类型。若想知道 Charles 出生时 Marie 的年龄，可以这样写：

Int diff = year[string("Charles Babbage")] - year[string("Marie Antoinette")];
cout << diff << '\n';

这将输出 "36"。

如果不使用 C++ 的 string 类，也可以（尽管不推荐）使用以空字符结尾的 C 字符串作为键，如下所示：

Hash<char*, Int, Hashfn2> yearc(1000, 1000);
yearc["Charles Babbage"] = 1791;

这样可行是因为 Hashfn2 对于看起来只有一字节长的键类型有特殊处理（如前所述）。注意这里不需要使用 [0]；实际上，"Charles Babbage" 作为 char* 传递，Hashfn2 会根据指针找到整个字符串。（语法 yearc["Charles Babbage"] 是等价的，同样传递的是字符串首地址。）

假设我们希望反向操作，即按出生年份为索引，将人名存储到哈希表中。由于同一年可能有多人出生，我们需要使用哈希多重映射（multimap）Mhash：

Mhash<Int, string, Hashfn2> person(1000, 1000);

person.store(1775, string("Jane Austen"));
person.store(1791, string("Charles Babbage"));
person.store(1767, string("Andrew Jackson"));
person.store(1791, string("James Buchanan"));
person.store(1767, string("John Quincy Adams"));
person.store(1770, string("Ludwig van Beethoven"));
person.store(1791, string("Samuel Morse"));
person.store(1755, string("Marie Antoinette"));

当然，将名字存入哈希表的顺序无关紧要。以下代码遍历年份，并打印该年出生的人物：

string str;
for (Int i = 1750; i < 1800; i++) {
    if (person.getinit(i)) {
        cout << '\n' << "born in " << i << ":\n";
        while (person.getnext(str)) cout << str.data() << '\n';
    }
}

输出结果为：

born in 1755:
Marie Antoinette

born in 1767:
Andrew Jackson
John Quincy Adams

born in 1770:
Ludwig van Beethoven

born in 1775:
Jane Austen

born in 1791:
Charles Babbage
James Buchanan
Samuel Morse

注意，在此示例中不能使用以空字符结尾的 C 字符串，因为 C++ 不将它们视为可作为 vector 元素存储的一等对象。当使用 Hash 或 Mhash 处理字符串时，通常最好使用 C++ 的 string 类。

在 §21.2 和 §21.8 中，我们将广泛使用 Hash 和 Mhash 类及其几乎所有成员函数；更多使用示例请参见这些章节。

顺便一提，Hacher 先生的名字源自法语 hacher，意为“剁碎”或“哈希”。

参考文献

Stroustrup, B. 1997, The C++ Programming Language, 3rd ed. (Reading, MA: Addison-Wesley), §17.6.2.[1]

Knuth, D.E. 1997, Sorting and Searching, 3rd ed., vol. 3 of The Art of Computer Programming (Reading, MA: Addison-Wesley), §6.4–§6.5.

Vitter, J.S., and Chen, W-C. 1987, Design and Analysis of Coalesced Hashing (New York: Oxford University Press).

7.7 简单的蒙特卡罗积分

数值方法的灵感可能来自意想不到的源头。“样条”（splines）最初是绘图员使用的柔性木条。“模拟退火”（我们将在§10.12中讨论）则源于热力学类比。而“蒙特卡罗方法”这一名称，谁又能不感受到一丝隐约的迷人魅力呢？

假设我们在一个多维体积 $V$ 中均匀地随机选取 $N$ 个点，记为 $x_0,\dots ,x_{N-1}$。那么蒙特卡罗积分的基本定理可估计函数 $f$ 在该多维体积上的积分：

$$\int fdV\approx V⟨f⟩\pm V\sqrt{\frac{⟨f^2⟩-{⟨f⟩}^2}{N}}$$

其中尖括号表示对 $N$ 个采样点取算术平均：

$$⟨f⟩\equiv \frac{1}{N}\sum _{i=0}^{N-1}f(x_i)⟨f^2⟩\equiv \frac{1}{N}\sum _{i=0}^{N-1}{f}^2(x_i)$$

公式（7.7.1）中的“正负”项是对积分的一倍标准差误差估计，并非严格的误差界；此外，也无法保证误差服从高斯分布，因此该误差项仅应视为可能误差的粗略指示。

假设你想对一个区域 $W$ 上的函数 $g$ 进行积分，但该区域难以直接进行随机采样。例如，$W$ 的形状可能非常复杂。这没有问题。只需找到一个包含 $W$ 且易于采样的区域 $V$，然后定义函数 $f$：当点位于 $W$ 内时，$f=g$；当点位于 $W$ 外（但仍处于采样区域 $V$ 内）时，$f=0$。你应尽量使 $V$ 紧密包围 $W$，因为 $f$ 的零值会增大公式（7.7.1）中的误差估计项。这理所当然：选在 $W$ 外的点不包含任何信息，因此有效采样点数 $N$ 实际上减少了。公式（7.7.1）中的误差估计已考虑了这一点。

图 7.7.1 展示了可用于对复杂区域 $W$ 进行采样的三种可能区域 $V$。其中第一个 $V_1$ 显然是一个糟糕的选择。较好的选择 $V_2$ 可通过选取一对均匀分布的随机数 $(s,t)$，再经线性变换映射为 $(x,y)$ 来实现采样。另一个较好的选择 $V_3$ 可通过以下方式采样：首先用一个均匀随机数按左右两个矩形子区域的面积比例选择其中一个，然后用另外两个随机数在所选矩形内选取一个点。

下面我们创建一个体现上述通用方案的对象。（具体实现代码稍后讨论。）基本思路是通过提供以下参数（作为构造函数参数）来创建一个 MCintegrate 对象：

• 向量 xlo：待采样矩形区域各坐标的下界；
• 向量 xhi：待采样矩形区域各坐标的上界；
• 向量值函数 funcs：返回一个或多个需同时积分的函数值；
• 布尔函数：判断某点是否位于我们希望积分的（可能很复杂的）区域 $W$ 内；该点已保证位于由 xlo 和 xhi 定义的区域 $V$ 内；
• 映射函数（稍后讨论），若无映射函数或你不想关注细节，可设为 NULL；
• 随机数生成器的种子。

MCintegrate 对象具有如下结构：

struct MCIntegrate {
    // 用于在 ndim 维区域中对一个或多个函数进行蒙特卡罗积分的对象。
    Int ndim, nfun, n;          // 维数、函数个数、采样点数。
    VecDoub ff, fferr;          // 结果：积分值及其标准误差。
    VecDoub xlo, xhi, x, xx, fn, sf, sferr;
    Doub vol;                   // 区域 V 的体积。

    // mcintegrate.h

    VecDoub (*funcsp)(const VecDoub &);   // 指向用户提供的函数。
    VecDoub (*xmapp)(const VecDoub &);    // 映射函数指针。
    Bool (*inregionp)(const VecDoub &);   // 判断点是否在区域 W 内的函数指针。
    Ran ran;                              // 随机数生成器。
};

随机数生成器。

MCintegrate(const VecDoub &xlow, const VecDoub &xhigh, 
            VecDoub funcs(const VecDoub &), 
            Bool irregion(const VecDoub &), 
            VecDoub xmap(const VecDoub &), 
            Int ranseed);

构造函数。参数顺序如上文所列。

成员函数 step 用于添加采样点，其参数指定采样点的数量。成员函数 calcanswers 用于更新向量 ff 和 fferr，它们分别包含函数的蒙特卡洛积分估计值及其估计误差。你可以检查这些值，如果需要，还可以再次调用 step 和 calcanswers 以进一步减小误差。

下面通过一个具体例子来说明该方法的基本原理。假设我们希望计算一个形状复杂的物体的质量及其质心位置，该物体是环面与一个大盒子的两个面相交的部分。具体来说，该物体由以下三个条件同时定义：

$$z^2+{\left(\sqrt{x^2+y^2}-3\right)}^2\le 1$$

（以原点为中心的环面，主半径为 3，次半径为 1）

$$x\ge 1y\ge -3$$

（盒子的两个面；见图 7.7.2）。暂时假设该物体具有恒定密度 $\rho =1$。

图 7.7.2. 蒙特卡洛积分示例（见正文）。感兴趣区域是环面的一部分，由两个平面的交线所界定。该区域的积分限难以用解析闭式表达，因此蒙特卡洛方法是一种有效的技术。

我们希望估计以下在该复杂物体内部的积分：

$$\int \rho dxdydz\int x\rho dxdydz\int y\rho dxdydz\int z\rho dxdydz$$

质心的坐标将是后三个积分（一次矩）与第一个积分（质量）的比值。

要使用 MCintegrate 对象，我们首先编写描述被积函数和积分区域 $W$（位于盒子 $V$ 内部）的函数：

VecDoub torusfuncs(const VecDoub &x) {
    // 返回公式 (7.7.5) 中的被积函数，其中 ρ = 1。
    Doub den = 1;
    VecDoub f(4);
    f[0] = den;
    for (Int i = 1; i < 4; i++) f[i] = x[i-1] * den;
    return f;
}

Bool torusregion(const VecDoub &x) {
    // 返回不等式 (7.7.3)。
    return SQR(x[2]) + SQR(sqrt(SQR(x[0]) + SQR(x[1])) - 3.) <= 1.;
}

实际执行积分的代码现在变得非常简单：

VecDoub xlo(3), xhi(3);
xlo[0] = 1.; xhi[0] = 4.;
xlo[1] = -3.; xhi[1] = 4.;
xlo[2] = -1.; xhi[2] = 1.;
MCIntegrate mymc(xlo, xhi, torusfuncs, torusregion, NULL, 10201);
mymc.step(1000000);
mymc.calcanswers();

这里我们通过 xlo 和 xhi 指定了盒子 $V$，创建了一个 MCintegrate 实例，并进行了 100 万次采样，然后更新了结果 mymc.ff 和 mymc.fferr，这些结果可用于打印或其他用途。

7.7.1 变量替换

在蒙特卡洛积分中，变量替换通常能带来极大的好处。例如，假设我们要计算相同的积分，但环面片段的密度是 $z$ 的强函数，具体按如下方式变化：

$$\rho (x,y,z)=e^{5z}$$

一种做法是在 torusfuncs 中简单地将语句

Doub den = 1.;

替换为

Doub den = exp(5.*x[2]);

这种方法可行，但并非最优。由于 (7.7.6) 随着 $z$ 减小（直至其下限 $-1$）迅速趋近于零，大多数采样点对质量和力矩的求和几乎毫无贡献。这些点实际上被浪费了，其效率损失几乎与落在区域 $W$ 之外的点一样严重。通过变量替换（与 §7.3 中的变换方法完全相同）可以解决这个问题。令

$$ds=e^{5z}dz\text{即}s=\frac{1}{5}{e}^{5z},z=\frac{1}{5}\ln (5s)$$

于是 $\rho dz=ds$，而 $z$ 的积分限 $-1<z<1$ 变为 $0.00135<s<29.682$。

MCintegrate 对象知道你可能希望执行这样的操作。如果它发现参数 xmap 不为 NULL，就会假定由 xlo 和 xhi 定义的采样区域并非处于物理空间中，而是需要先映射到物理空间，之后才能计算函数值或区域边界。因此，为了实现变量替换，我们无需修改 torusfuncs 或 torusregion，但需要修改 xlo 和 xhi，并为参数 xmap 提供如下函数：

返回由公式 (7.7.7) 中最后一个方程定义的从 $s$ 到 $z$ 的映射，其余坐标则使用恒等映射。

实际积分的代码现在如下所示：

VecDoub slo(3), shi(3);
slo[0] = 1.; shi[0] = 4.;
slo[1] = -3.; shi[1] = 4.;
slo[2] = 0.2*exp(5.*(-1.)); shi[2] = 0.2*exp(5.*(1.));
MCintegrate mymc2(slo, shi, torusfunctions, torusregion, torusmap, 10201);
mymc2.step(1000000);
mymc2.calcanswers();

稍加思考你就会意识到，公式 (7.7.7) 之所以有用，仅仅是因为我们希望消除的被积函数部分（$e^{5z}$）既可解析积分，其积分结果又可解析求逆。（参见 §7.3.2。）一般情况下，这些性质并不成立。那么问题来了：此时该怎么办？答案是：从被积函数中提取出“最佳”的因子，使其既可积分又可求逆。“最佳”的标准是尽可能使剩余的被积函数接近常数。

极限情况颇具启发性：如果你成功地将被积函数 $f$ 变为严格常数，并且已知体积的区域 $V$ 恰好完全包围了目标区域 $W$，那么你计算出的 $f$ 的平均值将严格等于该常数值，而公式 (7.7.1) 中的误差估计将严格为零。实际上，你已经精确完成了积分，蒙特卡罗数值计算变得多余。因此，从这一极端情况退一步讲，只要你能通过变量替换使 $f$ 近似为常数，并且采样区域仅略大于 $W$，就能提高蒙特卡罗积分的精度。文献中通常将这种技术称为方差缩减（variance reduction）。

简单蒙特卡罗积分的根本缺点在于，其精度仅随采样点数 $N$ 的平方根增长。如果你的精度要求不高，或者计算资源充足，那么该方法因其普适性强而非常值得推荐。在 §7.8 和 §7.9 中，我们将看到一些能够“突破 $N$ 的平方根障碍”的技术，在某些情况下可以用更少的函数求值获得更高的精度。

MCintegrate 的实现并无令人意外之处。构造函数存储用户函数的指针，对 funcs 进行一次（原本多余的）调用以确定返回向量的大小，然后据此调整 sum 和 answer 向量的尺寸。step 和 calcanswer 方法分别精确实现了公式 (7.7.1) 和 (7.7.2)。

MCintegrate::MCintegrate(const VecDoub &xlow, const VecDoub &xhigh,
                         VecDoub funcs(const VecDoub),
                         Bool inregion(const VecDoub &),
                         VecDoub xmap(const VecDoub &), Int ranseed)
    : ndim(xlow.size()), n(0), xlo(xlow), xhi(xhigh), x(ndim), xx(ndim),
      funcsp(funcs), xmapp(xmap), inregionp(inregion), vol(1.), ran(ranseed) {
    if (xmapp) nfun = funcs(xmap(xlo)).size();
    else nfun = funcs(xlo).size();
    ff.resize(nfun);
    fferr.resize(nfun);
    fn.resize(nfun);
    sf.assign(nfun, 0.);
    sferr.assign(nfun, 0.);
    for (Int j = 0; j < ndim; j++) vol *= abs(xhi[j] - xlo[j]);
}

void MCintegrate::step(Int nstep) {
    Int i, j;
    for (i = 0; i < nstep; i++) {
        for (j = 0; j < ndim; j++)
            x[j] = xlo[j] + (xhi[j] - xlo[j]) * ran.doub();
        if (xmapp) xx = (*xmapp)(x);
        else xx = x;
        if ((*inregionp)(xx)) {
            fn = (*funcsp)(xx);
            for (j = 0; j < nfun; j++) {
                sf[j] += fn[j];
                sferr[j] += SQR(fn[j]);
            }
        }
    }
    n += nstep;
}

void MCintegrate::calcanswers() {
    for (Int j = 0; j < nfun; j++) {
        ff[j] = vol * sf[j] / n;
        fferr[j] = vol * sqrt((sferr[j] / n - SQR(sf[j] / n)));
    }
}

参考文献

Robert, C.P., and Casella, G. 2006, Monte Carlo Statistical Methods, 2nd ed. (New York: Springer).
Sobol', I.M. 1994, A Primer for the Monte Carlo Method (Boca Raton, FL: CRC Press).
Hammersley, J.M., and Handscomb, D.C. 1964, Monte Carlo Methods (London: Methuen).
Gentle, J.E. 2003, Random Number Generation and Monte Carlo Methods, 2nd ed. (New York: Springer), Chapter 7.
Shreider, Yu. A. (ed.) 1966, The Monte Carlo Method (Oxford: Pergamon).
Kalos, M.H., and Whitlock, P.A. 1986, Monte Carlo Methods: Volume 1: Basics (New York: Wiley).

7.8 拟（即次）随机序列

我们刚刚看到，在 $n$ 维空间中均匀随机地选择 $N$ 个点，会导致蒙特卡罗积分的误差项以 $1/\sqrt{N}$ 的速度减小。本质上，每个新采样点线性地累加到将成为函数平均值的总和中，同时也线性地累加到将成为方差的平方和中（公式 7.7.2）。估计误差来自该方差的平方根，因此具有 $N^{-1/2}$ 的幂律。

然而，这种平方根收敛虽然常见，却并非不可避免。一个简单的反例是：选择位于笛卡尔网格上的采样点，并且每个网格点恰好采样一次（顺序任意）。

此时，蒙特卡罗方法就变成了一种确定性的求积方案——尽管很简单——其相对误差至少以 $N^{-1}$ 的速度减小（如果函数在采样区域边界处光滑地趋于零，或是区域内的周期函数，则收敛更快）。

网格的问题在于，你必须提前决定网格的精细程度，并且必须完成所有网格点的采样。使用网格时，不方便实现“持续采样直到满足某种收敛或终止条件”。人们可能会问，是否存在某种中间方案，即以某种“随机”的方式选择采样点，但又能以某种自回避的方式分布，避免均匀随机点可能出现的聚集现象。

类似的问题也出现在蒙特卡罗积分以外的任务中。例如，我们可能希望在 $n$ 维空间中搜索满足某种（局部可计算）条件的点。当然，为了使该任务在计算上是有意义的，条件必须具有连续性，使得满足条件的点存在于某个有限的 $n$ 维邻域内。然而，我们可能事先并不知道该邻域的大小。我们希望“持续采样直到”找到目标点，并随着采样数量的增加平滑地过渡到更精细的尺度。是否存在比无关联的随机采样更好的方法？

上述问题的答案是“肯定的”。那些比无关联随机点更均匀地填充 $n$ 维空间的 $n$ 元组序列，被称为拟随机序列（quasi-random sequences）。这一术语有些用词不当，因为拟随机序列实际上毫无“随机”可言：它们是精心构造的，实际上是次随机（subrandom）的。拟随机序列中的采样点在某种精确意义上是“最大程度相互回避”的。

一个概念上简单的例子是 Halton 序列 [1]。在一维情况下，序列中的第 $j$ 个数 $H_j$ 通过以下步骤获得：(i) 将 $j$ 表示为某个素数 $b$ 为底的数。（例如，$j=17$ 在底数 $b=3$ 下为 122。）(ii) 将数字反转，并在序列前加上小数点（即底数为 $b$ 的小数点）。（在该例子中，得到 0.221（底数 3）。）结果即为 $H_j$。为了在 $n$ 维空间中生成 $n$ 元组序列，每个分量使用不同素数底数的 Halton 序列。通常使用前 $n$ 个素数。

不难理解 Halton 序列的工作原理：每当 $j$ 的位数增加一位时，$j$ 的数字反转分数的网格就会细化 $b$ 倍。因此，该过程是在一系列越来越精细的笛卡尔网格上填充所有点——并且在每个网格上以一种“最大程度分散”的顺序进行（例如，$j$ 中变化最快的数字控制分数的最高有效位）。

Faure、Sobol'、Niederreiter 等人提出了其他生成拟随机序列的方法。Bratley 和 Fox [2] 提供了很好的综述和参考文献，并讨论了 Antonov 和 Saleev [4] 提出的 Sobol' [3] 序列的一种特别高效的变体。我们现在讨论的正是这一 Antonov-Saleev 变体的实现。

Sobol' 序列直接生成介于 0 和 1 之间的数，表示为长度为 $w$ 位的二进制小数，这些数来自一组 $w$ 个特殊的二进制小数 $V_i$（$i=1,2,\dots ,w$），称为方向数（direction numbers）。在 Sobol 的原始方法中，第 $j$ 个数 $X_j$ 通过对满足“$j$ 的第 $i$ 位非零”这一条件的所有 $V_i$ 进行异或（XOR，按位异或）运算得到。换句话说，随着 $j$ 递增，不同的 $V_i$ 以不同的时间尺度在 $X_j$ 中“闪烁”进出。$V_1$ 的出现与消失最为频繁，而

图 7.8.1. 二维 Sobol' 序列的前 1024 个点。该序列是通过数论方法而非随机方法生成的，因此在任何阶段，后续点都能“知道”如何填补先前生成分布中的空隙。

$V_k$ 仅在每 $2^{k-1}$ 步才从出现变为消失（或反之）。

Antonov 和 Saleev 的贡献在于指出：与其使用整数 $j$ 的二进制位来选择方向数，不如使用 $j$ 的格雷码（Gray code）$G(j)$ 的二进制位。（关于格雷码的快速回顾，请参见 §22.3。）

现在 $G(j)$ 和 $G(j+1)$ 恰好在一个比特位上不同，即在 $j$ 的二进制表示中最右边的零比特位的位置上（必要时可在 $j$ 前面补一个前导零）。由此可得，第 $j+1$ 个 Sobol'-Antonov-Saleev 数可以通过将第 $j$ 个数与单个 $V_i$ 进行异或（XOR）运算得到，其中 $i$ 是 $j$ 的二进制表示中最右边零比特的位置。这使得序列的计算非常高效，我们将在下文看到这一点。

图 7.8.1 绘出了二维 Sobol' 序列生成的前 1024 个点。可以看出，后续的点确实“知道”之前留下的空隙，并以分层的方式不断填充这些空隙。

我们一直推迟到此处才讨论方向数 $V_i$ 是如何生成的。其中涉及一些非平凡的数学内容，因此我们仅提供一个“菜谱式”的总结：每一个不同的 Sobol' 序列（或 $n$ 维序列中的每个分量）都基于一个不同的模 2 域上的本原多项式，即系数为 0 或 1 的多项式，并且能生成最大长度的移位寄存器序列。（模 2 的本原多项式已在 §7.5 中使用，并将在 §22.4 中进一步讨论。）假设 $P$ 是这样一个次数为 $q$ 的多项式：

$$P=x^q+a_1{x}^{q-1}+a_2{x}^{q-2}+\cdots +a_{q-1}x+1$$

我们通过如下 $q$ 项递推关系定义一个整数序列 $M_i$：

$$M_i=2a_1{M}_{i-1}\oplus 2^2{a}_2{M}_{i-2}\oplus \cdots \oplus 2^{q-1}{a}_{q-1}{M}_{i-q+1}\oplus (2^q{M}_{i-q}\oplus M_{i-q})$$

其中 $\oplus$ 表示按位异或（XOR）运算。该递推关系的初始值为：$M_1,\dots ,M_q$ 可以分别任意选取小于 $2,\dots ,2^q$ 的奇整数。然后，方向数 $V_i$ 由下式给出：

$$V_i=M_i/2^{i}i=1,\dots ,w$$

上表列出了所有次数 $q\le 10$ 的模 2 本原多项式。由于系数只能是 0 或 1，且最高次项和常数项的系数恒为 1，因此用一个二进制数来表示中间系数非常方便（高次幂对应更高位）。上表正是采用这种约定。

现在转向 Sobol' 序列的具体实现。在一次初始化调用之后，对函数 sobseq 的连续调用将返回基于表中前 $n$ 个本原多项式的 $n$ 维 Sobol' 序列中的连续点。当前给出的例程最多支持 6 维（即 $n=6$），字长 $w$ 为 30 位。通过修改 MAXBIT（$\equiv w$）和 MAXDIM，并为数组 ip（来自上表的本原多项式）、mdeg（它们的次数）和 iv（递推关系的初始值，见公式 7.8.2）添加更多初始化数据，即可调整这些参数。下一页的第二张表进一步阐明了例程中使用的初始化数据。

subseq 中使用的初始化值
次数	多项式	初始值
1	0	1	(3)	(5)	(15)...
2	1	1	1	(7)	(11)...
3	1	1	3	7	(5)...
3	2	1	3	3	(15)...
4	1	1	1	3	13...
4	4	1	1	5	9...
括号中的值不是自由指定的，而是由该次数对应的递推关系所决定。

const Int MAXBIT=30, MAXDIM=6;
Int j,k,l;
Unt i,im,ipp;
static Int mdeg[MAXDIM]={1,2,3,3,4,4};
static Uint in;
static VecUint ix(MAXDIM);
static NRvector<Unt*> iu(MAXBIT);
static Uint ip[MAXDIM]={0,1,1,2,1,4};
static Uint iv[MAXDIM*MAXBIT]={
    1,1,1,1,1,1,3,1,3,1,1,5,7,7,3,3,5,15,11,5,15,13,9
};
static Doub fac;

if (n < 0) {
    for (k=0; k<MAXDIM; k++) ix[k]=0;
    in=0;
    if (iv[0] != 1) return;
    fac=1.0/(1 << MAXBIT);
    for (j=0, k=0; j<MAXBIT; j++, k=(k+1)%MAXDIM)
        iu[j] = &iv[k]; // 允许一维和二维寻址
    for (k=0; k<MAXDIM; k++) {
        for (j=0; j<mdeg[k]; j++)
            iu[j][k] <<= (MAXBIT-1-j); // 存储的值只需归一化
        for (j=mdeg[k]; j<MAXBIT; j++) {
            ipp=ip[k];
            i=iu[j-mdeg[k]][k];
            i ^= (i >> mdeg[k]);
            for (l=mdeg[k]-1; l>=1; l--) {
                if (ipp & 1) i ^= iu[j-l][k];
                ipp >>= 1;
            }
            iu[j][k]=i;
        }
    }
} else {
    im=in++;
    for (j=0; j<MAXBIT; j++) {
        if (!(im & 1)) break;
        im >>= 1;
    }
}

if (j >= MAXBIT) throw("MAXBIT too small in sobseq");
im=j*MAXDIM;
for (k=0; k<MIN(n,MAXDIM); k++) {
    ix[k] ^= iv[im+k];
    x[k]=ix[k]*fac;
}

那么，Sobol' 序列到底有多好呢？对于 $n$ 维空间中光滑函数的蒙特卡罗积分，答案是：相对误差随样本数 $N$ 的减小速度为 $(\ln N{)}^n/N$，即几乎与 $1/N$ 一样快。举个例子，我们来积分一个在三维空间中环面（甜甜圈）内部非零的函数。若环面的大半径为 $R_0$，小半径为 $r_0$，则小径向坐标 $r$ 定义为：

$$r={\left([(x^2+y^2{)}^{1/2}-R_0{]}^2+z^2\right)}^{1/2}$$

我们尝试如下函数：

$$f(x,y,z)=\{\begin{array}{ll}1+\cos \left(\frac{\pi r^2}{r_0^2}\right)& r<r_0\\ 0& r\ge r_0\end{array}$$

该函数在柱坐标下可解析积分，结果为：

$$\iiint dxdydzf(x,y,z)=2{\pi }^2{r}_0^2{R}_0$$

取参数 $R_0=0.6$，$r_0=0.3$，我们在区域 $-1<x,y,z<1$ 内均匀采样，对公式 (7.8.4) 进行了 100 次连续的蒙特卡罗积分，分别使用不相关的随机点和由 sobseq 生成的 Sobol' 序列。图 7.8.2 显示了结果，绘制了 100 次积分的均方根（r.m.s.）平均误差随采样点数的变化。（对于任何单次积分，误差当然会在正负之间波动，因此对相对误差取对数绘图并不十分直观。）细虚线对应不相关的随机点，显示出熟悉的 $N^{-1/2}$ 渐近行为。细实灰线显示 Sobol' 序列的结果。预期的 $(\ln N{)}^3/N$ 中的对数项在曲线上表现为明显的弯曲，但其渐近的 $N^{-1}$ 趋势是毋庸置疑的。

为了理解图 7.8.2 的重要性，假设我们希望对函数 $f$ 进行蒙特卡罗积分并达到 1% 的精度。Sobol' 序列仅需几千个样本即可达到该精度，而伪随机采样则需要近 10 万个样本。若要求更高精度，这一比例差距还会更大。

另一种情况则不那么有利：当被积函数在采样区域内具有硬边界（不连续）时，例如环面内部为 1、外部为 0 的函数：

$$f(x,y,z)=\{\begin{array}{ll}1& r<r_0\\ 0& r\ge r_0\end{array}$$

其中 $r$ 由公式 (7.8.4) 定义。巧合的是，该函数的解析积分与公式 (7.8.5) 中的函数相同，即 $2{\pi }^2{r}_0^2{R}_0$。

精心设计的分层 Sobol' 序列基于一组笛卡尔网格，但环面的边界与这些网格并无特定关系。结果是，在环面表面附近厚度为 $N^{-1/3}$ 量级的薄层中（包含约 $N^{2/3}$ 个采样点），这些点落在环面内部还是外部本质上是随机的。对该薄层应用平方根定律，可得总和的涨落为 $N^{1/3}$，即蒙特卡罗积分的相对误差为 $N^{-2/3}$。图 7.8.2 中较粗的灰线验证了这一行为。

图 7.8.2. 对于两种不同的被积函数和两种不同的随机点选取方法，蒙特卡罗积分的相对精度随采样点数的变化。拟随机 Sobol' 序列比传统的伪随机序列收敛快得多。当被积函数光滑（“软边界”）时，拟随机采样的效果优于存在阶跃不连续（“硬边界”）的情况。图中曲线为 100 次试验的均方根平均值。

图 7.8.2 中较粗的虚线是使用独立随机点对公式 (7.8.7) 的函数进行积分的结果。虽然 Sobol' 序列在此情形下的优势不如光滑函数情形那么显著，但在 1% 这样的中等精度下仍可带来显著的改进（约 5 倍），且在更高精度下优势更大。

请注意，我们尚未为 sobseq 例程提供从序列中除开头以外的任意点开始的功能，但添加这一特性非常容易。一旦方向数 iv 的初始化完成，第 $j$ 个点就可以直接通过对 $j$ 的格雷码中所有非零位所对应的方向数进行异或（XOR）运算得到，如上所述。

7.8.1 拉丁超立方体

这里我们提到一种与前述方法无关的技术——拉丁方（Latin square）或拉丁超立方体（Latin hypercube）抽样。当你必须在 $N$ 维空间中极其稀疏地抽取 $M$ 个样本点时，这种方法非常有用。例如，你可能希望同时测试汽车在四个不同设计参数下的抗撞性能，但预算只允许使用三辆可报废的汽车。（问题不在于这个计划是否合理——它显然不合理——而在于如何在这种情况下做到最好！）

其基本思想是将每个设计参数（维度）划分为 $M$ 个区间，从而将整个空间划分为 $M^N$ 个单元格。（你可以根据需要选择每个维度中的区间是等长还是不等长。）例如，对于四个参数和三辆汽车的情形，最终会得到 $3\times 3\times 3\times 3=81$ 个单元格。

接下来，按照以下算法选择 $M$ 个单元格作为样本点：首先从 $M^N$ 个单元格中随机选择一个作为第一个样本点。然后剔除所有在任意一个参数上与该点相同的单元格（即划掉同一行、列等的所有单元格），剩下 $(M-1{)}^N$ 个候选单元格。从中随机选择一个，再剔除新的行和列，重复此过程，直到只剩下一个单元格，该单元格即为最后一个样本点。

这种构造的结果是，每个设计参数在其所有子区间中都得到了测试。如果被测系统的响应主要由某个设计参数（主效应）主导，那么通过这种抽样技术就能发现该参数的影响。然而，如果不同设计参数之间存在重要的交互效应，那么拉丁超立方体抽样就不再具有特别的优势。请谨慎使用。

统计学中有一个庞大的领域专门研究实验设计（design of experiments）。关于该主题的一个简明教学性介绍见文献 [5]。

参考文献

Halton, J.H. 1960, "On the Efficiency of Certain Quasi-Random Sequences of Points in Evaluating Multi-dimensional Integrals," Numerische Mathematik, vol. 2, pp. 84–90.[1]

Bratley P., and Fox, B.L. 1988, "Implementing Sobol's Quasirandom Sequence Generator," ACM Transactions on Mathematical Software, vol. 14, pp. 88–100.[2]

Lambert, J.P. 1988, "Quasi-Random Sequences in Numerical Practice," in Numerical Mathematics — Singapore 1988, ISNM vol. 86, R.P. Agarwal, Y.M. Chow, and S.J. Wilson, eds. (Basel: Birkhäuser), pp. 273–284.

Niederreiter, H. 1988, "Quasi-Monte Carlo Methods for Multidimensional Numerical Integration." in Numerical Integration III, ISNM vol. 85, H. Brass and G. Hämmerlin, eds. (Basel: Birkhäuser), pp. 157–171.

Sobol', I.M. 1967, "On the Distribution of Points in a Cube and the Approximate Evaluation of Integrals," USSR Computational Mathematics and Mathematical Physics, vol. 7, no. 4, pp. 86–112.[3]

Antonov, I.A., and Saleev, V.M 1979, "An Economic Method of Computing $l{p}_t$ Sequences," USSR Computational Mathematics and Mathematical Physics, vol. 19, no. 1, pp. 252–256.[4]

Dunn, O.J., and Clark, V.A. 1974, Applied Statistics: Analysis of Variance and Regression (New York, Wiley) [discusses Latin Square].

Czitrom, V. 1999, "One-Factor-at-a-Time Versus Designed Experiments," The American Statistician, vol. 53, pp. 126–131, online at http://www.amstat.org/publications/tas/czitrom.pdf.[5]

7.9 自适应与递归蒙特卡罗方法

本节讨论更高级的蒙特卡罗积分技术。作为这些技术应用的示例，我们介绍了两个相当不同且较为复杂的多维蒙特卡罗程序：vegas [1,2] 和 miser [4]。我们所讨论的技术都属于方差缩减（§7.7）这一大类，但彼此之间又有显著差异。

7.9.1 重要性抽样

重要性抽样（importance sampling）的思想在公式 (7.7.6) 和 (7.7.7) 中已经隐含体现。我们现在以稍加形式化的方式重新讨论它。假设被积函数 $f$ 可以写成一个近似常数的函数 $h$ 与另一个正函数 $g$ 的乘积。那么它在多维区域 $V$ 上的积分为：

$$\int fdV=\int (f/g)gdV=\int hgdV$$

在公式 (7.7.7) 中，我们将公式 (7.9.1) 解释为建议进行变量变换，令 $G$ 为 $g$ 的不定积分。这样 $gdV$ 就成为一个全微分。随后我们应用了蒙特卡罗积分的基本定理，即公式 (7.7.1)。对公式 (7.9.1) 更一般的解释是：我们可以通过对 $h$ 进行抽样来积分 $f$——但不是以均匀概率密度 $dV$ 抽样，而是以非均匀密度 $gdV$ 抽样。在第二种解释下，第一种解释成为特例，即通过变换方法（使用不定积分 $G$，见 §7.3）来生成非均匀抽样 $gdV$。

更直接地，我们可以回过头来将基本定理 (7.7.1) 推广到非均匀抽样的情形：假设点 $x_i$ 在区域 $V$ 内按概率密度 $p$ 选取，且满足

$$\int pdV=1$$

则推广后的基本定理指出：任意函数 $f$ 的积分可利用 $N$ 个样本点 $x_0,\dots ,x_{N-1}$ 估计为

$$I\equiv \int fdV=\int \frac{f}{p}pdV\approx ⟨\frac{f}{p}⟩\pm \sqrt{\frac{⟨f^2/p^2⟩-{⟨f/p⟩}^2}{N}}$$

其中尖括号表示对 $N$ 个点的算术平均，与公式 (7.7.2) 完全相同。与公式 (7.7.1) 类似，这里的“±”项是一个标准差误差估计。注意，公式 (7.7.1) 实际上是公式 (7.9.3) 在 $p=\text{常数}=1/V$ 时的特例。

那么，抽样密度 $p$ 的最佳选择是什么？直观上，我们已经看到目标是使 $h=f/p$ 尽可能接近常数。我们可以通过关注公式 (7.9.3) 中平方根内的分子（即每个样本点的方差）来更严格地分析。两个尖括号本身都是积分的蒙特卡罗估计量，因此我们可以写成：

$$S\equiv ⟨\frac{f^2}{p^2}⟩-{⟨\frac{f}{p}⟩}^2\approx \int \frac{f^2}{p^2}pdV-{\left[\int \frac{f}{p}pdV\right]}^2=\int \frac{f^2}{p}dV-{\left[\int fdV\right]}^2$$

我们现在通过泛函变分法，在约束方程（7.9.2）下寻找最优的 $p$：

$$0=\frac{\delta }{\delta p}\left(\int \frac{f^2}{p}dV-{\left[\int fdV\right]}^2+\lambda \int pdV\right)$$

其中 $\lambda$ 是拉格朗日乘子。注意中间项不依赖于 $p$。变分（在积分号内进行）给出 $0=-f^2/p^2+\lambda$，即

$$p=\frac{|f|}{\sqrt{\lambda }}=\frac{|f|}{\int |f|dV}$$

这里 $\lambda$ 的取值是为了满足约束条件（7.9.2）。

如果在积分区域内 $f$ 保持同一符号，那么显然最优的 $p$ 选择（前提是能找到一种实际可行的抽样方法）应与 $|f|$ 成正比。此时方差降为零。不太明显但可以验证为真的一点是，即使 $f$ 同时取正负值，$p\propto |f|$ 仍然是最优的。在这种情况下，每个样本点的方差（由方程 7.9.4 和 7.9.6 给出）为

$$S=S_{\mathrm{optimal}}={\left(\int |f|dV\right)}^2-{\left(\int fdV\right)}^2$$

一个有趣的现象是，我们可以给被积函数加上一个常数，使其在整个区域保持同一符号，因为这样只会使积分改变一个已知量，即常数 $\times V$。此时，最优的 $p$ 选择总能使方差为零，从而得到完全精确的积分！这个看似悖论的问题（已在 §7.7 末尾提及）的解释在于：要完美地知道方程（7.9.6）中的 $p$，就需要完美地知道 $\int |f|dV$，而这实际上等价于已经知道了你试图计算的积分！

如果你的函数 $f$ 在大部分体积 $V$ 中取一个已知的常数值，那么加上一个常数使其变为零显然是个好主意。完成这一步后，重要性抽样所能达到的精度实际上并不取决于方程（7.9.7）有多小，而是取决于对于一个可实现的 $p$（通常只是理想情况的粗略近似），方程（7.9.4）有多小。

7.9.2 分层抽样

分层抽样的思想与重要性抽样截然不同。我们稍微扩展一下记号：令 $⟨⟨f⟩⟩$ 表示函数 $f$ 在体积 $V$ 上的真实平均值（即积分除以 $V$），而 $⟨f⟩$ 仍表示该平均值最简单的（均匀抽样）蒙特卡罗估计量：

$$⟨⟨f⟩⟩\equiv \frac{1}{V}\int fdV⟨f⟩\equiv \frac{1}{N}\sum _{i}f(x_i)$$

估计量的方差 $\mathrm{Var}(⟨f⟩)$ 衡量了蒙特卡罗积分误差的平方，它与函数本身的方差 $\mathrm{Var}(f)\equiv ⟨⟨f^2⟩⟩-⟨⟨f⟩{⟩}^2$ 之间满足如下渐近关系：

$$\mathrm{Var}\left(⟨f⟩\right)=\frac{\mathrm{Var}\left(f\right)}{N}$$

（参见方程 7.7.1）。

假设我们将体积 $V$ 分成两个相等且互不重叠的子体积，记为 $a$ 和 $b$，并在每个子体积中抽取 $N/2$ 个点。那么对于 $⟨⟨f⟩⟩$，除了方程（7.9.8）之外，还有另一个不同的估计量，记为 $⟨f{⟩}^{\prime }$：

$$⟨f{⟩}^{\prime }\equiv \frac{1}{2}\left(⟨f{⟩}_a+⟨f{⟩}_b\right)$$

换句话说，即两个半区域样本平均值的均值。该估计量（7.9.10）的方差为

$$\begin{array}{rl}& \mathrm{Var}\left(⟨f{⟩}^{\prime }\right)=\frac{1}{4}\left[\mathrm{Var}\left(⟨f{⟩}_a\right)+\mathrm{Var}\left(⟨f{⟩}_b\right)\right]\\ & =\frac{1}{4}\left[\frac{{\mathrm{Var}}_a\left(f\right)}{N/2}+\frac{{\mathrm{Var}}_b\left(f\right)}{N/2}\right]\\ & =\frac{1}{2N}\left[{\mathrm{Var}}_a\left(f\right)+{\mathrm{Var}}_b\left(f\right)\right]\end{array}$$

这里 ${\mathrm{Var}}_a\left(f\right)$ 表示 $f$ 在子区域 $a$ 中的方差，即 ${⟨⟨f^2⟩⟩}_a-{⟨⟨f⟩⟩}_a^2$，对 $b$ 同理。

根据已给出的定义，不难证明如下关系：

$$\mathrm{Var}\left(f\right)=\frac{1}{2}\left[{\mathrm{Var}}_a\left(f\right)+{\mathrm{Var}}_b\left(f\right)\right]+\frac{1}{4}{\left(⟨⟨f⟩{⟩}_a-⟨⟨f⟩{⟩}_b\right)}^2$$

（在物理学中，这个组合二阶矩的公式称为“平行轴定理”。）比较方程（7.9.9）、（7.9.11）和（7.9.12）可以看出，分层（分成两个子体积）抽样的方差永远不会大于简单蒙特卡罗方法的方差——并且当分层样本的均值 $⟨⟨f⟩{⟩}_a$ 和 $⟨⟨f⟩{⟩}_b$ 不同时，方差会更小。

我们尚未考虑在两个子体积中使用不同数量的样本点进行抽样的可能性，例如在子区域 $a$ 中抽取 $N_a$ 个点，在子区域 $b$ 中抽取 $N_b\equiv N-N_a$ 个点。现在我们来考虑这种情况。此时估计量的方差为

$$\mathrm{Var}\left(⟨f{⟩}^{\prime }\right)=\frac{1}{4}\left[\frac{{\mathrm{Var}}_a\left(f\right)}{N_a}+\frac{{\mathrm{Var}}_b\left(f\right)}{N-N_a}\right]$$

（很容易验证）当

$$\frac{N_a}{N}=\frac{{\sigma }_a}{{\sigma }_a+{\sigma }_b}$$

时，该方差取得最小值。

这里我们采用了简写记号 ${\sigma }_a\equiv [{\mathrm{Var}}_a(f){]}^{1/2}$，对 $b$ 同理。

如果 $N_a$ 满足方程（7.9.14），那么方程（7.9.13）可简化为

$$\mathrm{Var}\left(⟨f{⟩}^{\prime }\right)=\frac{({\sigma }_a+{\sigma }_b{)}^2}{4N}$$

若 $\mathrm{Var}(f)={\mathrm{Var}}_a(f)={\mathrm{Var}}_b(f)$，则式 (7.9.15) 退化为式 (7.9.9)，此时对样本进行分层不会产生任何差异。

将上述结果推广的一种标准方法是考虑将体积 $V$ 划分为多于两个的等大小子区域。可以很容易地得出结论：在各区域之间最优的样本点分配方式是使每个区域 $j$ 中的点数与 ${\sigma }_j$（即函数 $f$ 在该子区域中方差的平方根）成正比。然而，在高维空间中（例如 $d\gtrsim 4$），这种方法在实践中用处不大。因为在每个维度上将体积划分为 $K$ 段会导致 $K^d$ 个子体积，通常这个数量太大，以至于难以估计所有对应的 ${\sigma }_j$。

7.9.3 混合策略

重要性抽样和分层抽样乍看之下似乎相互矛盾。前者将样本点集中在被积函数 $|f|$ 最大的区域，而后者则集中在 $f$ 的方差最大的区域。两者怎么可能都是正确的呢？

答案是（就像生活中的许多事情一样），这完全取决于你掌握的信息以及掌握的程度。重要性抽样依赖于你已经对积分有了某种近似，从而能够按照所需的概率密度 $p$ 生成随机点 $x_i$。如果你的 $p$ 并不理想，那么误差只会以 $N^{-1/2}$ 的速率减小。如果在被积函数 $f$ 变化剧烈的区域中，你的 $p$ 与理想情况相差甚远，情况会尤其糟糕，因为此时被抽样的函数 $h=f/p$ 将具有很大的方差。重要性抽样通过平滑被抽样函数 $h$ 的取值来工作，其效果仅取决于你能否成功实现这种平滑。

相比之下，分层抽样不一定要求你对 $f$ 有任何了解。分层抽样通过平滑各子区域中样本点数量的波动来工作，而不是通过平滑点的取值。最简单的分层策略是将 $V$ 划分为 $N$ 个相等的子区域，并在每个子区域中随机选择一个点，这种方法的误差渐近地以 $N^{-1}$ 的速率减小，远快于 $N^{-1/2}$。（注意：拟随机数，见 §7.8，是另一种平滑点密度波动的方法，其效果几乎与这种“盲目”的分层策略一样好。）

然而，“渐近地”是一个重要的限制条件：例如，如果被积函数在除一个子区域外的所有区域中都可忽略不计，那么这种单点抽样的积分几乎毫无用处。此时，即使是非常粗略的信息，只要能通过重要性抽样将大量点集中在活跃子区域中，也会远优于盲目的分层抽样。

当你有办法估计各子区域的方差时，分层抽样才能真正发挥其优势，这样你就可以根据式 (7.9.14) 或其推广形式，在不同子区域中分配不等数量的样本点；同时，你还需要找到一种方法，将区域划分为合理数量的子区域（尤其是避免在高维 $d$ 下产生 $K^d$ 个子区域），并显著降低每个子区域中函数的方差（相比于整个体积中的方差）。要做到这一点，需要对 $f$ 有大量了解，尽管这种了解与重要性抽样所需的知识不同。

在实践中，重要性抽样和分层抽样并非互不相容。在许多（如果不是大多数）感兴趣的情形中，被积函数 $f$ 在体积 $V$ 的大部分区域都很小，仅在一小部分“活跃区域”中显著。在这些区域中，$|f|$ 的大小与标准差 $\sigma =[\text{Var}(f){]}^{1/2}$ 相当，因此两种技术都会导致样本点的大致相同集中。在更复杂的实现中，还可以将这两种技术“嵌套”使用，例如先在粗糙网格上进行重要性抽样，然后在每个网格单元内进行分层抽样。

7.9.4 自适应蒙特卡罗方法：VEGAS

由 Peter Lepage [1,2] 发明的 VEGAS 算法被广泛用于粒子物理中的多维积分。VEGAS 主要基于重要性抽样，但如果维度 $d$ 足够小（具体来说，当 $(K/2{)}^d<N/2$，其中 $N$ 为样本点数），它也会进行一些分层抽样以避免 $K^d$ 的爆炸式增长。VEGAS 中重要性抽样的基本技术是自适应地构造一个可分离的多维权重函数 $g$：

$$p\propto g(x,y,z,\dots )=g_x(x)g_y(y)g_z(z)\dots$$

这种函数通过两种方式避免了 $K^d$ 的爆炸：(i) 它可以存储为 $d$ 个独立的一维函数，每个函数由 $K$ 个表格值定义，因此存储量从 $K^d$ 减少到 $K\times d$；(ii) 它可以通过依次对 $d$ 个一维函数进行抽样，从而得到坐标向量分量 $(x,y,z,\dots )$。

可以证明，最优的可分离权重函数为 [1]：

$$g_x(x)\propto {\left[\int dy\int dz\dots \frac{f^2(x,y,z,\dots )}{g_y(y)g_z(z)\dots }\right]}^{1/2}$$

（对 $y,z,\dots$ 有类似表达式）。注意，在一维情况下，这退化为 $g\propto |f|$（式 7.9.6）。式 (7.9.17) 立即提示了 VEGAS 的自适应策略：给定一组 $g$ 函数（初始时可设为常数），对函数 $f$ 进行抽样，不仅累积积分的整体估计值，还累积式 (7.9.17) 右侧的 $Kd$ 个估计值（即在 $d$ 个维度上每个维度划分为 $K$ 个子区间）。这些估计值随后用于确定下一次迭代中改进的 $g$ 函数。

当被积函数 $f$ 在 $d$ 维空间中集中于一个或少数几个区域时，权重函数 $g$ 会迅速在这些区域向坐标轴投影的位置上变得很大。此时，蒙特卡罗积分的精度将远高于简单蒙特卡罗方法。

VEGAS 的弱点显而易见：如果函数 $f$ 在各个坐标方向上的投影是均匀的，那么 VEGAS 在这些维度上就不会集中样本点。VEGAS 最糟糕的情况是被积函数集中在接近体对角线的区域，例如从 (0, 0, 0, ...) 到 (1, 1, 1, ...) 的直线附近。由于这种几何结构完全不可分离，VEGAS 将毫无优势。更一般地，当被积函数集中在一维（或更高维）的弯曲轨迹（或超曲面）上时，除非这些轨迹恰好接近坐标方向，否则 VEGAS 的表现可能不佳。

下面的 vegas 例程本质上是 Lepage 的标准版本，仅做了最小修改以符合我们的约定。（我们感谢 Lepage 允许在此处重现该程序。）为了与其他版本的 VEGAS 算法保持一致，我们保留了原始变量名。参数 NDMX 即我们所说的 $K$，表示每个轴上的最大分段数；MXDIM 是 $d$ 的最大值；其他一些参数在注释中有说明。

vegas 例程执行 $m=\text{itmx}$ 次统计独立的积分评估，每次使用 $N=\text{ncall}$ 次函数求值。尽管这些迭代在统计上是独立的，但它们会相互辅助，因为每次迭代的结果都用于优化下一次迭代的抽样网格。所有迭代的结果通过以下关系合并为一个最佳答案及其估计误差：

$$I_{\mathrm{best}}=\sum _{i=0}^{m-1}\frac{I_i}{{\sigma }_i^2}/\sum _{i=0}^{m-1}\frac{1}{{\sigma }_i^2}{\sigma }_{\mathrm{best}}={\left(\sum _{i=0}^{m-1}\frac{1}{{\sigma }_i^2}\right)}^{-1/2}$$

此外，还返回以下量：

$${\chi }^2/m\equiv \frac{1}{m-1}\sum _{i=0}^{m-1}\frac{(I_i-I_{\mathrm{best}}{)}^2}{{\sigma }_i^2}$$

如果该值显著大于 1，则表明各次迭代的结果在统计上不一致，所得答案值得怀疑。

以下是 vegas 的接口。（完整代码见 [3]。）

void vegas(VecDoub_I &regn, Doub fxn(VecDoub_I &), const Int init, const Int ncall, const Int itmx, const Int nprn, Doub &tgral, Doub &sd, Doub &chi2a)

对用户提供的 ndim 维函数 fxn 在由 regn[0..2*ndim-1] 指定的矩形体积上执行蒙特卡罗积分，其中 regn 是一个向量，前 ndim 个元素为区域的“左下”坐标，后 ndim 个元素为“右上”坐标。积分包含 itmx 次迭代，每次迭代大约调用 ncall 次函数。每次迭代后都会优化网格；通常超过 5 或 10 次迭代并无益处。输入标志 init 指示本次调用是全新开始还是后续调用以进行更多迭代（参见代码中的注释）。输入标志 nprn（通常为 0）控制诊断输出的量。返回的结果包括 tgral（积分的最佳估计值）、sd（其标准差）和 chi2a（自由度为 $m-1$ 的 ${\chi }^2$ 值，用于指示结果是否一致）。更多细节见正文。

输入标志 init 可以巧妙使用。例如，可以先调用 init=0, ncall=1000, itmx=5，紧接着调用 init=1, ncall=100000, itmx=1。这样做的效果是：先用少量样本进行五次迭代以生成抽样网格，然后在优化后的网格上进行一次高精度积分。

要将缪vas用于§7.7中讨论的环面（torus）示例（仅考虑密度被积函数），函数 $f_{xn}$ 应为：

Doub torusfunc(const VecDoub &x, const Doub wgt) {
    Doub den = exp(5.*x[2]);
    if (SQR(x[2]) + SQR(sqrt(SQR(x[0]) + SQR(x[1])) - 3.) <= 1.) return den;
    else return 0.;
}

主程序代码则为：

Doub tgral, sd, chi2a;
VecDoub regn(6);
regn[0] = 1.; regn[3] = 4.;
regn[1] = -3.; regn[4] = 4.;
regn[2] = -1.; regn[5] = 1.; 
vegas(regn, torusfunc, 0, 10000, 10, 0, tgral, sd, chi2a);
vegas(regn, torusfunc, 1, 900000, 1, 0, tgral, sd, chi2a);

注意，用户提供的被积函数 fxn 除了预期的求值点 x 外，还有一个参数 wgt。在大多数应用中，你可以在函数内部忽略 wgt。然而，有时你可能希望在对主函数 $f$ 积分的同时，也对某些附加函数 $g$ 进行积分。任何此类函数 $g$ 的积分可由下式估计：

$$I_g=\sum _i{w}_{i}g(\mathbf{x})$$

其中 $w_i$ 和 $\mathbf{x}$ 分别对应参数 wgt 和 x。在你的函数 fxn 内部累积该求和项，并通过全局变量将结果传回主程序是直接可行的。当然，$g(\mathbf{x})$ 最好在某种程度上与主函数 $f$ 相似，因为采样是针对 $f$ 进行优化的。

缪vas 的完整代码见网络附注 [3]。

7.9.5 递归分层抽样

我们已经看到，分层抽样的问题在于它可能无法避免在 $d$ 维体积的直观笛卡尔剖分中固有的 $K^d$ 爆炸问题。一种称为递归分层抽样（recursive stratified sampling）[4] 的技术试图通过逐次对体积进行二分来解决这一问题，但并非沿所有 $d$ 个维度同时进行，而是每次仅沿一个维度进行。

起点是公式 (7.9.10) 和 (7.9.13)，它们被应用于对越来越小的子区域进行逐次二分。

假设我们有 $N$ 次函数 $f$ 的求值配额，并希望在矩形平行六面体区域 $R=({\mathbf{x}}_a,{\mathbf{x}}_b)$ 中计算 $⟨f{⟩}^{\prime }$。（我们用其对角顶点的两个坐标向量来表示这样一个区域。）首先，我们将 $N$ 的一个比例 $p$ 分配用于探索 $R$ 中 $f$ 的方差：我们在 $R$ 中均匀采样 $pN$ 个函数值，并累积求和，以得到 $d$ 个不同坐标方向上对应的 $d$ 对方差估计值（这些方向对应于 $R$ 可能被二分的方向）。换句话说，在 $pN$ 个样本中，我们估计在 $R$ 可能被二分后各子区域中的 $\mathrm{Var}(f)$：

$$\begin{array}{r}{R}_{ai}\equiv ({\mathbf{x}}_a,{\mathbf{x}}_b-\frac{1}{2}{\mathbf{e}}_i\cdot ({\mathbf{x}}_b-{\mathbf{x}}_a){\mathbf{e}}_i)\\ R_{bi}\equiv ({\mathbf{x}}_a+\frac{1}{2}{\mathbf{e}}_i\cdot ({\mathbf{x}}_b-{\mathbf{x}}_a){\mathbf{e}}_i,{\mathbf{x}}_b)\end{array}$$

此处 ${\mathbf{e}}_i$ 是第 $i$ 个坐标方向的单位向量，$i=1,2,\dots ,d$。

其次，我们检查这些方差，找出最有利于进行二分的维度 $i$。例如，根据公式 (7.9.15)，我们可以选择使子区域 $R_{ai}$ 和 $R_{bi}$ 中方差估计值的平方根之和最小的那个 $i$。（实际上，如我们将解释的，我们做了一点不同的处理。）

第三，我们将剩余的 $(1-p)N$ 次函数求值分配给区域 $R_{ai}$ 和 $R_{bi}$。如果我们使用公式 (7.9.15) 来选择 $i$，则应根据公式 (7.9.14) 进行该分配。

现在我们有两个平行六面体，每个都有其自己的函数求值配额，用于估计 $f$ 的均值。我们的“RSS”算法此时显现出其递归性：为了计算每个区域中的均值，我们回到公式 (7.9.21) 上方段落中以“首先，...”开头的句子。（当然，当分配给某个区域的点数低于某个阈值时，我们就改用简单的蒙特卡罗方法，而不再继续递归。）

最后，我们利用公式 (7.9.10) 和公式 (7.9.11) 的第一行，将两个子体积的均值和估计方差结合起来。

至此，RSS 算法的最简形式就完成了。在描述一些统称为“实现细节”的额外技巧之前，我们需要简要回到公式 (7.9.13)–(7.9.15)，并推导出我们实际使用的替代公式。

公式 (7.9.13) 的右侧两次应用了熟悉的公式 (7.9.9) 的缩放律，一次用于 $a$，一次用于 $b$。如果估计值 $⟨f{⟩}_a$ 和 $⟨f{⟩}_b$ 都是通过简单蒙特卡罗方法（使用均匀随机采样点）得到的，那么这是正确的。然而，这两个均值估计实际上是递归得到的。因此，没有理由期望公式 (7.9.9) 成立。相反，我们可以用以下关系替代公式 (7.9.13)：

$$\mathrm{Var}\left(⟨f{⟩}^{\prime }\right)=\frac{1}{4}\left[\frac{{\mathrm{Var}}_a\left(f\right)}{N_a^{\alpha }}+\frac{{\mathrm{Var}}_b\left(f\right)}{(N-N_a{)}^{\alpha }}\right]$$

其中 $\alpha$ 是一个未知常数且 $\alpha \ge 1$（等号情况对应于简单蒙特卡罗）。在这种情况下，一个简短的计算表明，当

$$\frac{N_a}{N}=\frac{{\mathrm{Var}}_a{\left(f\right)}^{1/\left(1+\alpha \right)}}{{\mathrm{Var}}_a{\left(f\right)}^{1/\left(1+\alpha \right)}+{\mathrm{Var}}_b{\left(f\right)}^{1/\left(1+\alpha \right)}}$$

时，$\mathrm{Var}\left(⟨f{⟩}^{\prime }\right)$ 达到最小值，且其最小值为

$$\mathrm{Var}\left(⟨f{⟩}^{\prime }\right)\propto {\left[{\mathrm{Var}}_a{\left(f\right)}^{1/\left(1+\alpha \right)}+{\mathrm{Var}}_b{\left(f\right)}^{1/\left(1+\alpha \right)}\right]}^{1+\alpha }$$

当 $\alpha =1$ 时，公式 (7.9.22)–(7.9.24) 退化为公式 (7.9.13)–(7.9.15)。通过数值实验寻找 $\alpha$ 的自洽值，发现 $\alpha \approx 2$。也就是说，当递归地使用 $\alpha =2$ 的公式 (7.9.23) 来分配采样机会时，RSS 算法的观测方差大致按 $N^{-2}$ 衰减，而公式 (7.9.23) 中使用其他 $\alpha$ 值则会导致更差的衰减。（不过，对方差对 $\alpha$ 的敏感性并不很大；目前尚不清楚 $\alpha =2$ 是否是一个可解析证明的结果，还是仅仅一个有用的启发式规则。）

miser 实现与上述算法的主要区别在于如何估计公式 (7.9.23) 右侧的方差。我们通过经验发现，使用采样函数值的最大值与最小值之差的平方，而不是样本的真实二阶矩，会更加稳健。当然，随着样本量增大，该估计量的偏差会越来越大；然而，公式 (7.9.23) 仅用于比较具有大致相同样本数的两个子体积（$a$ 和 $b$）。当初步采样在被积函数的活跃区域仅得到一个点或少量点时，“最大值减最小值”估计量显示出其价值。在许多实际情况下，这些点预示着附近存在更重要的区域，而“最大值减最小值”所提供的更大采样权重有助于更好地捕捉这些区域。

代码中包含的第二个修改是引入了一个“抖动参数”（dithering parameter）dith，当其非零时，子体积的划分不再严格位于中点，而是划分为 $0.5\pm$ dith 的比例，其中 $\pm$ 的符号由内置的随机数例程随机选择。通常可将 dith 设为零。然而，如果被积函数的某种特殊对称性使其活跃区域恰好位于区域的中点，或位于区域某个 $2$ 的幂次子区域的中心，则将 dith 设为非零值会带来很大优势。我们希望避免极端情况，即活跃区域被均匀地划分到 $d$ 维空间中 $2^d$ 个相邻的角落中。在有用的情况下，dith 的典型非零值可能是 0.1。当然，当抖动参数非零时，我们必须考虑子体积大小的差异；代码通过变量 fract 来处理这一点。

代码中还有一个最终特性值得一提。RSS 算法使用同一组采样点在所有 $d$ 个方向上计算公式 (7.9.23)。在递归的底层，采样点的数量可能非常少。虽然罕见，但可能发生某个方向上所有样本都位于体积的一半中；在这种情况下，该方向将被忽略，不作为二分的候选方向。更罕见的情况是，所有方向上的样本都位于体积的一半中。此时，将随机选择一个方向。如果你的应用中这种情况发生过于频繁，则应增加 MNPT（参见代码中的 if (jb == -1)... 行）。

注意，如所给出的 miser 返回的是函数平均值 $⟨⟨f⟩⟩$ 的估计值 ave，而非 $f$ 在区域上的积分。而缪vas 例程采用另一种约定，返回积分值 tgral。这两种约定当然通过公式 (7.9.8) 简单关联，因为矩形区域的体积 $V$ 是已知的。

miser 例程的接口如下：

蒙特卡罗方法在由 regn[0..2*ndim-1] 指定的矩形体积内对用户提供的 ndim 维函数 func 进行采样，该向量由区域的 ndim 个“左下”坐标和随后的 ndim 个“右上”坐标组成。函数总共被采样 npts 次，采样位置由递归分层抽样方法确定。函数在区域内的平均值以 ave 返回；ave 的统计不确定性估计值（标准差的平方）以 var 返回。输入参数 dith 通常应设为零，但如果 func 的活跃区域恰好落在区域的 $2$ 的幂次子区域的边界上，则可设为（例如）0.1。

在 §7.7 中实现环面问题的代码如下：

Doub torusfunc(const VecDoub &x) {
    Doub den = exp(5.*x[2]);
    if (SQR(x[2]) + SQR(SQR(x[0]) + SQR(x[1])) - 3.) <= 1.) return den;
    else return 0;
}

主程序代码为：

Doub ave, var, tgral, sd, vol = 3.*7.*2.;
regn[0] = 1.; regn[3] = 4.;
regn[1] = -3.; regn[4] = 4.;
regn[2] = -1.; regn[5] = 1.;
miser(torusfunc, regn, 1000000, 0., ave, var);
tgral = ave * vol;
sd = sqrt(var) * vol;

（实际上，miser 并不太适合解决此问题。）

miser 的完整代码见网络附注 [5]。miser 例程调用一个简短的函数 ranpt，以在指定的 $d$ 维区域内生成一个随机点。网络附注中的 ranpt 版本通过连续调用均匀随机数生成器并进行显而易见的缩放来实现。我们可以很容易地修改 ranpt，使其通过拟随机序列例程 sobseq（§7.8）生成点。我们发现，使用 sobseq 的 miser 比使用均匀随机数的 miser 精度显著更高。然而，由于 RSS（递归分层抽样）与拟随机数的使用是完全可分离的，因此此处给出的代码并未依赖于 sobseq。类似地，重要性抽样原则上也可与 RSS 结合使用。（原则上也可以将 VEGAS 与 MISER 结合，尽管编程会相当复杂。））

参考文献

Hammersley, J.M. and Handscomb, D.C. 1964, Monte Carlo Methods (London: METHuen).

Kalos, M.H. and Whitlock, P.A. 1986, Monte Carlo Methods (New York: Wiley).

Bratley, P., Fox, B.L., and Schrage, E.L. 1983, A Guide to Simulation, 2nd ed. (New York: Springer).

Lepage, G.P. 1978, "A New Algorithm for Adaptive Multidimensional Integration," Journal of Computational Physics, vol. 27, pp. 192–203.[1]

Lepage, G.P. 1980, "VEGAS: An Adaptive Multidimensional Integration Program," Publication CLNS-80/447, Cornell University.[2]

Numerical Recipes Software 2007, "Complete VEGAS Code Listing," Numerical Recipes Webnote No. 9, at http://www.nr.com/webnotes?9 [3]

Press, W.H., and Farrar, G.R. 1990, "Recursive Stratified Sampling for Multidimensional Monte Carlo Integration," Computers in Physics, vol. 4, pp. 190–195.[4]

Numerical Recipes Software 2007, "Complete Miser Code Listing," Numerical Recipes Webnote No. 10, at http://www.nr.com/webnotes?10 [5]

本章解析PDF