随机数产生器¶

Monte Carlo 的输入看起来是哈密顿量、温度和系统尺寸，但程序真正反复调用的底层基础设施是随机数产生器。一次 Metropolis 接受或拒绝、一次 cluster bond 是否激活、一个随机自旋方向、一个随机初始构型、一次 bootstrap 重采样，都依赖随机数。

因此随机数产生器不是附属工具，而是 Monte Carlo 程序的基础部件。一个差的随机数产生器可能不会让程序崩溃，却会悄悄改变统计误差、引入相关性，甚至让临界点和指数估计产生系统偏差。

伪随机数¶

计算机通常不能直接产生真正的随机序列。Monte Carlo 程序使用的是伪随机数产生器：

\[ x_0\to x_1\to x_2\to\cdots. \]

给定初始种子 seed 后，整个序列完全确定。因此伪随机数有两个重要特点。

第一，它必须看起来像独立同分布随机变量。例如均匀随机数应满足

\[ u_i\in[0,1), \qquad u_i\sim {\rm Uniform}(0,1). \]

第二，它必须可复现。只要记录 seed、程序版本和参数，就应能重新得到同一条 Markov chain。这对调试、复现实验和定位异常 seed 很重要。

所以 Monte Carlo 中的随机数要同时满足两个要求：

对统计估计来说，序列要足够像独立随机样本。
对程序调试来说，序列又必须是确定可复现的。

从 seed 到随机数¶

程序里产生一个均匀随机数通常分成三步。

第一步，用 seed 初始化 RNG 的内部状态。这个状态可能是一个整数，也可能是一组很长的数组。抽象地写就是

\[ {\rm seed} \longrightarrow {\rm state}_0. \]

同一个 seed 必须给出同一个初始状态，所以同一份代码、同一组参数和同一个 seed 可以复现同一条随机序列。

第二步，用确定性递推更新内部状态：

\[ {\rm state}_{n+1}=F({\rm state}_n). \]

这里的 \(F\) 由具体 RNG 决定。线性同余、Mersenne Twister、PCG、xoshiro、lagged XOR 等方法，本质上都是不同的状态更新函数。虽然递推完全确定，但若 \(F\) 设计得好，输出序列在统计测试中会像独立随机数。

第三步，从内部状态提取整数，再映射到区间 \([0,1)\)。例如一个 32 位无符号整数 \(r_n\) 满足

\[ 0\le r_n<2^{32}, \]

就可以转换成

\[ u_n=\frac{r_n}{2^{32}}, \qquad 0\le u_n<1. \]

实际库函数会更仔细地处理整数范围、浮点精度和端点问题。对使用者来说，最重要的是理解这条链：

\[ {\rm seed} \to {\rm internal\ state} \to {\rm integer\ random\ bits} \to u\in[0,1). \]

不要在每次需要随机数时重新 seed。正确做法是在程序开始或每条独立 Markov chain 开始时初始化一次，然后连续调用 RNG，让状态自然向前推进。

常见接口¶

一个实际模拟通常至少需要以下随机数接口：

接口	用途
`uniform01()`	产生 \([0,1)\) 上的均匀随机数
`randint(n)`	产生 \(0,1,\cdots,n-1\) 的整数
`uniform(a,b)`	产生 \([a,b)\) 上的均匀实数
`gaussian()`	产生标准正态随机数
`random_unit_vector(N)`	产生 \(N\) 维单位球面上的均匀方向

物理算法中的很多随机步骤都可以由这些接口组合出来。例如 Metropolis 中接受率为 \(p\)，就比较

\[ u<p,\qquad u\sim{\rm Uniform}(0,1). \]

cluster 更新中一条 bond 的激活概率为 \(p_{ij}\)，也是同样比较

\[ u<p_{ij}. \]

传统自定义生成器的做法¶

许多 Monte Carlo 程序会先封装一层随机数接口，例如：

int getRandomNum();
int getRandomNum(int MaxNum);
long getRandomNum(long MaxNum);
double getRandomDouble();
double getRandomDouble(double a, double b);
double getGaussian();

一种传统做法是使用两个 lagged XOR 序列混合产生整数随机数。设两组状态数组参数为：

\[ {\rm len}_1=9689,\qquad {\rm ifd}_1=471, \]

和

\[ {\rm len}_2=127,\qquad {\rm ifd}_2=30. \]

每次刷新缓冲区时，可按下面方式生成两个滞后异或值：

\[ \ell=ir_1[i]\oplus ir_1[i+{\rm ifd}_1], \]

\[ k=ir_2[j]\oplus ir_2[j+{\rm ifd}_2], \]

再组合成

\[ {\rm rand}=k\oplus \ell. \]

这里 \(\oplus\) 是按位异或。生成出的整数再映射到 \([0,1)\)：

\[ u={\rm rand}\times 2^{-32}+\frac12. \]

