你没听过的梅森旋转算法
(标准开头)
如果单独提梅森旋转算法可能大家都很陌生,但如果说到C++11的random可能大家就都熟悉多了。事实上,C++,python等多种计算机语言的随机数都是通过梅森旋转算法产生的。(也有一个称呼是梅森缠绕算法)
那,本文就着重介绍这个梅森螺旋旋转算法
(算法本身挺学术的,我努力写得轻松点)
先在这里感谢一下@dgklr大佬的引导。如果没有他提及,笔者可能还不知道这个算法。
旋转算法简介
梅森旋转算法,也可以写作MT19937。是有由松本真和西村拓士在1997年开发的一种能快速产生优质随机数的算法。
其实这个算法跟梅森没有什么关系,它之所以叫做是梅森旋转算法是因为它的循环节是2^19937-1,这个叫做梅森素数。
如果看过我的那篇随机数的文章应该知道关于伪随机的一些知识。这个随机算法之所以说是产生“优质“”随机数,特点就是它的循环节特别长。而且产生的数分布是比较平均的。
可能有的同学对这个循环节有点质疑。可能觉得2^19937-1有点短?
我在这里大概给一个概念:
银河系中的恒星数量级10^11
撒哈拉沙漠中的沙子数数量级是10^26
宇宙中目前可观察的粒子数量级是10^87
219937数量级是106001
这个比较大概心里有数了吧
相差的已经不止是一个数量级了
同时他在623维中的分布都十分的均匀(这个不用理解)
知道分布平均就好了
(梅森镇楼)
->continue
前置知识
分析这个算法的原理需要的前置知识在网上讲的都比较绕,我在这里就通俗的科普一下,主要是认识这几个名词。
(用词不准确轻喷)
线性反馈移位寄存器(LFSR)
这个,就当它是随机数发生器就完事了,不要太去纠结定义。后面会讲。
本原多项式
简单的说来就是没法化简的多项式
比如 y=x4+x2 就可以化简
也是知道就好
级
计算机的一个二进制单位(0或1)就是一级
这个应该比较好理解
反馈函数
这个应该是网上看别的博客最绕的知识点
简单地理解成告诉你你要对这个寄存器干什么的一个函数就好了
(看到这里应该还没懵吧)
异或
这个...
还要我科普吗?
就是两个数,如果都是0或都是1就输出0,一个1一个0输出1.
->continue
原理分析
这个旋转算法实际上是对一个19937级的二进制序列作变换。
首先我们达成一个共识:
一个长度为n的二进制序列,它的排列长度最长为2^n。
当然这个也是理论上的,实际上可能因为某些操作不当,没挺到2^n个就开始循环了。
那么如何将这个序列的排列撑满2^n个,就是这个旋转算法的精髓。
如果反馈函数的本身+1是一个本原多项式,那么它的循环节达到最长,即2^n-1
这个数学证明本文不作过多论述,有兴趣者可以自己查阅资料
个人感觉单讲知识点挺难懂的(笔者就是这么被坑的)
我们就拿一个4级的寄存器模拟一下:
我们这里使用的反馈函数是 y=x4+x2+x+1(这个不是本原多项式,只是拿来好理解) 这个式子中x4,x2,x的意思就是我们每次对这个二进制序列的从后往前数第4位和第2位做异或运算 ,然后再拿结果和最后一位做异或运算。把最后的结果放到序列的开头,整个序列后移一位,最后一位舍弃(或者输出)
- 初始数组 { 1 , 0 , 0 , 0 } (为什么不是 0,0,0,0 你们可以自己想想,文章末尾揭晓)
- 将它的第四位和第二位抓出来做异或运算
- 把刚刚的运算结果和最后一位再做一次运算
- 把最后的运算结果放到第一位,序列后移。最后一位被无情的抛弃
这就是一次运算,然后这个算法就是不断循环这几步,从而不断伪随机改变这个序列。
上图是一个网上找的一个4级寄存器的模拟过程
大家可以推一下,它所使用的反馈函数(y=x^4+x+1)
因为这个是本原多项式
所以他最后的循环节是2^4-1=15
运算结果如下:
(图片摘自原文链接)
关于旋转
可能有人到这里还没看出“旋转”在哪里。因为我们每次计算出来的结果会放在开头,序列后移一位。看起来就像数组在向后旋转...
(本来想做gif的,后来不知道怎么做出旋转)
大家自行脑补
->continue
代码实现
(笔者很懒,直接搬原代码出处的代码)
#include <iostream>
#include <string.h>
#include <stdio.h>
#include <time.h>
using namespace std;
bool isInit;
int index;
int MT[624]; //624 * 32 - 31 = 19937
void srand(int seed)
{
index = 0;
isInit = 1;
MT[0] = seed;
for(int i=1; i<624; i++)
{
int t = 1812433253 * (MT[i-1] ^ (MT[i-1] >> 30)) + i;
MT[i] = t & 0xffffffff; //取最后的32位
}
}
void generate()
{
for(int i=0; i<624; i++)
{
// 2^31 = 0x80000000
// 2^31-1 = 0x7fffffff
int y = (MT[i] & 0x80000000) + (MT[(i+1) % 624] & 0x7fffffff);
MT[i] = MT[(i + 397) % 624] ^ (y >> 1);
if (y & 1)
MT[i] ^= 2567483615;
}
}
int rand()
{
if(!isInit)
srand((int)time(NULL));
if(index == 0)
generate();
int y = MT[index];
y = y ^ (y >> 11);
y = y ^ ((y << 7) & 2636928640);
y = y ^ ((y << 15) & 4022730752);
y = y ^ (y >> 18);
index = (index + 1) % 624;
return y; //笔者注:y即为产生的随机数
}
int main()
{
srand(0); //设置随机种子
int cnt = 0;
for(int i=0; i<1000000; i++) //下面的循环是用来判断随机数的奇偶概率的
{
if(rand() & 1)
cnt++;
}
cout<<cnt / 10000.0<<"%"<<endl;
return 0;
}
->continue
填一下前面的坑
这里回答一下前面的那个问题:
为什么循环节是2n-1而不是2n
这个问题的答案和为什么初始序列不能是 { 0 , 0 , 0 , 0 }是一样的,因为如果全是0的话,无论怎么异或运算都不能产生循环。那么还怎么伪随机啊。
因为不能是全0,所以循环节要-1
(* o *)
( ⊕ o ⊕ )
最后非常感谢你能有耐心读到这里。
大家都很强,可与之共勉。