跳转至

Mersenne Twister 19937

简介

MT19937 即梅森旋转算法（Mersenne twister）由松本眞和西村拓士在 1997 年开发，基于二进制有限域 \(\mathbb{F}_2\) ​上的矩阵线性递归，可以快速产生高质量的伪随机数。

该算法具有以下优点：

  • 周期很长，为 \(P = 2^{19937} - 1\);
  • \(623\) - 维 \(32\) - 比特正确

\(k\)- 维 \(v\)- 比特正确性

假设某一 PRNG 能够产生周期为 \(P\) 的 \(w\)- 比特的随机数序列 \(\{\vec{x}_i\}\)；同时，将 \(w\)- 比特的随机数 \(\vec{x}\) 的最高 \(v\) 位组成的数记作 \(\operatorname{trunc}_v (\vec{x})\)。现构造如下长度为 \(kv\)- 比特的二进制数

\[ \operatorname{PRNG}_{k, v} (i) := (\operatorname{trunc}_v (\vec{x}_i), \operatorname{trunc}_v (\vec{x}_{i + 1}), \ldots, \operatorname{trunc}_v (\vec{x}_{i + k - 1})) \]

显然 \(\operatorname{PRNG}_{k, v} (i)\) 有 \(2^{kv}\) 种不同的取值。当 \(i\) 遍历 \([0, P)\) 时，若 \(\operatorname{PRNG}_{k, v} (i)\) 的取值在这上面均匀分布，具体来说，就是非零的取值次数为 \(\dfrac{P + 1}{2^{kv}}\)，全为 \(0\) 的取值次数为 \(\dfrac{P + 1}{2^{kv}} - 1\)，则称该 PRNG 是 \(k\)- 维 \(v\)- 比特准确的。

对于固定的 \(v\)，如果该 PRNG 是 \(k\)- 维 \(v\)- 比特准确的，那么其也是 \((k - 1)\)- 维 \(v\)- 比特准确的，但不一定是 \((k + 1)\)- 维 \(v\)- 比特准确的。\(k\) 显然是有上限的，因此对于固定的 \(v\)，必然存在最大的 \(k = k(v)\) 使得 PRNG 是 \(k(v)\)- 维 \(v\)- 比特准确的。依据定义，有 \(2^{k(v)v} - 1 \leqslant P\)。

如果 PRNG 是 \(k\)- 维 \(v\)- 比特准确的，那么即使已知 PRNG 生成的 \(l < k\) 个伪随机数的最高 \(v\) 位，也无法推出 PRNG 生成的第 \(l + 1\) 个伪随机数的最高 \(v\) 位。

算法描述

旋转

MT19937 能够产生周期为 \(P\) 的 \(w\)- 比特的伪随机数序列 \(\{\vec{x}_i\}\)，其中 \(w = 32\)。定义如下记号来描述 MT19937 是如何进行旋转（线性移位）的：

  • \(n\): 参与 Mersenne Twist 的随机数个数；
  • \(r\): \([0, w)\) 之间的整数；
  • \(m\): \((0, n]\) 之间的整数；
  • \(\mathbf{A}\): \(w \times w\) 的常矩阵；
  • \(\vec{x}^{(u)}\): \(\vec{x}\) 的最高 \(w - r\) 比特组成的数（低位补零）；
  • \(\vec{x}^{(l)}\): \(\vec{x}\) 的最低 \(r\) 比特组成的数（高位补零）。

首先根据随机数种子初始化 \(n\) 个行向量：

\[ \vec{x}_0, \vec{x}_1, \ldots, \vec{x}_{n - 1} \]

而后根据下面这个式子，从 \(k = 0\) 开始依次计算 \(\vec{x}_n\)：

\[ \vec{x}_{k + n} = \vec{x}_{k + m} \oplus (\vec{x}_k^{(u)} \mid \vec{x}_{k + 1}^{(l)})\mathbf{A} \]

其中 \(\vec{x} \mid \vec{x}'\) 表示两个二进制数按位或，\(\vec{x} \oplus \vec{x}'\) 表示两个二进制数按位半加（不进位，也就是按位异或）；\(\vec{x}\mathbf{A}\) 表示按位半加的矩阵乘法。\(\mathbf{A}\) 的定义如下：

\[ \begin{pmatrix} & & 1 & & & & \\ & & & 1 & & & \\ & & & & \ddots & \\ & & & & & 1 & \\ & a_{w - 1} & a_{w - 2} & a_{w - 3} & \ldots & a_0 \end{pmatrix} \]

所以

\[ \vec{x}\mathbf{A} = \begin{cases} \vec{x} >> 1 & \text{if } x_0 = 0 \\ (\vec{x} >> 1) \oplus \vec{a} & \text{if } x_0 = 1 \end{cases} \]

所以整个旋转的计算过程完全有位运算组成（移位、按位与、按位或、按位异或）。

Linear Feedback Shift Register

但上面的运算过程怎么看也不像“旋转”，所以这里先介绍线性反馈移位寄存器（Linear Feedback Shift Register，LFSR）的“旋转”，再将其与 Mersenne Twister 相结合。

反馈移位寄存器是一种对二进制序列进行等长变换的函数，包括两个部分：

  1. 级：等长变换的长度；
  2. 反馈函数：若反馈函数是线性的，则称线性反馈移位寄存器；否则是非线性反馈移位寄存器。

一般的 LFSR 都是基于异或运算的，工作流程如下：

  1. 将原寄存器状态的最低位 \(x_0\) 作为输出；
  2. 执行线性反馈函数，即选取寄存器中的若干位，从高位到低位迭代异或；
  3. 将原寄存器状态向低位移位 1 位，并以上述迭代异或的结果作为填充值填入最高位。

评论