メルセンヌツイスター

メルセンヌツイスターは、 1997 年に松本 真(松本 眞)と西村 拓司(西村 拓士)によって開発された汎用擬似乱数生成器(PRNG) です。^{[1] [2]}その名前は、周期長として メルセンヌ素数を選択したことに由来しています。

Mersenne Twister は、古い PRNG に見られる欠陥のほとんどを修正するために特別に設計されました。

メルセンヌツイスターアルゴリズムの最も一般的に使用されているバージョンは、メルセンヌ素数に基づいています。その標準実装であるMT19937は、32ビットのワード長を使用します。64ビットのワード長を使用する別の実装（5つのバリエーション^[3] ）MT19937-64もあり、異なるシーケンスを生成します。 ${\displaystyle 2^{19937}-1}$

周期Pのwビット整数の疑似乱数シーケンスは、次の条件が成立する場合、 vビットの精度でk 分布していると言われます。${\displaystyle x_{i}}$

trunc _v ( x ) をxの先頭のvビットで構成される数とし、k vビットベクトル のPを考える。

${\displaystyle (\operatorname {trunc} _{v}(x_{i}),\operatorname {trunc} _{v}(x_{i+1}),\,\ldots ,\operatorname {trunc} _{v}(x_{i+k-1}))\quad (0\leq i<P)}$。

次に、すべてのゼロの組み合わせが 1 回少ない頻度で発生することを除き、ビットの可能な各組み合わせが 1 周期内に同じ回数発生します。${\displaystyle 2^{kv}}$

wビットのワード長の場合、メルセンヌツイスタは範囲 の整数を生成します。 ${\displaystyle [0,2^{w}-1]}$

メルセンヌツイスターアルゴリズムは、有限の2進体上の行列線形回帰に基づいています。このアルゴリズムは、状態ビットの反映と調整を備えた、有理正規形(TGFSR(R))のツイスト一般化フィードバックシフトレジスタ^[4] (ツイストGFSR、またはTGFSR)です。基本的な考え方は、単純な再帰関係を通じて級数を定義し、形式の数値を出力することです。ここで、Tは調整行列と呼ばれる可逆な -行列です。 ${\displaystyle {\textbf {F}}_{2}}$${\displaystyle x_{i}}$${\displaystyle x_{i}^{T}}$${\displaystyle {\textbf {F}}_{2}}$

一般的なアルゴリズムは、次の量によって特徴付けられます。

• w : ワードサイズ（ビット数）
• n : 再発の度合い
• m : 中間語、系列を定義する再帰関係で使用されるオフセット、${\displaystyle x}$${\displaystyle 1\leq m<n}$
• r : 1ワードの区切り点、または下位ビットマスクのビット数、${\displaystyle 0\leq r\leq w-1}$
• a : 有理正規形ツイスト行列の係数
• b、c : TGFSR(R) テンパリング ビットマスク
• s , t : TGFSR(R) 調整ビットシフト
• u、d、l : 追加のメルセンヌツイスター調整ビットシフト/マスク

はメルセンヌ素数であるという制限があります。この選択により、パラメータ検索に必要な 原始性テストとk分布テストが簡素化されます。${\displaystyle 2^{nw-r}-1}$

この級数は、再帰関係を持つ wビット量の級数として定義されます。${\displaystyle x}$

${\displaystyle x_{k+n}:=x_{k+m}\oplus \left(({x_{k}}^{u}\mid {x_{k+1}}^{l})A\right)\qquad k=0,1,2,\ldots }$

ここで、 はビットベクトルの連結（上位ビットが左側）を表し、ビットごとの排他的論理和（XOR）は の上位w − rビットを意味し、は の下位rビットを意味します。 ${\displaystyle \mid}$${\displaystyle \oplus}$${\displaystyle x_{k}^{u}}$${\displaystyle x_{k}}$${\displaystyle x_{k+1}^{l}}$${\displaystyle x_{k+1}}$

添え字はすべて-nでオフセットされる可能性がある

${\displaystyle x_{k}:=x_{k-(nm)}\oplus \left(({x_{kn}}^{u}\mid {x_{k-(n-1)}}^{l})A\right)\qquad k=n,n+1,n+2,\ldots }$

ここで、LHS は、RHS にある過去に生成された値に基づいて、シリーズ内で次に生成される値です。 ${\displaystyle x_{k}}$

