random()、線形合同法


線形合同法アルゴリズムの本を読んで気になっていたので、random()を調べてみた。

コンピュータでは通常、完全にランダムな値を作ることはできない。実用上問題無いくらいにバラバラの値を上手い計算方法で取り出して使っている。これを疑似乱数(pseudo-random number)という。

疑似乱数の計算方法はいくつも提案されているが、その中でも特に使われているのは、線形合同法のようだ。randomにも良し悪しがある、というのは聞いてはいたがActionScriptのMath.random()が悪いほうなのかどうかが気になっていた。

今回はさらにSFMTも並べて計算速度、偏り具合を見てみた。
 
 

●Math.random()と線形合同法とSFMTの三つをを比較してみた。

 
それぞれを実行時間(sec)を計るために100万回実行し、取り出した値の偏りを視覚的に確認するために、40万の点の座標をx,yで作り画像に埋めていった。計算結果に偏りがある場合はムラやパターン、筋が現れるはず。

◆結果

MacOSX10.5.7/2.4GHzCore2Duo、FlashPlayer10,0,22,87では
0.587sec:Math.random
0.037sec:線形合同法
0.515sec:SFMT
確認画像では、三つともほぼ同じ程度のムラとなった。

◆結論

理屈で言えば、完全に均一な面ができればいいのだろうが、そうはなっていない。しかしこれくらいであれば実用上問題になることはまず無いだろう。

・Math.randomについて
速度が遅くて困るとかは無いと思う。ほとんどの場合これで足りる、ということを確認できた。

・線形合同法について
randomSeedを与えられる、計算が速いという理由で、ひと手間かけてもいい場合には、この線形合同法が良い。

・SFMTについて
SFMTは今回くらいのテスト内容だと、特にメリットが見いだせなかった。それとも実装の問題なのだろうか。
 
 

●さらに詳しい話。

 

◆線形合同法(LCG:Linear congruential generator)

線形合同法の一般式は以下。

初期値numに設定値A、C、Mをかけたり足したりして、新たなnum値を作り出す。次回乱数が必要なときは、この「新たなnum値」を元に次の乱数を作る。

今回は以下のようなコード(extends SpriteはFlashCS4のコンパイラバグ対策のため)。

randomSeedで初期値を設定できる(get/setで作ったり、randomの呼び出し回数をカウントしたほうが使いやすいとかあるかも)。randomSeedを与えた場合には100回実行して偏りを抑える。

A= 1664525,C= 1013904223,M=0x100000000を設定値としている。計算の結果を32Bitの幅に値を収めるために、0x100000000(0xFFFFFFFFの長さ)であまりを出す。返り値として0以上1未満に収めるために0x100000000で割っている。

線形合同法は設定値の組み合わせによって振る舞いが変わる。言語によってバラバラなので、結果、良い線形合同法と悪い線形合同法ができてしまっており、また、randomSeedに対する擬似乱数が言語依存になってしまいがちだ。「線形合同法は駄目」という文章はネット上にたくさんあるが、文脈、どの線形合同法について述べているかを注意したほうが良さそうだ。

今回採用した値の組み合わせ(1664525、1013904223、0x100000000)は、Knuth先生が線形合同法の中では最上とお墨付きを与えたものらしい。どれくらいそれが確かなのかはわからないけど、、(笑)。

numをNumber型にしているが、uintにすると定義上32bit内に収まるので%0x100000000の必要が無くなる。ただし、テストした中ではuintにすると速度が遅かったのでNumber型にしてある。

◆randomSeed

多くの疑似乱数では再現可能なしくみ、randomSeed(乱数の種)を初期値として与えることができる。疑似乱数はあくまでも計算結果なため、初期値が同じなら同じ結果が得られるためだ。

例えば疑似乱数を使っての動作チェック時に、「同じ組み合わせの疑似乱数で再テストをしたい」ということはよくある。この場合、randomSeedをメモしておいて、再テスト時には、このrandomSeedを再び与えれば良い。ゲームのリプレイではrandomSeedを保っておけば、敵の出現率、モブの動きなどrandomに頼った部分も再現できる。

ナムコの「ドルアーガの塔」というゲームでは、階数をrandomSeedとした擬似乱数で迷路を作ったという。よって、異なるマシンへの移植の際には異なった迷路が生成されたらしい。
http://druaga.to/qanda_tod.html

◆AS3のMath.random()

ActionScriptでrandomの計算方法が非公開、randomSeed設定不能なのは恐らく、将来も含めていろいろなCPUで動作させること前提にしているので、CPU依存の仕様は避けているためだろう。ちなみに、AS3でMath.random()*0x100000000とやると結果は偶数の整数になる。疑似乱数の生成方法に由来しているはずだ。

◆SFMT(SIMD-oriented Fast Mersenne Twister)

メルセンヌ・ツイスタというかっこいい名前の擬似乱数用の計算式があって、それの高速版がSFMT。
SFMTの性質、要点を解説できるほど知らないので、興味がある人は解説サイトを見てほしい。

