你没听过的梅森旋转算法

（标准开头）

若是单独提梅森旋转算法可能你们都很陌生，但若是说到C++11的random可能你们就都熟悉多了。事实上，C++，python等多种计算机语言的随机数都是经过梅森旋转算法产生的。（也有一个称呼是梅森缠绕算法）

那，本文就着重介绍这个梅森螺旋旋转算法

（算法自己挺学术的，我努力写得轻松点）

先在这里感谢一下@dgklr大佬的引导。若是没有他说起，笔者可能还不知道这个算法。

旋转算法简介

梅森旋转算法，也能够写做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 大家能够本身想一想，文章末尾揭晓）

1. 将它的第四位和第二位抓出来作异或运算

1. 把刚刚的运算结果和最后一位再作一次运算

1. 把最后的运算结果放到第一位，序列后移。最后一位被无情的抛弃

这就是一次运算，而后这个算法就是不断循环这几步，从而不断伪随机改变这个序列。

上图是一个网上找的一个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

填一下前面的坑

这里回答一下前面的那个问题：

为何循环节是2^n-1而不是2^n

这个问题的答案和为何初始序列不能是 { 0 , 0 , 0 , 0 }是同样的，由于若是全是0的话，不管怎么异或运算都不能产生循环。那么还怎么伪随机啊。

由于不能是全0，因此循环节要-1

（* o *）

（ ⊕ o ⊕ )

最后很是感谢你能有耐心读到这里。

你们都很强，可与之共勉。