梅森旋转(Mersenne Twister/MT)算法

posted at 2022.3.25 06:04 by 风信子

   梅森旋转算法,也可以写作MT19937。是有由松本真和西村拓士在1997年开发的一种能快速产生优质伪随机数的算法。

   之所以叫做是梅森旋转算法是因为它的循环节是2^19937-1,这个叫做梅森素数。这个随机算法之所以说“优质”,是因为它的周期长,对于一个k位2进制数,梅森旋转算法可在[0,2^k-1]的范围内生成离散型均匀分布的伪随机数。在623维中的分布也十分的均匀。

一、优点

•许可免费,而且对所有它的变体专利免费(除CryptMT外) 

•几乎无处不在:它被包含在大多数编程语言和库中 

•通过了包括Diehard测试在内的大多数统计随机性测试(除TestU01测试外) 

•在应用最广泛的MT19937变体中,周期长达2^19937-1 

•在MT19937-32的情况下对1 ≤ k ≤ 623,满足k-分布 

•比其他大多数随机数发生算法要快

二、缺点

•需要大量的缓冲器(2.5kib),但在TinyMT版本中修正 

•吞吐量中等,但在SFMT版本中修正 

•产生的随机数与seed相关,不能用于蒙特卡洛模拟 

•由相同的初始序列产生的随机状态几乎相同,但在2002年的更新中对MT算法的初始化进行了改进,使得对于相同的初始序列也能产生不同的随机状态 

•非加密安全的,除CryptMT外

三、原理分析

   这个旋转算法实际上是对一个19937级的二进制序列作变换。

   首先一个长度为n的二进制序列,它的循环排列长度最长为2^n。当然这个是理论上的,实际上可能没到2^n个就开始循环了。

   那么如何将这个序列的排列撑满2^n个,就是这个旋转算法的精髓。如果反馈函数的本身+1是一个没法化简的多项式,那么它的循环节达到最长,即2^n-1。

   算法基于标准(线性)旋转反馈移位寄存器(twisted generalised feedback shift register/TGFSR)产生随机数。

四、算法三个阶段 

第一阶段:初始化,获得基础的梅森旋转链; 

第二阶段:对于旋转链进行旋转算法; 

第三阶段:对于旋转算法所得的结果进行处理。

五、C++代码

#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 it=0; it<33; it++)  //下面的循环是用来列示随机数的 
    {
        int y=rand();
        cout<<y<<endl;
    }
 
    return 0;
}
 

六、效果图

Tags: , , ,

IT技术

添加评论

  Country flag

biuquote
  • 评论
  • 在线预览
Loading