x_{n + 1} = (a x_{n} + c) mod m

計算の裏側

ゲームの乱数はどう作られているのか

How games make pseudo-random numbers

続けて出た 2 つの数を点に打つと、斜めの線が現れる。

計算機の乱数は決まった式で作る。線形合同法の点が線や平面に並ぶ理由、メルセンヌ・ツイスタと xorshift の工夫、正規分布への直し方、シードが重なる人数の見積もりまで。

種類: 計算の裏側
公式: 3 本
人物: 松本眞、西村拓士、ジョージ・マーサグリア、ジョージ・ボックス、マーヴィン・ミュラー

正規分布誕生日のパラドックスガウス積分

場面

37 を掛けて 1 を足す

電卓で、次の計算をくり返してみる。今の数に 37 を掛けて 1 を足し、256 で割った余りを次の数にする。0 から始めると、1、38、127、92、77、34、235、248、217、94 と続く。大きくなったかと思えば小さくなり、次の数の見当がつかない。256 回くり返すと 0 に戻り、それまでに 0 から 255 までの整数が 1 回ずつ出る。

ゲームのさいころも、宝箱の中身を決めるくじも、敵が右と左のどちらへ動くかも、計算機の中では多くの場合、こうした決まった式で作った数で決めている。式と始めの数が同じなら、何度やり直しても同じ数が同じ順に出る。それでも遊ぶ人には、でたらめに見える。

ところが、この並びには隠れた規則がある。続けて出た 2 つの数を組にして、(0, 1)、(1, 38)、(38, 127) と方眼紙に点を打っていく。256 個の点を打ち終えると、点は散らばらずに、9 本の平行な斜めの線の上に並ぶ。

決まった式から出る数は、どこまででたらめに見えるのか。線が現れる理由を式で確かめ、それを避ける xorshift とメルセンヌ・ツイスタの工夫、一様な乱数を正規分布に直す方法、シードの重なりを誕生日のパラドックスで見積もる方法までを追う。

考え方

決まった式が、でたらめに見えるまで

計算機の乱数は、内部に「状態」を持つ。1 回呼ぶたびに、決まった規則で状態を次へ進め、状態から数を 1 つ取り出す。いちばん初めの状態を決める数をシード(乱数の種)と呼ぶ。偶然から生まれた数ではないので、こうして作った数を疑似乱数という。

状態を表すビットの数には限りがあるので、状態は有限の通りしかない。いつか必ず前と同じ状態に戻り、そこから先は同じ並びをくり返す。くり返しの長さが周期で、状態が $b$ ビットなら、周期は $2^{b}$ を超えない。

冒頭の計算は線形合同法と呼ばれ、状態と出てくる数が同じ 1 つの整数である。37 を掛けると数は大きく増え、256 で割った余りをとると、あふれた分が捨てられる。 $37 \times 38 + 1 = 1407$ は、余りをとると 127 になる。掛け算で広げて余りで折り返すので、小さな違いが大きく散る。

周期の長さは係数で変わる。ハル=ドーベルの定理(1962 年)によれば、周期が法 $m$ に等しくなるのは、足す数 $c$ と $m$ が互いに素で、 $m$ のどの素因数も $a - 1$ を割り切り、 $m$ が 4 の倍数なら $a - 1$ も 4 の倍数のときに限られる。256 のように $m$ が 2 の累乗なら、 $c$ が奇数で、掛ける数 $a$ を 4 で割った余りが 1 という条件になる。 $a = 37$ は 4 で割って 1 余るので周期は 256、 $a = 3$ にすると周期は 128 に半減する。

よい疑似乱数には、1 つずつの数が一様に出るだけでなく、続けて出た数の組にも偏りが見えないことが求められる。メルセンヌ・ツイスタを作った松本眞と西村拓士は、これを次の形で表した。続けて出た $k$ 個の数を $k$ 次元の立方体の点とみなし、各辺を $2^{v}$ 等分して小部屋に分ける。1 周期のあいだに、どの小部屋にも同じ数の点が入るとき、 $v$ ビットの精度で $k$ 次元に均等に分布するという(すべてのビットが 0 の小部屋だけは 1 つ少ない)。冒頭の 256 個の 2 つ組では、正方形を縦横 2 等分した 4 つの小部屋でさえ、入る点の数が 59 から 69 とばらつく。9 本の線は、この均等さの崩れが目に見える形で表れたものである。