今回は「SFMT 簡易版クラス」から抜き出した。
http://www.eqliquid.com/blog/2009/02/sfmt-actionscript-30.html
wonderflで使いやすいように一部修正した。
この修正や使い方によって、重大な問題が発生している可能性があるので、「SFMT 簡易版クラス」やSFMT自体の評価は別途行うこと。
また、このコードは修正BSDライセンスなので、注意。

◆確認画像について

元のbitmapDataの色は0x000000で生成している。
計算結果により座標が得られると、そこの点が若干明るくなる。次に同じ座標が得られた場合は、さらに明るくなる。0x00から始まり、0xFFまで青が段々強くなっていく。
最終段階は目立ちやすいように赤0xFF0000にした。

◆お勧めリンク

Knuth先生お墨付きの話はこちらから
http://hp.vector.co.jp/authors/VA004808/random1.html

線形合同法 – wikipedia
http://ja.wikipedia.org/wiki/線形合同法

擬似乱数 – wikipedia
http://ja.wikipedia.org/wiki/擬似乱数

線形合同法とは – はてなキーワー
http://d.hatena.ne.jp/keyword/%C0%FE%B7%C1%B9%E7%C6%B1%CB%A1
http://en.wikipedia.org/wiki/Linear_congruential_generator

良い乱数・悪い乱数
http://www001.upp.so-net.ne.jp/isaku/rand.html

flash.display.BitmapDataのnoise()には乱数の種を指定
http://d.hatena.ne.jp/flashrod/20070223

SFMT
http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/index-jp.html

SFMT 簡易版クラス
http://www.eqliquid.com/blog/2009/02/sfmt-actionscript-30.html

メルセンヌ・ツイスタ
http://ja.wikipedia.org/wiki/メルセンヌ・ツイスタ

メルセンヌ・ツイスタ ActionScript 3.0に移植
http://onegame.bona.jp/tips/mersennetwister.html

Math.random() – ActionScript 3.0 言語およびコンポーネントリファレンス
http://help.adobe.com/ja_JP/AS3LCR/Flash_10.0/Math.html#random()

◆ソースコード

>>> 8) ^ (p[d + 3] << 18);
y = p[a + 2] ^ (p[a + 2] << 8) ^ (p[a + 1] >>> 24) ^ ((p[b + 2] >>> 11) & 0xbffaffff);
p[a + 2] = y ^ ((p

>>> 8) | (p

<< 24)) ^ (p[d + 2] << 18);
y = p[a + 1] ^ (p[a + 1] << 8) ^ (p[a] >>> 24) ^ ((p[b + 1] >>> 11) & 0xddfecb7f);
p[a + 1] = y ^ ((p

>>> 8) | (p

<< 24)) ^ (p[d + 1] << 18);
y = p[a] ^ (p[a] << 8) ^ ((p[b] >>> 11) & 0xdfffffef);
p[a] = y ^ ((p

>>> 8) | (p

[/c] << 24)) ^ (p[d] << 18);
c = d;
d = a;
a += 4;
b += 4;
if (b == 624) b = 0;
}
while (a != 624);
}

private function periodCertification():void
{
var work:int;
var inner:int = 0;
var i:int;
var j:int;
var parity:Vector.<int> = new Vector.<int>(4);
parity.push(0x00000001);
parity.push(0x00000000);
parity.push(0x00000000);
parity.push(0x13c9e684);
index = 624;

for (i = 0; i < 4; i++)
{
inner ^= x[i] & parity[i];
}
for (i = 16; i > 0; i >>>= 1)
{
inner ^= inner >>> i;
}
inner &= 1;
if (inner == 1) return;
for (i = 0; i < 4; i++)
{
for (j = 0, work = 1; j < 32; j++, work <<= 1 )
{
if ((work & parity[i]) != 0)
{
x[i] ^= work;
return;
}
}
}

}

/**
* 初期化
* @param s 整数のシード
*/
private function initMt(s:int):void
{
x[0] = s;
for (var p:int; p < 624; p++)
{
s = 1812433253 * (s ^ (s >>> 30)) + p;
x[p] = s;
}
periodCertification();
}

/**
* 32ビット符号あり整数乱数を返す
* @return 32ビット符号あり整数乱数
*/
public function nextMt():int
{
if (index == 624)
{
genRandAll();
index = 0;
}
return x[index++];
}

/**
* 0 以上 n 未満の整数乱数を返す
* @param n 乱数の上限値にする整数
* @return 0 以上 n 未満の整数乱数
*/
public function nextInt(n:int):int
{
var z:Number = nextMt();
if (z < 0) z += 4294967296.0;
return int(n * (1.0 / 4294967296.0) * z);
}

/**
* umhrによる追加
* 0 以上 1 未満の乱数(32ビット精度)を返す
* @return 0 以上 1 未満の乱数(32ビット精度)
*/
public function random():Number
{
var z:Number = nextMt();
if (z < 0){ z += 4294967296};
return z / 4294967296;
}

/**
* 0 以上 1 未満の乱数(53ビット精度)を返す
* @return 0 以上 1 未満の乱数(53ビット精度)
*/
public function nextUnif():Number
{
var z:Number = nextMt() >>> 11;
var y:Number = nextMt();
if (y < 0) y += 4294967296.0;
return (y * 2097152.0 + z) * (1.0 / 9007199254740992.0);
}

}