梅森旋轉算法
此條目可參照英語維基百科相應條目來擴充。 (2021年11月1日) |
梅森旋轉演算法(Mersenne twister)是一個偽隨機數發生算法。由松本眞和西村拓士[1]在1997年開發,基於有限二進制字段上的矩陣線性遞歸。可以快速產生高質量的偽隨機數,修正了古典隨機數發生算法的很多缺陷。
Mersenne Twister這個名字來自周期長度取自梅森質數的這樣一個事實。這個算法通常使用兩個相近的變體,不同之處在於使用了不同的梅森素數。一個更新的和更常用的是MT19937, 32位字長。還有一個變種是64位版的MT19937-64。對於一個k位的長度,Mersenne Twister會在的區間之間生成離散型均勻分布的隨機數。
應用
[編輯]梅森旋轉算法是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位整數序列。具有以下的優點:
- 周期非常長,達到219937−1。儘管如此長的周期並不必然意味着高質量的偽隨機數,但短周期(比如許多舊版本軟件包提供的232)確實會帶來許多問題。[2]
- 在1 ≤ k ≤ 623的維度之間都可以均等分布(參見定義)。
- 除了在統計學意義上的不正確的隨機數生成器以外,在所有偽隨機數生成器法中是最快的(當時)[3]
缺點
[編輯]為了性能,這個算法付出了巨大的空間成本(當時而言):需要 2.5 KiB 的緩存空間。2011年,松本真和西村拓士針對這一問題提出了一個更小的版本,僅占127 bits的 TinyMT (Tiny Mersenne Twister)。[4]
k-分布
[編輯]暫無
其他選擇
[編輯]算法詳細
[編輯]整個算法主要分為三個階段:
第一階段:獲得基礎的梅森旋轉鏈;
第二階段:對於旋轉鏈進行旋轉算法;
第三階段:對於旋轉算法所得的結果進行處理;
算法實現的過程中,參數的選取取決於梅森素數,故此得名。
相關代碼
[編輯]下面的一段偽代碼使用MT19937算法生成範圍在[0, 232 − 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
self.mti = 0
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):
if self.mti == 0:
self.twist()
y = self.mt[self.mti]
y = y ^ y >> 11
y = y ^ y << 7 & 2636928640
y = y ^ y << 15 & 4022730752
y = y ^ y >> 18
self.mti = (self.mti + 1) % 624
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 >> 1) ^ self.mt[(i + 397) % 624]
if y % 2 != 0:
self.mt[i] = self.mt[i] ^ 0x9908b0df
調用函數 MT19937(seed).extract_number()
將會返回隨機數,其中 seed
是已確定的種子。
SFMT
[編輯]實現
[編輯]- XLL Excel Addin Implementation of mersenne twister random number generator (頁面存檔備份,存於網際網路檔案館)
- Two implementations of Mersenne Twister in Java: one is the fastest known, and the other is a drop-in replacement for java.util.Random (頁面存檔備份,存於網際網路檔案館)
- The GNU Scientific Library (GSL), containing an implementation of the Mersenne Twister (頁面存檔備份,存於網際網路檔案館)
- C++ and binary function libraries for several platforms. Multithreaded. Includes Mersenne Twister and SFMT (頁面存檔備份,存於網際網路檔案館)
- Implementations of the Mersenne Twister in C and C++ (頁面存檔備份,存於網際網路檔案館)
- Implementation of the Mersenne Twister in C++
- Implementation of the Mersenne Twister in C++ within the Mobile Robot Programming Toolkit
- Implementation of Mersenne Twister as an add-in for Microsoft Excel (頁面存檔備份,存於網際網路檔案館)
- Implementation of Mersenne Twister as a free module for Visual Basic (Microsoft Excel, Microsoft Access and VB compilers) and for other Basic versions in the official site of the Mersenne Twister (頁面存檔備份,存於網際網路檔案館)
- Implementation of Mersenne Twister for REALbasic (requires REALbasic 2006r1 or greater)
- Implementation of Mersenne Twister for Lisp
- Implementation of Mersenne Twister in Euphoria
- Implementation of Mersenne Twister for C# (newer, System.Random drop-in replacement) (頁面存檔備份,存於網際網路檔案館) (Older implementation)
- Implementation of Mersenne Twister for Ada
- Implementation of Mersenne Twister for Fortran 95
- Implementation of Mersenne Twister for Mathematica
- Implementation of Mersenne Twister for MATLAB (頁面存檔備份,存於網際網路檔案館)
- Implementation of Mersenne Twister for Mitrion-C (頁面存檔備份,存於網際網路檔案館)
- Implementation of Mersenne Twister for Clean
- High-speed Implementation of Mersenne Twister (頁面存檔備份,存於網際網路檔案館) in Linoleum (a cross-platform Assembler), by Herbert Glarner
- CPAN module implementing the Mersenne Twister for use with Perl (頁面存檔備份,存於網際網路檔案館)
- Implementation of Mersenne Twister for Haskell (頁面存檔備份,存於網際網路檔案館)
- Implementation of Mersenne Twister for Standard ML
- Implementation of Mersenne Twister in F#
- It also is implemented in gLib and the standard libraries of at least PHP, Python and Ruby.
- RandomLib (頁面存檔備份,存於網際網路檔案館) C++ library implementing Mersenne Twister and SFMT.
- C++ implementation of Mersenne Twister for the IBM/Sony Cell Broadband Engine (Cell BE) specialized processing units
- Mersenne Twister ported to ActionScript
- Mersenne Twister in Clojure (頁面存檔備份,存於網際網路檔案館)
- Implementation of Mersenne Twister for ABAP (頁面存檔備份,存於網際網路檔案館)
參考列表
[編輯]- ^ Makoto Matsumoto, Takuji Nishimura. Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator. ACM Transactions on Modeling and Computer Simulation (TOMACS). 1998-01-01, 8 (1): 3–30 [2018-04-02]. ISSN 1049-3301. doi:10.1145/272991.272995.
- ^ 註:219937約等於4.3×106001,這個值比可觀測宇宙內粒子總數的估計值(1087)還要高出上千個數量級。
- ^ P. L'Ecuyer and R. Simard, ``TestU01: A C Library for Empirical Testing of Random Number Generators, ACM Transactions on Mathematical Software, 33, 4, Article 22, August 2007.
- ^ Tiny Mersenne Twister (TinyMT). hiroshima-u.ac.jp. [4 October 2015]. (原始內容存檔於2020-07-11).