おちゃめくらぶ掲示板 - 1659321896

975：御茶目菜子：2012/04/27(金) 15:14:27: プチコンで疑似乱数を作る
プチコンには標準でRND関数が備わっておりこれによってA=RND (10)で0～9の乱数を発生
させることが可能にょ。
普通に使う場合には全く問題ないけどゲームで使用する場合には同じ乱数列の乱数を発生
させることができないという問題があるにょ。
つまり、そのゲームと同じものを再現してプレイができないということにょ。
多くのBASICではRNDOMIZEの引数によって同じ乱数列の乱数を発生可能だし、ポケコンでは
シャープの機種においてはRND（負の数）で同じ乱数列の乱数を発生させることが可能になって
いるにょ。
プチコンの乱数でそういうものがないならば作ればいいだけの話にょ。

一般的なコンピュータによって生成される乱数は計算式によって導かれる疑似乱数になって
いるにょ。
その中で最もポピュラーなのが線形合同法によるものにょ。

　X(N+1)=(A*X(N)+B) mod M

このA、B、Mの組み合わせによって乱数の発生の仕方が最大値や周期が変わってくるにょ。
ここでMが乱数列の最大周期と最大値になるにょ。（実際は剰余であるため最大値はM-1）
例えば X=(X*3+2)%64 という式で考えてみるにょ。
まず使うためにはこの乱数の初期値を決めないといけないにょ。
この初期値によって乱数の発生の仕方が変わってくるけど逆にいえばこの初期値が同じならば
同じ値となるためゲームなどにおいて再現性を持たせることが可能になるにょ。

《　リスト1　》
　ACLS
　X=0
　FOR I=1 TO 64
　　GOSUB @RND
　　? X;
　NEXT
　END
@RND
　X=(X*3+2)%64
　RETURN

上記公式におけるMの値が64ということで64通りの乱数の発生を期待してみたけどよく見ると
同じものが4回タブっているにょ。
つまり、この式は周期16の乱数ということにょ。

次にこの乱数を使って画面上にドットを表示してみるにょ。

《　リスト2　》
　ACLS
　X=0
　FOR I=1 TO 64
　　GOSUB @RND:A=X
　　GOSUB @RND:B=X
　　GPSET A,B,15
　NEXT
　END
@RND
　X=(X*3+2)%64
　RETURN

1周期16の乱数であるため16個のドットが表示されることを期待してみるけど画面上には8個しか
点が打たれていないにょ。

実は線形合同法には以下の2つの問題点があるにょ。
(1) 上位の桁はランダム性が高く下位の桁はランダム性が低い
(2) いくつか組にして使えばランダムではなくなる

(1)リスト1において発生した乱数をよく見ると全部偶数になっているにょ。
線形合同法における乱数は下位の桁は全部奇数、全部偶数、奇数と偶数を交互のいずれかに
なっているにょ。
2進数で表した時には下位Kビットの最大周期は2^Kになるにょ。
下位1ビットだと最大周期が2であるため良くて偶数と奇数が交互になり、悪ければ偶数のみ
しか出なくなってしまうにょ。
下位3ビットを使用した場合は最大周期が8になってしまうため実用性はほとんどない乱数に
なってしまうにょ。
したがって、基本的に線形合同法の乱数は上位から使うのが望ましいにょ。

この線形合同法の問題は分かっていれば回避が可能（発生した乱数を生の状態で下位の
ビットから使うことを避ければ良い）なのだけど市販のゲームでも回避されなかったもの
さえあるからね。
Xbox360用ソフトのカルドセプトサーガはダイスの目が偶数と奇数が交互に出るという
この線形合同法の問題をそのまま抱えた状態で発売されているにょ。（2人でプレイする
場合は片方のプレイヤーはダイスの目が必ず偶数、もう片方のプレイヤーは必ず奇数になって
しまっていた模様）
これは実際にダイスの目という表面上ですぐに分かるものであるため十分なテストプレイが
行われていたら簡単に発覚可能だったけど内部計算に乱数を用いる場合はこのように奇数と
偶数が交互に出ても気づかないにょ。
乱数による自動マップ生成などならばそれほど問題にはならないけれど周期が短いことで
パターンが固定化されてしまうため下位ビットから使うというのはあまり望ましくはないにょ。

