擬似乱数とは何か

擬似乱数とは何か

主にコンピューターゲームにおいて「擬似乱数」というものがよく使われます。
この「擬似乱数」というのは一体どういうものなのか、気になったことはあり
ませんか？そんな身近な「擬似乱数」の仕組みをざっくりと解説します。

なるべく技術的な要素は減らして分かりやすく、それでいて本質的な部分も知
れるようにしました。雑学気分で気楽にどうぞ。

• 擬似乱数の原理
• 二進数
• 乱数列のサイズ
• 乱数列を計算で作る
• 擬似乱数のアルゴリズム
• 周期を最大にする意味
• 線形合同法の条件
• 乱数を使うには注意がいる
• 乱数は二度振るな
• 乱数は二度振るな、再び
• 組み合わせの数
• シャッフル問題
• k次元の均等分布
• 大きな擬似乱数の初期化は難しい
• 間違った初期化
• 乱数の確率分布
  • 一様分布
  • 重み付け分布
  • 正規分布
  • サイコロの和
• 様々な擬似乱数アルゴリズム
  • 線形合同法
  • 線形帰還シフトレジスタ (LFSR)
    • ドルアーガの塔のLFSR
  • GFSR
  • メルセンヌ・ツイスタ
  • Xorshift
  • PCG
• 最後に
• 改訂履歴
• 利用許諾

擬似乱数の原理

次のような数列があるとします。0から7までの連番です。

0	1	2	3	4	5	6	7

数列をランダムに並び替えます。

5	2	3	0	1	6	7	4

そして、先頭から順に数値を取り出していくと、擬似乱数になります。これが
擬似乱数の基本的な原理です。

この乱数列の長さを周期と言い、乱数列の開始地点をシード値と言います。取
り出された値の範囲は、出力幅・出力ビット・ビット幅などと言います。

周期=8	乱数列	シード値=3
⏫️	5
⏫️	2
⏫️	3	⏪️シード値
周	0
期	1
⏬️	6
⏬️	7
⏬️	4

この例では、周期が8の乱数で、出力幅が0～7、シード値が3の擬似乱数である。
と言えます。

マイクラをやっていればご存知かもしれませんが、擬似乱数のシード値を変え
ることで、違った地形を生み出すといった事ができます。実際の擬似乱数はもっ
と巨大な乱数列になるため、乱数列上の開始位置を変えるだけで大きな変化が
起こるように見えるのです。

二進数

擬似乱数は、コンピューター上で二進数を区切りとした取り扱いをするので、
説明を進める前に二進数について少し知っておきましょう。

コンピューターではアナログな数値は扱わず、0と1だけで表現した数値を扱い
ます。いわゆるデジタル処理というものです。

0から7の十進数の数値の場合、二進数では次のように表現されます。

十進数	二進数
0	000
1	001
2	010
3	011
4	100
5	101
6	110
7	111

0～7の数値を表現するには、0と1だけを組み合わせて3桁の二進数で表現でき
ます。この場合は3ビットの数値であると言えます。

コンピューターでは8ビットを1バイトという単位でメモリの容量を表します。
1GB（ギガバイト）といった単位には聞き覚えがあると思います。

ビット数が多ければより大きな範囲の数値を扱える。と、取りあえずは覚えて
おいてください。

乱数列のサイズ

先程の例では、3ビットの擬似乱数を説明しましたが、実際にはもっと大きな
擬似乱数を使用します。1999年頃までは32ビットの擬似乱数が主流でした。こ
れに当てはめて乱数列がどのくらいのサイズになるかを考えてみましょう。

32ビットの乱数という事は、2の32乗の周期を持ち、乱数列もこの長さになる
という事です。数値の出力幅が32ビットだとすると1つあたり4バイトの容量に
なります。

周期	2の32乗	4,294,967,296（約43億）
出力幅	32ビット	4 バイト
合計サイズ		17,179,869,184 バイト（16GB）

これらを掛け算すると、32ビットの乱数列は16GBという巨大なサイズになりま
す。平均的なゲーミングPCでメモリの容量が16GB前後なので、この乱数列をメ
モリ上に保持するのはかなり難しいでしょう。

2024年時点では、128ビットや256ビットの擬似乱数が主流になってきており、
乱数列は想像を絶する巨大なサイズになります。

乱数列を計算で作る

乱数列をメモリに保持するのは難しい。ではどうするのか。計算によって乱数
列を生成していきます。値を一つだけ保持しておき、乱数を取り出す時に計算
で次の値を求め、値を更新していきます。

振り数	値の変化	計算式	説明
	5		最初の値がシード値となる
1回目	2	= (5 * 5 + 1) mod 8	(前回の値 * 5 + 1) mod 8
2回目	3	= (2 * 5 + 1) mod 8
3回目	0	= (3 * 5 + 1) mod 8
4回目	1	= (0 * 5 + 1) mod 8
5回目	6	= (1 * 5 + 1) mod 8
6回目	7	= (6 * 5 + 1) mod 8
7回目	4	= (7 * 5 + 1) mod 8
8回目	5	= (4 * 5 + 1) mod 8	一周すると最初の値に戻る

• 「*」は掛け算。プログラミング言語や表計算ソフトで使われる記号です。
• 「mod」は剰余。この例では8で割った余りとなります。Windowsの電卓を関
  数電卓モードにして26 Mod 8=と入力すれば計算できます。表計算ソフト
  ならば=MOD(26, 8)と入力すれば計算できます。

前回の値に同じ計算を繰り返していくことで、乱数列を再現することができま
す。1つの値だけ保持しておけば良いので、必要な容量は3ビットで済みます。

しかし、適当に計算しているように見えるのに、なぜ0～7までの値がそれぞれ
1回ずつ現れるのか。綺麗に一周して元に戻るのか。不思議ですよね。

適当な数で乱数列を生成するとどうなるのか、試してみましょう。掛ける値を
2に変更します。

振り数	値の変化	計算式	説明
	5		シード値
1回目	3	= (5 * 2 + 1) mod 8	(前回の値 * 2 + 1) mod 8
2回目	7	= (3 * 2 + 1) mod 8
3回目	7	= (7 * 2 + 1) mod 8
4回目	7	= (7 * 2 + 1) mod 8	値が7で固定されてしまう

5,3,7と変化し、以降はずっと7しか出力されなくなってしまいました。これで
は擬似乱数としては使えません。擬似乱数として成立させるには、擬似乱数の
アルゴリズム（演算法）に基づいて計算のパラメーターを選ぶ必要があります。

擬似乱数のアルゴリズム

擬似乱数の演算は、数学的な理論と証明によって成り立っています。これまで
の例に出てきた計算方法は「線形合同法」という擬似乱数のアルゴリズムに基
づいています。

擬似乱数のアルゴリズムを使うことで、周期の中で全ての値が1回だけ現れる。
という事が保証されます。例えば周期が8の場合。8回乱数を振ると
0,1,2,3,4,5,6,7がそれぞれ1回だけ現れます。この状態を「周期を最大にする」
と言います。

周期を最大にする意味

なぜ、周期を最大にする必要があるのか。例えばサイコロを振った時。1～6の
値が同じ確率で出るのを期待するでしょう。もし4だけが出ない、となったら
とても困ります。

擬似乱数ではもっと大きな範囲の値を出力しますが、実際に利用する場合には、
この値を加工して1～6までの値を出力できるように計算して調整します。

サイコロの出目（1～6）とするには、次のように計算します。

サイコロの出目 = (乱数値 mod 6) + 1

周期が12で、3と9が出現しない乱数、という前提で計算します。分かりやすく
するため、順序はランダムではなく昇順にしてあります。

乱数値	計算式	出目
0	(0 mod 6) + 1 =	1
1	(1 mod 6) + 1 =	2
2	(2 mod 6) + 1 =	3
4	(4 mod 6) + 1 =	5
5	(5 mod 6) + 1 =	6
6	(6 mod 6) + 1 =	1
7	(7 mod 6) + 1 =	2
8	(8 mod 6) + 1 =	3
10	(10 mod 6) + 1 =	5
11	(11 mod 6) + 1 =	6
12	(12 mod 6) + 1 =	1
13	(13 mod 6) + 1 =	2

出目を集計すると次のようになります。

出目	出現回数
1	3	⚠️
2	3	⚠️
3	2
4	0	⚠️
5	2
6	2

4が出現せず、1と2が多く出現する。といった結果になりました。周期が12で
あれば、本来はそれぞれ2回ずつ出現するはずです。