式

点が線と平面に並ぶ計算

線形合同法。 $a$ は乗数、 $c$ は増分、 $m$ は法で、 $x_{n}$ は 0 以上 $m$ 未満の整数である。0 から 1 の数にするときは $u_{n} = x_{n} / m$ とする。
$x_{n + 1} = (a x_{n} + c) mod m$
余りをとる操作を、ある整数 $k_{n}$ を引く形に書き直し、両辺を $m$ で割る。
$u_{n + 1} = a u_{n} + \frac{c}{m} - k_{n}$
$u_{n}$ は 0 以上 1 未満なので、 $k_{n}$ がとれる値は 0 から $a$ 程度までの整数に限られる。点 $(u_{n}, u_{n + 1})$ は、傾き $a$ の平行な線のどれかに乗る。 $a = 5$ なら 5 本である。
$u_{n + 1} = a u_{n} + \frac{c}{m} - k, k = 0, 1, 2, \dots$
ほかの向きの線もある。整数 $h_{1}, h_{2}$ が左の条件を満たすと、2 つ組はいつも右の関係を満たす。
$h_{1} + a h_{2} \equiv 0 (mod m) ⟹ h_{1} x_{n} + h_{2} x_{n + 1} \equiv h_{2} c (mod m)$
両辺を $m$ で割ると、点は向き $(h_{1}, h_{2})$ を法線とする平行な線に乗る。線の数はおよそ $∣ h_{1} ∣ + ∣ h_{2} ∣$ 、間隔は $1/ h_{1}^{2} + h_{2}^{2}$ である。 $a = 37$ ・ $m = 256$ では $h_{1} = - 3$ 、 $h_{2} = 7$ が条件を満たし( $- 3 + 37 \times 7 = 256$ )、点は間隔 0.131 の 9 本の線に乗る。
$h_{1} u_{n} + h_{2} u_{n + 1} = (整数) + \frac{h _{2} c}{m}$
RANDU は $c = 0$ 、 $m = 2^{31}$ で、乗数は $65539 = 2^{16} + 3$ である。乗数の 2 乗を $2^{31}$ で割った余りを求める。
$(2^{16} + 3)^{2} = 2^{32} + 6 \cdot 2^{16} + 9 \equiv 6 (2^{16} + 3) - 9 (mod 2^{31})$
2 つ先の数は $a^{2} x_{n}$ を $2^{31}$ で割った余りなので、続けて出た 3 つの数には、いつも次の関係がある。
$x_{n + 2} \equiv 6 x_{n + 1} - 9 x_{n} (mod 2^{31})$
$2^{31}$ で割ると、 $9 u_{n} - 6 u_{n + 1} + u_{n + 2}$ は整数になる。この値は $- 6$ より大きく $10$ より小さいので、とれる整数は $- 5$ から $9$ までの 15 個しかない。3 つ組の点は、法線 $(9, - 6, 1)$ の 15 枚の平行な平面に乗る。
$9 u_{n} - 6 u_{n + 1} + u_{n + 2} \in {- 5, - 4, \dots, 9}$

手順 4 で $h_{1}^{2} + h_{2}^{2}$ がいちばん小さい組を選ぶと、線の間隔がいちばん広い向きになる。3 つ組では、 $h_{1} + a h_{2} + a^{2} h_{3} \equiv 0 (mod m)$ を満たす $(h_{1}, h_{2}, h_{3})$ が平面の法線になる。下のシミュレーションは、この条件を満たす組を探して線と平面を引いている。

シミュレーション

続けて出た数を点に打つ

方式「線形合同法」では増分を $c = 1$ に固定し、乗数 $a$ と法 $m$ をつまみで変える。はじめの $a = 37$ ・ $m = 256$ では 9 本の線が現れ、 $a$ を 129 にすると傾き 1 の 3 本、5 にすると傾き 5 の 5 本になる。 $a = 37$ のまま $m$ を 4096 まで上げても、傾き 37 の線が 37 本並び、並びは消えない。見方を「3 つ組を立方体に」にすると立方体が回り、真横を向いたところで平面が縦の線に見える(横にドラッグすると自分で回せる)。RANDU では 15 枚になり、xorshift128+ とメルセンヌ・ツイスタでは、どちらの見方でも並びは見えない。