(2)これはリスト2の実行結果を見ればよく分かるにょ。
周期16の乱数でも2つ組み合わせれば実際の周期は半減してしまうからね。（三次元座標に
用いた場合にはさらに周期が短くなってしまう）
そのため必要な周期の数倍の周期の乱数が求められてくるにょ。
では、どの程度必要かというと具体的にはジャンルによって変わるので難しいにょ。
ゲームではなく科学技術シミュレーションならば億単位、兆単位の周期でないと使い物には
ならないけどプチコンではそのような用途には用いることはほぼないにょ。
というのも扱える数字の上限が524287だからね。
これは32ビット固定小数点を採用しており、符号1ビット、小数部12ビットで残りの19ビットが
整数部であるためにょ。

メインルーチンがVSYNC1（60fps）で実行するゲームだと30分のプレイで最低でも108000回の
乱数発生が行われるにょ。（敵キャラ8体に乱数を使えばその8倍になる）
そうすると少なくとも周期は10万単位くらいは欲しいところにょ。
また、プチコンにおいては標準のRND関数によってランダムに表示したドットで画面を埋め
尽くすことが可能になっているにょ。

《　リスト3　》
　ACLS
　FOR I=1 TO 524287
　　GPSET RND(256),RND(192),15
　NEXT

プチコン（というかDS）の液晶画面解像度は256x192の49152ドットで構成されているため
すべてのドットを埋め尽くす期待値を計算すると559345回のループが必要にょ。
X座標、Y座標は独立しているため期待値から考えると最低でも100万回程度の周期が必要と
いえそうにょ。
つまり、「RNDの変わりに使えるレベル」となると周期は100万回以上が合格ラインといえる
けど現実問題からすれば線形合同法を使用する限りはそれは不可能にょ。
それは上記のようにプチコンで扱える上限が524287だからにょ。
これが32ビット符号無し整数が扱えるならば下記のものが最大数になるにょ。

　X=(1664525*X+1013904223) MOD 2^32

桁あふれの場合のことを考えなくてもいいために実際はMOD以下は不要にょ。
しかし、上限を考えるとプチコンだと下記のものが最大になりそうにょ。

　X=(X*5+3)%65536

この式だと周期は65536になり、上記の十分とされる100万には及ばないけどかなり使える
レベルにはなりそうにょ。
これを用いて上位8ビットのみを使って画面上にドットを表示してみるにょ。
8ビットの最大数は256であるため座標は0から255になるにょ。
プチコンの画面解像度においてX座標は問題ないけどY座標は収まりきらないため下画面も
使って表示をすることにしたにょ。
524287回のループで埋め尽くせる可能性は低いため2重ループにして最大1億回実行できる
ようにしたにょ。（もっともこの疑似乱数は周期が65536だから無意味だけど）

《　リスト4　》
　ACLS
　X=0:Y=0
　　FOR I=0 TO 9999
　　　FOR J=0 TO 9999
　　　　GOSUB @RND:A=X/256
　　　　GOSUB @RND:B=X/256
　　　　GPAGE 0:GPSET A,B,15
　　　　GPAGE 1:GPSET A,B-192,15
　　　NEXT
　　NEXT
　END
@RND
　X=(X*5+3)%65536
　RETURN

これを実行すると画面には5本の直線が引かれるだけにょ。
線形合同法による乱数はパターン化されておりそのパターンが可視化されてしまったと
いうわけにょ。
この乱数列の周期は65536で上位8ビットだけを見た場合に0～255の値が一様に発生している
ため画面の多くの部分にドットが表示されるべきなのだけどそれなのに256x5＝1280個のドット
しか描かれないというのはあまりに少なすぎるにょ。
このような状態になっているのは0～255の数字がほぼ均等に出ているといってもその数字の
組み合わせそのものが少ないと言えそうにょ。（線形合同法のにおける問題の(2)が明確に
現れた形になったということ）
そういった規則性は乱数を一次元的に見ても分からないので2次元もしくは多次元において
検証する必要があるにょ。
その検証は今回やったように2つの数字を組み合わせて使い2次元的な描画を行うことで
視覚的に可能になっているにょ。
このプログラムではドットが表示されているか否かということだけに注目いしているため
ドットの色は白で固定なのだけど色を随時変えることで同じ場所にどの程度の頻度で
ドットを表示しているかというのが視覚的に分かるようになるので下記のように変更して
みるのもいいかもしれないにょ。

　GPAGE 0:GPSET A,B,GSPOIT (A,B)+1
　GPAGE 1:GPSET A,B-192,GSPOIT (A,B-192)+1

このことから疑似乱数を使い下位ビットを捨てて使用するという使用方法をとる場合で
あっても頻度やばらつきは十分にあっても線形合同法の場合は2つの数字を組み合わせて
使用する場合において明確な規則性が出てしまい実用性の面において問題が発生する場合が
あるにょ。
しかし、プチコン標準のRNDでできていることができないというのはさすがに気になるため
これを何とかして改善したいところにょ。