元の乱数列に特定の値が出ない状態があると、どのように加工しても特定の値
が出なくなったり、出現の確率に偏りが出てしまう。という事が起こるのです。

擬似乱数においては、全ての値が1回ずつ現れ、その順序だけがランダムであ
る。ということが重要な前提条件となります。

線形合同法の条件

線形合同法は1999年頃までよく使われたアルゴリズムで古いのですが、シンプ
ルで分かりやすいので、少し解説をします。

線形合同法は次のような式で表します。

x = (x * a + c) mod m

プログラミング言語に近い書き方で少しややこしいのですが、現在のxの値か
ら計算してxの値を書き換える、といった意味になります。この中のa,c,mとい
うパラメーターを調整することで、周期を最大にする事ができます。

周期を最大にするには以下の条件があります。

• mは2の累乗である。
• cは0ではない。
• mとcは互いに素である。（最大公約数が1である）
• a - 1は、mの全ての素因数で割り切れる。
• a - 1は、4で割り切れる。

これに当てはめて、これまでの例では a=5,c=1,m=8というパラメーターを選択
して擬似乱数となるようにしていました。

これは線形合同法のいくつかある条件の一つです。他にも成立する条件があり、
そのパラメーターの選び方によって乱数としての品質が変わります。

乱数を使うには注意がいる

実際に擬似乱数を使ってサイコロの出目（1～6）を作る場合を考えてみましょ
う。

サイコロの出目を作るには以下のように計算します。

サイコロの出目 = (乱数値 mod 6) + 1

乱数の周期は8です。分かりやすくするため、順序はランダムではなく昇順に
してあります。

乱数値	計算式	出目
0	(0 mod 6) + 1 =	1
1	(1 mod 6) + 1 =	2
2	(2 mod 6) + 1 =	3
3	(3 mod 6) + 1 =	4
4	(4 mod 6) + 1 =	5
5	(5 mod 6) + 1 =	6
6	(6 mod 6) + 1 =	1	⚠️
7	(7 mod 6) + 1 =	2	⚠️

出目を集計すると以下のようになります。

出目	出現回数
1	2	⚠️
2	2	⚠️
3	1
4	1
5	1
6	1

1と2の出現回数に偏りが出ています。使う擬似乱数の出力幅によっては6では
割り切れないため、出現回数に差が出てしまいます。単純に剰余を使っただけ
では正しい確率のサイコロは作れません。

では、どうするのか？乱数を振って6と7が出た場合には、もう一度乱数を振っ
てやり直します。再び6,7が出ても、出なくなるまで振り直しをします。0～5
が出れば計算し、出目として返します。

乱数値	計算式	出目
0	(0 mod 6) + 1 =	1
1	(1 mod 6) + 1 =	2
2	(2 mod 6) + 1 =	3
3	(3 mod 6) + 1 =	4
4	(4 mod 6) + 1 =	5
5	(5 mod 6) + 1 =	6
6	(6 mod 6) + 1 =	1	🚫
7	(7 mod 6) + 1 =	2	🚫

出目	出現回数
1	1
2	1
3	1
4	1
5	1
6	1

こうする事で出目の重複を防ぎ、均等な確率のサイコロを再現できます。

乱数値が非常に大きな値で、出目の値が小さい場合には無視できるほど偏りが
小さくなることもありますが、出目の値の範囲を大きくしていくと偏りも大き
くなっていきます。

一見簡単な問題のように見えても、擬似乱数を使って乱数を作る場合には色々
と注意すべき点があり、油断していると思わぬバグの原因となります。

C++の標準ライブラリにはstd::uniform_int_distributionというクラスがあ
り、これを使うと均等な確率の乱数値を得られます。

乱数は二度振るな

平面上にランダムに点を打つ、という使い方を考えてみましょう。x,yの座標
があり、x,yそれぞれに乱数を生成すれば実現できます。

まず、rand()という関数を定義します。使う毎に乱数を生成して返す。そうい
うものだと考えてください。前提条件は次のようにします。

• rand()は0～15の値を返す擬似乱数
• rand()の周期は16
• x,yの範囲は0～3

この条件でランダムな点を生成するなら、次のように書けるでしょう。

x = rand() mod 4
y = rand() mod 4

明快ですね。あとは点を打つプログラムにx,yを渡せば完了です。

しかし、これには問題があります。どういった問題が起こるのか考えていきま
す。まずは全ての出現パターンを書き出します。分かりやすくするため乱数の
順序は昇順にします。

乱数値	x	y	(x, y)
0	0 mod 4 = 0
1		1 mod 4 = 1	(0, 1)	⚠️
2	2 mod 4 = 2
3		3 mod 4 = 3	(2, 3)	⚠️
4	4 mod 4 = 0
5		5 mod 4 = 1	(0, 1)	⚠️
6	6 mod 4 = 2
7		7 mod 4 = 3	(2, 3)	⚠️
8	8 mod 4 = 0
9		9 mod 4 = 1	(0, 1)	⚠️
10	10 mod 4 = 2
11		11 mod 4 = 3	(2, 3)	⚠️
12	12 mod 4 = 0
13		13 mod 4 = 1	(0, 1)	⚠️
14	14 mod 4 = 2
15		15 mod 4 = 3	(2, 3)	⚠️

何かがおかしいですね。xには0と2しか出現せず、yには1と3しか出現しません。
(0,0) (0,2) (0,3)といった一部の組み合わせが1度も出現しません。これが全
てのパターンなので、何度乱数を回しても結果は同じです。これはとても乱数
と呼べる振る舞いではありません。

「周期を最大にする意味」で説明しましたが、擬似乱数は全ての値が1回ずつ
出現するように設計されています。今回のように2回乱数を振ってx,yに割り当
てると、それぞれが乱数値を1つ飛ばしで取ることになってしまいます。折角
用意された全ての値のうち半分しか使えていないのです。

擬似乱数のアルゴリズムでは、隣り合う乱数値の組み合わせが、全て出る事は
保証されません。周期を大きくして、取り出す値を小さくしていけば、特定の
組み合わせが出ない確率は低くなります。しかし、どこまで小さくすれば安全
になるかは、はっきりとは判らないのです。

では、どうすればよいのか？1つの乱数値からx,y座標を生成するように計算し
ます。xは4で割った余りを。yは4で割って小数点以下を切り捨てます。

r = rand()
x = r mod 4
y = floor(r / 4)

• floor()は小数点以下を切り捨てる、という意味の関数です。

結果は以下のようになります。

乱数値	x	y	(x, y)
0	0 mod 4 = 0	floor(0 / 4) = 0	(0, 0)
1	1 mod 4 = 1	floor(1 / 4) = 0	(1, 0)
2	2 mod 4 = 2	floor(2 / 4) = 0	(2, 0)
3	3 mod 4 = 3	floor(3 / 4) = 0	(3, 0)
4	4 mod 4 = 0	floor(4 / 4) = 1	(0, 1)
5	5 mod 4 = 1	floor(5 / 4) = 1	(1, 1)
6	6 mod 4 = 2	floor(6 / 4) = 1	(2, 1)
7	7 mod 4 = 3	floor(7 / 4) = 1	(3, 1)
8	8 mod 4 = 0	floor(8 / 4) = 2	(0, 2)
9	9 mod 4 = 1	floor(9 / 4) = 2	(1, 2)
10	10 mod 4 = 2	floor(10 / 4) = 2	(2, 2)
11	11 mod 4 = 3	floor(11 / 4) = 2	(3, 2)
12	12 mod 4 = 0	floor(12 / 4) = 3	(0, 3)
13	13 mod 4 = 1	floor(13 / 4) = 3	(1, 3)
14	14 mod 4 = 2	floor(14 / 4) = 3	(2, 3)
15	15 mod 4 = 3	floor(15 / 4) = 3	(3, 3)

全てのx,yの組み合わせが生成されるようになりました。実際の乱数値は順序
がランダムであるので、これでランダムに点を打つ、という目的を達成できま
す。

x,yの求め方が少し難しいですが、二進数で考えると分かりやすくなります。
0～15の値は4ビットで表せます。4ビットのうち上位2ビットをyに、下位2ビッ
トをxに割り当てる。というのが上記の計算の仕組みです。2ビットの値の範囲
は0～3になるのでx,yの範囲と丁度一致します。

例えば、乱数値が14の場合は以下のようになります。

十進数	二進数		上位	下位
14	1110	⏩️	11	10
			⏬️	⏬️
			y=3	x=2

この例では、周期とx,yの組み合わせが丁度同じになるようにしました。しか
し、この周期では充分ではありません。真の乱数であれば、16回乱数を回した
場合。x,y=(0,0)の出る確率が64%になります。

