你沒聽過的梅森旋轉算法

（標準開頭）

若是單獨提梅森旋轉算法可能你們都很陌生，但若是說到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 ⊕ )

最後很是感謝你能有耐心讀到這裏。

你們都很強，可與之共勉。