ツイスト変換A は、有理正規形で次のように定義されます。 は単位行列 です。有理正規形の利点は、Aによる乗算を次のように効率的に表現できることです。(ここでは行列の乗算が で行われているため、加算の代わりにビット単位の XOR が使用されることに注意してください)。はの最下位ビットです。 $${\displaystyle A={\begin{pmatrix}0&I_{w-1}\\a_{w-1}&(a_{w-2},\ldots ,a_{0})\end{pmatrix}}}$$${\displaystyle I_{w-1}}$${\displaystyle (w-1)(w-1)}$${\displaystyle {\textbf {F}}_{2}}$$${\displaystyle {\boldsymbol {x}}A={\begin{cases}{\boldsymbol {x}}\gg 1&x_{0}=0\\({\boldsymbol {x}}\gg 1)\oplus {\boldsymbol {a}}&x_{0}=1\end{cases}}}$$${\displaystyle x_{0}}$${\displaystyle x}$

TGFSR(R) と同様に、メルセンヌツイスターは、等分布の次元の減少を補正するために、調整変換とカスケード接続されます ( Aが有理正規形であることを選択したため)。これは、 T が可逆行列である行列A を使用することと同等であり、したがって、以下で説明する特性多項式の分析が依然として有効であることに注意してください。 ${\displaystyle A=T^{-1}*AT}$

Aと同様に、簡単に計算できるような調整変換を選択するため、T自体は実際には構築しません。この調整は、メルセンヌツイスターの場合、次のように定義されます。

${\displaystyle {\begin{aligned}y&\equiv x\oplus ((x\gg u)~\And ~d)\\y&\equiv y\oplus ((y\ll s)~\And ~b)\\y&\equiv y\oplus ((y\ll t)~\And ~c)\\z&\equiv y\oplus (y\gg l)\end{aligned}}}$

ここで、は系列からの次の値、は一時的な中間値、 はアルゴリズムから返される値です。およびはビット単位の左シフトと右シフト、 はビット単位のANDです。 最初と最後の変換は、下位ビットの均等配分を改善するために追加されます。 TGFSR の特性から、上位ビットの均等配分の上限に到達するには が必要です。 ${\displaystyle x}$${\displaystyle y}$${\displaystyle z}$${\displaystyle \ll }$${\displaystyle \gg }$${\displaystyle \&}$${\displaystyle s+t\geq \left\lfloor {\frac {w}{2}}\right\rfloor -1}$

MT19937 の係数は次のとおりです。

${\displaystyle {\begin{aligned}(w,n,m,r)&=(32,624,397,31)\\a&={\textrm {9908B0DF}}_{16}\\(u,d)&=(11,{\textrm {FFFFFFFF}}_{16})\\(s,b)&=(7,{\textrm {9D2C5680}}_{16})\\(t,c)&=(15,{\textrm {EFC60000}}_{16})\\l&=18\\\end{aligned}}}$

メルセンヌツイスタの 32 ビット実装では、通常、d  = FFFFFFFF ₁₆であることに注意してください。その結果、dを含むビット単位のand は効果を持たないため、アルゴリズムの説明からdが省略されることがあります。

MT19937-64の係数は以下の通りである: ^[5]

${\displaystyle {\begin{aligned}(w,n,m,r)=(64,312,156,31)\\a={\textrm {B5026F5AA96619E9}}_{16}\\(u,d)=(29,{\textrm {5555555555555555}}_{16})\\(s,b)=(17,{\textrm {71D67FFFEDA60000}}_{16})\\(t,c)=(37,{\textrm {FFF7EEE000000000}}_{16})\\l=43\\\end{aligned}}}$

メルセンヌツイスターの実装に必要な状態は、それぞれwビットのn個の値の配列です。配列を初期化するには、wビットのシード値を使用してシード値を 設定し、その後に${\displaystyle x_{0}}$${\displaystyle x_{n-1}}$${\displaystyle x_{0}}$

${\displaystyle x_{i}=f\times (x_{i-1}\oplus (x_{i-1}\gg (w-2)))+i}$

からまでの。 ${\displaystyle i}$${\displaystyle 1}$${\displaystyle n-1}$

