XorShift 跳跃 101，第 1 部分：矩阵乘法
5.00/5 (5投票s)
XorShift RNG 的前后跳跃过程分步解释
引言
本文是System.Random 和无限猴子定理的续篇，我们探索了 .NET Framework 标准随机数生成器 (RNG) 的内部结构。利用线性 RNG 的理论，我们学会了如何以对数复杂度跳跃 \(N\) 步向前或向后导航随机数序列。这次，我们将为另一类线性 RNG（称为 xorshift）开发类似的功能。
我们将从构建此 RNG 的转换矩阵开始。与滞后斐波那契生成器不同，这些矩阵具有更复杂的结构，因此我们更喜欢自动计算它们。此计算由一个用强大但不流行的语言 Haskell 编写的程序执行。尝试 Haskell 是我这项研究的目标之一，事实证明它值得付出努力。之后，我们将考虑另一种基于多项式算术 [2] 的跳跃方法，该方法在上一篇文章中提到但尚未尝试。最后，我们将比较这两种方法。
将所有材料放在一篇文章中会使它有点太长，无法一次性阅读，所以我决定将其分成两部分。第 1 部分专门介绍基于矩阵的方法，而第 2 部分将考虑基于多项式的算法和性能测试。出于同样的原因，我将侧重于数学背景和一般思想，而不是软件实现。Haskell 中的代码生成器是个例外：假设这种语言不被广泛使用，我逐行解释了代码的每个部分，更详细地阐述了它们。
为什么我们需要再来一个 RNG 算法？不幸的是，System.Random 并非完美无缺。首先，统计测试暴露出滞后斐波那契生成器整个家族固有的问题 [3, §3.2.2]，[4, 第 7 章]。除此之外，.NET 实现中还发现了一个关键错误：代码中指定的滞后值不能保证最长可能的周期。事实上，我不知道这些滞后的周期长度的任何理论估计。这个问题已经被开源 .NET 社区多次讨论过 (#5974, #12746, #23298)，但为了向后兼容，代码保持原样。（另请参阅 Søren Fuglede Jørgensen 做的精彩研究。）
考虑到这些事实，如果您不仅需要随机数，还需要对其属性进行某种经过验证的估计，我不会推荐使用 System.Random。
现在自然会问一个问题：既然我们知道 System.Random 的缺点，为什么还要探索它呢？答案很简单：虽然它不适合实际使用，但这个 RNG 非常适合解释理论。事实上，滞后斐波那契算法只使用简单的整数算术，其转换矩阵足够清晰，无需复杂计算即可直接编写。此外，这个 RNG 的状态向量足够大，可以存储文本消息，这是我的上一篇文章得以实现的小技巧。
背景
假定读者熟悉System.Random 和无限猴子定理中解释的线性 RNG 理论基础。与那篇文章不同，本文的演示代码是用 C++ 11 编写的，因此读者应该熟练掌握这种方言的中级水平。
基于多项式的方法需要一些额外的代数基础知识，例如多项式的概念、其性质以及涉及多项式的运算。软件实现基于 Victor Shoup 开发的数论方法库 NTL。该库又需要 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 Platform。不熟悉 Haskell 的读者还需要一些入门级指南。我使用了 Hal Daumé III 撰写的 Yet Another Haskell Tutorial。
Google Test 1.8.1 用于单元测试基础结构。
入门
让我们从实现 xorshift 算法的基本 RNG 功能开始。
首先，我们应该选择该族中的某个特定算法用于我们的研究。让我们实现 George Marsaglia 在 [5, p. 5] 中定义的 xor128。
所有与基本 RNG 函数相关的代码都将放置在文件 xs.cc 中，并在头文件 xs.h 中对定义的实体进行前向声明。这段代码中做出了一个在研究编程中流行的重要设计权衡：为了实现最大的灵活性，我们不会遵循面向对象编程的封装原则，而是将 RNG 和其他对象的所有内部状态公开到公共范围。此外，我更喜欢使用“裸”数据结构，由非成员函数处理。同样，这是一个研究编程权衡：当您进行研究时，您经常需要定义一些关于数据的替代操作，并且在这些操作被使用的地方直接进行定义比在一些公共类定义中添加越来越多的类成员或为继承而烦恼要方便得多。
Marsaglia 将 xor128 的内部状态定义为四个独立的 32 位变量：x、y、z 和 w。由于我们将应用基于矩阵的方法，因此此状态应被视为单个 128 位向量。为了统一执行矩阵运算，将它们打包成 4 元素 32 位无符号整数向量会很方便，但有时我们将使用 x、y、z 和 w 作为相应向量项的别名，当它们被区别对待时。
using state_t = std::array<uint32_t, 4>;
此状态向量可以用任意值初始化，除了全零，但原始论文定义了一些“规范”初始种子，我们也将使用它。
void init(state_t &s)
{
s[0] = 123456789;
s[1] = 362436069;
s[2] = 521288629;
s[3] = 88675123;
}
下一步是实现从该状态生成随机数的函数。这个过程可以分为两个子问题：我们应该将 RNG 状态变异为下一个值，然后使用这个新状态生成下一个输出项。通常的做法是将这两个过程合并到一个代码块中（Marsaglia 的论文也遵循了这一点），但在我们的研究中，我们需要能够将这些操作分开。实际上，我们必须能够像往常一样执行下一步，或者使用某种跳转算法计算下一个状态，然后从该状态继续生成随机数，就像在前一次访问生成器时获得的那样。因此，我们将定义两个独立的函数
void step_forward(state_t &s)
{
uint32_t &x = s[0], &y = s[1], &z = s[2], &w = s[3];
uint32_t t = x ^ (x << 11);
x = y; y = z; z = w;
w = w ^ (w >> 19) ^ (t ^ (t >> 8));
}
uint32_t gen_u32(state_t &s)
{
step_forward(s);
return s[3];
}
基于矩阵的方法
定义 step_forward() 之后，现在是讨论 xorshift RNG 为何属于线性家族的好时机。部分 检查线性度 涉及一些数学知识，因此对编程更感兴趣的读者可以跳过它直到最终结论。下一节 计算转换矩阵 提供了如何以编程方式生成所需转换矩阵的分步解释。
检查线性度
我们的目标将是展示 xorshift 方法的单步操作可以表示为 \(S^{(i+1)} = TS^{(i)}\)，其中 \(T\) 是一个称为*转换矩阵*的常数矩阵，\(S^{(i)}\) 是在步骤 \(i\) 的状态向量。对于 xorshift，\(S\) 和 \(T\) 都定义在 \(\mathbb{F}_2\) 上：\(S \in \mathbb{F}_2^{N}\)，和 \(T \in \mathbb{F}_2^{N\times N}\)，其中 \(N\) 是状态向量的大小（本文其余部分 \(N=128\)）。用程序员的说法，这意味着这些向量或矩阵的项是单个位，允许的算术运算是乘法（按位“与”）和模 \(2\) 的加法（按位“异或”）。
为了简化我们的公式，让我们省略步骤索引，并用 \(S'=S^{(i+1)}\) 和 \(S=S^{(i)}\) 表示。那么我们的转换公式将是 \(S' = TS\)。此外，让我们使用带索引的小写字母表示向量和矩阵项，例如 \(s_i\) 或 \(t_{i,j}\)，\(0\le i < N\)，\(0\le j < N\)。
给定 \(S'=TS\)，\(S'\) 的每个分量为
让我们从更新状态向量中名为 x、y、z 的部分开始。它们覆盖状态向量的位 \(s_0\) 到 \(s_{95}\)，并通过直接复制其他 32 位块获得。
因此，我们转换矩阵的 \(75\%\) 行将是这样的
满足要求。
另一方面，块 w 的计算方式更为复杂：我们可以看到位运算符“异或”和位移操作的组合。为每个组件编写精确的公式将非常耗时，因此我们将采用归纳法。
- 公式最内层的元素是状态向量分量，因此它们与 (\(\ref{eq:component}\)) 的结构匹配。
- 给定匹配 (\(\ref{eq:component}\)) 的子表达式，我们可以应用模 2 加法并得到满足相同条件的和（这遵循加法的性质）。
- 给定一个匹配 (\(\ref{eq:component}\)) 的子表达式向量，我们可以看到移位操作要么用另一个分量的子表达式替换一个分量的子表达式，保持属性 (\(\ref{eq:component}\)) 不变，要么将其替换为零，这也符合该条件。
综上所述，我们可以得出结论，step_forward() 中定义的变换确实适用于我上一篇文章中描述的基于矩阵的跳转算法。
计算转换矩阵
如前所述，可以从实现 xorshift RNG 的代码中计算转换矩阵 \(T\)。为此，我们需要逐个组件地为代码中的每个操作编写公式，即为状态向量 \(S\) 的每个位。这既繁琐又容易出错。更糟糕的是，由于 xorshift 是一系列 RNG，我们必须从头开始为每种 xorshift 计算矩阵 \(T\)。解决方案很简单：让我们编写一个程序来自动完成这项工作。
该程序应获取描述给定 xorshift 类型的公式，并使用算术规则将这些公式转换为 (\(\ref{eq:component}\)) 形式。对状态向量 \(S\) 的每个 \(N\) 位执行此操作，我们将得到矩阵 \(T\) 的 \(N\times N\) 个项。
公式的符号操作需要合适的语言来实现。主要要处理的数据结构将是一个树状结构，其节点表示公式中的每个常量、变量或操作及其相互作用。函数式语言 Haskell 似乎是一个合适的工具，因为它对递归数据结构、代数数据类型和模式匹配提供了一流的支持。
此外，Haskell 提供了一个强大的功能，有助于进行实验和快速原型设计，即 REPL（Read-eval-print-loop）。使用 REPL，我们可以交互式地实现我们的矩阵生成器，逐个指定类型和函数定义，并立即从交互式 Haskell shell 中使用它们，以查看它们如何工作。
文章提供的源代码档案包含文件 xsTgen.hs。我将逐步解释它的工作原理。为此，我们将把这个文件加载到 Haskell 解释器中，如果您使用的是推荐的 Haskell Platform，则该解释器名为 *ghci*。之后，我们将在 REPL 中输入一些 Haskell 表达式并检查它们的输出。请注意，开发过程非常相似：我逐个编写函数定义并立即对其进行测试。（如果您使用 Haskell 感知编辑器，例如 Emacs，这个过程会容易得多。）当您进行探索性编程时，这种反馈非常有用。
要启动解释器，请打开命令行提示符窗口，将当前目录更改为包含我们的源文件 xsTgen.hs 的目录，然后运行 ghci。如果一切顺利，您应该会看到这样的启动消息
GHCi, version 8.6.3: http://www.haskell.org/ghc/ :? for help
Prelude>
使用命令 :l 加载源代码
Prelude> :l xsT
[1 of 1] Compiling Main ( xsTgen.hs, interpreted )
Ok, one module loaded.
*Main>
Windows 用户也可以启动 WinGHCi 作为 GUI 应用程序，并通过菜单加载源文件。
让我们从表示 xorshift 转换公式的数据类型开始我们的旅程
data Expr =
Zero
| Var String Integer
| Xor [Expr]
deriving (Eq, Show)
在这里，我们有一个带标签的联合，它可以存储常量零，或者由名称和索引引用的变量，或者一组任意表达式异或在一起。deriving 子句告诉 Haskell 自动为我们的数据类型创建一些辅助代码。
我们现在可以看到关于已定义数据类型的信息
*Main> :i Expr
data Expr = Zero | Var String Integer | Xor [Expr]
-- Defined at xsTgen.hs:1:1
instance [safe] Show Expr -- Defined at xsTgen.hs:5:17
instance [safe] Eq Expr -- Defined at xsTgen.hs:5:13
并输入一些表达式以查看其工作原理
*Main> Zero
Zero
*Main> Var "x" 1
Var "x" 1
*Main> Xor [Var "x" 1, Var "y" 2]
Xor [Var "x" 1,Var "y" 2]
*Main> Xor [Var "x" 1, Xor[Zero, Var "y" 2]]
Xor [Var "x" 1,Xor [Zero,Var "y" 2]]
这些表达式分别定义了诸如 \(0\), \(x_1\), \(x_1\oplus y_2\), \(x_1\oplus (0\oplus y_2)\) 等公式的内部表示。请注意，我们的 Expr 不足以表示完整的 \(\mathbb{F}_2\) 算术：我们没有常量 \(1\) 或乘法来写诸如 \(y_1\oplus 1\) 或 \(x_1\cdot y_2\) 之类的表达式。这有点限制，但我们不需要这种能力来解决我们当前的问题。
现在是时候教我们的代码处理表达式向量了，这样我们就不用为 RNG 状态的每个位编写公式了。如果我们查看 xsTgen.hs，数据类型 Expr 后面跟着这组定义
bus name cntr = map (Var name) [0 .. cntr - 1]
zeroes = Zero : zeroes
xor x y = zipWith (\x y -> Xor [x, y]) x y
crop x y = zipWith (\ x y -> x) x y
shl x n = crop (take n zeroes ++ x) x
shr x n = crop (drop n x ++ zeroes) x
这些方程定义了一些简单的函数。让我们看看它们是如何工作的。函数 bus 用于根据给定名称和大小创建索引变量的向量。
*Main> bus "x" 3
[Var "x" 0,Var "x" 1,Var "x" 2]
*Main> bus "y" 4
[Var "y" 0,Var "y" 1,Var "y" 2,Var "y" 3]
辅助函数 zeroes 生成一个无限的零常数序列。因此，我们无法直接评估它，但我们可以在需要时消耗一定数量的零。
*Main> take 5 zeroes
[Zero,Zero,Zero,Zero,Zero]
*Main> take 7 zeroes
[Zero,Zero,Zero,Zero,Zero,Zero,Zero]
函数 Xor 接收两个表达式向量并逐分量进行异或操作，生成一个新向量。
*Main> xor (bus "x" 3) (bus "y" 3)
[Xor [Var "x" 0,Var "y" 0],Xor [Var "x" 1,Var "y" 1],Xor [Var "x" 2,Var "y" 2]]
另一个辅助函数 crop 接收两个表达式向量，并返回第一个向量中长度与第二个向量相同的分量。
*Main> crop (bus "x" 4) (bus "y" 3)
[Var "x" 0,Var "x" 1,Var "x" 2]
利用这个辅助函数，我们现在能够定义两个在表达式向量上操作并移动其分量的移位函数
*Main> shl (bus "x" 4) 2
[Zero,Zero,Var "x" 0,Var "x" 1]
*Main> shr (bus "x" 4) 2
[Var "x" 2,Var "x" 3,Zero,Zero]
请注意，如果我们将函数结果看作表达式向量本身，函数名称似乎有些反直觉，但如果我们将向量分量解释为字中的位，那么向左移位确实意味着向最高有效位移位，这样就更清晰了。
利用这些基本元素，我们可以定义类似于 xorshift 的变换，其公式与原始 C 或 C++ 代码非常接近
*Main> let x = bus "x" 3
*Main> let y = bus "y" 3
*Main> xor x (shl y 2) ++ y
[Xor [Var "x" 0,Zero],Xor [Var "x" 1,Zero],Xor [Var "x" 2,Var "y" 0],
Var "y" 0,Var "y" 1,Var "y" 2]
我们还可以利用 Haskell 的另一个特性，即通过将函数名用反引号括起来，将其视为二进制操作
*Main> x `xor` (y `shl` 2) ++ y
[Xor [Var "x" 0,Zero],Xor [Var "x" 1,Zero],Xor [Var "x" 2,Var "y" 0],
Var "y" 0,Var "y" 1,Var "y" 2]
事实上，我们为我们的问题创建了一种迷你 DSL。
现在是时候学习如何转换我们的表达式了。最终目标将是线性形式 (\(\ref{eq:component}\))，它可以直接转换为转换矩阵。
查看类型 Expr 的定义，我们可以看到它的值是树，其中叶子可以是零或变量组件，所有内部节点都是应用于子树的 Xor 操作。在数学中，我们使用嵌套子表达式来表达这一点：\(x_1 \oplus (x_2 \oplus (y_3 \oplus \ldots))\) 或 \(x_1 \oplus ((x_2 \oplus y_3) \oplus \ldots))\)。
考虑到模 2 加法的性质，我们可以将这些表达式转换为线性形式，这两种表达式的线性形式是相同的：\(x_1 \oplus x_2 \oplus y_3 \oplus \ldots\)。函数 flatten 正是这样做的。
flatten::Expr -> [Expr]
flatten (Xor (x:xs)) = flatten x ++ (concatMap flatten xs)
flatten (Xor []) = []
flatten x = [x]
它接收一个表达式树并返回一个隐式异或在一起的原始元素序列（常量或变量）。为此，我们递归地遍历表达式树，对于每个 Xor 操作节点，我们扁平化每个子树并将它们连接在一起。其他类型的节点照原样收集到此列表中。
让我们将此函数应用于上面描述的表达式树
*Main> flatten (Xor [Var "x" 1, Xor [Var "x" 2, Xor [Var "y" 3, Zero]]])
[Var "x" 1,Var "x" 2,Var "y" 3,Zero]
*Main> flatten (Xor [Var "x" 1, Xor [Xor[Var "x" 2, Var "y" 3], Zero]])
[Var "x" 1,Var "x" 2,Var "y" 3,Zero]
（由于我们的表达式树不适合描述不完整的公式，我用零常量替换了省略号，使其格式正确。）另外请注意，扁平化后可能会出现重复项。
*Main> flatten (Xor [Var "x" 1, Xor [Xor[Var "x" 2, Var "x" 1], Zero]])
[Var "x" 1,Var "x" 2,Var "x" 1,Zero]
我们在实现最终目标方面取得了重大进展！现在我们能够以简洁的形式定义 xorshift 转换，并将公式转换为涉及变量的线性序列。我们的下一步是将这些变量列表转换为系数向量，而系数向量又将成为转换矩阵的行。
让我们再定义一个函数
matrix e vars = zipWith (\e v -> (v, makeRow e)) e vars where
makeRow e = map (\v -> length (filter (== v) e) `mod` 2) vars
参数 e 是一个表达式列表的列表 [[Expr]]。每个列表项对应状态向量的单个位，并且此项是 flatten 生成的 Var 或 Zero 项的列表。第二个参数 vars 是 Var 的列表，其中第 \(i\) 个项表示与状态向量的第 \(i\) 个项对应的变量。我们已经引入 x 和 y 作为我们的构建块，我们可以在 REPL 中评估它们以回忆它们的值。现在是时候定义我们要处理的公式了。
*Main> x
[Var "x" 0,Var "x" 1,Var "x" 2]
*Main> y
[Var "y" 0,Var "y" 1,Var "y" 2]
*Main> let e = x `xor` (y `shl` 2) ++ y
*Main> let v = x ++ y
*Main> e
[Xor [Var "x" 0,Zero],Xor [Var "x" 1,Zero],Xor [Var "x" 2,Var "y" 0],Var "y" 0,
Var "y" 1,Var "y" 2]
*Main> v
[Var "x" 0,Var "x" 1,Var "x" 2,Var "y" 0,Var "y" 1,Var "y" 2]
这里，我们有两个向量变量，x 和 y。这些三位向量组合成状态向量 \(S=\langle x_0, x_1, x_2, y_0, y_1, y_2\rangle\)。表达式 e 是这些向量上的类似 xorshift 的变换。为了获取 matrix 的第一个参数，我们应该逐分量地将 flatten 应用于表达式 e。变量 v 只是 x 和 y 中所有分量的集合，按照 matrix 的第二个参数的要求，按顺序排列。
*Main> matrix (map flatten e) v
[(Var "x" 0,[1,0,0,0,0,0]),(Var "x" 1,[0,1,0,0,0,0]),(Var "x" 2,[0,0,1,1,0,0]),
(Var "y" 0,[0,0,0,1,0,0]),(Var "y" 1,[0,0,0,0,1,0]),(Var "y" 2,[0,0,0,0,0,1])]
太棒了！输出结果与 (\(\ref{eq:component}\)) 完全一致：对于状态向量的每个位，我们都得到了一个系数向量。
为了了解 matrix 的工作原理，我们选择一个特定的状态向量项进行检查。就选择 \(x_2\) 吧，因为它的系数向量有 2 个非零项，而其他项只有一个。我们应该先准备好它
*Main> let e2 = flatten (e!!2)
*Main> e2
[Var "x" 2,Var "y" 0]
由于我们将只考虑单个状态位，计算 matrix 实际上与使用适当参数计算 makeRow 相同，而 makeRow 又可以用其定义替换。
*Main> map (\v' -> length (filter (== v') e2) `mod` 2) v
[0,0,1,1,0,0]
（请注意，我不得不将函数参数名 v 替换为 v'，以避免与用于保存涉及变量列表的顶级名称冲突。）
现在让我们看看如果我们应用最内层表达式会发生什么
*Main> v
[Var "x" 0,Var "x" 1,Var "x" 2,Var "y" 0,Var "y" 1,Var "y" 2]
*Main> map (\v' -> filter (== v') e2) v
[[],[],[Var "x" 2],[Var "y" 0],[],[]]
正如你所看到的，v 中提到的每个变量都映射到一个包含它自身和其余变量的列表，而其余变量则映射到空列表。
再增加一层，我们将最内层的 filter 用 length 包装起来
*Main> map (\v' -> length (filter (== v') e2)) v
[0,0,1,1,0,0]
看起来我们得到了想要的结果！但我们还在 matrix 定义中最后一步取了模 \(2\)。我们必须这样做才能覆盖在同一表达式中多次引用某个变量的情况。出现偶数次的变量由于这些项的相互抵消而应该消失，而取模 \(2\) 正是这样做的：所有偶数计数变为零，所有奇数计数变为一。
警告：查看我们的代码，我们可以看到 \(N\) 个表达式列表中的每个列表（其中 \(N\) 是状态向量长度）都与变量列表进行匹配，而变量列表的长度也是 \(N\)。这导致了 \(O(N^2)\) 的复杂性。我们也许可以使用排序表达式列表来改进它，但由于这段代码不常运行，我们牺牲效率以换取简洁性。
获取给定 xorshift 变换的转换矩阵问题已解决。然而，一些技术问题仍需解决。首先，我们需要将我们的位矩阵打包成整数，以便将来在 C++ 代码中使用。此目标通过以下函数实现
packMatrix m vars = map (\(v, b) -> (v, tail (packRow b vars "" 0))) m where
packRow xb@(b:bs) xv@((Var n i):vs) vname acc | n == vname = packRow bs vs vname (acc + b*2^i)
| otherwise = acc : packRow xb xv n 0
packRow [] _ _ acc = [acc]
packRow _ [] _ acc = [acc]
让我们测试一下
*Main> packMatrix (matrix (map flatten e) v) v
[(Var "x" 0,[1,0]),(Var "x" 1,[2,0]),(Var "x" 2,[4,1]),(Var "y" 0,[0,1]),
(Var "y" 1,[0,2]),(Var "y" 2,[0,4])]
函数 packMatrix 取用函数 matrix 构造的转换矩阵的每一行，并将对应相同变量族的位分组。也就是说，如果我们的表达式中有两个变量族 \(x_i\) 和 \(y_i\)，那么每行矩阵将有两个整数单词。在同一个单词内，位权重对应于变量索引，例如，\(x_i\) 的位行项在表示所有 \(x\) 分量的单词中占据第 \(i\) 位。例如，在 Var "x" 2 的矩阵行中，我们看到 [4, 1]。第一个单词是针对 \(x\)，其唯一的非零位位置为 \(2\)（从零开始计数）。实际上，\(x_2\) 的表达式包含 \(x_2\) 本身。对第二个单词执行相同的操作，我们可以看到另一个项将是 \(y_0\)。
下一个函数将接收一个打包的转换矩阵并生成 C++ 代码用于数组初始化器。
matrixCode pm = header ++ (concat $ intersperse ",\n"
(map (\(v, r) -> " {" ++ hexRow r ++ "}") pm)) ++ "\n};\n" where
header = printf "#include <cstdint>\n\nuint32_t T[%d][%d] = \
\{\n" (length pm) (length $ snd $ head pm)
hexRow r = concat $ intersperse ", " (map (printf "0x%08XU") r)
其结果如下所示
*Main> matrixCode $ packMatrix (matrix (map flatten e) v) v
"#include <cstdint>\n\nuint32_t T[6][2] = {\n {0x00000001U, 0x00000000U},\n {0
x00000002U, 0x00000000U},\n {0x00000004U, 0x00000001U},\n {0x00000000U, 0x0000
0001U},\n {0x00000000U, 0x00000002U},\n {0x00000000U, 0x00000004U}\n};\n"
我们快完成了！让我们定义 main 函数
main =
do
putStr (matrixCode $ packMatrix (matrix (map flatten xs128) vars) vars)
where x = bus "x" 32
y = bus "y" 32
z = bus "z" 32
w = bus "w" 32
vars = (x++y++z++w)
t = x `xor` (x `shl` 11)
xs128 = y ++ z ++ w ++ (w `xor` (w `shr` 19) `xor` t `xor` (t `shr` 8))
这里，我们有一个 xor128 定义，它与原始 C 代码非常相似。我们应用上述转换，并将结果字符串输出到标准输出。就这些！
还应向 Makefile 添加特殊规则，以实现此生成器的自动化
xsT.cc: xsTgen.hs
$(HS) xsTgen.hs > xsT.cc
其中 makefile 变量 HS 指向 Haskell Platform 中的 runhaskell 程序，或者您选择的其他 Haskell 实现中的等效程序。这意味着 xsT.cc 只有在文件不存在或生成器本身被修改时才会生成。本文的源代码存档包含预生成的 xsT.cc，因此您无需安装 Haskell 即可构建 C++ 代码。
向前一步
如果您知道转换矩阵 \(T\)，从 RNG 状态 \(S\) 向前跳跃 \(k\) 步就像计算 \(T^kS\) 一样简单。（有关详细信息，请参阅我的上一篇文章。）因此，我们首先需要实现这些算术原语。让我们在 xs.h 中将转换矩阵表示定义为
using tr_t = uint32_t[STATE_SIZE_EXP][STATE_SIZE];
而我们在 xs.cc 中的下一段代码将是
inline uint32_t get_bit(const tr_t A, size_t row, size_t col)
...
inline void set_bit(tr_t A, size_t row, size_t col, uint32_t val)
...
void mat_mul(tr_t C, const tr_t A, const tr_t B)
...
void mat_pow(tr_t B, const tr_t A, uint64_t n)
只要您熟悉矩阵算术、平方求幂和 C++ 中的位操作，这段代码或多或少都是直接的。因此，我将简要描述一些为了使其更清晰而不是更高效而做出的实现权衡
- 我们的矩阵 \(T\) 使用紧凑表示，其中每个整数词对应多个矩阵项，每个项一个位。因此，在 C++ 中使用按位逻辑运算，我们可以并行处理许多项，但在我为本文编写的代码中，我为了清晰度而牺牲了效率。也就是说，我简单地定义了单比特矩阵项的访问函数
get_bit和set_bit，并在这些访问器之上实现了简单的 \(O(N^3)\) 矩阵乘法。 - 另一个权衡是在矩阵乘法中使用内部临时存储。如果没有这些临时变量，调用者将负责避免别名，即，像
mat_mul(B, A, B)这样的函数调用，其中乘法结果覆盖了其中一个参数。让调用者负责分配无别名缓冲区可能更高效，但清晰度会降低。 - 计算矩阵幂的函数
mat_pow以 64 位整数作为指数。这不足以覆盖xor128的所有可能步长，但处理多字数字会使代码不够清晰。
后退一步
使用转换矩阵向后移动需要另一个特殊的矩阵，它是转换矩阵 \(T\) 的逆。也就是说，我们需要计算 \(T^{-1}\)，使得 \(TT^{-1}=I\)，其中 \(I\) 是单位矩阵。
虽然可以通过求解矩阵方程或使用其他矩阵求逆算法（这在 \(\mathbb{F}_2\) 中可能很困难）从 \(T\) 获得 \(T^{-1}\)，但有一种基于 xorshift 特性的更简单的方法。我们知道这个 RNG 的周期是 \(2^N-1\)，其中 \(N\) 是状态位数。这意味着从某个状态 \(S\) 开始并执行 \(2^N-1\) 步，我们应该回到该状态 \(S\)。在我们的矩阵表示法中，这可以写成 \(T^PS = S\)，其中 \(P = 2^N-1\)。由于根据单位矩阵的定义 \(IS = S\)，我们可以得出结论 \(T^P = I\)。我们也可以将后一个方程写成 \(TT^{P-1} = I\)，最后得出结论 \(T^{-1} = T^{P-1}\)。
无需数学也可以得到相同的结果。确实，如果执行 \(P\) 步将我们带回初始状态 \(S\)，那么如果少执行一步会发生什么？我们应该得到一个状态 \(S^*\)，使得下一步将使其变为 \(S\)，即 \(TS^*=S\)。因此，这个状态 \(S^*\) 是 \(S\) 的直接前驱，正如所需要的那样，我们可以使用 \(T^{P-1}\) 作为转换矩阵立即跳到它。
弄清楚如何计算 \(T^{-1}\) 后，实现就显而易见了。我们可以使用 mat_pow 函数来获得 \(T\) 的适当幂，但正如您所记得的，为了简化起见，我们将指数值限制为 64 位。因此，我实现了单独的函数 calculate_bwd_transition_matrix，它有效地计算矩阵幂，但硬编码的指数值为 \(2^N-2\)。这个值具有简单的二进制表示：\(N-1\) 个 1 后跟一个 0。因此，mat_pow 中检查下一个位是否为 1 的代码可以替换为位位置测试。
在进行任何向后步骤之前，应以某种方式执行此初始化。在生产代码中，我们应该使用预先计算的常量值，就像我们对矩阵 \(T\) 所做的那样，但在我们的概念验证代码中，我是在第一次使用之前计算它。这使得第一个向后步骤比向前步骤慢得多。
通用基于矩阵的转换函数
为了方便起见，有两个函数将所有矩阵操作封装成一个简单的接口。
enum class Direction {FWD, BWD};
void prepare_transition(tr_t T_k, uint64_t k, Direction dir);
void do_transition(state_t &s, const tr_t T);
第一个函数根据作为参数指定的方向，使用适当的单步矩阵计算 \(k\) 步的转换矩阵。将参数 \(k\) 设置为有符号整数并传递负值以向后跳跃看起来很有吸引力，但这个想法有其自身的缺点：即使使用 128 位整数，我们也需要完整的无符号精度来表示所有可能的跳跃大小，因此最好将方向作为单独的参数传递。
第二个函数 do_transition 只是接收 prepare_transition 准备好的矩阵，并将其应用于给定的 RNG 状态。
如果您要执行多次相同步数的跳跃，则需要将这两个转换阶段分开。由于矩阵幂运算比将转换矩阵应用于 RNG 状态更昂贵，我们可以预先计算转换矩阵一次，然后在以后多次重用它。
中途停留
文章的第一部分到此结束。我们回顾了基于矩阵的跳跃算法的背景知识，并将其应用于另一个 RNG 家族。使用函数式语言 Haskell 生成转换矩阵是一次非常不寻常但有趣的冒险。借助 Haskell，可以用几行代码开发出方便的类 DSL 符号来描述 xorshift 公式。这段代码可能看起来有些深奥，但使用 read-eval-print loop（Haskell 实现中流行（如果不是标准）的功能），您可以分部分、分步骤地开发和尝试它。
在第二部分中，我们将深入代数的奥秘，并学习如何以完全不同的方式使用多项式解决这个问题。
参考文献
[1] Hiroshi Haramoto, Makoto Matsumoto, 和 Pierre L’Ecuyer. 多项式空间中线性递推的快速跳跃算法. 见 第五届序列及其应用国际会议论文集, SETA ’08, 290–298 页, 柏林, 海德堡, 2008. Springer-Verlag. (可在线获取.)
[2] Hiroshi Haramoto, Makoto Matsumoto, Takuji Nishimura, François Panneton, 和 Pierre L’Ecuyer. \(\mathbb{F}_2\)-线性随机数生成器的高效跳跃. INFORMS Journal on Computing, 20(3):385–390, 2008. (可在线获取.)
[3] Donald E. Knuth. 计算机程序设计艺术，第二卷（第三版）：半数值算法. Addison-Wesley Longman Publishing Co., Inc., 波士顿, 马萨诸塞州, 美国, 1997.
[4] P. L’Ecuyer 和 R. Simard. TestU01: 用于随机数生成器经验测试的 C 库. ACM Transactions on Mathematical Software, 33(4): Article 22, 2007 年 8 月.
[5] George Marsaglia. Xorshift RNGs. Journal of Statistical Software, Articles, 8(14):1–6, 2003. (可在线获取.)
历史
- 2020 年 4 月 12 日：初始版本