これは線形合同法をやめて非線形合同にする（1次式ではなく2次式にする）とかM系列乱数を
使用するとかいう方法によって解決は可能にょ。
しかし、2次式を使えば扱える上限数値が低いプチコンの場合はOverflowの可能性がさらに
高まるためそれによって公式のMの値（つまり、乱数の周期）を小さくするしかないにょ。
最大に見積もっても周期256回が限界だからね。（X、Y完全に独立した乱数が発生可能で
あっても理論上最大で画面上の128カ所にしかドットを表示することができない）
その点ではM系列乱数ならば非常に優秀にょ。
どんな乱数かを掻い摘んで説明すると予めP個の乱数を配列変数などに確保しておいてその
中で最も古い乱数（新しいものから数えてP番目の乱数）と最新のものから数えてQ番目の
乱数のXORを取りそれを新しい乱数として生成するというものにょ。
リングバッファで確保しておけば処理そのものは単純になるにょ。
この場合はPとQの値で周期は変わるものの最大周期は2＾Pｰ1になるにょ。
さらに線形合同法とは異なり下位ビットも周期の小さい乱数にならないというメリットは
あるもののこれはいくら単純とはいえプログラムに組み込んで使うならばできるだけ短くて
処理の速いものが望ましいからね。
単純で高速な線形合同法で十分な周期やばらつきが得られるならばそれがベターにょ。

では、線形合同法ではこれ以上は無理なのかというとそうではないにょ。
さまざまな方法によって改善は可能だからね。
その1つの方法として私が考えたのはカウンタを併用することにょ。
1回ごとにインクリメントするカウンタを用意すれば座標が1つずつずれていくので理論上
画面をドットで埋め尽くすことが可能になるにょ。（単純なアイデアだけど組み合わせて使用
した場合にランダム性がほぼ無くなってしまう線形合同法においてはかなり効果的だと思う）
カウンタの周期を65536にしてしまえば乱数の周期と重なってしまうため周期は65536になって
しまうけどこれをずらすことで理論上の周期をほぼ2乗にすることが可能になるにょ。（この
カウンタの上限をもっと大きな値にした場合は周期そのものは伸びると考える人もいるかも
しれないけど同じ乱数列をダブって発生させるだけであり実質的な周期は変わらないので
MとM-1が互いに素ならばM-1となる値にするのがベストだと思われる）

《　リスト5　》
　ACLS
　X=0:Y=0
　　FOR I=0 TO 9999
　　　FOR J=0 TO 9999
　　　　GOSUB @RND:A=X/256
　　　　GOSUB @RND:B=X/256
　　　　GPAGE 0:GPSET A,B,15
　　　　GPAGE 1:GPSET A,B-192,15
　　　NEXT
　　NEXT
　END
@RND
　Y=(Y+1)%65535
　X=(X*5+Y)%65536
　RETURN

これを実行すると確かにどんどん埋め尽くされていくもののランダムという感じが全く
せずに明らかにパターン性を感じてしまうにょ。
これは上記の線形合同法の公式におけるAに相当する5が65536に対して小さすぎるためにょ。
これを大きくすれば解決する（たとえば12345とかに変えればかなり自然にドットが描画
されていく）問題だけどそれはプチコンで扱える数字に上限があるため無理にょ。（最大値を
計算すると809095109になってしまい余裕でOverflowになる）
したがって、公式のMの値を小さくするしかないにょ。
そこでMの値は4096にしたにょ。（この4096の理由は後述）
それによって最大100少々までAの値を上げることが可能になり、かなり見た目のばらつきを
持たせることが可能になったにょ。
周期が4096になるAの値はたくさんあるけど実際にドットを表示してみて発生した乱数に
見た目の規則性が最も見られなかった「117」を採用したにょ。

《　リスト6　》
　ACLS
　X=0:Y=0
　　FOR I=0 TO 9999
　　　FOR J=0 TO 9999
　　　　GOSUB @RND:A=X/16
　　　　GOSUB @RND:B=X/16
　　　　GPAGE 0:GPSET A,B,15
　　　　GPAGE 1:GPSET A,B-192,15
　　　NEXT
　　NEXT
　END
@RND
　Y=(Y+1)%4095
　X=(X*117+Y)%4096
　RETURN？

