梅森旋转算法

概念

梅森旋转算法（Mersenne twister），可以快速产生高质量的伪随机数，修正了古典随机数发生算法的很多缺陷。

常见的两种为基于32位的 MT19937和基于64位的 MT19937-64。

由于梅森旋转算法是利用线性反馈移位寄存器（LFSR）产生随机数的，

对于LFRS有结论：一个 $k$ 位的移位寄存器，选取合适的特征多项式（即加1为本原多项式）时，取得最大周期 $2^k-1$.

而MT19937梅森旋转算法的周期为 $2^{19937}-1$（正如算法名，这是一个梅森素数），说明它是一个19937级的线性反馈移位寄存器，梅森旋转算法是利用线性反馈寄存器一直进行移位旋转，因此实际上 MT19937-32只需要用32位就能做到。它能做到在 $1 leq kleq 623$ 都可以均匀分布。

注：梅森素数：如果形如 $2^n-1$ 的数是一个素数，则称之为梅森素数。

算法步骤

整个算法主要分为三个阶段：

• 第一阶段：获得基础的梅森旋转链；
• 第二阶段：对于旋转链进行旋转算法；
• 第三阶段：对于旋转算法所得的结果进行处理；

算法实现的过程中，参数的选取取决于梅森素数，故此得名。

实现

使用MT19937算法生成范围在 $[2^{32}-1]$ 的均匀分布的32位整数。

#include <iostream>
#include <string.h>
#include <stdio.h>
#include <time.h>
using namespace std;

bool isInit;
int index;
unsigned int MT[624];  //624 * 32 - 31 = 19937

void srand(int seed)
{
    printf("seed:%d
", seed);
    index = 0;
    isInit = 1;
    MT[0] = seed;
    //对数组的其它元素进行初始化
    for(int i = 1; i < 624; i++)
    {
        unsigned int t = 1812433253 * (MT[i - 1] ^ (MT[i - 1] >> 30)) + i;
        MT[i] = t & 0xffffffff;   //取最后的32位赋给MT[i]
    }
}

void generate()
{
    for(int i = 0; i < 624; i++)
    {
        // 2^31 = 0x80000000
        // 2^31-1 = 0x7fffffff
        unsigned int y = (MT[i] & 0x80000000) + (MT[(i + 1) % 624] & 0x7fffffff);
        MT[i] = MT[(i + 397) % 624] ^ (y >> 1);
        if (y & 1)
            MT[i] ^= 2567483615;
    }
}

unsigned int Rand()
{
    if(!isInit)
        srand((int)time(NULL));
    if(index == 0)
        generate();
    unsigned int y = MT[index];
    y = y ^ (y >> 11);                 //y右移11个bit
    y = y ^ ((y << 7) & 2636928640);   //y左移7个bit与2636928640相与，再与y进行异或
    y = y ^ ((y << 15) & 4022730752);  //y左移15个bit与4022730752相与，再与y进行异或
    y = y ^ (y >> 18);                 //y右移18个bit再与y进行异或
    index = (index + 1) % 624;
    return y;
}

int main()
{
    int cnt = 0;
    for(int i = 0; i < 1000000; i++)
    {
        if(Rand() & 1)  cnt++;
    }
    cout<<cnt / 1000000.0<<"%"<<endl;
    return 0;
}

应用

梅森旋转算法是R、Python、Ruby等的默认伪随机数产生器。从C++11开始，C++也可以使用这个算法。

例如C++11中：

#include <iostream>
#include <random>
using namespace std;

int main()
{
  random_device rd;  //可视为真随机数，但Windows下调用的rand_s就不是了
  mt19937 mt(rd());
  for(int i = 0; i < 10; i++)
    cout << mt() << endl;
  return 0;
}

参考链接：

1. 维基百科——梅森旋转算法

2. 维基百科——线性反馈移位寄存器

3. 维基百科——梅森素数

4. 线性反馈移位寄存器和梅森旋转算法——acdreamers

5.C++11带来的随机数生成器——egmkang