基本介紹
- 中文名：梅森旋轉算法
- 外文名：Mersenne twister
套用[編輯],優點,缺點,算法詳細,相關代碼,偽代碼,Python 代碼,
套用[編輯]
梅森旋轉算法是R、Python、Ruby、IDL、Free Pascal、PHP、Maple、Matlab、GNU多重精度運算庫和GSL的默認偽隨機數產生器。從C++11開始，C++也可以使用這種算法。在Boost C++,Glib和NAG數值庫中，作為外掛程式提供。
在SPSS中，梅森旋轉算法是兩個PRNG中的一個：另一個是產生器僅僅為保證舊程式的兼容性，梅森旋轉被描述為“更加可靠”。梅森旋轉在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 << 7 & 2636928640 y = y ^ y << 15 & 4022730752 y = y ^ y >> 18 return _int32(y) def twist(self): for i in range(0, 624): y = _int32((self.mt[i] & 0x80000000) + (self.mt[(i + 1) % 624] & 0x7fffffff)) self.mt[i] = y ^ self.mt[(i + 397) % 624] >> 1 if y % 2 != 0: self.mt[i] = self.mt[i] ^ 0x9908b0df
調用函式MT19937(seed).extract_number()將會返回隨機數，其中seed是已確定的種子。