規則性はほぼ無くなり上位8ビットを使う限りは全く問題はないように見えるにょ。
ちなみに771582回のループ（1543164回の乱数発生）で256x256の範囲をドットで埋め尽くす
ことができたにょ。（期待値は764646回なのでほぼ期待値通り）
これならばRNDの置き換えは十分にできそうにょ。
線形合同法による乱数が苦手としている二次元分布にも十分使えると思うのでプチコンで
普通に使う分には特に困ることはないと思うにょ。（再現性が要求されるシミュレーション
とかにも活用できると思う）

これで終わってもいいけどやはり0～4095の乱数というのは気に入らないため4096で割った
数字を発生させることにしたにょ。
つまり、0以上1未満の乱数ということにょ。

《　リスト7　》
　ACLS
　X=0:Y=0
　　FOR I=0 TO 9999
　　　FOR J=0 TO 9999
　　　　GOSUB @RND:A=X*256
　　　　GOSUB @RND:B=X*256
　　　　GPAGE 0:GPSET A,B,15
　　　　GPAGE 1:GPSET A,B-192,15
　　　NEXT
　　NEXT
　END

@RND
　Y=(Y+1)%4095
　X=(X*479232+Y)%4096/4096
　RETURN

これによって下記の2つのメリットが生まれてくるにょ。

　(1)使いやすい
　(2)上位ビットのみが使用される

(1)は言うまでもないことだけど0～1の乱数というのは一般的なBASIC言語では非常に
ポピュラーであるため使用に慣れている人も多いということにょ。
これが0～4095の乱数だとダイスで1～6の目を表示したいというだけでも4096で一端割ってから
6をかけて1を足して整数化するという手間が加わるにょ。
つまり、4096の約数の乱数を発生させない限りは必ず4096で割るという手間が入ってしまう
ためにそれならば事前に4096で割ったものを発生した方がいいということにょ。

これが4096ではなく8192だったらそういうわけにはいかないにょ。
それはプチコンの小数部は12ビットであるため1/4096単位の数に丸められてしまうためにょ。
http://ww5.tiki.ne.jp/~ochame/petitcom/p005.htm#column
丸められたら4096の場合よりも周期が短くなってしまうにょ。
つまり、0～1の乱数を発生させるという時点で0～4095の乱数である必要が出てくるわけにょ。

(2)出てくる数字が小数であれば0～3の乱数が欲しい場合はAに乱数の返り値が入っている
場合はFLOOR (A*4)もしく0 OR A*4 いう方法をとるのが普通だけどこれが0～4095であれば
A%4という方法をとるのが普通だろうからね。
つまり、下位ビットを優先して使う人が出てくるということにょ。
カウンタを併用することで乱数のパターン性はかなり減ったとはいえ下位ビットにおいては
上記のような小さい間隔の周期性が残されているために下位ビットのみを使うというのは
推奨できないにょ。
わざわざ推奨できない方法をコメントで書くのも何なので普通にかつお手軽に使えるものに
するならばこれがベストだと思うにょ。
C言語の場合<stdlib.h>に含まれるrand関数は線形合同法によるものだけど下位ビットを
捨てて0～32767の値で返しているので上記のような偶数と奇数が交互出るという問題は発生
しないもののそれでもRAND_MAXの値で割って「事実上0～1の乱数の状態」にして使用する
ことが推奨されているくらいだからね。

あとプチコンで扱える数字が1/4096単位ならば上限2048や1024の乱数でもいいという考えの
人もいるかもしれないにょ。
確かにそれはそうだけど少しでも周期が長い方がいい（1024にしてしまうと周期は1/16に
なってしまう）わけだし、それだけではなくプチコンで扱える数字が1/4096単位という
ことはそれ以下の数字は無視されるということにょ。
その結果が上記リンク先のコラムに書いている演算誤差に繋がっているわけだけど4096で
事前に割っておくことで普通に使うだけで下位ビットの周期性の短さは演算誤差で打ち
消されてしまう可能性が高くなるということにょ。
つまり、線形合同法の問題点が自然消滅するというわけにょ。
もちろん、わざわざ4096倍して使えば下位ビットの周期性の短さは健在になるけどこれが
100倍とか1000倍とかの「4096の倍数でない」ならば演算誤差で消えて無くなる可能性が
高いくなるにょ。
ちなみにこの乱数の周期は16773120であるため普通にゲームとかで使用していても問題
ないと思われるにょ。（1フレームで8回乱数を発生させても10時間の周期になる）

では、数値的な偏りがどの程度あるのかをモンテカルロ法によって検証してみたにょ。

《　リスト8　》
　X=0:Y=0
　FOR I=1 TO 100000
　　GOSUB @RND
　　A=X
　　GOSUB @RND
　　B=X
　　IF A*A+B*B<1 THEN C=C+1
　NEXT
　? C*4/100000
　END
@RND
　Y=(Y+1)%4095
　X=(X*479232+Y)%4096/4096
　RETURN

