2015年7月17日金曜日

乱数を作ろう

最近のマイブーム

最近、SFCのすーぱーぷよぷよ通Remixにハマってます。今更ながら、ですが(笑)。


さて、これ遊びながら思ったんですが、あまりにもぷよの出現が偏りがあって、

「あんま良い乱数使ってないんじゃないの?」

とかぶつくさボヤいてたんですが(笑)。
しかし待てよ、と。どういう風に乱数プログラミングしてたんだろ?

※ もっとも、良い乱数使ったから、と言って、このテの落ち物パズルのゲーム性が上がる、とは限らないんですがね(笑)。

機械語には乱数は無い

ここんとこ、色々古いCPU、例えば6502とかZ80の資料を眺めていたんですが。
基本的には(当たり前かもしれないですが)機械語では「擬似乱数列」ってのは用意されてないんですね。
スーファミのCPUは65816と言うApple II GSで使用されていた、あまりメジャーじゃないCPUが使われています。
90年代初頭のスーファミの時代、ゲームプログラマはアセンブリ言語でプログラムする事が多かった模様です。そして当時は(ある程度は用意されてたみたいですが)開発ツールなんてのも今のようにゴージャスじゃなかった模様です。
と言う事は、かなりの確率で、例えばアクションゲームとかシューティングならさておき、RPG等、自作のアセンブリで書いた乱数ライブラリを使ってたりしたんじゃないでしょうか。そうじゃないと16bit CPU程度ではスピードはあんま稼げなかった筈なんですよね。
擬似乱数。果たしてどうやってプログラムするんでしょう。

SICPの乱数

話は唐突に飛びますが。
計算機プログラムの構造と解釈でも中途半端に(笑)、擬似乱数の作り方が載っています。


   「ランダムに選ばれた」の意味は明瞭ではない. おそらく望んでいるのは randの次々の呼び出しは一様分布という統計の性質を持つ数列を生成することである. 適切な数列の発生法はここでは論じない. そうではなく, 手続きrand-updateがあり, ある数x1から出発し,
x2 = (rand-update x1)x3 = (rand-update x2)
を作ると, 値の列x1, x2, x3, ... は望みの統計的性質を持つと仮定する.
   randをある固定した値random-initに初期化される局所状態変数xを持った手続きとして実装出来る. randを呼び出す度にxの現在値のrand-updateを計算し, これを乱数として返し, またそれをxの新しい値とする.
(define rand
  (let ((x random-init))
    (lambda ()
      (set! x (rand-update x))
      x)))

   もちろん同じ乱数の列を単にrand-updateを直接呼び出すことで, 代入なしに生成することが出来る. しかしプログラムの乱数を使う部分は rand-updateに引数として渡すxの値をしっかり覚えていなければならないことを意味する. 
そして注釈には次のような記述があります。

rand-updateを実装する通常のやり方は, a, b, mを適切に選んだ整数とし, mを法としてxをax + bで更新する規則を使うことである. 

要するに具体的な例は全く紹介されてないんですね(笑)。
そこでWikipediaを頼ったコードは次のようなモノです。
;;; Wikipediaの作例によるrand-update
(define (rand-update x)
(let ((a 3) (b 5) (m 13))
(modulo (+ (* a x) b) m)))
;;; 乱数の初期値も Wikipedia に倣って8とする
(define random-init 8)
;;; SICP の rand
(define rand
(let ((x random-init))
(lambda ()
(set! x (rand-update x))
x)))
view raw rand-update.rkt hosted with ❤ by GitHub


これは線形合同法と言う極めて単純なやり方なんですが。ところが精度がメチャクチャ悪いんですね。
SICPではこの後、円周率πをモンテカルロ・シミュレーションで求めよう、ってトピックになるわけですが。

(define (estimate-pi trials)
(sqrt (/ 6 (monte-carlo trials cesaro-test))))
(define (cesaro-test)
(= (gcd (rand) (rand)) 1))
(define (monte-carlo trials experiment)
(define (iter trials-remaining trials-passed)
(cond ((= trials-remaining 0)
(/ trials-passed trials))
((experiment)
(iter (- trials-remaining 1) (+ trials-passed 1)))
(else
(iter (- trials-remaining 1) trials-passed))))
(iter trials 0))
view raw monte-carlo.rkt hosted with ❤ by GitHub


ところが、これやってみりゃ分かるんですが、円周率πに届かないんですよ(笑)。

> (estimate-pi 10000)
2.449489742783178
view raw monte-carlo-fact.rkt hosted with ❤ by GitHub