• アルゴリズムが生成する最初の値はではなくに基づきます。${\displaystyle x_{n}}$${\displaystyle x_{0}}$
• 定数f は、アルゴリズム自体の一部ではありませんが、ジェネレーターの別のパラメーターを形成します。
• MT19937 のfの値は1812433253 です。
• MT19937-64のfの値は6364136223846793005である。 ^[5]

#include <stdint.h> 

#define n 624 
#define m 397 
#define w 32 
#define r 31 
#define UMASK (0xffffffffUL << r) 
#define LMASK (0xffffffffUL >> (wr)) 
#define a 0x9908b0dfUL 
#define u 11 
#define s 7 
#define t 15 
#define l 18 
#define b 0x9d2c5680UL 
#define c 0xefc60000UL 
#define f 1812433253UL

typedef struct { uint32_t state_array [ n ]; // 状態ベクトルの配列int state_index ; // 状態ベクトル配列のインデックス、常に 0 <= state_index <= n-1 } mt_state ; 

              
                      
 


void initial_state ( mt_state * state , uint32_t seed ) { uint32_t * state_array = & ( state -> state_array [ 0 ]); state_array [ 0 ] = seed ; // 推奨される初期シード = 19650218UL for ( int i = 1 ; i < n ; i ++ ) { seed = f * ( seed ^ ( seed >> ( w -2 ))) + i ; // Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. state_array [ i ] = seed ; } state -> state_index = 0 ; }     

       
    
                                
    
        
    
                      
           
    
    
      



uint32_t random_uint32 ( mt_state * state ) { uint32_t * state_array = & ( state -> state_array [ 0 ]); int k = state -> state_index ; // 現在の状態の位置を指す// 常に 0 <= state_index <= n-1 // int k = k - n; // n 反復前の状態を指す// if (k < 0) k += n; // n を法とする循環インデックス// 前の 2 行は実際には何もしません// 説明のみint j = k - ( n -1 ); // n-1 反復前の状態を指すif ( j < 0 ) j += n ; // n を法とする循環インデックス  

       
    
             
                                     
    


                                     
                                     
    
                        
                         

    uint32_t x = ( state_array [ k ] & UMASK ) | ( state_array [ j ] & LMASK ); uint32_t xA = x >> 1 ; if ( x & 0x00000001UL ) xA ^= a ; j = k - ( n - m ); // nm 反復前の状態を指すif ( j < 0 ) j += n ; // n を法とする循環インデックスx = state_array [ j ] ^ xA ; // 状態の次の値を計算state_array [ k ++ ] = x ; // 新しい状態値を更新if ( k >= n ) k = 0 ; // n を法とする循環インデックスstate -> state_index = k ; uint32_t y = x ^ ( x >> u ); // 調整y = y ^ (( y << s ) & b ); y = y ^ (( y << t ) & c ); uint32_t z = y ^ ( y >> l ); zを返す; }         
    
         
          
    
                           
                         
    
                 
                  
    
                         
      
    
                  
                     
                     
           
    
      

T GFSRの周期の理論的な上限を達成するためには、は原始多項式でなければならず、これは次の 特性多項式である。${\displaystyle 2^{nw-r}-1}$${\displaystyle \phi _{B}(t)}$${\displaystyle \phi _{B}(t)}$

${\displaystyle B={\begin{pmatrix}0&I_{w}&\cdots &0&0\\\vdots &&&&\\I_{w}&\vdots &\ddots &\vdots &\vdots \\\vdots &&&&\\0&0&\cdots &I_{w}&0\\0&0&\cdots &0&I_{w-r}\\S&0&\cdots &0&0\end{pmatrix}}{\begin{matrix}\\\\\leftarrow m{\text{-th row}}\\\\\\\\\end{matrix}}}$

${\displaystyle S={\begin{pmatrix}0&I_{r}\\I_{w-r}&0\end{pmatrix}}A}$

ツイスト変換により、従来のGFSR が以下の重要な特性によって改善されます。

• 周期が理論上の上限に達する（0で初期化された場合を除く）${\displaystyle 2^{nw-r}-1}$
• n次元での均等分布（例えば、線形合同型ジェネレータはせいぜい5次元での合理的な分布を管理できる）

CryptMTはストリーム暗号であり、内部でメルセンヌツイスターを使用する暗号学的に安全な疑似乱数生成器です。 ^{[6] [7]}これは、松本と西村が萩田真理子と斉藤睦夫と共同で開発しました。eCRYPTネットワークのeSTREAMプロジェクトに提出されました。^[6]メルセンヌツイスターや他の派生とは異なり、CryptMTは特許を取得しています。

MTGPは、斉藤睦夫と松本誠によって発表された、グラフィックス処理装置向けに最適化されたメルセンヌツイスターの変種です。 ^[8]基本的な線形再帰演算はMTから拡張され、多くのスレッドが状態空間を共有してメモリ負荷を軽減しながら並行して再帰を計算できるようにパラメータが選択されています。論文では、 MTよりも均等配分が改善され、古い（2008年頃）GPU（ 192コアのNvidia GTX260）で5×10 ⁷のランダムな32ビット整数 に対して4.7ミリ秒のパフォーマンスが得られたと主張しています。

SFMT（SIMD指向の高速メルセンヌツイスター）は、2006年に導入されたメルセンヌツイスターの変種であり、^[9] 128ビットSIMDで実行した場合に高速になるように設計されています。

• これはメルセンヌツイスターの約2倍の速度です。^[10]
• これは、MT よりも v ビット精度の等分布特性が優れていますが、 WELL ("Well Equidistributed Long-period Linear")よりも劣っています。
• ゼロ過剰初期状態からの回復は MT よりも速いですが、WELL よりも遅いです。
• ^{2 607}  − 1 から 2 ²¹⁶⁰⁹¹ − 1 までのさまざまな周期をサポートします 。

SFMTはIntel SSE2とPowerPC AltiVecをサポートしています。また、PlayStation 3のCell BEを使用したゲームにも使用されています。^[11]

TinyMTは、2011年に斎藤と松本によって提案されたメルセンヌツイスタの変種です。^[12] TinyMTはわずか127ビットの状態空間を使用し、オリジナルの2.5KiBの状態空間と比較すると大幅に減少しています。しかし、その周期は であり、オリジナルよりもはるかに短いため、メモリが貴重である場合にのみ著者によって推奨されています。 ${\displaystyle 2^{127}-1}$

利点:

• CryptMT を除くすべてのバリアントは、許可ライセンスおよび特許フリーです。
• DiehardテストやTestU01テストのほとんど（すべてではない）を含む、統計的ランダム性に関する多数のテストに合格します。^[13]
• 非常に長い周期。長い周期は乱数生成器の品質を保証するものではないが、多くの古いソフトウェアパッケージで一般的な短い周期は問題になる可能性があることに注意してください。^[14]${\displaystyle 2^{19937}-1}$${\displaystyle 2^{32}}$
• 32ビット精度でk分布する（ k分布の定義については以下を参照）${\displaystyle 1\leq k\leq 623}$
• 実装は一般に、ハードウェア実装の方法よりも速く乱数を生成します。ある研究によると、メルセンヌツイスターは、ハードウェア実装のプロセッサベースのRDRAND命令セットよりも約20倍速く64ビット浮動小数点乱数を生成します。^[15]

デメリット:

• TinyMT バリアントを使用しない限り、状態バッファは約 2.5 kBと比較的大きくなります 。
• SFMTバリアント（後述）を使用しない限り、現代の基準では平凡なスループットです。^[16]
• TestU01スイートのCrushとBigCrushの両方で2つの明らかな失敗（線形複雑性）を示しています。このテストは、メルセンヌツイスターと同様に、-代数に基づいています。^[13]${\displaystyle {\textbf {F}}_{2}}$
• シード値のみが異なる（他のパラメータは異なる）複数のインスタンスは、独立した乱数ジェネレータを必要とするモンテカルロシミュレーションには一般的に適していませんが、複数のパラメータ値セットを選択する方法は存在します。^{[17] [18]}
• 拡散が悪い: 初期状態が非常に非ランダムな場合、特に初期状態に多くのゼロがある場合、ランダム性テストに合格する出力を生成するのに長い時間がかかることがあります。この結果、ほぼ同じ初期状態で開始されたジェネレータの2つのインスタンスは、最終的に発散する前に、多くの反復でほぼ同じシーケンスを出力します。 MTアルゴリズムの2002年の更新では初期化が改善されたため、このような状態から開始することはほとんどありません。^[19] GPUバージョン（MTGP）はさらに優れていると言われています。^[20]
• 1 より 0 が多いサブシーケンスが含まれています。これにより拡散特性が悪くなり、多数のゼロがある状態からの回復が困難になります。
• CryptMTバリアント (後述) を使用しない限り、暗号的に安全ではありません。その理由は、十分な数の反復 (MT19937 の場合は 624。これは将来の反復が生成される状態ベクトルのサイズであるため) を観察することで、将来のすべての反復を予測できるためです。

Mersenne Twister は、次のソフトウェアによってデフォルトの PRNG として使用されます。

• プログラミング言語: Dyalog APL、^[21] IDL、^[22] R、^[23] Ruby、^[24] Free Pascal、^[25] PHP、^[26] Python （ NumPyでも利用可能、ただしバージョン1.17 ^[27]以降ではデフォルトがPCG64に変更された）、^{[28] [29] [30]} CMU Common Lisp、^[31] Embeddable Common Lisp、^[32] Steel Bank Common Lisp、^[33] Julia（Julia 1.6 LTSまで、それ以降でも利用可能、1.7以降ではより優れた/高速なRNGがデフォルトで使用される）^[34]
• Unix系ライブラリとソフトウェア: GLib、^[35] GNU Multiple Precision Arithmetic Library、^[36] GNU Octave、^[37] GNU Scientific Library ^[38]
• その他: Microsoft Excel、^[39] GAUSS、^[40] gretl、^[41] Stata、^[42] SageMath、^[43] Scilab、^[44] Maple、^[45] MATLAB ^[46]

また、Apache Commons ^[47]、標準C++ライブラリ（C++11以降）^{[48] [49]}、Mathematica ^[50]でも利用可能です。アドオン実装は、Boost C++ライブラリ^[51]、CUDAライブラリ^[52]、NAG Numerical Library ^[53]など、多くのプログラムライブラリで提供されています。

メルセンヌツイスターは、 SPSSの2つのPRNGのうちの1つです。もう1つのジェネレーターは、古いプログラムとの互換性のためだけに残されており、メルセンヌツイスターは「より信頼性が高い」と言われています。^{[54]メルセンヌツイスターは、同様に}SASのPRNGの1つです。他のジェネレーターは古く、非推奨です。^[55]メルセンヌツイスターはStataのデフォルトのPRNGであり、もう1つはKISSで、古いバージョンのStataとの互換性のためです。^[56]

代替ジェネレータであるWELL（「Well Equidistributed Long-period Linear」）は、より迅速な回復、均等なランダム性、ほぼ同等の速度を提供します。^[57]

マルサリアのXORシフトジェネレータとその変種はLFSRクラスの中で最速である。^[58]

64ビットMELG（「メルセンヌ素数周期を持つ64ビット最大等分布線形ジェネレータ」）はk分布特性に関して完全に最適化されています。^[59]${\displaystyle {\textbf {F}}_{2}}$

ACORNファミリー(1989 年公開) は、別のk分散 PRNG であり、MT と同様の計算速度と、現在の (2019 年) TestU01 基準をすべて満たすため、より優れた統計特性を示します。適切なパラメーターを選択して使用すると、ACORN は任意の長さの周期と精度を持つことができます。

PCGファミリーは、より現代的な長周期ジェネレータであり、キャッシュの局所性が向上し、最新の分析方法を使用して検出可能なバイアスが少なくなっています。^[60]

• Harase, S. (2014)、「メルセンヌツイスター擬似乱数生成器の -線形関係について」、 Mathematics and Computers in Simulation、100 : 103–113、arXiv : 1301.5435、doi :10.1016/j.matcom.2014.02.002、S2CID  6984431${\displaystyle \mathbb {F} _{2}}$。
• Harase, S. (2019)、「メルセンヌツイスターの倍精度浮動小数点数への変換」、Mathematics and Computers in Simulation、161 : 76–83、arXiv : 1708.06018、doi :10.1016/j.matcom.2018.08.006、S2CID  19777310。

• MTの学術論文と松本誠氏の関連記事
• Mersenne Twister のホームページ。C、Fortran、Java、Lisp などの言語のコードが掲載されています。
• Mersenne Twister の例 — 複数のプログラミング言語での Mersenne Twister 実装のコレクション - GitHubで
• SFMT の実践: パート I - SSE2 サポートを含む DLL の生成 - Code Projectで