$$1 - \bigg(1 - \frac{1}{16}\bigg)^{16} = 0.643925\ldots \approx 64％$$

ですが、この例では、16回乱数を回した場合。x,y=(0,0)が100%出てしまいま
す。乱数らしい振る舞いをするためには、もっと長い周期の擬似乱数を使う必
要があります。

また、擬似乱数の出力幅が、x,yの組み合わせの数で割り切れない場合には、
出現率に偏りが発生します「乱数を使うには注意がいる」で説明したように、
この偏りも考慮する必要があります。

今回は2次元で極狭い範囲の乱数生成でしたが、もっと広い範囲の値や多次元
になってくると、より大きな周期や出力幅をもつ擬似乱数が必要になってきま
す。

多数の要素を持つランダムな物を生成する時には、その組み合わせの数と、擬
似乱数の周期や出力値の範囲を、慎重に計算しなければなりません。その許容
範囲を超えた場合には、特定の組み合わせが絶対に出現しない、といったバグ
も起こり得るのです。

全てのパラメーターが最高値のアイテムを求めて敵を狩りまくる。というのは
ゲームではよくある事です。もし、絶対にでない組み合わせの中に、全てが最
高値のアイテムが含まれていたとしたら、無駄な努力になってしまいます。プ
レイヤーにとってこれほど恐ろしい事はありません。

ただ、擬似乱数のアルゴリズムによっては、複数回乱数を振っても組み合わせ
の安全性が保証されるものも一部あります。これについては後述します。

乱数は二度振るな、再び

乱数を二度振ると、組み合わせの問題が起こる事を学びました。しかし、気を
付けていても意図しないところで、乱数を二度振っている事もあります。次の
ような例を考えてみましょう。

• 攻撃力、防御力のパラメーターを持つアイテムを生成する。
  攻撃力、防御力の範囲は0～3のランダム。
• アイテム生成後、画面のランダムな位置x,yに花火を表示する。
  x,yの範囲は0～3のランダム。
• 周期が16、出力値の範囲が0～15の擬似乱数を使う。

結果は以下のようになります。

乱数値	下位	上位	アイテム (攻, 防)	花火 (x, y)
0	0 mod 4 = 0	floor(0 / 4) = 0	(0, 0)		⚠️
1	1 mod 4 = 1	floor(1 / 4) = 0		(1, 0)
2	2 mod 4 = 2	floor(2 / 4) = 0	(2, 0)		⚠️
3	3 mod 4 = 3	floor(3 / 4) = 0		(3, 0)
4	4 mod 4 = 0	floor(4 / 4) = 1	(0, 1)		⚠️
5	5 mod 4 = 1	floor(5 / 4) = 1		(1, 1)
6	6 mod 4 = 2	floor(6 / 4) = 1	(2, 1)		⚠️
7	7 mod 4 = 3	floor(7 / 4) = 1		(3, 1)
8	8 mod 4 = 0	floor(8 / 4) = 2	(0, 2)		⚠️
9	9 mod 4 = 1	floor(9 / 4) = 2		(1, 2)
10	10 mod 4 = 2	floor(10 / 4) = 2	(2, 2)		⚠️
11	11 mod 4 = 3	floor(11 / 4) = 2		(3, 2)
12	12 mod 4 = 0	floor(12 / 4) = 3	(0, 3)		⚠️
13	13 mod 4 = 1	floor(13 / 4) = 3		(1, 3)
14	14 mod 4 = 2	floor(14 / 4) = 3	(2, 3)		⚠️
15	15 mod 4 = 3	floor(15 / 4) = 3		(3, 3)

アイテムの生成と花火の演出で、同じ乱数を使っているため、アイテムの一部
の組み合わせが出現しないようになってしまいました。それぞれの機能を別々
の担当者が制作していると、このように意図せず乱数を二度振った状態になる
事もあります。

この場合は、アイテムの生成と花火の演出で使う擬似乱数を、二つに分ける事
が必要です。組み合わせの欠落は絶対にあってはならない。という用途の場合
には、他に流用されない独立した擬似乱数を持つ必要があります。

この例は単純なものですが、規模が大きくなるとより深刻になってきます。

例えば、1秒間に2回攻撃するキャラクターがサーバー内に30人居て、攻撃の都
度乱数を引くとすると、1秒間に60もの乱数を消費します。共通の擬似乱数を
使っていた場合、アイテム生成で使える乱数はかなり飛び飛びとなってしまい
ます。これでは擬似乱数本来の性能を発揮することは難しいでしょう。

オンラインのシステムでは沢山のユーザーが同時アクセスするため、システム
やユーザー毎の乱数に、どのくらいの独立性を持たせるのかも考慮する必要が
あります。

組み合わせの数

擬似乱数を使ってランダムな物を生成する場合には、その要素の組み合わせ数
と擬似乱数の周期や出力幅を考慮しなければなりません。実際に計算してみま
しょう。次のような条件とします。

• 攻撃力、防御力のパラメーターを持つアイテムを生成する。
• 攻撃力、防御力の範囲は0～3のランダム。

要素	範囲	組み合わせ数
攻撃力	0～3	4
防御力	0～3	4
合計		4 * 4 = 16

攻撃力が4種類、防御力が4種類の組み合わせで合計で16種類の組み合わせとな
ります。擬似乱数に必要な性能は周期が16以上、出力幅も16以上が必要と分か
ります。

これは小さい範囲の例なので簡単に扱えますが、実際にはもっと大きな数にな
るため数値では分かりにくくなってきます。なので、二進数で何ビット必要な
のか、というのを求めていきます。

次のような例を考えます。

• STR,DEX,VIT,AGI,INT,MND,CHRという7種類のパラメーターを持つアイテム生
  成する。
• パラメーターの値の範囲は0～15のランダム

7種類のパラメーターに、数値の組み合わせが16。なので、16の7乗が組み合わ
せ数となります。

16^7 = 268,435,456

• 「^」はべき乗を意味します。数学ではなくプログラミング言語や表計算ソ
  フトで使われる記号です。
• Windowsの電卓を関数電卓モードにして16^7=と入力すれば同様の計算がで
  きます。表計算ソフトならば=16^7と入力すれば計算する事ができます。

かなり大きな組み合わせ数になりました。これを表現するには二進数で何ビッ
ト必要なのかを計算します。

log(268,435,456) / log(2) = 28

• log()は対数関数です。難しい数学関数が出てきたな、と思うかもしれませ
  んが大丈夫です。このように計算すれば何ビットになるのかを計算できる。
  そういうものだと覚えておけばよいです。
• Windowsの電卓を関数電卓モードにして268435456 log / 2 log =と入力す
  れば同様の計算ができます。表計算ソフトならば=LOG(268435456, 2)と入
  力すれば計算する事ができます。

28とでました。28ビット以上の周期と出力幅を持つ擬似乱数が必要だと分かり
ます。

「乱数は二度振るな」でも説明しましたが、ビット数の要件を満たしていて
も、乱数を複数回振って生成した場合には、組み合わせ問題は解決しない事
があるので注意してください。

パラメーターが増えていくと、どんどんと組み合わせ数が増えていきます。パ
ラメーターが32種類、値の範囲を0～15としてみましょう。

log(16^32) / log(2) = 128

128ビットも必要になります。擬似乱数のアルゴリズムによっては256ビット以
上のものもありますが、周期は256ビットでも、出力幅は64ビットまで。とい
うのが一般的です。組み合わせ数が巨大になると、高性能な擬似乱数を持って
しても実現するのが困難になってきます。

シャッフル問題

一見簡単そうに見える処理でも、大きな問題を含んでいる事があります。典型
的な例がカードをシャッフルする問題です。新品のトランプのカードを使った
時、充分に混ぜ合わせるまで苦労した覚えはないでしょうか。擬似乱数を使っ
てカードをシャッフルするには、どのくらいの組み合わせが必要になるのか考
えてみましょう。

トランプのカードの場合。ジョーカー1枚を加えると、53枚のカードになりま
す。これをシャッフルするには次のような操作を行います。

1. 53枚の中からランダムに1枚選び、取り出す。
2. 52枚の中からランダムに1枚選び、取り出す。
   ...
3. 上記の操作を最後の1枚まで繰り返す。

とても簡単ですね。これは「フィッシャー・イェーツのシャッフル」と呼ばれ
るアルゴリズムです。この組み合わせを計算すると以下のようになります。

組み合わせ数 = 53 * 52 * 51 * ... * 1

