XorShift Jump 101，第 2 部分：多项式算术
5.00/5 (6投票s)
为 XorShift 随机数生成器 (RNG) 逐步讲解向前/向后跳转过程
引言
这是我关于 xorshift 系列随机数生成器跳转过程教程的第二部分。第一部分阐述了本文的动机，并探讨了一种基于矩阵乘法的算法。这一次，我们将学习另一种使用多项式的方法。最后，我们将分析这两种算法的理论性能和实证性能。
请注意，本文的这两个部分共享相同的方程式和参考文献编号以及源代码示例。
背景
为方便读者，本节完全复制自第一部分。
本文假设读者熟悉 System.Random 和无限猴子定理 中解释的线性随机数生成器理论基础。与那篇文章不同，本文的演示代码是用 C++ 11 编写的，因此读者应能熟练掌握该语言的中级水平。
基于多项式的方法需要一些额外的代数基础知识，例如多项式的概念、其性质以及涉及多项式的运算。软件实现基于 NTL，这是一个由 Victor Shoup 开发的数论方法库。该库又需要 GNU GMP 和 GF2X。本文使用的版本分别为 NTL 11.0.0、GMP 6.1.2 和 GF2X 1.2。
虽然 NTL 文档提到可以使用 Visual C++ 构建此库，但推荐的环境应为 Unix 风格，例如 Linux 或 Cygwin 下的 G++（及其配套工具链）（我使用的是后者）。
本文的一些 C++ 代码是由一个用 Haskell 编写的程序自动生成的。源代码存档中包含了该生成器已生成好的输出，因此读者不必一定运行它。否则，就需要一个可用的 Haskell 系统，例如 Haskell 平台。不熟悉 Haskell 的读者还需要一些入门级指南。我使用的是 Hal Daumé III 编写的 《Yet Another Haskell Tutorial》。
单元测试基础架构使用了 Google Test 1.8.1。
基于多项式算术的方法
第一部分中描述的基于矩阵的跳转存在一个显著的缺点：效率不高。确实，如果矩阵很大，矩阵乘法运算会变得非常昂贵：朴素实现的复杂度为 \(O(N^3)\)，而可能的最优实现的复杂度下界为 \(\Omega(N^2)\)。
文献 [2] 中描述了一种基于多项式算术的更高效的跳转过程。
在本节中，我们将从该算法背后的理论开始，然后继续讨论向前和向后跳转的实现细节。由于这些算法几乎相同，并行描述它们而不是分开描述似乎是个好主意。
理论背景
让我们从特征多项式的定义开始。给定矩阵 \(A\)，多项式
被称为该矩阵的特征多项式。（请注意，维基百科中的定义使用了另一个表达式 \(\det(zI - A)\)，但在 \(\mathbb{F}_2\) 中加法和减法是相同的。）该多项式的一个基本性质是
多项式算术定义了一个多项式对另一个多项式取模的概念
我们将考虑一个特例
可以用相对于 \(k\) 的对数复杂度来计算。
利用上述关系，我们可以得出结论
并且，如果我们代入 \(A\) 作为参数并重新排列各项，
最后，由于 \(p(A) = 0\)，
让我们用这个等价关系来计算我们的跳转变换 \(T^kX\)
这也可以使用霍纳法则写成
注意这个公式的结构：首先，我们取状态向量 \(X\)，计算 \(\alpha_1 X\) 并将其乘以 \(T\)。然后我们再次加上我们的状态向量乘以 \(\alpha_2\)，得到一个新的状态向量，它又通过乘以 \(T\) 进行变换，依此类推。现在，某个状态向量乘以 \(T\) 意味着什么？是的，这意味着我们取某个状态下的 xorshift 随机数生成器并向前推进一步。我们已经知道一个快速完成这一步的简单算法，即 step_forward()。
这个算法会比基于矩阵的算法更好吗？ 首先，它在内存消耗方面是显而易见的赢家：我们的“转移多项式”可以存储在 \(N\) 位的内存中，而转移矩阵则需要 \(N\times N\) 位。这对于具有大状态向量的随机数生成器尤其重要，例如梅森旋转算法，其状态向量由 624 个 32 位字组成。在这种情况下，转移矩阵将占用 \((624*32)^2=398721024\) 位（约 47 兆字节）的内存！在我们 128 位 xorshift 的情况下，增益不会那么显著，但也很可观（\(128/8=16\) 字节对 \(128^2/8=2048\) 字节）。
关于计算速度，如果我们已经准备好了转移矩阵或多项式，渐进复杂度是相同的。确实，矩阵与向量的乘法需要将 \(N\) 位的矩阵行与相同长度的状态向量进行异或运算，并且该运算对输出的每一位重复进行，这使得我们得到的最终复杂度为 \(O(N^2)\)。对于基于多项式的转移，我们对每个非零系数（最坏情况下为 \(N\) 个）执行一步 RNG 算法和状态向量加法（\(O(N)\)）。所以，我们也会得到 \(O(N^2)\) 的总复杂度，但这些估计隐藏了常数因子，这可能使它们在我们的特定情况下有显著不同。我们将通过实验来观察会发生什么。
[1] 中提到了对此算法的另一个改进。矩阵 \(A\) 的特征多项式是唯一具有 \(p(A)=0\) 属性的多项式吗？答案是“否”。有一整类零化多项式，它们可以有各种各样的次数，而次数最小的那个被称为最小多项式。如果我们使用最小多项式而不是特征多项式，我们可以预期它的次数可能会小于 \(N\)，因此即使在最坏的情况下，我们计算新 RNG 状态所需的操作也会更少。
向前跳转
本节及之后讨论的代码基于 jump_ahead.tar.gz，这是梅森旋转算法作者提供的跳转算法参考实现。请注意，目前有该代码的更新版本，但我的实现使用的正是那个版本。
首先，让我们计算最小多项式。为此，我们将使用 NTL 库中的 MinPolySeq 函数。该函数接受由线性变换（例如 xorshift）产生的值序列和要计算的最小多项式的上界。我们知道最小多项式的次数不能大于特征多项式的次数，在我们的例子中即为 STATE_SIZE_EXP。序列长度必须至少是该上界的两倍。
此代码实现为函数 init_transition_polynomials。我们必须在首次使用基于多项式的转移算法之前显式调用它。这种方法在研究中足够好，但对于生产用途，我们最好预先计算所需的多项式，并将其系数硬编码到源代码中，作为一个静态字节数组，可以按需用于初始化模多项式。
执行转移的其他步骤是计算 \(x^k\) 对最小多项式取模，并将其应用于给定的 RNG 状态。这些函数的代码对于向前和向后跳转非常相似，因此将在下一节之后一并讨论，届时我们将了解向后跳转的实现细节。
向后跳转
实现向后跳转的一种可能方法是利用 xorshift 算法的回绕特性，即将向后 \(k\) 步表示为向前 \(P-k = 2^N-(k+1)\) 步，其中 \(P=2^N-1\) 是状态大小为 \(N\) 的 xorshift 的周期长度。我将使用另一种方法，并为反向序列计算所需的原语。换句话说，我将创建另一个 RNG 算法，它产生与 xorshift 相同的数字，但方向相反。
首先，让我们计算最小多项式。我们可以在不了解反向 xorshift 的情况下做到这一点。只需获取 xorshift 生成的序列，但在将其传递给 MinPolySeq 之前将其反转。
下一步会更棘手。将多项式应用于初始状态（\(\ref{eq:horner}\)）需要 RNG 的单步操作，在我们的例子中，是 xorshift 向后一步。我们能从 xorshift 的定义中推导出这个公式吗？
首先，让我们回顾一下 xorshift128 的定义
uint32_t t = x ^ (x << 11);
x = y; y = z; z = w;
w = w ^ (w >> 19) ^ (t ^ (t >> 8));
其中 x, ..., w 代表 128 位状态的不同 32 位块。我们可以看到，更新后状态的 x、y 和 z 分量等于初始状态的 y、z 和 w。因此，恢复初始状态的这三个分量是微不足道的。
要恢复 x，让我们使用异或运算的性质
正如你从上面的代码中看到的，w 分量的新值是使用其先前的值（w ^ (w >> 19)，我们简称为 w-子表达式）和某个中间表达式的值（t ^ (t >> 8)，相应地称为 t-子表达式）计算的，而后者又是从 x 计算得出的。
我们的第一步是恢复 t-子表达式的值。为此，我们可以使用更新后状态的 z-分量。实际上，这个 z 的值直接来自过去的 w。因此，如果我们计算 z ^ (z >> 19)，它将等于从过去状态计算出的 w ^ (w >> 19)，如果我们将其与 w 的新值进行异或，我们就会“剥离”掉 w-子表达式，得到纯粹的 t-子表达式的值。
uint32_t tt = w ^ (z ^ z >> 19); // (t ^ (t >> 8))
注意，这个 tt 不是 t 本身，而是关于 t 的一个表达式。我们能恢复 t 吗？
让我们考虑一个更简单的例子，用一个更短的位向量来理解基本思想。假设有一个 8 位量
例如，将该向量右移三位意味着如下操作
将这些值进行异或运算会得到
请注意，最高有效的三位保留了它们原始的值，我们可以从位向量中提取它们并右移三位
如果我们用这个位向量与 \(v \oplus (v \gg 3)\) 进行异或，我们会得到
可以简化为
就是这样！利用变换的结果，我们可以提取参数的某些位，然后恢复该参数的更多位。对最近恢复的位（\(v_4\) 和 \(v_3\)）进行同样的操作，我们可以更进一步，计算出 \(v\) 的剩余两位。
我们可以使用这种方法来恢复 t 的值
uint32_t t_3 = tt & 0xFF000000U;
uint32_t t_2 = (tt & 0x00FF0000U) ^ (t_3 >> 8);
uint32_t t_1 = (tt & 0x0000FF00U) ^ (t_2 >> 8);
uint32_t t_0 = (tt & 0x000000FFU) ^ (t_1 >> 8);
uint32_t t = t_3 | t_2 | t_1 | t_0; // x ^ (x << 11)
最后，是 x 的值
uint32_t x_10_00 = t & 0x000007FFU;
uint32_t x_21_11 = (t & 0x003FF800U) ^ (x_10_00 << 11);
uint32_t x_31_22 = (t & 0xFFC00000U) ^ (x_21_11 << 11);
// ...
x = x_31_22 | x_21_11 | x_10_00;
现在我们准备好实现两个方向的跳转过程了。
通用实现
让我们从初始化代码开始，它必须在任何其他跳转函数之前调用
NTL::GF2XModulus fwd_step_mod;
NTL::GF2XModulus bwd_step_mod;
void init_transition_polynomials()
{
state_t s;
init(s);
const size_t N = 2*STATE_SIZE_EXP;
NTL::vec_GF2 vf(NTL::INIT_SIZE, N);
NTL::vec_GF2 vb(NTL::INIT_SIZE, N);
for(long i = 0; i < N; i++)
{
step_forward(s);
vf[i] = vb[N - 1 - i] = s[3] & 0x01ul;
}
NTL::GF2X fwd_step_poly;
NTL::GF2X bwd_step_poly;
NTL::MinPolySeq(fwd_step_poly, vf, STATE_SIZE_EXP);
NTL::MinPolySeq(bwd_step_poly, vb, STATE_SIZE_EXP);
NTL::build(fwd_step_mod, fwd_step_poly);
NTL::build(bwd_step_mod, bwd_step_poly);
}
这个函数无疑是基于多项式的跳转实现中最神秘的部分。首先，它创建了一个 xorshift 随机数生成器实例，并生成一个长度为 STATE_SIZE_EXP (RNG 状态的位数) 两倍的随机序列。该序列中每个数的最低有效位被存储在两个向量 vf 和 vb 中，这两个向量逐项填充，但顺序不同。
然后，这些序列被传递给 NTL 函数 MinPolySeq，该函数计算与我们的 RNG 及其对应物（以相反方向生成相同数字）指定的线性变换相对应的最小多项式。这个 NTL 函数真是神奇：通过观察 RNG 输出中单个位的值，它就能提取出关于整个 RNG 的信息，足以构建其最小多项式。
最后，此函数初始化两个全局对象，fwd_step_mod 和 bwd_step_mod。这是 NTL 提供的一种优化。如果你需要多次将不同的多项式对同一个多项式取模后求幂，预处理模多项式、保存结果并反复重用它可能会更有效率。
void prepare_transition(tr_poly_t &tr_k, uint64_t k, Direction dir)
{
tr_k.dir = dir;
NTL::GF2X x(1, 1);
switch(dir)
{
case Direction::FWD:
NTL::PowerMod(tr_k.poly, x, k, fwd_step_mod);
break;
case Direction::BWD:
NTL::PowerMod(tr_k.poly, x, k, bwd_step_mod);
break;
}
}
void add_state(state_t &x, const state_t &y)
{
for(size_t i = 0; i < x.size(); i++)
x[i] = x[i] ^ y[i];
}
void horner(state_t &s, const tr_poly_t &tr)
{
state_t tmp_state = s;
tmp_state.fill(0);
int i = NTL::deg(tr.poly);
if(i > 0)
{
for( ; i > 0; i--)
{
if(NTL::coeff(tr.poly, i) != 0)
add_state(tmp_state, s);
if(tr.dir == Direction::FWD)
step_forward(tmp_state);
else
step_backward(tmp_state);
}
if(NTL::coeff(tr.poly, 0) != 0)
add_state(tmp_state, s);
}
s = tmp_state;
}
void do_transition(state_t &s, const tr_poly_t &tr)
{
horner(s, tr);
}
这些函数不言自明：prepare_transition 简单地用 \(x^k \bmod f(x)\) 填充 tr_k，其中 \(f(x)\) 根据 dir 的不同，是 fwd_step_mod 或 bwd_step_mod。最后一个函数 do_transition 将多项式 tr 应用于 RNG 状态 s，实现了霍纳法则。在这种情况下，dir 参数决定应执行哪个单个 RNG 步骤，向前还是向后。
性能对比
我之前已经提供了一些关于这两种方法性能的理论上的考虑。关于内存复杂度，没有什么可补充的，我们只能重复说基于多项式的方法是明确的赢家。因此，我们下面的分析将只关注时间复杂度。
我们之前的估计与 \(N\)（RNG 状态的大小）有关。如果你比较不同 RNG 的相同跳转算法，这些估计可能有用，但请谨慎使用：大 \(O\) 估计是渐近的，对于对应于状态大小的有限 \(N\) 值可能效果不佳。在我们的例子中，我们将两种方法应用于同一个 RNG，所以 \(N\) 是相同的。因此，我们是否可以得出结论说我们的跳转算法具有常数复杂度？不，我们当然还不能，因为我们还有一个参数需要考虑：跳转大小 \(k\)。
请注意，我们的跳转算法有两个独立的阶段：准备和转移。它们的复杂度可能对参数的依赖性不同。
对于基于矩阵的方法，将 \(k\) 加入到估计中更简单：该参数决定了准备阶段的复杂度，因为我们将矩阵求 \(k\) 次幂。在 \(N\) 为常数的情况下，求幂是在 \(O(\log k)\) 步内计算的。转移阶段甚至更好：它涉及一次矩阵到向量的乘法，这根本不依赖于 \(k\)，所以转移确实是*关于 \(k\)* 的 \(O(1)\)。
对于基于多项式的方法，分析不太清晰，因为我们使用了第三方库，所以我们不太了解其性能。从常识来看，将多项式求 \(k\) 次幂也可以通过 \(O(\log k)\) 的平方求幂方法实现，与矩阵类似。转移阶段更难估计，因为 \(k\) 与结果多项式的次数之间的关系不太明确，而结果多项式的次数最终决定了转移复杂度。幸运的是，这个次数有上界，即最小多项式的次数，而最小多项式的次数又受 \(N\) 的限制。在 \(N\) 为常数的情况下，我们可以再次得出结论，基于多项式的转移是关于 \(k\) 的 \(O(1)\)。
综合来看，我们可以认为这两种方法在渐近上是相等的。如果你在写一篇计算机科学的论文，这是个不错的结果，但对于软件开发来说还不够：在 \(k\) 的值有限的情况下，那些被大 \(O\) 忽略的隐藏因子可能会超过 \(k\) 的影响。为了估计这些因子，我们来做一个实验。
我将测量准备转移所需的时间以及对一些随机跳转长度应用转移所需的时间。考虑到我们关于对数跳转复杂度的结论，在每个数量级上探测 \(k\) 似乎是合理的。为了将其与随机性结合起来，我们可以实现一种分层抽样：鉴于 \(k\) 可以保存在 \(\left\lceil\log_2 k\right\rceil\) 位中，对于这个序列的每一位，我们都将其视为最高有效位（因此将其设置为 1），并随机填充其余的位。通过这种方式，我们将创建一个随机样本，其中每个数量级的数字都得到平等地表示。
此抽样过程在 xstime.cc 中实现。函数 generate_test_cases 为 msb_pow 内的每个层级生成 n_trials 个样本。此外，此函数还用作转移代码的最终健全性检查：对于每个生成的样本跳转大小，都会用基于矩阵和基于多项式的两种方法准备并执行向前跳转，以确保最终的 RNG 状态相同。
模板函数 run_test 很直接：它对每个指定的测试用例进行带时间测量的转移，并将结果以五列制表符分隔的格式输出到标准输出。这些列是：使用的算法名称、跳转方向、数量级（以 \(\left\lfloor\log_2 k\right\rfloor\) 表示），以及最后的准备时间和跳转时间（单位为微秒）。
我们测量的目标是获得运行时间的一个粗略估计，牺牲了准确性以换取代码的可移植性。因此，我使用了标准 C++ 库中的高分辨率时钟，而不是一些特定于系统的计时器或硬件计数器。类似地，xstime.cc 中没有添加任何特殊代码来控制进程优先级或线程亲和性等。
在我的笔记本电脑（i7-4702MQ，G++ 5.4.0）上运行测试程序后，我得到了一个数据样本，需要对其进行适当的后处理和可视化。让我们使用 R 编程语言来完成这项工作。
我的 R 会话记录如下。首先，让我们加载两个流行的数据处理和可视化库
library(dplyr)
library(ggplot2)
加载数据
timings <- read.table("timing.1.txt",
col.names = c("Alg", "Dir", "logK", "Prepare", "Jump"))
从数据中获取初步见解
timings %>% group_by(Alg, Dir) %>% summarise(Prepare = mean(Prepare), Jump = mean(Jump))
每种方法和跳转方向在准备和执行跳转上花费的平均时间
# A tibble: 5 x 4
# Groups: Alg [3]
Alg Dir Prepare Jump
<fct> <fct> <dbl> <dbl>
1 matr BWD 175080. 3.16
2 matr FWD 174040. 3.14
3 matr_ntl FWD 4202. 0.0569
4 poly BWD 9.97 1.42
5 poly FWD 9.30 1.04
我们的第一个见解是：在准备阶段，我们用“matr”表示的基于矩阵的转移方法比基于多项式的方法慢了 \(10000\) 倍以上！鉴于这种差异，转移时间之间 \(2\) 倍的差距看起来微不足道。
结果的第三行包含了一种我们尚未讨论的跳转方法，matr_ntl。在获得了其他方法的初步性能数据后，很明显我的矩阵算术实现确实效率极低。因此，自然而然地出现了一个问题：我们能做得更好吗？NTL 库包含了它自己的矩阵算术实现，我相信这应该是最先进的。所以，matr_ntl 是使用 NTL 的基于矩阵的跳转的另一种实现。由于我只对检查其在跳转准备阶段的性能感兴趣，所以该方法的 do_transition 函数是空的。另外，由于此代码并非用于计时实验之外，它被放在了 xstime.cc 文件中。
基于 NTL 的矩阵方法性能要好得多，比我朴素的实现快了 41 倍，但它仍然比多项式方法慢了几个数量级。
让我们绘制一些漂亮的图表来展示性能对 \(k\) 的依赖关系
sparse.labels <- function(breaks) { ifelse(1:length(breaks) %% 10 == 0, breaks, "") }
ggplot(data = timings %>% filter(Alg == "matr")) +
geom_boxplot(aes(x = as.factor(logK), y = Prepare, color = Dir)) +
xlab("logK") + scale_x_discrete(labels = sparse.labels)
ggplot(data = timings %>% filter(Alg == "matr_ntl")) +
geom_boxplot(aes(x = as.factor(logK), y = Prepare, color = Dir)) +
xlab("logK") + scale_x_discrete(labels = sparse.labels)
ggplot(data = timings %>% filter(Alg == "poly")) +
geom_boxplot(aes(x = as.factor(logK), y = Prepare, color = Dir)) +
xlab("logK") + scale_x_discrete(labels = sparse.labels)
观察图 1、2 和 3，我们可以看到跳转准备时间与 \(\log_2 k\) 之间存在相当明显的线性依赖关系。这证实了我们关于准备时间为 \(O(\log k)\) 的理论估计。
为了给我们的假设提供一些定性的证据，我们可以对我们的数据拟合一个线性模型。为了使文章更短，我只对向前跳转的测量数据进行拟合。
matr_prep.lm <- lm(Prepare ~ logK, data = timings %>% filter(Alg == "matr" & Dir == "FWD"))
summary(matr_prep.lm)
Call:
lm(formula = Prepare ~ logK, data = timings %>% filter(Alg ==
"matr" & Dir == "FWD"))
Residuals:
Min 1Q Median 3Q Max
-42474 -5425 367 6143 45880
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6308.19 841.64 7.495 2.31e-13 ***
logK 5324.81 23.23 229.207 < 2e-16 ***
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 10350 on 618 degrees of freedom
Multiple R-squared: 0.9884, Adjusted R-squared: 0.9884
F-statistic: 5.254e+04 on 1 and 618 DF, p-value: < 2.2e-16
matr_ntl_prep.lm <- lm(Prepare ~ logK, data = timings %>% filter(Alg == "matr_ntl" & Dir == "FWD"))
summary(matr_ntl_prep.lm)
Call:
lm(formula = Prepare ~ logK, data = timings %>% filter(Alg ==
"matr_ntl" & Dir == "FWD"))
Residuals:
Min 1Q Median 3Q Max
-1253.0 -152.3 8.1 147.2 1385.0
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -320.2577 22.0018 -14.56 <2e-16 ***
logK 143.5789 0.6073 236.42 <2e-16 ***
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 270.6 on 618 degrees of freedom
Multiple R-squared: 0.9891, Adjusted R-squared: 0.989
F-statistic: 5.589e+04 on 1 and 618 DF, p-value: < 2.2e-16
对于基于矩阵的跳转，准备过程与线性模型拟合得相当好：logK 变量是显著的，并且 \(R^2\) 分数非常接近 \(1\)。
poly_prep.lm <- lm(Prepare ~ logK, data = timings %>% filter(Alg == "poly" & Dir == "FWD"))
summary(poly_prep.lm)
Call:
lm(formula = Prepare ~ logK, data = timings %>% filter(Alg ==
"poly" & Dir == "FWD"))
Residuals:
Min 1Q Median 3Q Max
-1.7293 -0.6262 -0.3608 -0.0511 18.9861
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.519085 0.171505 3.027 0.00258 **
logK 0.278856 0.004734 58.905 < 2e-16 ***
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 2.109 on 618 degrees of freedom
Multiple R-squared: 0.8488, Adjusted R-squared: 0.8486
F-statistic: 3470 on 1 and 618 DF, p-value: < 2.2e-16
请注意，对于基于多项式的跳转，准备过程的拟合效果稍差：\(R^2\) 仅约为 \(84\%\)。这是怎么发生的呢？看图 3，你可以注意到离群值相对离主序列更远，因此我们的线性模型在解释数据变异性方面表现较差。
ggplot(data = timings %>% filter(Dir == "FWD")) +
geom_boxplot(aes(x = as.factor(logK), y = Jump, color = Alg)) +
xlab("logK") + scale_x_discrete(labels = sparse.labels)
ggplot(data = timings %>% filter(Dir == "BWD")) +
geom_boxplot(aes(x = as.factor(logK), y = Jump, color = Alg)) +
xlab("logK") + scale_x_discrete(labels = sparse.labels)
转移过程的线性模型
matr_jump.lm <- lm(Jump ~ logK, data = timings %>% filter(Alg == "matr" & Dir == "FWD"))
summary(matr_jump.lm)
Call:
lm(formula = Jump ~ logK, data = timings %>% filter(Alg == "matr" &
Dir == "FWD"))
Residuals:
Min 1Q Median 3Q Max
-0.4843 -0.3436 -0.3029 0.1088 15.0334
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.085267 0.102040 30.236 <2e-16 ***
logK 0.001768 0.002817 0.628 0.53
---
Signif. codes: 0 `***' 0.001 `*' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 1.255 on 618 degrees of freedom
Multiple R-squared: 0.0006374, Adjusted R-squared: -0.0009797
F-statistic: 0.3941 on 1 and 618 DF, p-value: 0.5304
poly_jump.lm <- lm(Jump ~ logK, data = timings %>% filter(Alg == "poly" & Dir == "FWD"))
summary(poly_jump.lm)
Call:
lm(formula = Jump ~ logK, data = timings %>% filter(Alg == "poly" &
Dir == "FWD"))
Residuals:
Min 1Q Median 3Q Max
-0.8818 -0.2126 -0.0505 0.2311 9.3183
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.844373 0.049360 17.106 < 2e-16 ***
logK 0.006241 0.001362 4.581 5.61e-06 ***
---
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 0.6071 on 618 degrees of freedom
Multiple R-squared: 0.03284, Adjusted R-squared: 0.03127
F-statistic: 20.98 on 1 and 618 DF, p-value: 5.607e-06
我们可以看到，对于基于矩阵的转移，模型证实了我们关于 \(O(1)\) 复杂度（相对于 \(k\)）的假设。然而，对于基于多项式的转移，模型表明存在一些对 \(k\) 的依赖，但系数非常小。同时请注意两种情况下 \(R^2\) 分数都非常小。这令人失望，但也可以解释。实际上，由于我们的回归线几乎是水平的，它无法很好地捕捉每个 \(\log k\) 值的测量变异性。这是有道理的，因为这种变异性不是由 \(\log k\) 的变化引起的，而是由性能测量的不确定性引起的。另外，请注意我们实际上是在建模对舍入值 \(\left\lfloor \log_2 k \right\rfloor\) 的依赖关系，丢失了一些关于 \(k\) 的信息，这反过来又导致了更多隐藏的不确定性。
结论
我们从本教程中学到了很多。首先，我们看到了如何根据 RNG 定义构建非平凡的转移矩阵。用 Haskell 进行符号转换展现了这种语言的强大之处：一个方便的 DSL 在十几行代码中就定义好了。之后，我们考虑了一种基于多项式算术的替代方法。看似不相关的数学对象能够一起使用，这真是太神奇了！最后，我们比较了这两种方法的性能。
从纯技术的角度来看，基于多项式的方法在内存和时间复杂度上都是明显的赢家。然而，基于矩阵的方法有一个重要的非技术优势：简单性。这种方法可以由一个主修数学或计算机科学的大二学生轻松地从头实现，而基于多项式的方法则需要更高级的技能。如果你因任何原因（无论是法律上还是技术上）无法使用像 NTL 这样的现成实现，这一点可能很重要。另外，如果你只需要对某个预定义的步长 \(k\) 进行跳转（例如，迭代随机子流），转移矩阵可以预先计算。就执行转移本身而言，两种方法的速度是相当的。当前基于矩阵的转移实现慢了大约两倍，但它可能可以通过使用矢量化和/或 SIMD 扩展来优化，以并行化矩阵行与向量的乘法，并通过在某些 CPU 架构上可用的特殊操作“population count”来加速最终的按位求和。
参考文献
为方便读者，本节完全复制自第一部分。
[1] Hiroshi Haramoto, Makoto Matsumoto, and Pierre L’Ecuyer. A fast jump ahead algorithm for linear recurrences in a polynomial space. In Proceedings of the 5th International Conference on Sequences and Their Applications, SETA ’08, pages 290–298, Berlin, Heidelberg, 2008. Springer-Verlag. (可在线获取。)
[2] Hiroshi Haramoto, Makoto Matsumoto, Takuji Nishimura, François Panneton, and Pierre L’Ecuyer. Efficient jump ahead for \(\mathbb{F}_2\)-linear random number generators. INFORMS Journal on Computing, 20(3):385–390, 2008. (可在线获取。)
[3] Donald E. Knuth. 计算机程序设计艺术，第2卷（第3版）：半数值算法。 Addison-Wesley Longman Publishing Co., Inc., Boston, MA, USA, 1997.
[4] P. L’Ecuyer and R. Simard. TestU01: A C library for empirical testing of random number generators. ACM Transactions on Mathematical Software, 33(4): Article 22, August 2007.
[5] George Marsaglia. Xorshift RNGs. Journal of Statistical Software, Articles, 8(14):1–6, 2003. (可在线获取。)