考え方

ずらして混ぜる

線形合同法では、状態が 1 つの整数で、2 つ組の点は全部で $m$ 個しかなく、それが格子の形に並ぶ。 $m$ を $2^{31}$ のように大きくすれば線の間隔は細かくなるが、RANDU のように 3 つ組で粗い平面が残ることもある。

ジョージ・マーサグリアが 2003 年に発表した xorshift は、状態のビットを左右にずらし(シフト)、もとの値との排他的論理和(XOR)をとる操作を重ねる。XOR は、ビットごとに繰り上がりを無視して足す計算で、整数の掛け算とは種類が違う。状態が全部 0 だと 0 のまま動かないので、それ以外の状態を 1 周でめぐるように、ずらす量を選ぶ。

Chrome の JavaScript エンジン V8 は、2015 年に Math.random をこの仲間の xorshift128+ に切り替えた。状態は 64 ビットの数 2 つである。1 回の更新では、片方の数に、それを左に 23 ビットずらした値を XOR し、次に右に 17 ビットずらした値を XOR する。さらに、もう片方の数と、それを右に 26 ビットずらした値を XOR する。名前の + は、更新した 2 つの数の和を出力にすることを表す。V8 のブログによれば周期は $2^{128} - 1$ で、乱数の検査をまとめた TestU01 の検査をすべて通る。ブログの追記にあるとおり、V8 は 7.1 から和をとるのをやめ、片方の数だけから 0 以上 1 未満の数を作っている。

メルセンヌ・ツイスタ(MT19937)は、32 ビットの数 624 個を状態に持つ。1 回ごとに、ある語の最上位の 1 ビットと次の語の下位 31 ビットをつなぐ。つないだ値を 1 ビット右にずらして 397 語先の語と XOR をとり、もとの語を書き換える(つないだ値の最下位が 1 なら、定数 9908B0DF も XOR する)。数を取り出すときには、右に 11、左に 7、左に 15、右に 18 ビットずらして XOR をとる 4 段の変換(テンパリング)をかけ、上位のビットで見た均等さを高める。状態のうち実際に効くのは $32 \times 624 - 31 = 19937$ ビットで、周期は $2^{19937} - 1$ 、6002 桁の数になる。この数は $2^{p} - 1$ の形の素数、つまりメルセンヌ素数である。周期をメルセンヌ素数にするために、配列の 1 語のうち 31 ビットは使わない。論文の題にある「623 次元」は、続けて出た 623 個までの数が 32 ビットの精度で均等に分布することを指す。

xorshift もメルセンヌ・ツイスタも、ビットの XOR とシフトだけでできた線形な方法である。メルセンヌ・ツイスタの作者は、ページのよくある質問で、線形な漸化式による方法は、十分な長さの出力を見ればその後の数列を完全に予言できるので、そのままでは暗号には使えないと注意している。この記事を書くときに試すと、MT19937 の続けて出た 624 個の出力からテンパリングを逆にたどって状態が求まり、その先の 1000 個をすべて当てられた。V8 のブログも、xorshift128+ は暗号用としては安全でないとしている。

式

0 から 1 の数を釣り鐘の形に