これは階乗というものです。1ずつ減らして1までの値を全て掛け算します。数
学の書き方では「!」を使って以下のよう書きます。

組み合わせ数 = 53!

• 階乗はWindowsの電卓を関数電卓モードにして53!と入力すると計算できま
  す。表計算ソフトならば=FACT(53)で計算できます。

この組み合わせに何ビット必要なのかを計算すると、以下のようになります。

log(53!) / log(2) = 231.3089

232ビット以上の擬似乱数が必要です。

別の例を考えてみましょう。麻雀牌の場合です。麻雀牌は136個あります。同
じ種類の牌を含んでいますが、それを考慮しても組み合わせ数はそれほど減り
ません。計算やプログラミングがとても複雑になるので、136種類として計算
をします。

log(136!) / log(2) = 772.5588

773ビット以上の擬似乱数が必要になります。

このようにシャッフルにはとても巨大な組み合わせ問題があり、1999年頃まで
は、擬似乱数を使って完全なシャッフルを行うのは不可能である、と考えられ
ていました。しかし「メルセンヌ・ツイスタ」という擬似乱数が登場し、状況
が変わっていきます。

k次元の均等分布

擬似乱数の性能の指標に「k次元の均等分布」というものがあります。難しい
言い方になっていますが、簡単に言うとk回乱数を振って1つの要素とした場合、
全ての組み合わせが同じ回数出現しますよ。ということです。k=2なら2回乱数
を振ってx,yにしても組み合わせの安全が保証されます。

いくつかの擬似乱数アルゴリズムを比較してみます。

擬似乱数アルゴリズム	k次元	周期	出力幅	主な実装
線形合同法	1	2^31	15ビット	MSVC
Xorshift128+	1	2^128-1	53ビット	Webブラウザ
Xoshiro128**	4	2^128-1	31ビット	.NET 6 (x86)
Xoshiro256**	4	2^256-1	31ビット	.NET 6 (x64)
メルセンヌ・ツイスタ	623	2^19937-1	32ビット	C++

kが少ない擬似乱数が多いですが、メルセンヌ・ツイスタだけが圧倒的です。k
が623もあります。これがあると何ができるのかというと、シャッフル問題を
解決できます。

麻雀では136個の牌をシャッフルしますが、メルセンヌ・ツイスタであれば単
純に135回乱数を振ってランダムに1つずつ牌を選ぶだけで、完全なシャッフル
を実現することができます。623のうち135ですから、これでもまだまだ余裕が
あります。

高次元の均等分布だけが擬似乱数の性能の良さ、という訳ではありませんが、
多次元の要素を乱数で生成するには、今までは慎重に慎重を重ねて設計する必
要がありました。メルセンヌ・ツイスタの登場によりその手間が大きく軽減さ
れたので、これは乱数の革命と言っても過言ではありません。

大きな擬似乱数の初期化は難しい

擬似乱数を使う前には初期化が必要です。線形合同法のような32ビット程度の
小さい乱数であれば、32ビットの数値を用意してセットすれば、それがそのま
まシード値となって初期化できます。256ビットの擬似乱数であれば、32ビッ
トの数値を8個用意する必要があります。

しかし、固定の数値で初期化を行うと、擬似乱数は必ず同じパターンの乱数を
返します。ゲームを起動した直後は、毎回同じカードが配られる。といった状
況は困りますよね。

起動直後でも擬似乱数をランダムに初期化するにはどうすればよいのか？それ
には「暗号論的擬似乱数生成器」というものを使います。また難しい言葉が出
てきました。これは、暗号通信の際に必要となる乱数を生成するために作られ
た専用の擬似乱数です。

暗号通信といっても今はとても身近です。httpsで始まるURLのWebサイトは全
て暗号通信になっています。暗号通信では一時的に乱数を交換することで、安
全に通信できるようなアルゴリズムになっていて、その乱数は予測できないも
のである必要があるのです。そのために暗号アルゴリズム等を使って乱数を生
成し、予測されないように工夫がされています。

これまで使おうとしてきた線形合同法等の擬似乱数は、高速に動作し、再現性
がある（シード値による再現性）事に重きが置かれています。再現性があれば
シミュレーション等で繰り返しテストができ便利ですが、再現性がある事で予
測可能性も高まってしまうのです。

「暗号論的擬似乱数生成器」は暗号通信に必要なので、OSの機能として必ず備
わっています。WindowsであればBCryptGenRandom()。Linuxであれば
getrandom()等といったAPIから暗号論的擬似乱数を得る事ができます。これ
らを使って擬似乱数を初期化すれば、毎回違った状態でスタートできる訳です。

しかし、メルセンヌ・ツイスタのような巨大な擬似乱数の登場で少し注意が必
要になりました。メルセンヌ・ツイスタは簡単に言うと19,937ビットの乱数で
す。メモリとしては2,496バイトもの容量を持ちます。これを「暗号論的擬似
乱数生成器」で完全に初期化ができるのか？という疑問があります。

「暗号論的擬似乱数生成器」はAESのような暗号アルゴリズムを使って実装さ
れています。AESは128～256ビットの鍵を使って、128ビットのデーターを暗号
化、復号ができます。「暗号論的擬似乱数生成器」の規格としてはNIST SP
800-90A/Bがありますが、OSによっては独自の手法を使っている事もあります。

一例として、以下のような仕組みで「暗号論的擬似乱数」が生成されます。

1. OS/CPU等のノイズデーターから256ビットのkeyを作る。
2. OS/CPU等のノイズデーターから64ビットのnonceを作る。
3. 64ビットのcounterを0で作る。
4. counterとnonceを繋げて128ビットのblockを作る。
5. keyを使ってblockをAESで暗号化する。
6. 暗号化されたblockを乱数として出力する。
7. counterを1増加させる。

規格に沿ったものではなく正確ではないですが、大まかにはこういった動作に
なると思います。少ないノイズデータを使い、それを暗号で引き伸ばして予測
のできない擬似乱数としています。

ここからAESによる擬似乱数はどれだけの組み合わせを出力できるのかを考え
ると、次のようになります。

256ビット(key) + 128ビット(block) = 384ビット

key,nonce,counterがどのタイミングで更新されるかにもよりますが、最大で
384ビット相当の組み合わせを持てると考えられます。メルセンヌ・ツイスタ
が19,937ビット相当として比較すると、充分とは言えません。これはシャッフ
ル問題と同じ構造です。

どういう事かというと、起動する度にメルセンヌ・ツイスタを初期化した場合。
「暗号論的擬似乱数生成器」を使っても完全に初期化することはできない。と
いう事です。

麻雀牌のシャッフルのように773ビット相当の組み合わせを持つ場合には「暗
号論的擬似乱数生成器」のビット数を超えることもあります。

組み合わせの数を器のように見立てて移し替えていく、と考えると分かりやす
いでしょう。暗号論的擬似乱数生成器で擬似乱数を初期化した後、1回目に出
てくる組み合わせは暗号論的擬似乱数生成器の範囲に限定されます。

暗号論的擬似乱数生成器   384/384 |########|
メルセンヌ・ツイスタ   384/19937 |########..............................|
麻雀牌のシャッフル       384/773 |########........|

「暗号論的擬似乱数生成器」も充分に大きな擬似乱数なので、実用上問題にな
る事はほとんどないと思います。ただ、完全に初期化はできない。という事は
頭の片隅に留めておくとよいでしょう。

知っていれば、毎回メルセンヌ・ツイスタを初期化するべきか、内部状態を保
存して使い続けるべきか、そういった選択をする事ができます。

間違った初期化

メルセンヌ・ツイスタのような大きな内部状態をもつ擬似乱数は、初期化にも
大きな乱数源を必要とします。完全に初期化するには、32ビットの数値で624
個必要になります。

これを用意するのは大変なので、メルセンヌ・ツイスタには少ない乱数源であっ
ても、それを引き伸ばして、なるべく偏りが出ないように内部状態を初期化す
る機能が備わっています。32ビットの数値であれば1つからでも初期化できま
す。

初期化が簡単で使い易いがために、間違った初期化をされてしまう事がありま
す。32ビットの数値1つだけを使ってメルセンヌ・ツイスタを初期化をすると、
どのような問題が起こるのでしょうか。

間違った初期化プログラムの例 (C++)

// 暗号論的擬似乱数生成器
// 呼び出すごとに32ビットの乱数を1つ返す
std::random_device rd;

// メルセンヌ・ツイスタ (mt19937)
// 32ビットの乱数を1つ渡し、初期化する
std::mt19937 mt(rd());

