梅森旋转算法

科技工作者之家 2020-11-17

梅森旋转算法Mersenne twister)是一个伪随机数发生算法。由松本真和西村拓士在1997年开发,基于有限二进制字段上的矩阵线性递归。可以快速产生高质量的伪随机数,修正了古典随机数发生算法的很多缺陷。

应用[编辑]梅森旋转算法是R、Python、Ruby、IDL、Free Pascal、PHP、Maple、Matlab、GNU多重精度运算库和GSL的默认伪随机数产生器。从C++11开始,C++也可以使用这种算法。在Boost C++,Glib和NAG数值库中,作为插件提供。

在SPSS中,梅森旋转算法是两个PRNG中的一个:另一个是产生器仅仅为保证旧程序的兼容性1,梅森旋转被描述为“更加可靠”。梅森旋转在SAS中同样是PRNG中的一个,另一个产生器是旧时的且已经被弃用。

优点最为广泛使用Mersenne Twister的一种变体是MT19937,可以产生32位整数序列。具有以下的优点:

周期非常长,达到2−1。尽管如此长的周期并不必然意味着高质量的伪随机数,但短周期(比如许多旧版本软件包提供的2)确实会带来许多问题。

在1 ≤k≤ 623的维度之间都可以均等分布(参见定义)。

除了在统计学意义上的不正确的随机数生成器以外,在所有伪随机数生成器法中是最快的(当时)

缺点为了性能,这个算法付出了巨大的空间成本(当时而言):需要 2.5KiB的缓存空间。2011年,松本真和西村拓士针对这一问题提出了一个更小的版本,仅占127 bits的 TinyMT (Tiny Mersenne Twister)。

算法详细整个算法主要分为三个阶段:

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

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

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

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

相关代码下面的一段伪代码使用MT19937算法生成范围在[0, 2−1]的均匀分布的32位整数:

伪代码//创建一个长度为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 } } }Python 代码def _int32(x): return int(0xFFFFFFFF & x)class MT19937: def __init__(self, seed): self.mt = [0] * 624 self.mt[0] = seed for i in range(1, 624): self.mt[i] = _int32(1812433253 * (self.mt[i - 1] ^ self.mt[i - 1] >> 30) + i) def extract_number(self): self.twist() y = self.mt[0] y = y ^ y >> 11 y = y ^ y 1 if y % 2 != 0: self.mt[i] = self.mt[i] ^ 0x9908b0df调用函数MT19937(seed).extract_number()将会返回随机数,其中seed是已确定的种子。

本词条内容贡献者为:

王伟 - 副教授 - 上海交通大学

科技工作者之家

科技工作者之家APP是专注科技人才,知识分享与人才交流的服务平台。