さては、この結果知ってて「乱数作成」をサラーっと流したんじゃねぇの(笑)。
もう、SICP、こんなんばっかだよ(苦笑)。

Schemeと乱数

ところで。
実はSchemeって仕様には乱数が含まれてないんですね。
まあ、もちろん、各実装には実装依存で乱数生成プログラムが入ってるんでしょうが、あくまで「仕様としては」乱数は持ってません。
そう言う意味では結構メンドくせぇプログラミング言語なんですよね(笑)。

最強の乱数生成法

統計処理言語Rやお馴染みのPythonではメルセンヌ・ツイスタと言われる乱数生成プログラムが組み込まれています。メルセンヌ・ツイスタは今のトコ、史上最強と言われる擬似乱数生成プログラムで、仕様がない言語で最強の乱数列が使い放題なんですよねぇ。
一方Lispはどうか、と言うと例えば、Common Lispなんかにはメルセンヌ・ツイスタのライブラリが公開されていて、共有されています。
Schemeは・・・・・・まあ、実装次第なんでしょうねぇ。Gaucheなんかはメルセンヌ・ツイスタを使ってたと思うんですが、他の実装だと保証の限りじゃないです。まあ、こう言う時、「ライブラリを共有しづらい」Schemeだとやりにくいですね(R7RSでだいぶ仕様上は改善されていますが)。
さて、Racketはどうでしょうか?

Racketの乱数

残念ながらRacketの乱数はメルセンヌ・ツイスタではございません(笑)。RacketのDocumentationによると、L’Ecuyer’s MRG32k3a algorithmと言うアルゴリズムを用いてるそうです。
何だそりゃ。初めて聞きました。
実はCのソースコードが公開されていました。

/*
32-bits Random number generator U(0,1): MRG32k3a
Author: Pierre L'Ecuyer,
Source: Good Parameter Sets for Combined Multiple Recursive Random
Number Generators,
Shorter version in Operations Research,
47, 1 (1999), 159--164.
---------------------------------------------------------
*/
#include "MRG32k3a.h"
#define norm 2.328306549295728e-10
#define m1 4294967087.0
#define m2 4294944443.0
#define a12 1403580.0
#define a13n 810728.0
#define a21 527612.0
#define a23n 1370589.0
/***
The seeds for s10, s11, s12 must be integers in [0, m1 - 1] and not all 0.
The seeds for s20, s21, s22 must be integers in [0, m2 - 1] and not all 0.
***/
#define SEED 12345
static double s10 = SEED, s11 = SEED, s12 = SEED,
s20 = SEED, s21 = SEED, s22 = SEED;
double MRG32k3a (void)
{
long k;
double p1, p2;
/* Component 1 */
p1 = a12 * s11 - a13n * s10;
k = p1 / m1;
p1 -= k * m1;
if (p1 < 0.0)
p1 += m1;
s10 = s11;
s11 = s12;
s12 = p1;
/* Component 2 */
p2 = a21 * s22 - a23n * s20;
k = p2 / m2;
p2 -= k * m2;
if (p2 < 0.0)
p2 += m2;
s20 = s21;
s21 = s22;
s22 = p2;
/* Combination */
if (p1 <= p2)
return ((p1 - p2 + m1) * norm);
else
return ((p1 - p2) * norm);
}
view raw MRG32k3a.c hosted with ❤ by GitHub

これって結構シンプルじゃね?
さぁ、やっと本題です(笑)。
今回はこのMRG32k3aと言うアルゴリズムをRacketで「遅延評価を使って」実装してみたいと思います。SICPの遅延評価の乱数ストリームの仇をここで取ります(笑)。

遅延評価の友

ところで。
Schemeでは仕様上、遅延評価に絡んだ機能はdelayとforceしか提供していません。
これだけだと甚だメンド臭く、「マトモに遅延評価をさせる」には、色々基本的な遅延評価用の手続きを実装しないといけません。事実、SICPの遅延評価の項目だとそこに結構な量が割かれているんですよね。
それじゃああんまりにもメンドくせえだろ、って言う貴方に朗報です(笑)。
実は分かりづらいんですが、SRFI-41ってのがあって、これを呼び出せば遅延評価の基本機構は使い放題、となります。
主なトコではstream-cons、stream-car、stream-cdr、stream-null、stream-map、stream-filter、stream-constant、stream-ref、stream-from等等等、SICP程度のお題ならこの程度で大体書けるようになってるライブラリです。
便利でこの辺考えなくって良いんで、これを使いましょう。