サンプルコードで稀に見かけますが、これは良くない使い方です。

メルセンヌ・ツイスタは19,937ビット相当の内部状態を持つので、32ビットの
乱数値1つで初期化しても、その組み合わせのうちほんの一部しか再現できま
せん。

32ビットの数値もそれなりに大きな数値ではあるのですが、意外と同じ数値が
出やすいのです。同じ数値が出ることを衝突と言います。どのくらいの頻度で
衝突が出るのか計算してみましょう。計算するのは難しいのですが、近似する
公式があります。

衝突までの平均回数 = sqrt(π / 2 * H)

• sqrt() 平方根を求める関数です。数学記号では√です。
• π 円周率です。
• H 集合の数です。ここでは32ビットの数値なので2^32です。

これに当てはめると、

sqrt(π / 2 * 2^32) = 82137.1953

およそ8万回に1回は衝突が起こるという事になります。意外と少なく感じます
よね。

擬似乱数の初期化で衝突が起こるという事は、例えばオンラインのカードゲー
ムで沢山の部屋が生成され、部屋毎に毎回擬似乱数を初期化しているような場
合。全く同じカード配置でゲームが開始される。という事が8万回に1回起こる
という事です。

組み合わせを器のように移し替えていくイメージで考えると分かりやすいでしょ
う。

暗号論的擬似乱数生成器   32/384 |#.......|
メルセンヌ・ツイスタ   32/19937 |#.....................................|
カードのシャッフル       32/232 |#....|

暗号論的擬似乱数生成器はもっと大きな擬似乱数ですが、そのうちの32ビット
しか使っていない状態です。

擬似乱数の初期化は、その擬似乱数の特性に合わせた乱数を充分に与えて初期
化する必要があります。

望ましい初期化プログラムの例 (C++)

std::mt19937 create_mt32()
{
    // メルセンヌ・ツイスタの内部状態分のサイズを確保した配列を用意する
    std::array<std::random_device::result_type,
               std::mt19937::state_size> a;

    // 暗号論的擬似乱数生成器
    std::random_device rd;

    // 配列に乱数を生成する
    std::generate(a.begin(), a.end(), std::ref(rd));

    // 乱数ライブラリ共通のシードオブジェクトに、
    // 生成した乱数をセットする
    std::seed_seq q(a.begin(), a.end());

    // メルセンヌ・ツイスタに、
    // シードオブジェクトを渡して初期化する
    return std::mt19937(q);
}

std::mt19937_64 create_mt64()
{
    // 64ビット版はstate_sizeが半分になる事に注意する
    std::array<std::random_device::result_type,
               std::mt19937_64::state_size * 2> a;

    std::random_device rd;
    std::generate(a.begin(), a.end(), std::ref(rd));
    std::seed_seq q(a.begin(), a.end());
    return std::mt19937_64(q);
}

C++標準のstd::random_deviceは2024時点ではほとんどの場合「暗号論的擬
似乱数生成器」になっています。ですが過去には、実際の実装はただの擬似
乱数であった。という問題がありました。使用しているコンパイラの実装が
どうなっているのか事前に確認しておくのがよいでしょう。

乱数の確率分布

一様分布

全ての値の出現率が等しい分布です。最もよく使われる分布です。

擬似乱数そのままの値は一様分布になっています。しかし、任意の範囲の値を
取り出す場合には偏りが出ます。「乱数を使うには注意がいる」で説明したよ
うに、余りが出た分の乱数値を弾く等をして、一様分布になるようにする必要
があります。

C++ではstd::uniform_int_distributionや
std::uniform_real_distributionというクラスを使うことで、任意の範囲の
値を一様分布にできます。

一様分布のグラフ（0～10の値）

    +---------------------------------------------------+
    |    +    +     +    +    +    +    +     +    +    |
    |                                                   |
    |                                                   |
    |                                                   |
    |                                                   |
    |                                                   |
    |                                                   |
 9% |***************************************************|
    |                                                   |
    |                                                   |
    |                                                   |
    |                                                   |
    |                                                   |
    |                                                   |
    |    +    +     +    +    +    +    +     +    +    |
    +---------------------------------------------------+
    0    1    2     3    4    5    6    7     8    9    10

重み付け分布

各要素毎に重み(Weight)という値を付けて、意図的に出現率に偏りを持たせる
方法です。分かりやすく、細かく設定でき、実装も比較的簡単であるため、ゲー
ムではよく使われる印象があります。

欠点としては。要素数が多くなってくると、管理が煩雑になりミスも起こりや
すくなります。また、設定によっては、プレイヤーからは予測し辛い動作をす
るために、不評を買うこともあります。

以下のような重み付けの表を作り、データーを設定します。

値	重み		出現率
0	3	3 / 22 =	13.636%
1	3	3 / 22 =	13.636%
2	3	3 / 22 =	13.636%
3	2	2 / 22 =	9.091%
4	2	2 / 22 =	9.091%
5	2	2 / 22 =	9.091%
6	2	2 / 22 =	9.091%
7	2	2 / 22 =	9.091%
8	1	1 / 22 =	4.545%
9	1	1 / 22 =	4.545%
10	1	1 / 22 =	4.545%
合計	22

C++ではstd::discrete_distributionというクラスで、Pythonでは
random.choices()という関数で、重み付けを持ったランダム選択をする事が
できます。

重み付け分布のグラフ

 14% +--------------------------------------------------+
     |***************    +     +    +    +    +    +    |
 13% |-+            *                                 +-|
 12% |-+            *                                 +-|
     |              *                                   |
 11% |-+            *                                 +-|
 10% |-+            *                                 +-|
     |              *                                   |
  9% |-+            ***************************       +-|
     |                                        *         |
  8% |-+                                      *       +-|
  7% |-+                                      *       +-|
     |                                        *         |
  6% |-+                                      *       +-|
  5% |-+                                      *       +-|
     |    +    +    +    +     +    +    +    **********|
  4% +--------------------------------------------------+
     0    1    2    3    4     5    6    7    8    9    10

正規分布

正規分布というのは、自然界の事象で現れる確率分布を表す数式です。自然界
の事象というと、野菜の重さや、テストの成績（偏差値、うっ…頭が…）等、
平均に近い値が最も現れやすい事象です。人間の感覚には馴染みやすい確率分
布と言えます。

RPGのようなゲームでは使用例はあまり見かけませんが、ガンシューティング
のようなゲームでは、銃の集弾性を再現するために使われる事もあるかと思い
ます。

正規分布の乱数では主に二つのパラメーターを使います。

• 平均 - 出現する値の平均値。
• 標準偏差 - 出現する値の幅。平均±標準偏差の値が、出現率68.27%になりま
  す。

C++であればstd::normal_distributionというクラスを使うことで、正規分
布乱数を得られます。

正規分布乱数の値の幅は理論的には無限の幅があるので、確率は低いですが、
大きな数値が出てくることがあります。これでは扱いにくいので、想定する最
小値や最大値の範囲を超えた場合には、範囲内の値が出るまで正規分布乱数を
振り直す必要があります。

この振り直しをすると、「正規分布」のグラフで両端が切り取られた状態にな
ります。この場合は「切断正規分布」と呼びます。「正規分布」とは計算が異
なってくるため、表計算ソフト等で試算する場合には注意が必要です。

設定においては、標準偏差を広げすぎたり、極端に最小値・最大値を狭くする
と乱数の振り直し回数が増え、意図せず負荷を増やす事になります。その場合
には「一様分布」や「重み付け分布」等を検討してください。無段階の乱数値
が必要でなく、決まった整数値の範囲であるならば「正規分布」を「重み付け
分布」に変換する事もできます。

正規分布乱数のアルゴリズムはいくつかありますが「ボックス＝ミュラー法」
の場合では、2つの乱数を必要とし、内部的に乱数を2回呼び出す事があります。

また、出力幅が32ビットの擬似乱数を使用し、倍精度浮動小数点（53ビット）
を得る場合には、さらに乱数を2回呼び出します。つまり、最大で4回の乱数を
呼び出すことになります。ですので、k次元均等分布の高い擬似乱数を選ぶの
が良いでしょう。

次のような例で、切断正規分布を試算する方法を説明します。

平均	5
標準偏差	2
最小値	0
最大値	10
数値丸め	四捨五入して整数とする

まずチェックすべきは「数値丸め」があるかです。正規分布乱数は整数ではな
く、浮動小数点数の連続性のある値です。丸めが四捨五入である場合は、1の
値が出てくる範囲は0.5～1.499...の範囲になるので、その範囲の出現率を求
める必要があります。計算上は目的の値に±0.5して最小値、最大値とすれば
良いです。