0 より大きく 1 以下の一様な乱数を、互いに独立に 2 つ作る。
$U_{1}, U_{2} \in (0, 1]$
半径と角を作る。
$R = - 2 ln U_{1}, Θ = 2 π U_{2}$
直交座標に戻すと、互いに独立に標準正規分布 $N (0, 1)$ に従う 2 つの数になる。これがボックス=ミュラー法である。
$Z_{1} = R cos Θ, Z_{2} = R sin Θ$
理由を確かめる。独立な 2 つの標準正規分布の同時の確率密度を、極座標 $z_{1} = r cos θ$ , $z_{2} = r sin θ$ に移す。面積の要素は $d z_{1} d z_{2} = r d r d θ$ になる。
$\frac{1}{2 π} e^{- \frac{z _{1}^{2} + z _{2}^{2}}{2}} d z_{1} d z_{2} = \frac{d θ}{2 π} \cdot e^{- \frac{r ^{2}}{2}} r d r$
右辺は $θ$ の部分と $r$ の部分の積なので、角は 0 から $2 π$ に一様で、半径とは独立である。半径が $r$ を超える確率は次のとおり。
$P (R > r) = \int_{r}^{\infty} e^{- \frac{s ^{2}}{2}} s d s = e^{- \frac{r ^{2}}{2}}$
$U = e^{- R^{2} /2}$ と置くと、 $P (U < t) = P (R > - 2 ln t) = t$ なので、 $U$ は 0 から 1 の一様分布に従う。逆に解くと手順 2 の $R$ になる。
$U = e^{- \frac{R ^{2}}{2}} ⟺ R = - 2 ln U$
平均 $μ$ 、標準偏差 $σ$ の正規分布にするには、標準化の式 $Z = (X - μ) / σ$ を逆に使う。
$X = μ + σ Z$
例として $U_{1} = 0.3$ , $U_{2} = 0.6$ なら、 $R = 1.552$ , $Θ = 3.770$ になる。もう 1 つの数は $Z_{2} = 1.552 sin 3.770 = - 0.912$ である。
$Z_{1} = 1.552 cos 3.770 = - 1.255$

手順 4 は、ガウス積分の証明で 2 乗して極座標に移るのと同じ計算である。プログラムでは $U_{1} = 0$ を避ける。JavaScript の Math.random は 0 以上 1 未満の数を返すので、1 − Math.random() を $U_{1}$ に使えば $ln 0$ にならない。

もっと素朴な方法もある。1970 年の IBM の手引きに載った GAUSS というサブルーチンは、一様な乱数を 12 個足して 6 を引いていた。一様分布の分散は $1/12$ なので、12 個の和の分散はちょうど 1 になり、中心極限定理で正規分布に近づく。ただし 12 個の和は 0 から 12 の間にしかないので、平均から 6 より外の値は決して出ない。本物の正規分布で平均から $6 σ$ より外に出るのは 10 億回に 2 回ほどで、ゲームではまず困らないが、まれな出来事の確率を調べる計算には向かない。

歴史

65539 から 2¹⁹⁹³⁷ − 1 まで

1958 年、ジョージ・ボックスとマーヴィン・ミュラーは、数理統計学の学術誌『The Annals of Mathematical Statistics』に、2 ページの論文「正規乱数の生成についての覚え書き」を載せた。上の手順 2 と 3 の方法は、2 人の名前で呼ばれている。

IBM は大型計算機 System/360 向けに、FORTRAN で書いた計算部品を集めた科学計算サブルーチンパッケージ(SSP)を配った。1970 年 8 月の第 5 版の手引きに、乱数を作る RANDU が載っている。本体は 8 行で、前の数に 65539 を掛け、結果が負なら $2^{31}$ を足し、 $2^{31}$ の逆数にあたる $0.4656613 \times 1 0^{- 9}$ を掛けて 0 から 1 の数にする。説明と注意書きには、最初のシードを 9 桁以下の奇数にすること、System/360 専用であること、 $2^{29}$ 個の数を出すとくり返すことが書かれている。この記事でシード 1 から数えると、536,870,912 回目、つまり $2^{29}$ 回目に 1 へ戻った。

1968 年 9 月、ジョージ・マーサグリアは米国科学アカデミー紀要(PNAS)に 4 ページの論文を載せた。題は「乱数はおもに平面に落ちる(Random numbers fall mainly in the planes)」である。

1997 年 10 月、慶應義塾大学の松本眞と西村拓士はメルセンヌ・ツイスタのページを開き、同じ月の 27 日に朝日新聞の夕刊が紹介した。論文は、1998 年 1 月の『ACM Transactions on Modeling and Computer Simulation』に載った。論文の謝辞によれば、前身の TGFSR という方法の着想は、東京大学の近くの KURUMAYA というパブでの米田信夫との議論から生まれた。米田は 1996 年 4 月に亡くなり、2 人は論文をその思い出に捧げた。2002 年 1 月には、シードの最上位のビットが状態にうまく伝わらない問題を直した初期化の手順を公開した。いまの C++ の mt19937 も、この初期化の係数 1812433253 を使う。

