rand() の問題点
これまで広く使用されてきた rand() による乱数生成には以下の問題点がある。
- 生成される範囲が [0, 32767] と狭い
- rand() % N は一様では無い
- 周期があまり長くない
- 乱数生成アルゴリズムが固定(通常は線形合同法)
- 正規分布など、一様でない乱数生成が面倒
rand() で生成される値の範囲は [0, 32767] と狭い。15ビットしかない。
(rand() << 15) + rand() とすれば、30ビット乱数にすることは可能だが、スマートではないし、
生成アルゴリズムの関係で、乱数に偏りが出る場合がある。
通常、乱数として欲しいのは [0, N) 範囲の値であることが多い。この場合は rand() % N とすることが一般的だ。
しかし、このようにして生成した乱数は一様でなく偏りがある。特に N が大きい時にその現象が顕著になる。
例えば、N が1万だとすると、rand() % N は [0, 2768) の範囲の値を返す確率は [2768, 10000) の範囲の値を返す確率より有意に大きくなる。
具体的には 4:3 の比率になる。100万回乱数を発生させるとすると、前者の期待値は 120回, 後者の期待値は 90回となる。
当然、この問題を回避するコードは記述可能なのだが、ライブラリがそれを用意するべきだ。
rand() の周期は 2^31 である。充分な長さだと思う読者もいるかもしれないが、適用分野によっては短い場合もある。
通常、rand() は線形合同法により乱数を生成する。このアルゴリズムは原始的であり、高速だが、あまりよい性質の乱数を生成してくれない。
場合によっては乱数系列を推測されてしまう危険性もある。
オンラインポーカーのシャフルにこの乱数を使い、もしユーザに次のカードを予測されてしまっては、大問題だ。
rand() は一様な乱数しか生成してくれない。正規分布の乱数を生成させたい場合は、そのような分布になるよう、プログラムを記述しなくていけない。
上記のような問題点を解消するために、C++11 で乱数ライブラリが追加された。
Visual Studio に於ける本格的なC++11対応は VS2013 からであるが、auto など一部の仕様は VS2010 から利用可能になっている。
嬉しいことに、この random ライブラリも利用可能なので、VS2010 以上を使用しているのであれば random ライブラリ使わなきゃ損だぜ。
演習問題:解答例
- rand() % 10000 により [0, 10000) の乱数を100万回発生させ、0, 1, 2, 9997, 9998, 9999 が出る回数をカウントし、表示しなさい。
- 同様に、rand() % 1000 の場合はどうなるか検証してみなさい。
C++ 乱数ライブラリ std::random の使い方
準備
乱数ライブラリを使用するには、以下のように <random> をインクルードする。
#include <random>
非決定的な乱数生成
非決定的な乱数生成器として「std::random_device」が用意されている。
「非決定的」とは、次の値が一意に決まらないという意味だ。内部的にはハードウェアの状態や時刻を用いて実現している。
使い方は「std::random_device rnd;」で乱数生成器(rnd は乱数生成器の名前、お好みの名前を付けてよい)を生成し、 「rnd()」をコールすると、乱数を生成してくれる。生成される乱数は 32bit で、[0, 0xffffffff] の範囲だ。
※ rnd は std::random_device クラスのオブジェクトなのに、それをあたかも関数のように「()」を付けてコール出来るの? と不思議に思った読者もいるかもしれない。 実はクラス定義で () 演算子(oprator()())をオーバーロードすることで、オブジェクト名() でその演算子をコールすることが出来るのだ。 これを「ファンクタ」と呼ぶ。これは初級者にはわかりづらいかもしれないので、よくわからない場合は、あまりこだわらず、そういうものだと思い、 先に進んで欲しい。
#include <random> int main() { std::random_device rnd; // 非決定的な乱数生成器 for (int i = 0; i < 10; ++i) { std::cout << std::hex << rnd() << "\n"; } getchar(); return 0; }
上記は10個の乱数を生成し、表示するサンプルだ。
演習問題:
- 上記ソースコードをビルド・実行し、結果を確かめなさい。
決定的な乱数(擬似乱数)生成
「std::random_device」は若干速度が遅いので、通常は決定的な乱数生成器を使用する。
「決定的」とは、現在の値により次の値が一意に決まるという意味。擬似乱数とも呼ばれ、数学的に完全な乱数ではない。
が、rand() の srand(int) 同様、乱数のシードを指定できるので、決定的でも実用上問題は無い
決定的な乱数生成アルゴリズムとしては、線形合同法、
メルセンヌ・ツイスタなどがあり、選択することができる。
線形合同法は大昔から使われていたアルゴリズムで、単純・高速ではあるが、乱数の質はあまり良くない。
メルセンヌ・ツイスタは1990年代後半に日本人により開発されたもので、周期が 2^219937-1 もあるなど、高品質の乱数を高速に生成するように設計されている。
かなりお薦めだ。
本稿では、この メルセンヌ・ツイスタ を用いた乱数生成について説明する。
まずは、シード生成のために「std::random_device rnd;」のようにして、非決定的乱数生成器を生成する。 これまで通り、time(0) を使ってもよい。が、お薦めはしない。
これを引数として指定し、「std::mt19937 mt(rnd());」の様にメルセンヌ・ツイスタ生成器を生成する(mt は生成器の名前。自由に命名してよい)。
後は「mt()」により乱数を生成するだけだ。
生成される乱数の範囲は 32ビット、[0, 0xffffffff] の範囲だ。
64ビット範囲の乱数を生成したい場合は std::mt19937_64 を使用する。
std::random_device rnd; // 非決定的な乱数生成器 std::mt19937 mt(rnd()); // メルセンヌ・ツイスタの32ビット版、引数は初期シード //std::mt19937 mt((int)time(0)); // メルセンヌ・ツイスタの32ビット版、引数は初期シード for (int i = 0; i < 10; ++i) { std::cout << std::hex << mt() << "\n"; }
シードの指定はコンストラクタだけでなく、seed(val) メソッドにより指定することも可能。
下記はその例。先にデフォルトコンストラクタでオブジェクトを構築し、後でシードを指定している。
自分で定義するクラスのメンバ変数として std::mt19937 を保持する場合は、この方法を使いコンストラクタなどでシードを指定するとよい。
std::mt19937 mt; // メルセンヌ・ツイスタの32ビット版 std::random_device rnd; // 非決定的な乱数生成器 mt.seed(rnd()); // シード指定 for (int i = 0; i < 10; ++i) { std::cout << std::hex << mt() << "\n"; }
演習問題:
- 上記ソースコードをビルド・実行し、結果を確かめなさい。
- 上記ソースコードで、初期シードを指定しないようにして、ビルド・実行し、結果を確かめなさい。
- mt19937 を mt19937_64 に変更・ビルド・実行し、結果を確かめなさい。
範囲指定乱数生成
これまで、std::random_device、std::mt19937 を生で使ってきたが、これらは範囲が固定だった。
rand() の時と同様に「% N」演算を行うことで [0, N) の範囲の(ある程度)一様乱数を得ることは可能だが、
もっと正確な一様乱数を得る仕組みが用意された。それが「uniform_int_distribution」だ。
「std::uniform_int_distribution<> rand100(0, 99);」のように宣言すると、[0, 99] の範囲の一様分布器を生成できる。
これを先の乱数生成器と組み合わせ、「rand100(mt)」のように記述することで一様分布乱数を生成できる。
※ 引数が「mt()」ではなく「mt」であることに注意されたい。生成された乱数値を渡しているではなく、乱数生成器自体を渡している。
乱数生成器が持つ情報を参照可能だから、正確に一様分布する乱数を発生可能なのである。
std::random_device rnd; // 非決定的な乱数生成器を生成 std::mt19937 mt(rnd()); // メルセンヌ・ツイスタの32ビット版、引数は初期シード値 std::uniform_int_distribution<> rand100(0, 99); // [0, 99] 範囲の一様乱数 for (int i = 0; i < 20; ++i) { std::cout << rand100(mt) << "\n"; }
上記は、初期シードに非決定的乱数生成器で生成した値を使ったメルセンヌ・ツイスタ乱数生成器と、 [0, 100) の範囲の一様分布器を組み合わせて、20個の乱数を生成する例。
整数ではなく、double の一様分布器は「uniform_real_distribution」で生成する。
一様乱数ではなく、「normal_distribution」を使って、指定された平均値、分散値の正規分布乱数を生成することもできる。
コードは以下の様になる。
std::random_device rnd; // 非決定的な乱数生成器でシード生成機を生成 std::mt19937 mt(rnd()); // メルセンヌツイスターの32ビット版、引数は初期シード //std::uniform_int_distribution<> rand100(0, 99); // [0, 99] 範囲の一様乱数 std::normal_distribution<> norm(50.0, 10.0); // 平均50, 分散値10の正規分布 for (int i = 0; i < 20; ++i) { //std::cout << rand100(mt) << "\n"; std::cout << norm(mt) << "\n"; }ma
演習問題:解答例
- 上記ソースコードをビルド・実行し、結果を確かめなさい。
- 上記ソースコードを、[1.0, 2.0] の範囲の一様なdouble型乱数を20個表示するよう修正・ビルド・実行してみなさい。
- 上記ソースコードを、平均 0.0、分散 0.2 の正規分布乱数を20個表示するように修正してみなさい。
- [0, 9999] の範囲の一様な整数乱数を100万回生成し、0, 1,2 と 9997, 9998, 9999 の値が出る回数がほぼ等しいことを確認しなさい。
- rand() % 100 と、本稿の両方の方法で [0, 99] の範囲の一様な乱数を1000回生成し、平均・標準偏差を計算し、それらを比較してみなさい。
まとめ
ここに紹介した以外にも、random ライブラリには便利なクラス・機能がたくさん用意されている。 詳細は cpprefjp - random(C++11) や、 <random> などを読んで自分で勉強してね。
これまでの、srand(), rand() に慣れている人には、ちょっと面倒に感じるかもしれないが、
慣れればたいしたことないし、パズルの問題生成・カードシャフル・敵機の動きなど、ちゃんとした乱数にした方がよいことは確かなので、
random ライブラリ使いこなせるようになっておくれやす。
演習問題解答例
- rand() % 10000 により [0, 10000) の乱数を100万回発生させ、0, 1, 2, 9997, 9998, 9999 が出る回数をカウントし、表示しなさい。
- 同様に、rand() % 1000 の場合はどうなるか検証してみなさい。
#include <iostream> #include <time.h> #include <vector> int main() { srand((int)time(0)); const int N = 10000; // [0, N) の乱数を発生 std::vector<int> vc(N, 0); // 回数を保存する vector for (int i = 0; i < 1000*1000; ++i) { // 百万回ループ int r = rand() % N; // 乱数生成 vc[r] += 1; // 生成された乱数の個数をインクリメント } std::cout << "vc[0] = " << vc[0] << "\n"; std::cout << "vc[1] = " << vc[1] << "\n"; std::cout << "vc[2] = " << vc[2] << "\n"; std::cout << "vc[" << N-3 << "] = " << vc[N-3] << "\n"; std::cout << "vc[" << N-2 << "] = " << vc[N-2] << "\n"; std::cout << "vc[" << N-1 << "] = " << vc[N-1] << "\n"; getchar(); return 0; }
#include <iostream> #include <time.h> #include <vector> int main() { srand((int)time(0)); const int N = 1000; // [0, N) の乱数を発生 std::vector<int> vc(N, 0); // 回数を保存する vector for (int i = 0; i < 1000*1000; ++i) { // 百万回ループ int r = rand() % N; // 乱数生成 vc[r] += 1; // 生成された乱数の個数をインクリメント } std::cout << "vc[0] = " << vc[0] << "\n"; std::cout << "vc[1] = " << vc[1] << "\n"; std::cout << "vc[2] = " << vc[2] << "\n"; std::cout << "vc[" << N-3 << "] = " << vc[N-3] << "\n"; std::cout << "vc[" << N-2 << "] = " << vc[N-2] << "\n"; std::cout << "vc[" << N-1 << "] = " << vc[N-1] << "\n"; getchar(); return 0; }
- 上記ソースコードをビルド・実行し、結果を確かめなさい。
- 上記ソースコードを、[1.0, 2.0] の範囲の一様なdouble型乱数を20個表示するよう修正・ビルド・実行してみなさい。
- 上記ソースコードを、平均 0.0、分散 0.2 の正規分布乱数を20個表示するように修正してみなさい。
- [0, 9999] の範囲の一様な整数乱数を100万回生成し、0, 1,2 と 9997, 9998, 9999 の値が出る回数がほぼ等しいことを確認しなさい。
- rand() % 100 と、本稿の両方の方法で [0, 99] の範囲の一様な乱数を1000回生成し、平均・標準偏差を計算し、それらを比較してみなさい。
解答例省略
std::random_device rnd; // 非決定的な乱数生成器を生成 std::mt19937 mt(rnd()); // メルセンヌ・ツイスタの32ビット版、引数は初期シード値 std::uniform_real_distribution<> rand12(1.0, 2.0); // [1.0, 2.0] 範囲の一様乱数 for (int i = 0; i < 20; ++i) { std::cout << rand12(mt) << "\n"; }
std::random_device rnd; // 非決定的な乱数生成器を生成 std::mt19937 mt(rnd()); // メルセンヌ・ツイスタの32ビット版、引数は初期シード値 std::normal_distribution<> norm(0.0, 0.2); // 平均0.0, 分散値0.2の正規分布 for (int i = 0; i < 20; ++i) { std::cout << norm(mt) << "\n"; }
#include <iostream> #include <vector> int main() { const int N = 10000; // [0, N) の乱数を発生 std::random_device rnd; // 非決定的な乱数生成器を生成 std::mt19937 mt(rnd()); // メルセンヌ・ツイスタの32ビット版、引数は初期シード値 std::uniform_int_distribution<> randN(0, N-1); // [0, N) 範囲の一様乱数 std::vectorvc(N, 0); // 回数を保存する vector for (int i = 0; i < 1000*1000; ++i) { // 百万回ループ int r = randN(mt); // 乱数生成 vc[r] += 1; // 生成された乱数の個数をインクリメント } std::cout << "vc[0] = " << vc[0] << "\n"; std::cout << "vc[1] = " << vc[1] << "\n"; std::cout << "vc[2] = " << vc[2] << "\n"; std::cout << "vc[" << N-3 << "] = " << vc[N-3] << "\n"; std::cout << "vc[" << N-2 << "] = " << vc[N-2] << "\n"; std::cout << "vc[" << N-1 << "] = " << vc[N-1] << "\n"; getchar(); return 0; }
#include <iostream> #include <time.h> #include <vector> #include <random> int main() { const int N = 1000; { srand((int)time(0)); int sum = 0; // Σx int sum2 = 0; // Σx^2 for (int i = 0; i < N; ++i) { int r = rand() % 100; sum += r; sum2 += r * r; } double ave = (double)sum / N; std::cout << "ave = " << ave << "\n"; std::cout << "σ = " << sqrt((double)sum2/N - ave * ave) << "\n"; } { int sum = 0; // Σx int sum2 = 0; // Σx^2 std::random_device rnd; // 非決定的な乱数生成器を生成 std::mt19937 mt(rnd()); // メルセンヌ・ツイスタの32ビット版、引数は初期シード値 std::uniform_int_distribution<> rand100(0, 99); // [0, 99] 範囲の一様乱数 for (int i = 0; i < N; ++i) { int r = rand100(mt); sum += r; sum2 += r * r; } double ave = (double)sum / N; std::cout << "ave = " << ave << "\n"; std::cout << "σ = " << sqrt((double)sum2/N - ave * ave) << "\n"; } getchar(); return; }