你可能没听说过的梅森旋转随机数算法

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