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 = (num*A+C)%M;
初期値numに設定値A、C、Mをかけたり足したりして、新たなnum値を作り出す。次回乱数が必要なときは、この「新たなnum値」を元に次の乱数を作る。
今回は以下のようなコード（extends SpriteはFlashCS4のコンパイラバグ対策のため）。
class Mas extends Sprite{
private static var num:Number = 1;
public static function randomSeed(s:Number):void{
num = s;
for (var i:int = 0; i < 100; i++) {
random();
}
}
public static function random():Number{
num = (num*1664525+1013904223)%0x100000000;
return num/0x100000000;
}
}
randomSeedで初期値を設定できる（get/setで作ったり、randomの呼び出し回数をカウントしたほうが使いやすいとかあるかも）。randomSeedを与えた場合には100回実行して偏りを抑える。
A= 1664525,C= 1013904223,M=0×100000000を設定値としている。計算の結果を32Bitの幅に値を収めるために、0×100000000（0xFFFFFFFFの長さ）であまりを出す。返り値として0以上1未満に収めるために0×100000000で割っている。
線形合同法は設定値の組み合わせによって振る舞いが変わる。言語によってバラバラなので、結果、良い線形合同法と悪い線形合同法ができてしまっており、また、randomSeedに対する擬似乱数が言語依存になってしまいがちだ。「線形合同法は駄目」という文章はネット上にたくさんあるが、文脈、どの線形合同法について述べているかを注意したほうが良さそうだ。
今回採用した値の組み合わせ（1664525、1013904223、0×100000000）は、Knuth先生が線形合同法の中では最上とお墨付きを与えたものらしい。どれくらいそれが確かなのかはわからないけど、、（笑）。
numをNumber型にしているが、uintにすると定義上32bit内に収まるので%0×100000000の必要が無くなる。ただし、テストした中では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()*0×100000000とやると結果は偶数の整数になる。疑似乱数の生成方法に由来しているはずだ。
◆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の色は0×000000で生成している。
計算結果により座標が得られると、そこの点が若干明るくなる。次に同じ座標が得られた場合は、さらに明るくなる。0×00から始まり、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()
◆ソースコード
/*
SFMTを含んでいるため、修正BSDライセンスです。
SFMTについては、260行目あたりを参照。
本当は↓からforkにしたかったけど、
http://wonderfl.kayac.com/code/feb94dc51f5dfdeb1aeef3d16b35c70661e6b963
ライセンス的にダメだったので方向転換。
◆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は今回くらいのテスト内容だと、
特にメリットが見いだせなかった。
それとも実装の問題なのだろうか、、、。
＝＝＝＝＝
最近、アルゴリズムの本を読んで気になっていたので、
random()を調べてみた。
コンピュータでは通常、完全にランダムな
値を作ることはできない。
実用上問題無いくらいにバラバラの値を
上手い計算方法で取り出して使っている。
これを疑似乱数（pseudo-random number）
という。
疑似乱数の計算方法はいくつも提案されているが、
その中でも特に使われているのは、
線形合同法のようだ。
今回は以下のような関数。
function random():Number{
num = (num*1664525+1013904223)%0x100000000;
return num/0x100000000;
}
初期値numに1664525、1013904223、0x100000000をかけたり足したりして、新たなnum値を作り出す。次回乱数が必要なときは、この「新たなnum値」を元に次の乱数を作る。
32Bitの幅に値を収めるので、0x100000000で
線形合同法は設定値の組み合わせによって、
振る舞いが変わる。言語によってバラバラなので、
結果、良い線形合同法と悪い線形合同法が
できてしまっている。注意が必要だ。
また、疑似乱数では再現可能なしくみ、
乱数の種（randomSeed）を与えることができる。
疑似乱数はあくまでも計算結果なため、
初期値が同じなら、同じ結果が得られるためだ。
例えば疑似乱数を使っての動作チェック時に、
「同じ組み合わせの疑似乱数で再テストをしたい」
ということはよくある。
この場合、randomSeedをメモしておいて、
再テスト時には、このrandomSeedを再び与えれば良い。
ASでは疑似乱数の生成方法が公開されておらず、
Math.random()にrandomSeedを与えることもできない。
注：
元のbitmapDataの色は0x000000であり、
計算結果により座標が得られると、そこの点が若干明るくなる。
次に同じ座標が得られた場合は、さらに明るくなる。
0x00から始まり、0xFFまで青が段々強くなっていく。
最終段階は目立ちやすいように赤0xFF0000にした。
SFMTについては、260行目あたりを参照。
参考
http://www.eqliquid.com/blog/2009/02/sfmt-actionscript-30.html
http://www001.upp.so-net.ne.jp/isaku/rand.html
http://hp.vector.co.jp/authors/VA004808/random1.html
http://d.hatena.ne.jp/keyword/%C0%FE%B7%C1%B9%E7%C6%B1%CB%A1
*/
package {
import flash.display.Sprite;
import flash.text.TextField;
public class Main extends Sprite {
public var mySfmt:Sfmt = new Sfmt(1);
public var my00DrawBitmap:DrawBitmap = new DrawBitmap(this,0,0);
public var my01DrawBitmap:DrawBitmap = new DrawBitmap(this,0,1);
public var my02DrawBitmap:DrawBitmap = new DrawBitmap(this,0,2);
public function Main():void {
var time:Number = new Date().getTime();
benchMarkj(_a0);
benchMarkj(_a1);
benchMarkj(_a2);
benchMarkj(_a99);
my00DrawBitmap.txt(""+benchMarkj(_a0)/1000+"sec：Math.random");
benchDraw(_b0);
my01DrawBitmap.txt(""+benchMarkj(_a1)/1000+"sec：線形合同法");
benchDraw(_b1);
my02DrawBitmap.txt(""+benchMarkj(_a2)/1000+"sec：SFMT");
benchDraw(_b2);
}
//100万回関数を実行して、かかった時間を計っている。
private function benchMarkj(_fn:Function):int {
var time:Number = (new Date()).getTime();
_fn(1000000);
return (new Date()).getTime() - time;
}
//40万回関数を実行して、ビットマップの色を変える。
private function benchDraw(_fn:Function):void {
_fn(400000);
}
private function _a0(n:uint):void {
for (var i:int = 0; i < n; i++) {
Math.random();
}
}
private function _a1(n:uint):void {
for (var i:int = 0; i < n; i++) {
Mas.random();
}
}
private function _a2(n:uint):void {
for (var i:int = 0; i < n; i++) {
mySfmt.random();
}
}
private function _b0(n:uint):void {
for (var i:int = 0; i < n; i++) {
my00DrawBitmap.pxy(Math.random(),Math.random());
}
}
private function _b1(n:uint):void {
for (var i:int = 0; i < n; i++) {
my01DrawBitmap.pxy(Mas.random(),Mas.random());
}
}
private function _b2(n:uint):void {
for (var i:int = 0; i < n; i++) {
my02DrawBitmap.pxy(mySfmt.random(),mySfmt.random());
}
}
private function _a99(n:uint):void {
for (var i:int = 0; i < n; i++) {
zero();
}
}
private function zero():Number {
return 0;
}
}
}
import flash.display.Sprite;
import flash.display.BitmapData;
import flash.display.Bitmap;
import flash.text.TextField;
class DrawBitmap extends Sprite{
private var w:Number = 465;
private var h:Number = 142;
private var bmp_data:BitmapData = new BitmapData(w,h, false , 0x000000);
private var bmp:Bitmap = new Bitmap(bmp_data);
private var text_field:TextField = new TextField();
function DrawBitmap(ta:Object,_x:Number,_y:Number){
bmp.x = _x*w;
bmp.y = 155*_y+12;
ta.addChild(bmp);
text_field.x = _x*w;
text_field.y = 155*_y-2;
text_field.height=50;
text_field.width = w;
ta.addChild(text_field);
}
public function txt(_str:String):void{
text_field.text=_str;
}
public function pxy(numX:Number,numY:Number):void{
var _x:int = int(numX*w);
var _y:int = int(numY*h);
if(bmp_data.getPixel(_x,_y) < 0x11){
bmp_data.setPixel(_x,_y ,0x11);
}else if(bmp_data.getPixel(_x,_y) < 0x22){
bmp_data.setPixel(_x,_y ,0x22);
}else if(bmp_data.getPixel(_x,_y) < 0x33){
bmp_data.setPixel(_x,_y ,0x33);
}else if(bmp_data.getPixel(_x,_y) < 0x44){
bmp_data.setPixel(_x,_y ,0x44);
}else if(bmp_data.getPixel(_x,_y) < 0x55){
bmp_data.setPixel(_x,_y ,0x55);
}else if(bmp_data.getPixel(_x,_y) < 0x66){
bmp_data.setPixel(_x,_y ,0x66);
}else if(bmp_data.getPixel(_x,_y) < 0x77){
bmp_data.setPixel(_x,_y ,0x77);
}else if(bmp_data.getPixel(_x,_y) < 0x88){
bmp_data.setPixel(_x,_y ,0x88);
}else if(bmp_data.getPixel(_x,_y) < 0x99){
bmp_data.setPixel(_x,_y ,0x99);
}else if(bmp_data.getPixel(_x,_y) < 0xAA){
bmp_data.setPixel(_x,_y ,0xAA);
}else if(bmp_data.getPixel(_x,_y) < 0xBB){
bmp_data.setPixel(_x,_y ,0xBB);
}else if(bmp_data.getPixel(_x,_y) < 0xCC){
bmp_data.setPixel(_x,_y ,0xCC);
}else if(bmp_data.getPixel(_x,_y) < 0xDD){
bmp_data.setPixel(_x,_y ,0xDD);
}else if(bmp_data.getPixel(_x,_y) < 0xEE){
bmp_data.setPixel(_x,_y ,0xEE);
}else if(bmp_data.getPixel(_x,_y) < 0xFF){
bmp_data.setPixel(_x,_y ,0xFF);
}else{
bmp_data.setPixel(_x,_y ,0xFF0000);
}
}
}
/*
線形合同法
randomSeedを与えた場合は、100回実行して偏りを抑える。
参考
http://hp.vector.co.jp/authors/VA004808/random1.html
http://d.hatena.ne.jp/keyword/%C0%FE%B7%C1%B9%E7%C6%B1%CB%A1
*/
class Mas extends Sprite{
private static var num:Number = 1;
public static function randomSeed(s:Number):void{
num = s;
for (var i:int = 0; i < 100; i++) {
random();
}
}
public static function random():Number{
num = (num*1664525+1013904223)%0x100000000;
return num/0x100000000;
}
}
/*
SFMT 簡易版クラス
author:kimura@sionnet.jp
http://www.eqliquid.com/blog/2009/02/sfmt-actionscript-30.html
から抜き出した。
umhrにより使いやすいようにごく一部修正した。
umhrによる修正や使い方によって、
重大な問題が発生している可能性があるので、
「SFMT 簡易版クラス」やSFMT自体の評価は別途行うこと。
以下、元のReadme_jp.txtよりコピペ
====================
SFMT 簡易版クラス
author:kimura@sionnet.jp
パッケージ:
jp.sionnet.utils.Sfmt
====================
概要:
要約すれば、高度な擬似乱数を発生するメルセンヌ・ツイスタの高速版がSFMTです。
このSFMTの各言語版を公開されているisakuさんの簡易版ライブラリをActionScript3.0に移植したものです。
内容はちっとも理解できてないので、間違いがあったらご指摘いただけると幸いです。
SFMT:
http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/index-jp.html
乱数ライブラリ:
http://www001.upp.so-net.ne.jp/isaku/rand2.html
ライセンス:
SFMTのライセンスに依ります。2009年2月現在、修正BSDライセンスです。
====================
以上、元のReadme_jp.txtよりコピペ
*/
class Sfmt
{
private var index:int;
private var x:Vector.<int> = new Vector.<int>(624);//状態テーブル
/**
* コンストラクタ<br />
* 初期化する
* @param s 整数のシード
*/
public function Sfmt(s:int)
{
initMt(s);
}
private function genRandAll():void
{
var a:int = 0;
var b:int = 488;
var c:int = 616;
var d:int = 620;
var y:int;
var p:Vector.<int> = x;//Javaの振る舞い、処理の内容からシャローコピーと判断
do
{
y = p[a + 3] ^ (p[a + 3] << 
^ (p[a + 2] >>> 24) ^ ((p[b + 3] >>> 11) & 0xbffffff6);
p[a + 3] = y ^ (p[c 1="3" language="+"][/c] >>> ^ (p[d + 3] << 18);
y = p[a + 2] ^ (p[a + 2] << ^ (p[a + 1] >>> 24) ^ ((p[b + 2] >>> 11) & 0xbffaffff);
p[a + 2] = y ^ ((p[c 1="2" language="+"][/c] >>> | (p[c 1="3" language="+"][/c] << 24)) ^ (p[d + 2] << 18);
y = p[a + 1] ^ (p[a + 1] << ^ (p[a] >>> 24) ^ ((p[b + 1] >>> 11) & 0xddfecb7f);
p[a + 1] = y ^ ((p[c 1="1" language="+"][/c] >>> | (p[c 1="2" language="+"][/c] << 24)) ^ (p[d + 1] << 18);
y = p[a] ^ (p[a] << ^ ((p[b] >>> 11) & 0xdfffffef);
p[a] = y ^ ((p[c][/c] >>> | (p[c 1="1" language="+"][/c][/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);
}
}