※ 何故SRFI-41が「分かりづらい」かと言うと、遅延評価、つまりlazy-evaluationとして題名が付いてなくって「stream」と言う分かりづらい用語になってるからに他ならない。これは歴史的要因で、今ではlazy-evaluationと言う呼び方の方が妥当になってるが、過去にはそのままstreamと呼ばれてた事に由来する。

実装方針

最終的には単一のプロシージャにまとめますが、途中までは色々バラバラのパーツとして組み立てて行こうと思います。じゃないと遅延評価って良くわかりませんからね。

大域変数


まずは大域変数をCコードに則って定義していきます。まあ、最終的にはこれらはletで局所変数化されるんですが、最初はまあ、Cのコードそのままで書いていきます。

(define norm 2.328306549295728e-10)
(define m1 4294967087.0)
(define m2 4294944443.0)
(define a12 1403580.0)
(define a13n 810728.0)
(define a21 527612.0)
(define a23n 1370589.0)
(define SEED 12345)
view raw Global Variables for MRG32k3a.rkt hosted with ❤ by GitHub

この辺はまあ、いいですよね。乱数の種、SEEDも12345と言うフザけた(笑)値が与えられています。
ちなみに、実はm1、m2と言う2つの値は2の32乗から下2つの素数の模様です。

Component

さて、Cソースのコメントで、取り敢えずComponent 1と書かれたところに着目してみます。

/* Component 1 */
p1 = a12 * s11 - a13n * s10;
k = p1 / m1;
p1 -= k * m1;
if (p1 < 0.0)
p1 += m1;
s10 = s11;
s11 = s12;
s12 = p1;
view raw Component 1.c hosted with ❤ by GitHub


さすがC、Lisp慣れしてると物凄くゴチャゴチャして見えますが(笑)。
原因の一つってのはここの箇所って実は「2つの事柄が1つにまとまってる」からなんですよね。

ここで行われてるのは次の2つです。


  1. p1 なる値の計算
  2. s10、s11、s12と言う3つの数の更新

当然Lisperとしては(笑)、「2つの計算目的は2つのプロシージャで書くべきだ」と言う原則を守ろうと思います。
さて、では手始めに2番から行きますか。
実はさっき定義しなかったんですが、原版ではs10、s11、s12と言う3つの値も大域変数として定義されています。
しかし「遅延評価を持ったLisp」であるSchemeではこれはやるべきじゃないトコなんですよね。
つまり、「各値を更新する」んじゃなくって、s10、s11、s12、そしてその後にはs13、s14・・・って「続いていくだろう」無限ストリームがある、って考えた方がスマートなわけです。
実際は、元のCコードを見ると、

(SEED SEED SEED p1 ....)

と言うようなカタチの無限ストリームを想定出来そうだ、って事ですよね。そしてp1以降は適当なプロシージャ、具体的には1の「計算によって」得られた数値がハマっていく。
そう言う想定から言うとまずは次のような「ストリーム生成」プロシージャをでっち上げられます。
(define (stream1 s10 s11 s12)
(stream-cons s10 (stream1 s11 s12 (Component1 s10 s11))))
view raw stream1.rkt hosted with ❤ by GitHub




そしてp1を計算するComponent1なるプロシージャを設計すれば良い。オリジナルのCコードだと次の部分ですね。
p1 = a12 * s11 - a13n * s10;
k = p1 / m1;
p1 -= k * m1;
if (p1 < 0.0)
p1 += m1;
view raw Part of Component1.c hosted with ❤ by GitHub



ちなみに、最初、この2行目と3行目で何やってんだかサッパリ分かんなかったんですけど、何とこれ剰余計算やってるみたいです(笑)。え"え"え"え"え"、とか思ったんですが(苦笑)。うげぇ(笑)。
何で%使わねぇのかな。不思議です。ま、いいや(苦笑)。
つまり、この部分をSchemeで書くとこうなる、って事ですね。

(define (Component1 s10 s11)
(let ((p1 (modulo (- (* a12 s11) (* a13n s10)) m1)))
(if (negative? p1)
(+ p1 m1)
p1)))
view raw Component1.rkt hosted with ❤ by GitHub



stream1とComponent1を組み合わせると、例えば実行結果は次のようになります。
> (define s1 (stream1 SEED SEED SEED))
> s1
#<stream>
> (stream->list 10 s1)
'(12345
12345
12345
3023790853.0
3023790853.0
3385359573.0
1322208174.0
2930192941.0
2462079208.0
2386811717.0)
view raw test of Component1.rkt hosted with ❤ by GitHub