2003 年、マーサグリアは、統計ソフトウェアの学術誌『Journal of Statistical Software』に「Xorshift RNGs」を発表した。Unity の文書によれば、Unity の乱数はこの論文にもとづく Xorshift 128 である。

2015 年 12 月、V8 の開発者ヤン・グオはブログで、Math.random を MWC1616 から xorshift128+ に替えたと報告した。MWC1616 は、出力の上位半分がほとんど片方の状態だけで決まり、短い周期がたくさんあって、悪いシードでは周期が 4000 万に届かないこともあった。TestU01 の検査にも数多く落ちた。新しい実装は V8 4.9.41.0 に入って Chrome 49 から使われ、Firefox と Safari も xorshift128+ に切り替えた。

使いどころ

シードは何人で重なるか

エンジンと言語の乱数。 公式の文書とソースで確かめた範囲(2026 年 10 月)をまとめる。

環境	方式	文書にある注意
Unity の Random	Xorshift 128	起動時に OS から得たシードで 1 回だけ初期化する。Random.InitState でシードを決められる
Godot の RandomNumberGenerator	PCG32	方式は実装の詳細で、頼らないように。randfn はボックス=ミュラー法
Lua 5.4 の math.random	xoshiro256**	起動のたびに違う並びになるよう初期化する
Python の random	MT19937	安全のための用途には secrets モジュールを使う
JavaScript の Math.random(V8)	xorshift128+	暗号用ではない。crypto.getRandomValues を使う
C++ の標準ライブラリ	mt19937・minstd_rand0 など	既定のシードからの 10000 番目の値を規格が決めている

シードで同じステージを作る。 Unity は Random.InitState の文書で、手続き的に作るステージをシードの整数 1 つで保存すれば、記憶領域を節約できると勧めている。Godot も文書で、近いシードから似た並びが出ることがあるので、外から受け取ったシードはハッシュ関数に通すよう注意している。

重なりの見積もり。 InitState のシードは 32 ビットの整数なので、 $N = 2^{32}$ 通りある。ステージのシードを $n$ 個でたらめに選ぶと、同じシードが 2 つ以上出る確率は、誕生日のパラドックスの近似で次のようになる。

P (n) \approx 1 - exp (- \frac{n ( n - 1 )}{2 N})

半分を超えるのは $n \approx 1.18 N \approx 77, 000$ のときである。1 万個なら 1.2 %、10 万個なら 69 % で、100 万個なら同じシードの組が平均 116 組できる。シードを秒単位の時刻から作ると、 $N$ はずっと小さい。1 時間のあいだに 72 人がばらばらの時刻に遊び始めると、同じ秒に始めた 2 人がいる確率は 51 % になる。周期が $2^{19937} - 1$ のメルセンヌ・ツイスタでも、32 ビットのシードから始めるなら、出てくる並びは $2^{32}$ 通りしかない。作者も、使い道によっては初期値が $2^{32}$ 通りでは足りないとして、長い配列から初期化する手順を用意している。64 ビットのシードなら、半分を超えるのは約 51 億個からである。

出力の値も重なる。 シミュレーションのメルセンヌ・ツイスタは、32 ビットの出力の値が初めて重なった番目を読み出しに出す。シード 5489 では 101,995 個目、5490 では 79,829 個目、5491 では 14,609 個目だった。でたらめな 32 ビットの数でも、7.7 万個までに重なる確率は半分で、1.5 万個までに重なることも 40 回に 1 回ほどある。反対に、 $m = 2^{32}$ で周期が $m$ の線形合同法は、 $2^{32}$ 個を出し終えるまで同じ値を 1 度も出さない。重なりがなさすぎることも、でたらめでない手がかりになる。

