借助梅森旋转原理设计 mt44497
LG_jyc
·
·
科技·工程
设计想法
众所周知，C++11 标准加入了基于梅森旋转的伪随机数生成器 mt19937，而它名字里的 19937 则是梅森素数 2^{19937}-1 的指数部分，下面是截至目前已经发现的所有梅森素数。
| 序号 |
指数 p |
梅森素数 M_p = 2^p - 1 |
位数 |
发现年份 |
发现者/验证者 |
| 1 |
2 |
2^2 - 1 = 3 |
1 |
古代 |
已知 |
| 2 |
3 |
2^3 - 1 = 7 |
1 |
古代 |
已知 |
| 3 |
5 |
2^5 - 1 = 31 |
2 |
古代 |
已知 |
| 4 |
7 |
2^7 - 1 = 127 |
3 |
古代 |
已知 |
| 5 |
13 |
2^{13} - 1 = 8191 |
4 |
1456 |
匿名 |
| 6 |
17 |
2^{17} - 1 = 131071 |
6 |
1588 |
Cataldi |
| 7 |
19 |
2^{19} - 1 = 524287 |
6 |
1588 |
Cataldi |
| 8 |
31 |
2^{31} - 1 = 2147483647 |
10 |
1772 |
Euler |
| 9 |
61 |
2^{61} - 1 |
19 |
1883 |
Pervushin |
| 10 |
89 |
2^{89} - 1 |
27 |
1911 |
Powers |
| 11 |
107 |
2^{107} - 1 |
33 |
1914 |
Powers |
| 12 |
127 |
2^{127} - 1 |
39 |
1876 |
Lucas |
| 13 |
521 |
2^{521} - 1 |
157 |
1952 |
Robinson |
| 14 |
607 |
2^{607} - 1 |
183 |
1952 |
Robinson |
| 15 |
1279 |
2^{1279} - 1 |
386 |
1952 |
Robinson |
| 16 |
2203 |
2^{2203} - 1 |
664 |
1952 |
Robinson |
| 17 |
2281 |
2^{2281} - 1 |
687 |
1952 |
Robinson |
| 18 |
3217 |
2^{3217} - 1 |
969 |
1957 |
Riesel |
| 19 |
4253 |
2^{4253} - 1 |
1281 |
1961 |
Hurwitz |
| 20 |
4423 |
2^{4423} - 1 |
1332 |
1961 |
Hurwitz |
| 21 |
9689 |
2^{9689} - 1 |
2917 |
1963 |
Gillies |
| 22 |
9941 |
2^{9941} - 1 |
2993 |
1963 |
Gillies |
| 23 |
11213 |
2^{11213} - 1 |
3376 |
1963 |
Gillies |
| 24 |
19937 |
2^{19937} - 1 |
6002 |
1971 |
Tuckerman |
| 25 |
21701 |
2^{21701} - 1 |
6533 |
1978 |
Noll、Nickel |
| 26 |
23209 |
2^{23209} - 1 |
6987 |
1979 |
Noll |
| 27 |
44497 |
2^{44497} - 1 |
13395 |
1979 |
Nelson、Slowinski |
| 28 |
86243 |
2^{86243} - 1 |
25962 |
1982 |
Slowinski |
| 29 |
110503 |
2^{110503} - 1 |
33265 |
1988 |
Colquitt、Welsh |
| 30 |
132049 |
2^{132049} - 1 |
39751 |
1983 |
Slowinski |
| 31 |
216091 |
2^{216091} - 1 |
65050 |
1985 |
Slowinski |
| 32 |
756839 |
2^{756839} - 1 |
227832 |
1992 |
Slowinski、Gage |
| 33 |
859433 |
2^{859433} - 1 |
258716 |
1994 |
Slowinski、Gage |
| 34 |
1257787 |
2^{1257787} - 1 |
378632 |
1996 |
Slowinski、Gage |
| 35 |
1398269 |
2^{1398269} - 1 |
420921 |
1996 |
Armengaud |
| 36 |
2976221 |
2^{2976221} - 1 |
895932 |
1997 |
Spence |
| 37 |
3021377 |
2^{3021377} - 1 |
909526 |
1998 |
Clarkson |
| 38 |
6972593 |
2^{6972593} - 1 |
2098960 |
1999 |
Hajratwala |
| 39 |
13466917 |
2^{13466917} - 1 |
4053946 |
2001 |
Cameron |
| 40 |
20996011 |
2^{20996011} - 1 |
6320430 |
2003 |
Shafer |
| 41 |
24036583 |
2^{24036583} - 1 |
7235733 |
2004 |
Findley |
| 42 |
25964951 |
2^{25964951} - 1 |
7816230 |
2005 |
Nowak |
| 43 |
30402457 |
2^{30402457} - 1 |
9152052 |
2005 |
Cooper、Boone |
| 44 |
32582657 |
2^{32582657} - 1 |
9808358 |
2006 |
Cooper、Boone |
| 45 |
37156667 |
2^{37156667} - 1 |
11185272 |
2008 |
Elvenich |
| 46 |
42643801 |
2^{42643801} - 1 |
12837064 |
2009 |
Strindmo |
| 47 |
43112609 |
2^{43112609} - 1 |
12978189 |
2008 |
Smith |
| 48 |
57885161 |
2^{57885161} - 1 |
17425170 |
2013 |
Cooper |
| 49 |
74207281 |
2^{74207281} - 1 |
22338618 |
2016 |
Cooper |
| 50 |
77232917 |
2^{77232917} - 1 |
23249425 |
2017 |
Pace |
| 51 |
82589933 |
2^{82589933} - 1 |
24862048 |
2018 |
Patrick Laroche |
| 52 |
136279841 |
2^{136279841} - 1 |
41024320 |
2024 |
Luke Durant |
可以看到 2^{19937}-1 的后面还有很多梅森素数，而其中有一个 2^{44497}-1 的指数 44497 不仅没有后面那么大，而且还好记。于是很容易想到设计一个原理类似的 mt44497。
思路推导
我们先看 mt19937 的原理（这样方便给出一些常数的值）。
1 内部状态
mt19937 内部有一个状态数组 state[624] 共 624 个整形，这个数也约等于 \dfrac{19937}{32}。
2 递推方式
线性递推
有以下公式：
x_{k+n}=x_{k+m} \oplus ((\operatorname{upper}(x_k) \parallel \operatorname{lower}(x_{k+1}))\cdot A)
其中，\operatorname{upper}(a) 指 a 的第 31 位，而 \operatorname{lower}(a) 指 a 的第 0 至第 30 位，符号 \parallel 表示位拼接，\oplus 是异或。
在 mt19937 中，n=624，m=397。
$$
(x \cdot A)=\begin{cases}
\lfloor \dfrac{x}{2} \rfloor & \operatorname{LSB}(x)=0\\
\lfloor \dfrac{x}{2} \rfloor\oplus 9908B0DF_{(16)} & \operatorname{LSB}(x)=1
\end{cases}
$$
其中 $\operatorname{LSB}(x)$ 表示 $x$ 的第 $0$ 位。
#### 输出变换
对生成的 $x_{k+n}$ 进行非线性变换提高质量：
$$
y=x_{k+n} \oplus \lfloor \dfrac{x_{k+n}}{2^u} \rfloor\\
y=y \oplus ((y \times 2^s) \operatorname{and} b)\\
y=y \oplus ((y \times 2^t) \operatorname{and} c)\\
y=y \oplus \lfloor \dfrac{y}{2^l} \rfloor
$$
在 `mt19937` 中，$u=11,s=7,t=15,l=18,b=9D2C5680_{(16)},c=EFC60000_{(16)}$。
$y$ 即是返回的结果，具体地，`mt19937` 在每生成 $624$ 个随机数后再统一生成 $624$ 个整数在 `state` 中。
### 3 让我们类推一下
我们对于 `mt44497` 发现参数基本上可以照搬，只有像 $n,m$ 这样的数需要修改。
于是我们就得到了 `mt44497` 的基本代码实现。
## 代码实现
这里使用无符号 $32$ 位整数。
```cpp
class mt44497{
public:
typedef unsigned int uint32;
enum{STATE_SIZE=1392,DEFAULT_SEED=5489U,MAX_VALUE=0xFFFFFFFFU};
explicit mt44497(uint32 Seed=DEFAULT_SEED){
seed(Seed);
}
void seed(uint32 seed=DEFAULT_SEED){
state[0]=seed;
for(pos=1;pos<STATE_SIZE;++pos){
state[pos]=1812433253U*(state[pos-1]^(state[pos-1]>>30))+pos;
}
}
uint32 operator()(){
if(pos>=STATE_SIZE){twistState();}
uint32 y=state[pos++];
y^=(y>>15)|(y<<17);
y^=(y>>11);
y^=(y<<7)&0x9D2C5680U;
y=(y&0xFFFFFF00)|((y^(y>>8))&0xFFFFFF);
y^=(y<<15)&0xEFC60000U;
y^=(y<<18);
return y;
}
uint32 min()const{return 0;}
uint32 max()const{return MAX_VALUE;}
private:
void twistState(){
const uint32 MATRIX_A=0x9908B0DFU;
const uint32 UPPER_MASK=0x80000000U;
const uint32 LOWER_MASK=0x7FFFFFFFU;
for(int k=0;k<STATE_SIZE;++k){
uint32 x=(state[k]&UPPER_MASK)|(state[(k+1)%STATE_SIZE]&LOWER_MASK);
state[k]=state[(k+397)%STATE_SIZE]^(x>>1);
if(x&0x1U){
state[k]^=MATRIX_A;
}
}
pos=0;
}
uint32 state[STATE_SIZE];
int pos;
};
```