出現率は表計算ソフトではNORM.DIST()という正規分布関数を使って求めら
れます。

1の値の確率は正規分布関数では以下のように求めますが、1という値の1点に
絞られた確率なので、今回の場合は使えません。

=NORM.DIST(1, 5, 2, FALSE)

1は0.5～1.499...の範囲であるので、範囲の確率を求めなければなりません。
0.5～1.499...の範囲の確率は、1.499...までの累積確率から、0.5までの累積
確率を引く事で求められます。計算上は1.5として構いません。

=NORM.DIST(1.5, 5, 2, TRUE) - NORM.DIST(0.5, 5, 2, TRUE)

• 正規分布関数の引数TRUEは累積確率を求める。という意味です。

0.5～1.499...の範囲の「正規分布」の確率が求められましたが、まだ充分で
はありません。ここから「切断正規分布」の確率を求める必要があります。範
囲外の値が出た場合には、乱数の振り直しを行うので、そのぶん全体の確率が
上がるのです。

「切断正規分布」求めるにはまず、範囲内の値が出る確率を求めます。最小値
は0、最大値は10で、四捨五入があるので、-0.5～10.499..になります。計算
上は10.5として構いません。

=NORM.DIST(10.5, 5, 2, TRUE) - NORM.DIST(-0.5, 5, 2, TRUE)

範囲内の値が出る確率で、正規分布の確率を割れば、切断正規分布の確率にな
ります。四捨五入された1の値が0～10の範囲で出る確率は以下のように求めら
れます。

=(NORM.DIST(1.5, 5, 2, TRUE) - NORM.DIST(0.5, 5, 2, TRUE))
 / (NORM.DIST(10.5, 5, 2,TRUE) - NORM.DIST(-0.5, 5, 2, TRUE))

これらを踏まえて計算すると以下のような表になります。

値	最小値	最大値	平均	標準偏差	出現率
0	-0.5	0.5	5	2	0.930%
1	0.5	1.5	5	2	2.800%
2	1.5	2.5	5	2	6.598%
3	2.5	3.5	5	2	12.170%
4	3.5	4.5	5	2	17.571%
5	4.5	5.5	5	2	19.860%
6	5.5	6.5	5	2	17.571%
7	6.5	7.5	5	2	12.170%
8	7.5	8.5	5	2	6.598%
9	8.5	9.5	5	2	2.800%
10	9.5	10.5	5	2	0.930%
範囲内	-0.5	10.5	5	2	99.404%
範囲外					0.596%

範囲外の値が出る確率が、振り直し発生の確率に影響します。意図しない負荷
を避けるため、この確率が上がり過ぎないように注意する必要があります。

色々と注意点があり、試算と実装されたプログラムの結果を合わせるのが難し
いですが、一度システム化できればしっかりと出現率を管理することができま
す。

切断正規分布のグラフ

 20% +--------------------------------------------------+
     |    +    +    +    + *** + ** +    +    +    +    |
 18% |-+                 **        **                 +-|
 16% |-+                *            *                +-|
     |                **              **                |
 14% |-+             *                  *             +-|
 12% |-+            *                    *            +-|
     |             *                      *             |
 10% |-+          *                        *          +-|
     |           *                          *           |
  8% |-+        *                            *        +-|
  6% |-+       *                              *       +-|
     |       **                                **       |
  4% |-+   **                                    **   +-|
  2% |-+***                                        ***+-|
     |**  +    +    +    +     +    +    +    +    +  **|
  0% +--------------------------------------------------+
     0    1    2    3    4     5    6    7    8    9    10

切断正規分布乱数のプログラム例 (C++)

#include <algorithm>
#include <array>
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <limits>
#include <random>
#include <vector>

static auto create_rng()
{
    std::array<std::random_device::result_type,
               std::mt19937_64::state_size * 2> a;
    std::random_device rd;
    std::generate(a.begin(), a.end(), std::ref(rd));
    std::seed_seq q(a.begin(), a.end());
    return std::mt19937_64(q);
}

class trunc_normal_distribution
{
    // distは乱数を内部にキャッシュする事がある
    std::normal_distribution<double> dist;
    int min;
    int max;

public:
    trunc_normal_distribution(double mean, double stddev, int min, int max)
        : dist(mean, stddev), min(min), max(max) {}

    template <class RNG>
    int operator()(RNG& rng) {
        while (true) {
            double n = std::round(dist(rng));
            if (min <= n && n <= max)
                return n;
        }
    }
};

int main()
{
    constexpr double mean = 5;   // 平均
    constexpr double stddev = 2; // 標準偏差
    constexpr int min = 0;
    constexpr int max = 10;

    auto rng = create_rng();
    trunc_normal_distribution dist(mean, stddev, min, max);

    constexpr int range = max - min + 1;
    std::vector<int> v(range);

    // 出現回数を集計する
    constexpr int times = 0x100'0000;
    for (int i = 0; i < times; i++) {
        int n = dist(rng);

        size_t index = n - min;
        v[index]++;
    }

    // 出現率を表示する
    for (int i = 0; i < range; i++) {
        printf("%2d %6.3f%%\n", i + min, 100.0 * v[i] / times);
    }
}

サイコロの和

2つのサイコロを振って値の和を取ると三角形状の確率分布になります。

古くはテーブルトークRPGで使用され、6面のサイコロを2回振った合計を2d6、
3回振る場合は3d6、といった表記で馴染みがあるかもしれません。

コンピューターゲームでも出現率の統計をとっていると、このような分布を見
ることがあります。簡潔で計算量が少ないため、正規分布の近似としてゲーム
に使われることがあるようです。

0～10の乱数の場合は、0～5の乱数2つを足して乱数とします。他には0～10の
乱数を複数回足して、平均を取る。といった方法もあるようです。足し算をす
る回数を増やしていくと、分布のグラフが丸みを帯びた形になっていきます。

「乱数は二度振るな」でも説明しましたが、足し算の回数を増やしすぎると組
み合わせ問題にも繋がるので注意が必要です。回数を増やすなら「正規分布」
か「重み付け分布」を使うことを検討してください。

また、値の範囲が大きくなったり、足し算の回数が増えると、組み合わせの計
算が難しくなり、出現率を試算し辛くなります。そのため調整も難しくなりま
す。

乱数の和のグラフ（(0～5) + (0～5) = (0～10)）

 18% +--------------------------------------------------+
     |    +    +    +    +    **    +    +    +    +    |
 16% |-+                    **  **                    +-|
     |                    **      **                    |
 14% |-+                 *          *                 +-|
     |                 **            **                 |
 12% |-+             **                **             +-|
     |              *                    *              |
 10% |-+          **                      **          +-|
     |          **                          **          |
  8% |-+      **                              **      +-|
     |      **                                  **      |
  6% |-+  **                                      **  +-|
     |  **                                          **  |
  4% |**                                              **|
     |    +    +    +    +     +    +    +    +    +    |
  2% +--------------------------------------------------+
     0    1    2    3    4     5    6    7    8    9    10

様々な擬似乱数アルゴリズム

線形合同法

名前	線形合同法
英語名	Linear congruential generator
登場	1951,1958年頃

古くからあり、最も長く使われてきた擬似乱数と言えます。2024年時点では、
より高性能な擬似乱数が登場し、使われる機会は減ってきています。

以下のような計算式で表されます。

x = (x * a + c) mod m

条件に合ったa,c,mのパラメーターを選ぶ事で周期が最大となり、擬似乱数と
して成立します。このパラメーターの選び方によって大きく性能が変わってき
ます。条件については英語版のWikipediaに詳しく書かれています。

32ビットの線形合同法で良いとされたのは、1993年にParkとMillerによって提
案された以下のパラメーターです。

x = (x * 48,271 + 0) mod (2^31 - 1)

cが0であるため、内部状態が0になりません。そのため乱数値として0が出現し
なくなるため、注意が必要です。また、mが素数となっているため剰余の計算
が遅くなり、計算の過程も64ビットである必要があります。1993年頃のコン
ピューターにとっては少々計算量の多い乱数であったと言えます。

この他にも、64ビット版と128ビット版が存在します。

1993年頃のC言語のライブラリでは以下のような実装のものがありました。理
由はよく分かりませんが、このパラメーターを使った線形合同法が流行してい
たのです。

x = (x * 1,103,515,245 + 12,345) mod 2^31