これを実行すると「3.14」となりこれは円周率にかなり近い値になっていることから
普通に使用して問題ないレベルだと推測されるにょ。

あと気になるのは速度かもしれないにょ。
10万回ループによって実行してみたところRNDは198フレーム、この疑似乱数は590フレーム
だったにょ。
遅いようなイメージがあるけど1フレームあたりだと169回発生できているためメインルーチン
1回につき8個の乱数を発生させても1フレームの1/20以下の時間で済むにょ。
またこの乱数で使用しているカウンタだけどメインルーチンで他にカウンタを使用している
ならばそのカウンタを使えばY=(Y+1)%4095が不要になるにょ。
その場合はそのカウンタの変数がPの場合は X=(X*479232+P%4095)%4096/4096 とする必要が
あるにょ。
リストが短くなっただけではなく速度も10万回で463フレーム（1フレームあたり216回）と
いうことで1.5倍も高速化可能になるにょ。
プチコンでできるだけ短くで実用になる疑似乱数を求めているならばこれはかなりおすすめ
だと思うにょ。
初期値を与えるのが面倒ならばY=MAINCNTL%4096を入れておけばいいにょ。
X=MAINCNTL%4096/4096を入れてもいいけど「/4096」の分だけリストが長くなるにょ。
まぁ再現性があるのがRNDよりもこの疑似乱数を使う意味があるわけだから初期値を
お任せにしたら全く意味が無くなってしまうにょ。

あとそこまで長期の周期性が求められてない場合でなおかつプレイヤーにとって気分の良い
ばらつきの疑似乱数を発生させる場合には乱数テーブルを用意してそれを読み取るという方法もあるにょ。
1/2の成功確率のイベントに10回連続で失敗とかいうのは線形合同法等の演算による疑似乱数
なら起きてしまう可能性は十分にあるけど予めテーブル化しておけばそれは防げるからね。（もっとも演算による疑似乱数でも完全に一定確率ではなく失敗が続けば確率が上がるように
して期待値が2回になるような設定にすれば良いだけのことだけど）
しかし、そのテーブル化の手間を考えると速度面さえ気にならなければ今回書いた疑似乱数の
方が短いリストで済むと思われるにょ。

今回いろいろと疑似乱数を作ってみて思ったのだけど画面上にドットを表示したときの
様子を見るとプチコンのRNDによる動作によく似ているにょ。
もしかするとプチコンのRNDは線形合同法＋乱数を呼び出すたびにインクリメントという
私が考えた方法と同じような方法によって乱数を発生させているのかもしれないにょ。
上記のように乱数のシードにMAINCNTLのカウンタ用いているために同じ乱数列の乱数を
発生できないと考えているにょ。

再現性のある疑似乱数ということで初期値さえ同じならば環境が違っても同じものを再現が
可能になるにょ。（例えば先日作った横スクロールゲーム「CAVE」でもプチコン版とPC版で
全く同じステージでのプレイが可能になる）
ということで、PCを使って自由な縦横サイズにドットを表示するプログラム（要は上記
リスト7に解像度可変とドットが埋められたかどうかの判定を付けたもの）を作ったにょ。
http://ww5.tiki.ne.jp/~ochame/rand.zip
これを見ればどの程度のバラツキがあるのかはプチコンを持って無くても分かると思うにょ。
また、プチコンを持っている場合は事前にPCでシミュレーションをしておくことで乱数の
初期値をいくらに設定するのが良いのかというのも分かるにょ。（上記プログラムでは
初期値変更の機能まではサポートけど）
ちなみに1024x1024を埋めるのに7099847回のループ（14199694回の乱数発生）が必要で
2048x2048ドットだとこの疑似乱数の最大周期の最後でちょうど埋められるにょ。
さすがに4096x4096だと線形合同法特有のラインが発生してしまう（2048本の線が描かれる）
もののプチコンで使うならば十分だと思うにょ。

《追記》
上記のPC版の疑似乱数プログラムをバージョンアップしたにょ。
オプション機能を付けてX、Yの初期値が設定できるようになったのに加えて疑似乱数の
計算式も3種類が選べるようになったにょ。
今回私が考案した線形合同法におけるカウンタの効果の確認やMAX4096ではなくMAX8192に
変更した場合の違いの確認もできるようになったにょ。
8192と4096の違いは表示ウェイトを10くらいに設定したら4096の方がいかにランダムに
描画されているかということが分かると思うにょ。