Skip to content

Commit

Permalink
乱数書き始め
Browse files Browse the repository at this point in the history
  • Loading branch information
EzoeRyou committed Mar 6, 2019
1 parent 340545b commit d9aad45
Show file tree
Hide file tree
Showing 3 changed files with 1,315 additions and 693 deletions.
211 changes: 211 additions & 0 deletions 043-random.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,211 @@
# 乱数

乱数はプログラミングにおいてよく使う。例えば6面ダイスをプログラムで実装するには、1,2,3,4,5,6までのいずれかの目を出す。

~~~
$ ./dice
1
5
$ ./dice
3
5 1 6
$ ./dice
10
5 1 6 6 1 6 6 2 4 2
~~~

このプログラム`dice`は標準入力から整数型の値`n`を取り、1,2,3,4,5,6のいずれかをそれぞれ$\frac{1}{6}$の確率で`n`個出力する。

まずこの`dice`プログラムを作ることを目標にC++の乱数アルゴリズムである`<random>`の使い方を学んでいく。


## 疑似乱数

コンピューターで使われる乱数のほとんどは疑似乱数と呼ばれる方法で生成されている。様々なアルゴリズムがあるが、とても簡単に理解できる疑似乱数のアルゴリズムに、線形合同法(Linear congruential generator)がある。

線形合同法では今の乱数を$X_n$、次の乱数$X_{n+1}$とすると、$X_{n+1}$は以下のように求められる。

$$
X_{n+1} = (a \times X_{n} + c) \bmod m
$$

たとえば$a = 3, c = 5, m = 2^{sizeof(std::uint32_t) \times 8}$の場合で、$X_0 = 0$のとき、

$$X_0 = 0$$

$$X_1 = 3 \times 0 + 5 \bmod 2^{32}-1 = 5$$

$$X_2 = 3 \times X_1 + 5 \bmod 2^{32}-1 = 20$$

$$X_3 = 3 \times X_2 + 5 \bmod 2^{32}-1 = 65$$

「これは全然乱数ではない。予測可能じゃないか」と考えるかも知れない。しかし中でどのように乱数が生成されているかわからなければ、外部からは乱数のように見える。これが擬似乱数の考え方だ。

## 乱数エンジン

`乱数エンジン`は生の乱数を生成するライブラリだ。クラスで実装されている。

乱数エンジンはメンバー関数`min()`で最小値を、メンバー関数`max()`で最大値を、`operator()`で最小値から最大値の間の乱数を返す

~~~c++
template < typename Engine >
void f( Engine & e )

// 最小値
auto a = e.min() ;
// 最大値
auto b = e.max() ;
// 乱数
auto r1 = e() ;
// 次の乱数
auto r2 = e() ;
}
~~~
乱数エンジンのオブジェクト`e`は`operator ()`を呼び出すたび、つまり`e()`をするたびに変更される。これは疑似乱数のための内部状態を更新するためだ。そのため、乱数エンジンはconstでは新しい乱数を作るのに使えない。
標準ライブラリはデフォルトの乱数エンジンとして`std::default_random_engine`を提供している。
以下のプログラムはデフォルトの乱数エンジンから乱数を10個出力する。
~~~cpp
int main()
{
// 乱数エンジン
std::default_random_engine e ;
for ( int i = 0 ; i != 10 ; ++i )
{
// 乱数を出力
std::cout << e() << "\n"sv ;
}
}
~~~

標準ライブラリの提供する乱数エンジンには様々なものがあるが、本書ではもう一つ、メルセンヌツイスターというアルゴリズムを実装した乱数エンジンを紹介する。`std::mt19937`だ。

`std::mt19937`を使うには、`st::defualt_random_engine`を置き換えるだけでいい。

~~~cpp
int main()
{
std::mt19937 e ;
for ( int i = 0 ; i != 10 ; ++i )
{
std::cout << e() << "\n"sv ;
}
}
~~~

メルセンヌツイスターはとても優秀な乱数エンジンだ。乱数が必要な多くの場面では、メルセンヌツイスターを使っておけばまず問題はない。

では乱数エンジンを使って、生の乱数を標準入力で得た個数だけ出力するプログラムを書いてみよう。

~~~cpp

int main()
{
// 乱数エンジン
std::mt19937 e ;

// 標準入力からnを得る
unsigned n {} ;
std::cin >> n ;
// n個出力
for ( unsigned int i = 0 ; i != n ; ++i )
{
std::cout << e() << " "sv ;
}
}
~~~