さいころの目を作る。 0 以上 1 未満の数 $u$ に 6 を掛けて切り捨てれば、0 から 5 の目になる。上限の 1 を含む関数で同じことをすると、まれに目が 6 になる。メルセンヌ・ツイスタの作者も、閉区間の版では約 42 億回に 1 回この誤りが起きると注意している。Unity の Random.value と小数の Random.Range は上限を含むので、整数の目は上限を含まない整数の版の Random.Range(0, 6) で作る。目の数 $N$ が 2 の累乗でないときは、掛けて切り捨てる方法に $N \times 2^{- 32}$ 程度の偏りが出る。作者はよくある質問で、必要なビット数だけ取り出し、範囲を超えた値は捨てて引き直す方法を勧めている。

予測されては困る値。 対戦の組み合わせの抽選や認証のトークンのように、外から当てられると困る値は、上の 624 個の実験のように割り出されるおそれがある。Python なら secrets モジュール、JavaScript なら crypto.getRandomValues で作る。敵の動きや地形のように当てられても困らない場面では、速くて再現できる疑似乱数のほうが向く。

確かめ方。 生成器を自分で書いたら、既知の出力と照らす。C++ の規格は、既定のシードからの 10000 番目の値を、mt19937 なら 4123659995、minstd_rand0( $a = 16807$ ・ $m = 2^{31} - 1$ の線形合同法)なら 1043618065 と決めている。2003 年の標準化の提案文書によれば、実装の小さな誤りは出力を眺めても気づきにくいので、この値を入れた。メルセンヌ・ツイスタの作者は、プログラム mt19937ar.c と出力例 mt19937ar.out を配っている。このページのシミュレーションの生成器は、この 2 つの 10000 番目の値と、出力例の 1000 個に一致した。Node.js 24 を node --random-seed でシードを決めて動かすと、Math.random の出力 256 個(シード 3 通り)が、シミュレーションの xorshift128+ と一致した。V8 は 64 個ずつまとめて作って後ろから返すので、比べるときはその順に並べ直す。

次にやること。

使っているエンジンの乱数の方式と、シードのビット数を文書で確かめる。
遊び手やステージの数を $n$ として、シードが重なる確率を $1 - exp (- n (n - 1) /2 N)$ で見積もる。足りなければ 64 ビット以上のシードにするか、ハッシュに通してから使う。
ボックス=ミュラー法を 10 行ほどで書き、1 万個作って、平均から $\pm 1 σ$ に 68.3 % が入るかを正規分布の値と比べる。
シミュレーションで $m$ を 4096 にして、線の数がいちばん多くなる $a$ を探す。

数値計算の誤差が積もる例はオイラー法とルンゲ=クッタ法に、乱数で撃つ位置をずらして抜ける割合を数える試験は衝突判定とすり抜けにある。

参考文献

M. Matsumoto, T. Nishimura, "Mersenne Twister: A 623-dimensionally equidistributed uniform pseudo-random number generator", ACM Transactions on Modeling and Computer Simulation 8(1) (1998), 3-30. DOI 10.1145/272991.272995(作者のページの論文、mt19937ar.c・mt19937ar.out・よくある質問)
G. Marsaglia, "Random numbers fall mainly in the planes", Proceedings of the National Academy of Sciences 61(1) (1968), 25-28. DOI 10.1073/pnas.61.1.25
G. Marsaglia, "Xorshift RNGs", Journal of Statistical Software 8(14) (2003). DOI 10.18637/jss.v008.i14
G. E. P. Box, M. E. Muller, "A Note on the Generation of Random Normal Deviates", The Annals of Mathematical Statistics 29(2) (1958), 610-611. DOI 10.1214/aoms/1177706645
IBM, System/360 Scientific Subroutine Package Version III, Programmer's Manual, GH20-0205-4 (1970), Subroutine RANDU・Subroutine GAUSS
Y. Guo, "There's Math.random(), and then there's Math.random()", V8 blog (2015 年 12 月 17 日)と、V8 13.6 のソースコード src/base/utils/random-number-generator.h・src/numbers/math-random.cc・src/builtins/math.tq
ISO C++ 規格の作業草案 [rand.predef]、J. Maurer, "A Proposal to Add an Extensible Random Number Facility to the Standard Library (Revision 2)", WG21 N1452 (2003)
Unity Scripting API の Random・Random.InitState・Random.Range、Godot の RandomNumberGenerator、Lua 5.4 Reference Manual の math.random、Python の random モジュールの文書