今回はちょっと趣向を変え、疑似乱数の作り方について話します。乱数については第１部第７５章でやりました。この時、「疑似乱数はある初期値を元に計算を行って生成されます」といいました。今回はその方法の１つ「線形合同法」について話します。
疑似乱数生成ルーチンを自作すると、いろいろ面白いこともできるでしょう。
では、今回の要点です。
では、いってみましょう。
一見でたらめな数を計算で出すにはどうすればいいでしょうか？
普通の四則演算などを行っていると、何か規則性の見える結果がでてきます。かといって、ややこしい計算をしてでたらめな数がでてきても、速度が欲しいときには実用的とは言えません。
しかし、簡単で、しかも規則性の見えにくい計算というのがあります。それは「割り算」です。といっても、求めるのは商ではなく「余り」の方です。線形合同法は、この余りを含めた数列を使うことによって、一見乱雑な数列を作ろうというものなのです。
しかし、余りを求めるという計算は必ず元の数以下になってしまう計算です。そこで、一次式を中に含めることにします。
x = (a * x + c) % M;
あとはｘに初期値を入れて上の計算を行うと、一見乱雑な値が得られます。これが線形合同法の基本です。
そして、同じ初期値からは同じ値がでてくるので、最大でもＭ回で元に戻ってしまいます。つまり、合同法乱数は周期性を持っているわけです。
では、ａやｃやＭには何を選んでもいいのでしょうか？ 答はＮＯです。
例えば、Ｍが１だと必ず０にしかなりませんし、Ｍが大きくてもａ＝１，ｃ＝０だと同じ値しかでてきません。
それは極端な例としても、ある程度の選び方というものがあります。今度はそのことについて話しましょう。
先ずは除数Ｍについてです。
なるべく面倒な計算はしたくありません。そして、掛け算でオーバーフロー（計算結果が大きすぎて表現できなくなること）がおこる可能性を考えると、一見線形合同法は実現しにくいような気がします。
しかしこのオーバーフローというコンピュータの特性を逆に利用してしまうことができます。
例えば３２ビット機の int 型で計算すると、計算結果が３２ビットで表現できないときは下位３２ビットのみを残して上位ビットを切り捨ててしまいます。
これは２３２の余りを求めることと同じになります。つまり、Ｍ＝２３２にすれば、いちいち割り算をしなくても勝手に余りを求めてくれるというわけです。
また、Ｍが大きい方が周期を長くすることができるので、Ｍは２３２にするといいでしょう。
次は乗数ａと定数ｃについてです。
ａ＝１，ｃ＝０のときは必ず毎回同じ値がでてくるというように、ａやｃの数が悪いと極端に周期が短くなります。乱数の周期性が長い方が隠れた規則性が目立たないので、周期が最大の２３２になるように選びたいものです。
先ず、Ｍは偶数なのでｃも偶数だと困ります。ａ×ｘが奇数の時はａ×ｘ＋ｃも奇数になり、Ｍで割った余りも奇数になります。そして、ａ×ｘが偶数の時はａ×ｘ＋ｃも偶数になり、Ｍで割った余りも偶数になります。
ということで、ｃが偶数だと周期は必ずＭになりません。ということで、ｃは奇数にするといいでしょう。
ただ、速度を求めるときはｃを０にしたいと思います。０も偶数なので、ｃを０にすると周期は短くなります。そのときは、周期が短くなるのだということを覚悟して使って下さい。
そしてａですが、詳しいことは抜きにしますが、これはａ ＭＯＤ ８が５になるようにすると周期が長くなります（ＭＯＤ は余りを求める一般的な演算子です）。ｃが奇数の時はその周期はＭになります。
ということで、ａを８ｎ＋５、ｃを奇数（ｎは整数）にするといいでしょう。
ちなみに、ｃを０にしたいときは、ａを上のようにしてｘの初期値を奇数にすれば周期が２３０になります（ただし、常に奇数になります）。速度が要求されるときはこういう設定にしましょう。
そして、周期２３２の０〜２３２−１の数を２つ組み合わせていくと、当然周期は２３１になります。０〜２３２−１の数を２つ組み合わせるには２６４の組み合わせがあるので、ほとんどの組み合わせはでてこないことになります。/p>
ａの取り方によってはさらにここでの質が悪くなるときがあります。つまり、２６４通りの組み合わせのうち、偏った分布になるということがあるわけです。
これは３つ、４つと組み合わせていったときにも起こる問題で、合同法乱数は１個ずつ使えばランダムだが、組にして使うとランダム性が低いということになります。
こういった質をあまり気にならないようにするには、ａを特によい数にする必要があります。その数には
69069 1664525 39894229 48828125 1566083941 1812433253 2100005341
などが知られています。
ａにはこれらのうちの１つを使うようにしましょう。
ちなみに、よく使われるが悪い数に６５５３９があります。これを奨めているところがあっても、それは信用しないようにしましょう。
以上のことをまとめると、
ということになります。
これで一応は周期２３２の０〜２３２−１の疑似乱数ができるわけですが、使うときにちょっと気を付けることがあります。
実は上のａを８ｎ＋５にすると周期がＭになるというのは、Ｍが２の累乗であるとき全てに共通することなのです。つまり、こうやって得た乱数を使って０〜２ｋ−１の乱数が欲しいときに「ｘ ＭＯＤ ２ｋ」としてしまうと、周期が２ｋになってしまいます。
これは下位ビットをとるとこうなるのであって、上位ビットを使えばこんなことにはなりません。０〜２ｋ−１の乱数が欲しいときは、必ず「ｘ ＞＞ （３２−ｋ）」とビットシフトを使うことになります。
つまり、合同法乱数は上位の桁はランダム性が高いが、下位の桁はランダム性が低いということです。
合同法乱数を生のまま使うとこういう問題が起きるので、合同法乱数を使うときは上位ビットを取り出してから使うとよいでしょう。
以上のことを踏まえて乱数生成関数を作ってみましょう。
| プログラム１ |
|---|
// Rand2.h
#include <stddef.h> // NULL 用
#include <time.h> // time 用
// 使用する上位ビットの数
const int LCRAND_BITS = 15;
// 乱数の上限
const unsigned int LCRAND_MAX = (1 << LCRAND_BITS) - 1;
// 初期化
void LCInitSeed(unsigned int nSeed);
// 疑似乱数生成
int LCRand();
// 現在の時間で初期化
inline void LCInitWithTime()
{
LCInitSeed((unsigned int)time(NULL));
}
// ある数未満の乱数を生成
inline int LCRand(int nLimit)
{
// % nLimit よりも上質の分布になります
// ただし、１６ビット機では long でやらないといけません
return LCRand() * nLimit / (LCRAND_MAX + 1);
} |
| プログラム２ |
// Rand2.cpp
#include "Rand2.h"
static unsigned int gs_nSeed; // 乱数種
// 初期化
void LCInitSeed(unsigned int nSeed)
{
gs_nSeed = nSeed;
}
// 疑似乱数生成
int LCRand()
{
// ａには 48828125 を、ｃには 1 を使用しました
gs_nSeed = gs_nSeed * 48828125 + 1;
// 上位ビットを返します
return gs_nSeed >> (32 - LCRAND_BITS);
} |
| プログラム３ |
// Rand3.cpp
#include <iostream.h>
#include <stdio.h>
#include <memory.h>
#include "Rand.h"
#define numof(array) (sizeof (array) / sizeof *(array))
// サイコロを何度も振ったときの分布を調べます
void CheckDist(int nTimes)
{
const int DICE_LIMIT = 6; // サイコロの目の上限
int nDist[DICE_LIMIT * 2 - 1]; // 結果の分布
int i;
memset(nDist, 0, sizeof nDist);
for(i = 0; i < nTimes; i++)
nDist[LCRand(DICE_LIMIT) + LCRand(DICE_LIMIT)]++;
for(i = 0; i < numof(nDist); i++)
printf("%2d : %6.2lf%% ",
i + 2, (double)nDist[i] * 100 / nTimes);
putchar('\n');
}
int main()
{
int nTimes; // 試行回数
LCInitWithTime(); // 乱数の初期化
while(true)
{
cout << "何回サイコロを振りますか > " << flush;
cin >> nTimes;
if(nTimes <= 0)
break;
CheckDist(nTimes);
}
return 0;
} |
| 実行結果例 |
何回サイコロを振りますか > 20 2 : 0.00% 3 : 10.00% 4 : 5.00% 5 : 5.00% 6 : 20.00% 7 : 15.00% 8 : 15.00% 9 : 15.00% 10 : 5.00% 11 : 0.00% 12 : 10.00% 何回サイコロを振りますか > 100 2 : 1.00% 3 : 4.00% 4 : 7.00% 5 : 9.00% 6 : 22.00% 7 : 22.00% 8 : 12.00% 9 : 9.00% 10 : 7.00% 11 : 6.00% 12 : 1.00% 何回サイコロを振りますか > 10000 2 : 2.72% 3 : 5.75% 4 : 8.62% 5 : 10.87% 6 : 13.56% 7 : 16.51% 8 : 14.01% 9 : 11.19% 10 : 8.47% 11 : 5.54% 12 : 2.76% 何回サイコロを振りますか > 10000000 2 : 2.78% 3 : 5.56% 4 : 8.34% 5 : 11.11% 6 : 13.89% 7 : 16.68% 8 : 13.88% 9 : 11.11% 10 : 8.33% 11 : 5.55% 12 : 2.78% 何回サイコロを振りますか > 0 |
サイコロを２回振ったとき、各合計のでる確率は
2 : 2.78% 3 : 5.56% 4 : 8.33% 5 : 11.11% 6 : 13.89% 7 : 16.67% 8 : 13.89% 9 : 11.11% 10 : 8.33% 11 : 5.56% 12 : 2.78%
です。試行回数を増やす毎に極限値に近づいていくということを見事に再現できていますね。
乱数の作成に使用する変数のことを乱数種と呼びます。実際には乱数種をそのまま返すわけではないのですが、これを返す関数を作っておくと便利なことができます。
それは、ゲームなどでリセットしても同じ乱数を生成することができるということです。某シミュレーションＲＰＧや某競走馬育成ゲームなどにみられるアレです。こういったゲームを作るときには、乱数種を保存するようにしてみてはどうでしょうか？
では、今回の要点です。長いです。
では、次回をお楽しみに。
Last update was done on 2000.11.3
この講座の著作権はロベールが保有しています