这种写法属于传统数值模拟中常见的自定义伪随机数方案。它的优点是接口简单、可复现、速度快，也方便在旧代码中移植。缺点是需要自己负责质量验证、并行 seed 设计和长期维护。

生成正态随机数时，可以使用 Box-Muller 变换：

\[ z = \sqrt{-2\ln(1-u_1)} \cos(2\pi u_2), \]

其中 \(u_1,u_2\) 是独立均匀随机数。这里写成 \(1-u_1\) 是为了避免 \(u_1=0\) 时出现 \(\ln0\)。

C++ 标准库的选择¶

现代 C++ 中更常用 <random>。例如 Mersenne Twister：

#include <random>

std::mt19937 rng(seed);
std::uniform_real_distribution<double> uniform01(0.0, 1.0);
std::normal_distribution<double> normal(0.0, 1.0);

double u = uniform01(rng);
double z = normal(rng);

常见选择包括：

产生器	特点
`std::mt19937`	32 位 Mersenne Twister，周期长，经典默认选择
`std::mt19937_64`	64 位版本，适合需要 64 位随机整数的场景
`std::minstd_rand`	简单线性同余，速度快但质量较弱，不建议做严肃大规模模拟
`std::random_device`	用于获取非确定性 seed，不适合作为主循环随机数引擎

mt19937 的名字来自 Mersenne Twister 的周期长度：

\[ 2^{19937}-1. \]

它不是加密安全随机数，但对许多普通 Monte Carlo 已经足够可靠。严肃的大规模模拟仍应做随机数质量检查，尤其是并行、多 seed、临界慢化或高精度指数估计的任务。

播种与复现¶

随机数程序最常见的错误不是公式写错，而是 seed 管理混乱。

一个推荐做法是：

每个独立任务使用明确记录的 seed。
输出文件中保存 seed、系统尺寸、温度、热化步数和测量步数。
多 seed 并行时，不要简单使用连续 seed 而不做检查。
不要在内层循环中反复重新初始化 RNG。
调试时固定 seed；正式统计时使用多个独立 seed。

在 C++ 中，如果希望从一个整数 seed 扩展出较好的初始状态，可以使用 std::seed_seq：

std::seed_seq seq{seed, L, D, Nspin, run_id};
std::mt19937 rng(seq);

这样比简单地把相邻任务设成 seed=1,2,3,... 更稳一些。并行模拟中更严格的做法是使用支持跳跃或分流的 RNG，或者为每条链使用经过测试的独立 seed 方案。

整数抽样与 modulo bias¶

若要从 \(0,\cdots,n-1\) 中均匀抽一个整数，最危险的写法是

int x = rng() % n;

如果 RNG 的整数范围不是 \(n\) 的整数倍，就会产生 modulo bias。偏差通常很小，但在长时间模拟和离散事件抽样中没有必要冒这个风险。

标准库推荐使用

std::uniform_int_distribution<int> dist(0, n - 1);
int x = dist(rng);

若自己实现，可以用拒绝采样：只接受落在最大可整除区间内的随机整数，超出部分重抽。

从均匀数生成常见分布¶

大多数 Monte Carlo 程序先产生均匀随机数，再变换成其他分布。

Bernoulli 事件。 概率为 \(p\) 的事件：

\[ {\rm event}= (u<p). \]

离散分布。 若概率为 \(p_1,\cdots,p_m\)，先构造累计概率

\[ P_k=\sum_{i=1}^k p_i, \]

然后找第一个满足

\[ u<P_k \]

的 \(k\)。若 \(m\) 很大，可用二分搜索、alias method 或树状数组。

指数分布。 若事件等待时间满足

\[ P(\tau)=\lambda e^{-\lambda\tau}, \]

则

\[ \tau=-\frac1\lambda\ln(1-u). \]

event-chain、连续时间 Monte Carlo 和 rejection-free 算法中常出现这种抽样。

高斯分布。 Box-Muller 变换为

\[ z_1=\sqrt{-2\ln u_1}\cos(2\pi u_2), \]

\[ z_2=\sqrt{-2\ln u_1}\sin(2\pi u_2). \]

实际程序中要避免 \(u_1=0\)。

单位球面方向。 对 \(O(N)\) 自旋，生成 \(N\) 个高斯数 \(g_a\)，再归一化：

\[ S_a=\frac{g_a}{\sqrt{\sum_b g_b^2}}. \]

这样得到的是 \(S^{N-1}\) 上的均匀方向。

随机数质量看什么¶

随机数质量不是一句“周期很长”就能保证。Monte Carlo 关心的是：用这些随机数驱动算法时，是否会产生可见偏差。

可以从几个层次检查。

一维均匀性¶

最基本要求是 \(u_i\) 在 \([0,1)\) 上均匀。可以把区间分成 \(m\) 个 bin，统计每个 bin 的频数 \(n_k\)。若总样本数为 \(N\)，期望频数为

