Mersenne Twister on Mathematica

前口上

Mathematica と "メルセンヌ　ツイスター" でググっても出てこなかった．

Cのソースをもとにした実装

とりあえず動くものを作ってみよう．けれども， Mathematica らしいプログラムにしたい気持ちも沸いてきた…

パラメータと配列

n = 624

m = 397

matrixA = 16^^9908b0df

Off[General::spell1]

upperMask = 16^^80000000

lowerMask = 16^^7fffffff

Array[mt, n];

mti の値で初期化の有無を調べるようだ

mti = n + 1;

種の初期化

上で宣言するのがめんどうなので，ひとつにまとめてしまった．

initGenrand[s_Integer] := 
      Module[{}, 
n = 624; 
m = 
          397; 
matrixA = 
          16^^9908b0df; 
upperMask = 
          16^^80000000; 
lowerMask = 
          16^^7fffffff; 
Array[mt, 
          n]; 
  
mt[0] = 
          BitAnd[s, 16^^ffffffff]; 
Table[
          mt[mti] = 
            1812433253*BitXor[mt[mti - 1], IntegerPart[mt[mti - 1]/2^30]] + 
              mti; 
mt[mti] = 
            BitAnd[mt[mti], 16^^ffffffff];, 
{mti, 
            1, (n - 1)}]; 
mti = 
          n; 
]

initGenrand[100]

配列による初期化（未実装）

initByArray[initKey_List] := 
      Module[{},
]

32ビット乱数の生成

genrandInt32[] := 
      Module[{y,  mag01, kk}, 
mag01[0] = 
          0; 
mag01[1] = matrixA; 
If[
          mti ≥ n, 
If[
            mti == (n + 1),  
            initGenrand[
              5489]]; 
Table[
y = 
              BitOr[BitAnd[mt[kk], upperMask], 
                BitAnd[mt[kk + 1], lowerMask]]; 
mt[kk] = 
              BitXor[mt[kk + m], IntegerPart[y/2],  
                mag01[BitAnd[y, 1]]];, {kk, 
              0, (n - m - 
                  1)}
]; 
Table[
y = 
              BitOr[BitAnd[mt[kk], upperMask], 
                BitAnd[mt[kk + 1], lowerMask]]; 
mt[kk] = 
              BitXor[mt[kk + (m - n)], IntegerPart[y/2], 
                mag01[BitAnd[y, 1]]];, {kk, n - m, 
              n - 2}
]; 
y = 
            BitOr[BitAnd[mt[n - 1], upperMask], 
              BitAnd[mt[0], lowerMask]]; 
mt[n - 1] = 
            BitXor[mt[m - 1], IntegerPart[y/2], 
              mag01[BitAnd[y, 1]]]; 
mti = 
            0;
]; 
y = 
          mt[mti];  
mti++; 
y = 
          BitXor[y, IntegerPart[y/(2^11)]]; 
        y = BitXor[y, BitAnd[y*2^7, 16^^9d2c5680]]; 
y = 
          BitXor[y, BitAnd[y*2^15, 16^^efc60000]]; 
y = 
          BitXor[y, IntegerPart[y/(2^18)]]
]

initGenrand[100]

Table[genrandInt32[], {1000}]

31ビット乱数の発生

32ビット乱数を加工するので速さは遅い．

genrandInt31[] := IntegerPart[genrandInt32[]/2]

Table[genrandInt31[], {10}]

区間 [0,1] での実数値

genrandReal1[] := genrandInt32[]/4294967295.0

Table[genrandReal1[], {30}]

区間 [0,1) での実数値

genrandReal2[] := genrandInt32[]/4294967296.0

Table[genrandReal2[], {30}]

区間 (0,1) での実数値

genrandReal3[] := (genrandInt32[] + 0.5)/4294967296.0

Table[genrandReal3[], {30}]

区間 (0,1] での実数値

genrandReal4[] := genrandInt32[]/4294967296.0

Table[genrandReal4[], {30}]

53ビット解像度で区間 [0,1) での実数値

genrandRes53[] := (IntegerPart[genrandInt32[]/(2^5)]*67108864.  + 
            IntegerPart[genrandInt32[]/(2^6)])/
        9007199254740992.0

Table[genrandRes53[], {30}]

標準正規乱数

\(\(\(GaussianRandomMT[] := 
      Module[{ran1, ran2}, 
        If[BMFlag, ran1 = RandomMT[]; 
          ran2 = RandomMT[]; 
BMRest = 
            Sqrt[\(-2\)*Log[ran1]]*Cos[2*Pi*ran2];  (*nokori\ wo\ simau*) 
            BMFlag = False; 
Sqrt[\(-2\)*Log[ran1]]*
            Sin[2*Pi*ran2], 
          BMFlag = True; 
BMRest]]\)\(
\)
    \)\)

Table[GaussianRandomMT[], {1000}]

パッケージ化

RandomMT 関数を作ろう．方針としては，確率密度関数を与えるとその分布に従う関数にすること．整数の範囲が与えられたらその数を表示するようにすること．つまり完全に Random 関数から置き換えられるようにすること．種を変えられることなど．

SetDirectory["z:\\mathematica"]

<< RandomMT.m

?initGenrand

initGenrand[100]

?RandomMT

恒常的に使う場合には，C:\Program Files\Wolfram Research\Mathematica\5.2\Configuration\Kernel
などに RandomMT.m ファイルを置いて， init.m の End[] の前に「<< "RandomMT.m"」という行を追加する