算法定义

伪随机数

伪随机数就是可以通过特定的初始值可以得出特定的伪随机数，这些伪随机数在统计上具有均匀性、独立性。所以在现实生活中，我们可以通过伪随机数生成算法来满足生成随机数的需求。

梅森素数旋转法

它是一种伪随机数生成算法，它能够根据某个初始种子生成随机串。由于该该算法的"旋转"特性，生成的随机串是有周期的。在周期外的数是可以通过上一周期产生的随机串推出来，故在周期外就不具备随机中的独立性了。但是在此算法中，我们不必担心，因为它的周期是梅森素数$2^{19937}-1$，这数非常大，大到足以让我们安心的使用此随机数生成算法。比如像c++、R、python、MATLAB的随机数产生都是运用此算法实现。比较常见的有两种算法，分别是基于32位的MT19937-32和基于64位的MT19937-64，在此我们介绍的是MT19937-32，MT19937-32算法能够生成$[0,2^{32}-1]$范围内的整数随机数.

k-distributed to v-bit accuracy（k分布v位准确度）

生成随机数算法有很多，所以需要一个标准来评估这些随机数生成算法的质量。其中一个度量算法就是k-distributed to v-bit accuracy（翻译不准确，故在此不翻译啦~）。接下来具体来讲下这个评估算法：
假设某一个随机数生成算法能够产生周期为$P$的$w$比特的随机串{$x_i$},同时将$w$比特的随机数$x$的最高$v$（前$v$）位组成的数记作为$trun{c}_v(x)$.因此可以构造出如下的二进制数，长度为$kv$比特：
$PRN{G}_{k,v}(i)=(trun{c}_v(x_i),trun{c}_v(x_{i+1}),...,trun{c}_v(x_{i+k-1}))$

显而易见，$PRN{G}_{k,v}(i)$的长度是kv，所以总共有$2^{kv}$种取值。如果我们说随机串{$x_i$}是满足k-distributed to v-bit accuracy的，则当$x_i$取遍周期P内的值时，$PRN{G}_{k,v}(i)$在$2^{kv}$种取值中是均匀的。

对于固定的v，我们易知，如果{$x_i$}满足k-distributed to v-bit accuracy，那么他一定是满足 (k-1)-distributed to v-bit accuracy.在原论文中我们证明了MT19937-32算法生成的随机串{$x_i$}是满足 623-distributed to 32-bit accuracy的。所以生成的随机数序列{$x_i$}一定满足1-distributed to 32-bit accuracy.所以可证明随机数序列{$x_i$}在周期长度为P($2^{19937}-1$)内是均匀分布的。

算法步骤

旋转

梅森旋转素数法总共分两步，第一步为旋转，第二步为提取。此节介绍下旋转过程，旋转过程主要有如下公式实现：
$$x_{n+k}:=x_{m+k}\oplus (x_k^u|x_{k+1}^l)\ast A$$

符号解释：

• $x_i$:随机数序列{$x_i$}中的第$i$个元素，在此计算中$x_i$可看成是$1\ast w$的行向量，即可写成向量$x=(x_{w-1},x_{w-2},...,x_0)$.$x_i$实际可以看成我们生成的随机整数的二进制形式。
• $n:$随机数列{$x_i$}的长度，在不同形式的随机算法，$n$取值各不相同，在MT19937-32中，我们取$n$为623.
• $m:m\in [1,n].$在不同形式的随机算法，$m$取值各不相同，在MT19937-32中，我们取$m$为397.
• $x^u:$指的是由$x$前$u$位组成$1\ast w$的行向量，向右补零。
• $x^l:$指的是由$x$后$l$位组成$1\ast w$的行向量，向左补零。
• 矩阵$A_{w\ast w}$的形式如下：
  $$A=\left[\begin{array}{cccc}0& 1& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & 1\\ a_{w-1}& a_{w-2}& \cdots & a_0\end{array}\right]$$