\[ E_k=\frac{N}{m}. \]

常用 \(\chi^2\) 统计量：

\[ \chi^2 = \sum_{k=1}^m \frac{(n_k-E_k)^2}{E_k}. \]

如果 \(\chi^2\) 远大于自由度 \(m-1\)，说明分布可能不均匀；如果远小于 \(m-1\)，也可疑，因为真实随机数不应过于平滑。

也可以检查样本矩：

\[ \mathbb E[u]=\frac12, \qquad {\rm Var}(u)=\frac1{12}. \]

实际样本均值应满足

\[ \bar u-\frac12=O(N^{-1/2}). \]

相邻相关性¶

均匀性只看单个随机数，不检查序列相关性。Monte Carlo 更怕相邻随机数存在隐藏结构。可以计算自相关：

\[ C(k) = \frac1{N-k} \sum_{i=1}^{N-k} (u_i-\bar u)(u_{i+k}-\bar u). \]

归一化后

\[ \rho(k)=\frac{C(k)}{C(0)}. \]

理想独立序列中，\(\rho(k)\) 应在统计误差范围内围绕 \(0\) 波动。若某些固定 lag 有明显峰值，就可能影响 Markov chain。

还可以画二维散点：

\[ (u_i,u_{i+1}), \]

或三维散点：

\[ (u_i,u_{i+1},u_{i+2}). \]

差的线性同余产生器常会在这些图上显示格子结构。

低位比特¶

某些 RNG 的低位比特质量明显差。如果程序用

rng() % 2

这类操作，就可能主要使用低位比特。检查方法包括：

统计每一位上 \(0\) 和 \(1\) 的比例。
检查 bit 序列的 runs，即连续相同 bit 的长度分布。
检查 bit 自相关。

使用 uniform_int_distribution 或基于高质量高位的整数转换，通常比手写取模稳。

多维均匀性和谱结构¶

Monte Carlo 更新常连续使用多个随机数，例如一个随机数选站点，一个随机数选方向，一个随机数决定接受。于是需要检查多维向量

\[ (u_i,u_{i+1},\cdots,u_{i+s-1}) \]

是否在 \(s\) 维单位超立方中均匀分布。

简单方法是做二维、三维网格 \(\chi^2\) 检验。更系统的方法包括 spectral test、birthday spacings、gap test、serial test 等。实际项目中可以直接使用 TestU01、PractRand、Dieharder 这类测试套件，而不是自己从头实现所有测试。

物理量回归测试¶

随机数最终服务于物理模拟，所以还应做物理层面的测试。

常见方法包括：

用不同 RNG 重复同一小系统，比较能量、磁化率、Binder ratio 是否在误差范围内一致。
对很小系统做精确枚举，把 Monte Carlo 结果与精确值比较。
对已知临界点的模型，例如二维 Ising，检查交点和临界指数是否合理。
用多个 seed 独立运行，检查不同 seed 的结果是否服从误差估计。
在高温极限或低温极限检查简单解析结果。

这类测试很重要，因为某个 RNG 可能通过简单均匀性测试，却仍与特定算法产生相关性。例如 cluster 算法、长程 bond 抽样、并行 tempering 都可能放大特定模式的随机数缺陷。

Monte Carlo 中的常见陷阱¶

反复播种。 在每次 sweep 或每次函数调用时重新 seed，会破坏随机序列。RNG 应在程序开始时初始化一次，然后持续使用。

只保存最终平均值。 若没有保存 seed 和基本运行信息，发现异常时很难复现。

连续 seed 的独立性未经检查。 seed=1,2,3,... 不一定灾难，但不能把它自动当成严格独立流。多 seed 结果需要做一致性检查。

混用不同随机数接口。 同一个程序里一部分用自定义 RNG，一部分用 std::rand()，容易造成不可复现和质量不可控。

用全局 RNG 做并行。 多线程共享同一个 RNG 会产生竞争和不可复现。更稳的做法是每个线程或每条链有自己的 RNG 状态。

误把随机数误差当成热化误差。 如果不同 seed 的平均值系统性偏离，可能是热化不足，也可能是 RNG 或 seed 流有问题，需要分别排查。

实用建议¶

对一般研究代码，建议遵循以下原则。

默认使用经过广泛测试的 RNG，例如 std::mt19937_64、PCG、xoshiro 系列等。
自定义 RNG 可以保留，但要配套质量测试和物理回归测试。
所有输出必须记录 seed 和 RNG 类型。
离散整数抽样使用标准分布或拒绝采样，避免取模偏差。
并行任务要有明确的 seed 生成方案。
高精度临界指数、长程模型和大规模并行模拟，应使用专门测试套件检查 RNG。

随机数产生器的目标不是“看起来随机”，而是在目标模拟的尺度上不引入可见系统偏差。对 Monte Carlo 来说，这和 detailed balance、热化、自关联时间一样，是结果可信度的一部分。