この乱数は品質が良くありませんでした。特に下位16ビットに規則性があり。
さらに奇数・偶数が必ず交互に現れる。という特性がありました。当時はC言
語のrand()を使う場合には気を付けろ。と、よく言われていました。

C言語の擬似乱数は色々と問題が多かったため、その後、下位16ビットを捨て
て上位15ビットを使ったり、サイズを拡張して48ビットや64ビットにしたりと、
C言語ライブラリにも少しずつ改良が加えられていき、1999年頃には落ち着き
を見せます。

具体的な例では、2008年。Final Fantasy XI(FF11)というゲームにて、モグボ
ナンザという宝くじのようなイベントがありました。5桁のくじを買って、当
たれば景品が貰えるというものです。くじの番号をランダムで買う、という機
能があり、PS2版のFF11では番号の桁毎に奇数・偶数が交互に出る。という事
が起こったのです。

これは全ての組み合わせのうち6%の程度の範囲しか生成できない状態です。当
選確率に大きな差がついてしまう事から問題になりました。

結果から見ると、以下のようなプログラムになっていたと考えられます。

1桁目 = rand() mod 10
2桁目 = rand() mod 10
3桁目 = rand() mod 10
4桁目 = rand() mod 10
5桁目 = rand() mod 10

推測ですが、PS2のC言語ライブラリが古かったか、独自の擬似乱数を使用して
いたのが原因ではないかと思われます。

他には、2006年にカルドセプトサーガというゲームでも似た事例があったよう
です。

線形帰還シフトレジスタ (LFSR)

名前	線形帰還シフトレジスタ
英語名	Linear-feedback shift register (LFSR)
登場	1960年頃

1ビットの乱数を出力する擬似乱数です。いつ誰が考案したのか、よく判って
いませんが、かなり古くから使われていたようです。あくまで1ビットの乱数
なので、数値として扱うのには向いていません。

古いアルゴリズムですが、デジタル通信を行うハードウェア等では現役で使わ
れています。デジタル通信の信号は0と1のみを扱いますが、0や1が長い期間連
続で現れるようなデーターがあると、信号のズレが生じてうまく通信できなく
なります。

実際のデーターには0が多く出現します。通信データーのビット列に、LFSRで
出力されたビット列を、排他的論理和(XOR)で掛け合わせます。こうする事で0
と1が均等になり、通信の安定化を実現しています。

受信側でも同じシード値のLFSR使い、再度排他的論理和(XOR)を掛ける事で、
元のデータに戻す事ができます。

排他的論理和(XOR)というのは、コンピューターでよく使われる演算法で、1ビッ
ト単位で処理を行います。情報量を減らす事なくビットを撹拌するのに向いて
いるので、様々なアルゴリズムで利用されます。

排他的論理和(XOR)演算の組み合わせは以下のようになります。

0 XOR 0 = 0
0 XOR 1 = 1
1 XOR 0 = 1
1 XOR 1 = 0

8ビットのLFSRの場合、以下のように演算します。

1. 8,6,5,4ビットを取り出しXORで結合する
2. 右に1ビットシフトする（ビット全体をずらし、8ビット目は消滅する）
3. XORで結合したビットを、先頭に入れる

    | --------> Shift --------> |
+---+---+---+---+---+---+---+---+ - +
| * | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 :
+-^-+---+---+---+---+---+---+---+ - +
  |               |   |   |       | 
  |               |   |   +--XOR--+
  |               |   |       |
  |               |   +--XOR--+
  |               |       |
  |               +--XOR--+
  |                   |
  +-------------------+

この取り出すビットの事を「タップ」と呼びます。このタップ位置によって最
大周期になるかが決まります。8ビットの場合はこの例のように、8,6,5,4が最
大周期になります。

最大周期になるタップ位置を求めるのはとても計算が難しいのですが、ハード
ウェアメーカーの資料や、LFSRに関する論文等から、ビット数別のタップ位置
を知ることができます。WikipediaのLFSRページからも探せます。

LFSRの特徴として内部状態が0にならないというのがあります。これはXOR演算
のみを使っている影響です。シード値を0として初期化すると永久に0しか出力
しない状態になるので、シード値には0以外を使う必要があります。

また、内部状態が0にならないという事は、連続して0が出力される回数は、8
ビットであれば7回まで、32ビットであれば31回まで、という事になります。

この例では「フィボナッチLFSR」というアルゴリズムを説明してきました。も
う一つ「ガロアLFSR」という派生アルゴリズムがあります。これはコンピュー
ターで高速に処理できるように改良されたものです。同じビット列を出力でき
ますが、シード値には互換性がありません。

ドルアーガの塔のLFSR

「ドルアーガの塔」というゲームがあります。このゲームではマップが擬似乱
数によって生成されている事で有名です。毎回異なるマップになるという事で
はありません。フロア番号をシード値として擬似乱数を初期化し、その乱数を
使ってマップ要素を生成するようになっています。毎回同じマップになります
が、マップデーターがシード値の1バイトだけで済む。という訳です。

擬似乱数には独自の8ビットLFSRが使われています。タップ位置は8,5で、新し
いビットにはNOT演算が加えられているのが特徴です。

NOT演算は単にビットを反転する処理です。

NOT 0 = 1
NOT 1 = 0

ドルアーガの塔のLFSRは以下のように演算します。ビットの順序とビットシフ
トが逆になっています。

    | <-------- Shift <-------- |
+ - +---+---+---+---+---+---+---+---+
: 8 | 7 | 6 | 5 | 4 | 3 | 2 | 1 | * |
+ - +---+---+---+---+---+---+---+-^-+
  |           |                   |
  +----XOR----+                   |
        |                         |
        +-----------NOT-----------+

LFSRは内部状態が0になる事がないのですが、NOTを加えることで内部状態が0
になっても動作が可能になります。逆に内部状態の全ビットが1、という状態
にならなくなります。初期値が全ビット1の場合は、内部状態は1のまま変化し
なくなります。

ドルアーガの塔ではLFSRを1ビットの乱数としてではなく、内部状態全部を参
照して8ビットの乱数値として使用しています。そのため乱数として0が出ない
のは都合が悪かったのではないかと思われます。

8ビットでタップ位置が8,5なので、このLFSRは最大周期になっていません。シー
ド値によって周期が4つのグループに分かれます。

周期	値
217	0,1,2,3,4,5,7,8,9, ... 254
31	6,13,16,27,32,43, ... 244
7	26,52,70,104,141,163,209
1	255

255は全ビットが1となる内部状態であるため、255しか出力されません。フロ
ア60の壁が一直線になっているマップは、シード値を255にする事で実現して
います。

こうして見ると品質の良い擬似乱数とは言えないのですが、これによって特徴
のある階段状の壁が生成されています。壁の生成は、乱数を引いて4方向のう
ちどちらに進むかを決め、他の壁に当たるまでを繰り返してマップが形成され
ます。

LFSRは1ビットの乱数なので、内部状態をそのまま乱数値として取り出すと1ビッ
ト目が新しい乱数ビット。2ビット目以降は過去の履歴となります。0と1が出
現する確率は50%なので、最初の2ビットを見ると、01,10,01,10 ...というパ
ターンが多くなるわけです。これを4方向として解釈すると右・下・右・下と
いうパターンになり、これによって階段状の壁が出来やすくなるのです。

改良するならば、LFSRを2度回してから2ビットを取り出す。といった処理にす
ればかなり乱数の品質は改善され、階段状の壁はなくなると思います。ですが、
あの階段状のマップの方が味があって良い。という意見もあるようですね。

また、この乱数の規則性には別の効果があり、閉鎖路が形成されにくいという
効果もあります。閉鎖路というのはロの字形の壁ですが、ドルアーガの塔には
現れません。本来ならば閉鎖路を検出して壁を生成しない、といった処理が必
要なのですが、乱数の規則性により閉鎖路が生成されないため、閉鎖路検出を
せずとも正しくマップを生成でき、プログラムが簡潔になります。

ここまではアーケード版の話でした。ファミリーコンピュータ(FC)版では擬似
乱数に改良が加えられており、LFSRと線形合同法の2つの擬似乱数を持ってい
ます。2つの擬似乱数の結果をXORしたものを乱数として出力しているため、擬
似乱数の品質が向上しています。

擬似乱数の品質向上に伴い、閉鎖路も生成されやすくなってしまったため。FC
版では閉鎖路の検出も行われているようです。

LFSRの話が長くなってしまいましたが、これもドルアーガの塔の乱数を語りた
いが為、順を追って説明してきました。