stream1の最初と2番目と3番目の値はSEEDの値、4番目、5番目の値はそれぞれ1番目と2番目、2番目と3番目の値を使って計算してるから同じ値、6番目からガラッと変わってきますね。
これはstream->listと言うプロシージャでアタマ10個分だけ「実体化」させたわけですが、これがずーっと「無限ストリーム」として続いていくわけです。

Component2

基本ロジックはComponent1と全く同じなんでサクッと書いていきます。
(define (stream2 s20 s21 s22)
(stream-cons s20 (stream2 s21 s22 (Component2 s20 s22))))
(define (Component2 s20 s22)
(let ((p2 (modulo (- (* a21 s22) (* a23n s20)) m2)))
(if (negative? p2)
(+ p2 m2)
p2)))
view raw Component2.rkt hosted with ❤ by GitHub



ただしComponent1とは使ってる引数の値と引数の順番が若干違うんで引っかからないように(笑)。

Combination

現時点ストリームが2つあります。p1、p2はそれぞれのストリームの4番目からスタートして、条件に従ってどっちかを選択してちょこちょこと計算して新しいストリームに組み込まれるわけです。
この辺をstream-mapを使って上手いこと実装します。
(define (Combination stream1 stream2)
(let ((stream1 (stream-cdr (stream-cdr (stream-cdr stream1))))
(stream2 (stream-cdr (stream-cdr (stream-cdr stream2)))))
(stream-map (lambda (p1 p2)
(* (- p1 (if (<= p1 p2)
(- p2 m1)
p2)) norm))
stream1 stream2)))
view raw Combination.rkt hosted with ❤ by GitHub


Combinationは新しいストリーム、ここでは乱数列を返します。
ちょっとまた、アタマ10個くらいstream->listで実体化させて見てみます。

> (define s (Combination (stream1 SEED SEED SEED)
(stream2 SEED SEED SEED)))
> s
#<stream>
> (stream->list 10 s)
'(0.12701112204657714
0.3185275653967945
0.3091860155832701
0.8258468629271136
0.2216299157820229
0.5333953879182788
0.4807742033156181
0.3555598794381262
0.13598841039594017
0.7558522371615436)
view raw test for Combination.rkt hosted with ❤ by GitHub


おお、上手く動いてるようですね。
例えば乱数列の100番目、1000番目、10000番目はそれぞれ

> (stream-ref s 100)
0.07967177768510994
> (stream-ref s 1000)
0.9014158012565428
> (stream-ref s 10000)
0.1110427223837204
view raw test for Combination 2.rkt hosted with ❤ by GitHub


となります。
おお、面白い(笑)。これで(基本的に)実装は完了です。

ソースコード

最後に各パーツを局所手続き化して全コードをまとめます。

#lang racket
(require srfi/41)
;; 32-bits Random number generator U(0,1): MRG32k3a
;; Author: Pierre L'Ecuyer,
;; Source: Good Parameter Sets for Combined Multiple Recursive Random
;; Number Generators,
;; Shorter version in Operations Research,
;; 47, 1 (1999), 159--164.
;; Ported to Scheme by Cametan
;; ---------------------------------------------------------
(define (make-MRG32k3a)
(let ((norm 2.328306549295728e-10)
(m1 4294967087.0)
(m2 4294944443.0)
(a12 1403580.0)
(a13n 810728.0)
(a21 527612.0)
(a23n 1370589.0)
(SEED 12345))
(letrec ((stream1
(lambda (s10 s11 s12)
(stream-cons s10 (stream1 s11 s12 (Component1 s10 s11)))))
(Component1
(lambda (s10 s11)
(let ((p1 (modulo (- (* a12 s11) (* a13n s10)) m1)))
(if (negative? p1)
(+ p1 m1)
p1))))
(stream2
(lambda (s20 s21 s22)
(stream-cons s20 (stream2 s21 s22 (Component2 s20 s22)))))
(Component2
(lambda (s20 s22)
(let ((p2 (modulo (- (* a21 s22) (* a23n s20)) m2)))
(if (negative? p2)
(+ p2 m2)
p2))))
(Combination
(lambda (stream1 stream2)
(let ((stream1 (stream-cdr (stream-cdr (stream-cdr stream1))))
(stream2 (stream-cdr (stream-cdr (stream-cdr stream2)))))
(stream-map (lambda (p1 p2)
(* (- p1 (if (<= p1 p2)
(- p2 m1)
p2)) norm))
stream1 stream2)))))
(Combination (stream1 SEED SEED SEED)
(stream2 SEED SEED SEED)))))
view raw make-MRG32k3a.rkt hosted with ❤ by GitHub


投稿者 cametan_001 時刻: 15:20 0 コメント Translate
ラベル: , , , ,
登録: 投稿 (Atom)