所以我们可知：
$$xA=\{\begin{array}{rlr}x>>1& & x_0=0\\ x>>1\oplus \hat{a}& & x_0=1\end{array}$$
在梅森旋转算法中，它是利用线性反馈移位寄存器来实现旋转的。所以我们先来介绍线性反馈移位寄存器。反馈移位寄存器主要有两部分组成：

• 级，参与变换的长度。
• 反馈函数，如果反馈函数为线性则称为线性反馈移位寄存器，若为非线性，则称为非线性反馈移位寄存器。

一般来说，LFSR是基于异或操作的，具体步骤如下：

• 将原寄存器中的最低位作为输出。
• 按照反馈函数，选择若干高位依次取抑或，然后与步骤一中的输出（最低位）取抑或作为输出。
• 将原寄存器中的元素向右移一位，然后将步骤二中的输出作为最高位。最终寄存器中的数即为此次旋转后得到的数。

举个例子，假设反馈函数是$f(x)=x^4+x^2+x+1$，即需要将最高位与第三高位取抑或，然后再与最低位取抑或得到的数作为移位后的最高位，假设初始值位1000，则变换如下：

• 1000
• 1100
• 1110
• 0111
• 0011
• 0001
• 1000

此变换周期为6，当然不同的反馈函数周期是不一样的，比如反馈函数$f(x)=x^4+x+1$，它的周期达到了15，它能够取遍4位的所有值（0除外，0不管怎么变换都是0），大家有兴趣可以去验证下。在此，MT19937-32在旋转中，周期可高达$2^{19937}-1$

提取

MT19937具备两个特性,一个是周期大，另一个是满足 623-distributed to 32-bit accuracy.上面的寄存器达到了我们周期大的目的。接下来的提取，则可以满足第二特性。
提取的方法就是在每次旋转后再乘以个k可逆矩阵T即可，步骤如下：

• $y\leftarrow x\oplus (x>>u)$
• $y\leftarrow y\oplus ((y<<s)$ AND $b)$
• $y\leftarrow y\oplus ((y<<t)$ AND $c)$
• $z\leftarrow y\oplus (y>>l)$

算法步骤

由此可见，梅森旋转法主要分为两步：旋转和提取。具体算法步骤如下：

参数说明：

• w：生成随机数的二进制长度
• n：参与旋转的随机数个数(旋转的深度)
• m：参与旋转的中间项
• r : $x_k^{(u)}$与$x_{k+1}^{(l)}$的切分位
• a: 矩阵A的最后
• u,s,t,l：整数参数，移位运算的移动距离;
• b,c ：比特遮罩。

于是算法步骤如下：

1.初始化

• $lower\leftarrow 0$
• $lower\leftarrow (lower<<1)$ AND 1 for r times
• $upper\leftarrow$ COMPL $lower$
• $a\leftarrow a_{w-1}{a}_{w-2}...a_1{a}_0$
• $i\leftarrow 0$

2.非零初始化
$x[0],x[1],...,x[n-1]$

3.旋转

• $t\leftarrow (x[i]$ $AND$ $upper)$ OR $(x[(i+1])$ mod $n]$ AND $lower)$
• $x[i]\leftarrow x[(i+m)$ MOD $n]$ XOR $(t>>1)$ XOR $\{\begin{array}{ll}0,& \text{if }{t}_0=0\\ a,& \text{otherwise}\end{array}$

4.提取输出

• $y\leftarrow x[i]$
• $y\leftarrow yXOR(y>>u)$
• $y\leftarrow yXOR((y<<s)$ AND $b)$
• $y\leftarrow yXOR((y<<t)$ AND $c)$
• $z\leftarrow yXOR(y>>l)$

5.更新循环变量$i\leftarrow (i+1)$ MOD $n$，而后跳转至步骤3 继续进行。

文献：

• https://liam0205.me/2018/01/12/Mersenne-twister/
• https://blog.csdn.net/Touch_Dream/article/details/68948708