実行結果は以下のようになる。

~~~
$ dice
10
3499211612 581869302 3890346734 3586334585 545404204 4161255391 3922919429 949333985 2715962298 1323567403
~~~

乱数エンジンで生成されるのは生の乱数だ。これは通常、32bit符号なし整数とか64bit符号なし整数で表現できる全範囲の値として生成される。これは実際に必要な乱数とは値の範囲が違う。実際に必要な乱数とは、例えば6面ダイスの場合は、`int`型で1,2,3,4,5,6のいずれかの値がそれぞれ$\frac{1}{6}$の確率で出てほしい。

## 乱数分布

`乱数分布`とは生の乱数を望みの範囲の乱数に加工するためのライブラリだ。クラスで実装されている。

乱数分布ライブラリにも様々なものがあるが、6面ダイスのプログラムを実装するのに使うのは`std::uniform_int_distribution<T>`だ。

この乱数文法ライブラリは、`T`にほしい乱数の整数型を指定する。コンストラクター引数を2つ取るので、1つ目の引数に最小値、2つ目の引数に最大値を指定する。

~~~c++
std::uniform_int_distribution<int> d(a, b) ;
~~~
この乱数分布クラスの変数`d`は、$a \leq r \leq b$までの範囲の乱数`r`を作り出す。
6面ダイスを作るには、`d(a, b)`を`d(1, 6)`にすればよい。
~~~c++
std::uniform_int_distribution<int> d(1, 6) ;
~~~

乱数分布クラスのオブジェクト`d`を作ったならば、`operator()`に乱数エンジンのオブジェクトを引数に渡すことで乱数が作れる。乱数エンジンのオブジェクトを`e`とすると、`d(e)`だ。

~~~cpp
template < typename Engine, typename Distribution >
void f( Engine & e, Distribution d)
{
// 乱数
auto r1 = d(e) ;
// 次の乱数
auto r2 = d(e) ;
// 次の乱数
auto r3 = d(e) ;
}
~~~
以上の知識を利用して、プログラム`dice`を作ってみよう。
~~~cpp
int main()
{
// 乱数エンジン
std::mt19937 e ;
// 乱数分布
std::uniform_int_distribution<int> d(1, 6) ;
// 入力を処理
unsigned n {} ;
std::cin >> n ;
for ( unsigned int i = 0 ; i != n ; ++i )
{
// 乱数出力
std::cout << d(e) << " "sv ;
}
}
~~~

早速実行してみよう。

~~~
$ ./dice
5
5 1 6 6 1
$ ./dice
10
5 1 6 6 1 6 6 2 4 2
$ ./dice
20
5 1 6 6 1 6 6 2 4 2 1 4 2 2 4 6 6 6 6 6
~~~

確かに動く。しかし毎回同じ出力になる。これでは実用的な6面ダイスプログラムとは言えない。

## シード

線形合同法を思い出してみよう。線形合同法で次の乱数$X_{n+1}$を計算するには、今の乱数$X_{n}$に対して$X_{n+1} = (a \times X_{n} + c) \bmod m$という計算をする。

線形合同法とは現在の乱数値を内部状態として持ち、そこに計算を加えることで次の乱数を作り出すのだ。

一般化すると、疑似乱数は内部状態$S_n$を持ち、計算を加える関数$f(x)$を適用することで、次の内部状態$S_{n+1}=f(S_n)$を作り出すのだ。単純な線形合同法の場合、内部状態がそのまま乱数の値になるが、複雑な疑似乱数アルゴリズムでは、内部状態から乱数を求めるのにさらに計算を加えるものもある。


2 changes: 1 addition & 1 deletion Makefile
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@ book : docs/index.html


docs/index.html : *.md style.css
pandoc -s --toc --toc-depth=6 --mathjax -o $@ -H style.css pandoc_title_block *-*.md
pandoc -s --toc --toc-depth=6 "--mathjax=https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/latest.js?config=TeX-MML-AM_CHTML" -o $@ -H style.css pandoc_title_block *-*.md

filter.json : *.md style.css
pandoc -t json -s --toc --toc-depth=6 --mathjax -o $@ -H style.css pandoc_title_block *-*.md
Expand Down
Loading

0 comments on commit d9aad45

Please sign in to comment.