ドルアーガの塔の乱数アルゴリズムについて詳しく知りたければ、拙著のサン
プルコードをご覧になってください。

https://github.com/funanz/todmap

GFSR

名前	GFSR
英語名	Generalized feedback shift register
登場	1973年頃

LFSRは1ビットの乱数であるため、数値として扱えませんでした。ならばLFSR
を32本束ねて32ビットの乱数にしよう。というのがGFSRです。長大な周期を持
ち主に科学分野で利用されたようです。

現存する実装では、GSL(GNU Scientific Library)のgfsr4というものがありま
す。16,384ビットのLFSRを32本持っており。周期が2^16,384-1。出力幅が32ビッ
ト。1次元の均等分布になっています。内部状態は巨大で、65,540バイトのメ
モリを必要とします。

非常に長い周期を持ちますが、メモリの消費量が多いため。高速に動作させる
には扱いが難しく、初期化に大量の乱数を必要とし、最大周期が出るように初
期化するには工夫が必要でした。このためあまり広くは利用されなかったよう
です。

GFSRの内部状態を、単純に暗号論的擬似乱数生成器で初期化した場合。運が悪
いと最大周期にならない事があります。最大周期にするには以下のようなビッ
トの並びを挿入して対策します。

8 7 6 5 4 3 2 1 bit
---------------
1 x x x x x x x
0 1 x x x x x x
0 0 1 x x x x x
0 0 0 1 x x x x
0 0 0 0 1 x x x
0 0 0 0 0 1 x x
0 0 0 0 0 0 1 x
0 0 0 0 0 0 0 1

• LFSR8本、出力幅8ビットの例。縦方向が1つのLFSR。xは任意のビットでよい。
  32ビットの場合は32x32になる。

かなり規則性のあるビット列なので、これらを近くにまとめて挿入すると乱数
の出力にも規則性が現れてしまいます。なので、これをどのように分散して挿
入するかという工夫も必要になってきます。

GSLのように16,384ビットもある大きなGFSRであれば、暗号論的擬似乱数生成
器で初期化した場合。最大周期にならない確率は低くなります。しかし、万が
一どれかの値が出ない。という事が起こっては困るのが擬似乱数です。許容す
るかどうかは判断に迷うところです。

また、ビット数の少ないGFSRの場合は、それぞれのLFSRの内部状態が0になら
ないように注意する必要があります。0のLFSRがあるとそのビットだけが0を出
力し続ける事になります。

メルセンヌ・ツイスタ

名前	メルセンヌ・ツイスタ
英語名	Mersenne Twister (MT)
登場	1998年

GFSRを改良して、少ないメモリと、初期化を容易にした「Twisted GFSR」とい
う擬似乱数が作られました。そして、それをさらに改良してできたのが「メル
センヌ・ツイスタ」です。

2^19,937-1という長大な周期と、623次元の均等分布、という高性能さ、日本
発の技術ということも相まって大変注目されました。そして、これによって乱
数の二度振り問題や、シャッフル問題が簡単に解決できるようになり、擬似乱
数の扱いがとても簡単になりました。

なぜここまでの高次元均等分布になるのか。正確ではないですが、簡単に説明
すると、19,937ビットの巨大な乱数を生成して、そこから32ビットずつ切り出
して乱数としているためです。1つの値をビット分割してx,yを作るのと同じ効
果があるわけです。

メモリの使用量は2,496バイトとGFSRに比べれば格段に少なくなっていますが、
高速に動作させるには、充分に小さいとはいえないため、2024年時点では利用
を控える傾向にあるようです。

とはいえ、623次元の均等分布という性能は唯一無二であるので、絶対に組み
合わせの欠落を出してはいけない。という用途には積極的に使うべきです。パ
フォーマンスの問題が出てから、他の擬似乱数アルゴリズムを検討しても遅く
はないと個人的には思います。

C++では標準ライブラリに含まれており、二種類の実装があります。

クラス	周期	k次元	出力幅
std::mt19937	2^19937-1	623	32ビット
std::mt19937_64	2^19937-1	311	64ビット

広い範囲の整数値や、浮動小数点数を多用する場合には64ビット版の方が高速
に動作します。しかし、k次元均等分布は半分になります。

Xorshift

名前	Xorshift
英語名	Xorshift
登場	2003年

LFSRに似た構造をしており、XORとビットシフトの組み合わせを工夫すること
で、内部状態をそのまま取り出して乱数として使えるようにした擬似乱数です。
非常にコンパクトで、高速に動作し、さらにdiehardテストという擬似乱数の
統計的テストに合格した事から、注目されるようになりました。

32ビットのXorshiftは以下のように演算します。

x = x XOR (x << 13)
x = x XOR (x >> 17)
x = x XOR (x << 5)

• 「<<」は左ビットシフト。「>>」は右ビットシフト。右辺の数はシフトする
  ビット数。

この例では、周期が2^32-1。出力幅が1～2^32-1。1次元に均等分布します。
LFSRと同様、0が出現しないので注意する必要があります。よくある実装では
下位1ビットを捨てる等をして対策をしています。

2024年時点では、改良された多くの派生アルゴリズムが考案されており、一部
の言語や実行環境等にも採用されています。

アルゴリズム	周期	k次元	出力幅	実装上の出力幅	実装
xorshift128+	2^128-1	1	64ビット	倍精度浮動小数点数 (53ビット)	Chrome, Edge, Firefox (JavaScriptエンジン)
xoshiro128**	2^128-1	4	32ビット	Int32 (31ビット)	.NET 6 (x86) (シード値未設定の場合のみ)
xoshiro256**	2^256-1	4	64ビット	Int32 (31ビット)	.NET 6 (x64) (シード値未設定の場合のみ)

k次元の均等分布はそれほど多くないので、多要素な物を乱数で生成する場合
には注意が必要です。

xoshiro256**はTestU01という擬似乱数の統計的テストにて、もっとも難しい
BigCrushというテストに合格しています。

TestU01はdiehardテストの後継として開発された擬似乱数テストです。
BigCrushはとても厳しいテストで、メルセンヌ・ツイスタでも合格できない程
のテストです。

PCG

名前	PCG
英語名	Permuted congruential generator
登場	2014年

線形合同法をベースとして、Xorshiftに似た手法を組み合わせた擬似乱数です。
TestU01という擬似乱数の統計的テストにて、もっとも難しいBigCrushという
テストに合格した事から注目されるようになりました。

PCG32の場合、以下のように演算します。

x = state
count = x >> 59
state = x * 6364136223846793005 + increment
x = x XOR (x >> 18)
r = rotr32(x >> 27, count)

• state 内部状態を保持する値。64ビット。
• increment 1ビット目は1固定で、上位63ビットが任意の値。
• rotr32() 32ビットで右ローテートする関数。右にビットシフトし、溢れた
  ビットを左に戻していくビット操作。
• count ローテートするビット数。
• r 出力する乱数値。32ビット。

他の擬似乱数と違う点として、複数の乱数系統（streamとも呼ばれる）を持っ
ているのが特徴です。incrementパラメーターがこれにあたり、値を変更する
ことで系統が変わります。

PCG32の場合、64ビットの線形合同法で内部状態を持っているので周期は2^64
になります。PCGのWebサイトでは、任意に周期を拡張できると説明されていま
す。

これは、周期の終わり毎に乱数系統を切り替えることで、周期を拡張できる。
という趣旨であると思われます。incrementパラメーターは上位63ビットが有
効に使えるので、合わせると127ビット相当の周期になります。

PCGは基本的に1次元の均等分布です。PCGのWebサイトでは、k次元均等分布は
任意に設定できると書いてあります。これは、別の系統の乱数を複数作成する
事で多次元の均等分布にも対応できる。という趣旨であると思われます。

周期とk次元均等分布の拡張については、具体的な使用例が見つからなかった
ので、実際に使うとなると実装は難しそうな印象です。

最後に

いかがでしたか？難しかったでしょうか。そう、擬似乱数は難しい。間違いあ
りません。そう思ってもらうのが趣旨です。

擬似乱数は、使うのも、検証するのも、とても難しいものです。32ビット程度
の擬似乱数であれば、全ての組み合わせを4秒でテストできます。しかし、64
ビットになると500年、128ビットでは不可能になります。

そういった出るかも分からないものを、理論的に出ると確信するためには、擬
似乱数の知識が必要なのです。

では、良い乱数を！

改訂履歴

https://gist.github.com/funanz/b226e0b2acd4a180212d3a0cfb0e6ba8/revisions

利用許諾

© 2024 funanz「擬似乱数とは何か」CC BY-SA 4.0