Mersenne Twister (梅森旋转算法）的MATLAB实现

想自己编写一个MATLAB的梅森旋转算法产生随机数的程序，但是找遍网络也没找到拿MATLAB写的，所以试着写一个~

适合学过相关统计知识的朋友

Mersenne Twister 简介

• 迄今为止最好的随机数发生器之一
• 被标准C++、 Python、 MATLAB 等流行编程工具采用
• MT19937-32和MT19937-64的使用最为广泛
• 周期为2¹⁹⁹³⁷-1

不了解移位寄存器产生随机数的小伙伴可以去看这几个博客：
伪随机数生成算法-梅森旋转（Mersenne Twister/MT）算法介绍
伪随机数生成——梅森旋转（Mersenne Twister/MT）算法笔记

Mersenne Twister 的使用

1. Mersenne Twister 递推式
2. Twist
3. A的特殊结构使得xA的计算得到简化：
   1. xA = x >> 1 ; (x0 = 0)
   2. xA = (x >> 1) xor a ; (x0 = 1)
4. Temper: 为了改善均匀性引入的计算：
   1. y = x xor (x >> u)
   2. y = y xor ((y << s) & b)
   3. y = y xor ((y << t) & c)
   4. z = z xor (y >> l)

Mersenne Twister 相关代码

1.伪代码

//创建一个长度为624的数组来存储发生器的状态
int[0..623] MT int index = 0  
//用一个种子初始化发生器
//

function initialize_generator(int seed) 
{
     i := 0  
     MT[0] := seed 
     for i from 1 to 623 
     { 
     // 遍历剩下的每个元素         
     MT[i] := last 32 bits of(1812433253 * (MT[i-1] xor (right shift by 30 bits(MT[i-1]))) + i) 
     // 1812433253 == 0x6c078965     
     } 
}  

// Extract a tempered pseudorandom number based on the index-th value,
// calling generate_numbers() every 624 numbers 

function extract_number() 
{     
	if index == 0 
	{    
	     generate_numbers()    
	 }      
	 int y := MT[index]     
	 y := y xor (right shift by 11 bits(y))    
	 y := y xor (left shift by 7 bits(y) and (2636928640)) 
	 // 2636928640 == 0x9d2c5680    
	 y := y xor (left shift by 15 bits(y) and (4022730752))
	 // 4022730752 == 0xefc60000     
	 y := y xor (right shift by 18 bits(y))     
	 index := (index + 1) mod 624     
	 return y 
}  

// Generate an array of 624 untempered numbers
function generate_numbers() 
{
     for i from 0 to 623 
     {
              int y := (MT[i] & 0x80000000)                       
              // bit 31 (32nd bit) of MT[i]                     + (MT[(i+1) mod 624] & 0x7fffffff)   
              // bits 0-30 (first 31 bits) of MT[...]         
              MT[i] := MT[(i + 397) mod 624] xor (right shift by 1 bit(y))        
              if (y mod 2) != 0 { // y is odd             
              MT[i] := MT[i] xor (2567483615) 
              // 2567483615 == 0x9908b0df         
              }     
              }
               }

2. Python 代码

#梅森旋转算法   -xlxw
#参考:mersenne twister from wikipedia

#import
from time import time
import numpy as np

#var
index = 624
MT = [0]*index
# MT[0] ->seed

def inter(t):
    return(0xFFFFFFFF & t) #取最后32位->t

def twister():
    global index
    for i in range(624):
        y = inter((MT[i] & 0x80000000) +(MT[(i + 1) % 624] & 0x7fffffff))
        MT[i] = MT[(i + 397) % 624] ^ y >> 1
        if y % 2 != 0:
            MT[i] = MT[i] ^ 0x9908b0df
    index = 0

def exnum():
    global index
    if index >= 624:
        twister()
    y = MT[index]
    y = y ^ y >> 11
    y = y ^ y << 7 & 2636928640
    y = y ^ y << 15 & 4022730752
    y = y ^ y >> 18
    index = index + 1
    return inter(y)

def mainset(seed):
    MT[0] = seed    #seed
    for i in range(1,624):
        MT[i] = inter(1812433253 * (MT[i - 1] ^ MT[i - 1] >> 30) + i)
    return exnum()

def main():
    br = input("请输入随机数产生的范围(用,隔开):")
    mi = eval(br.split(',')[0])
    ma = eval(br.split(',')[1])    
    so = mainset(int(time())) / (2**32-1)
    rd = mi + int((ma-mi)*so)
    print("产生的随机整数为：",rd)
main()

3. Matlab 代码
根据Python就很容易可以写出MATLAB实现了