あなたの使っている乱数、大丈夫？
-危ない標準乱数と、メルセンヌ・ツイスター開発秘話-

松本眞

広島大学理学研究科数学専攻

2014/11/18, 第 50 回市村学術賞記念 先端技術講演会

m-mat “at mark と呼ばれるもの” math.sci.hiroshima-u.ac.jp1昔から、乱数発生は（ばくちに）必要だ。

• サイコロ：{1, 2, 3, 4, 5, 6} の中からでたらめに数を発生させる装置• 「一様かつ独立に次々とでたらめな数を発生させる方法」
を乱数発生法という。

• コイン投げ、は表と裏で {0, 1} の要素を発生
ルーレットは {0, 1, . . . , 36} の要素を発生• くじ引きも、乱数発生法と言える。

2

Q. シミュレーションに乱数はなぜ必要？
A. 世界は、ばくちだから。
確率的現象のシミュレーションは、すごろくそのもの。
例１：連鎖核分裂• 一つ一つの原子核が、一定の確率で分裂• 中性子が放出される（放出される方向も確率的）• 中性子がぶつかった核は、また分裂
動画は atomicarchive.com, (c) Copyright 1998-2013 AJ Software & Multimedia All Rights Reserved

3一つ一つの確率現象をシミュレートするごとに乱数が消費される。
核分裂 ：たとえば、各原子につきサイコロを 10 回振って、全部１の目なら分裂、それ以外は非分裂。
中性子の方向 ：中性子の飛ぶ方向・速度もサイコロにより決める。それにより別の原子核にぶつかるか決まる。
例２：人生ゲーム：ルーレットを回して、就職や結婚などのイベントが確率的に決まっていく人生シミュレータ。

4• モンテカルロ法
確率的現象を乱数を用いてシミュレーションすることをモンテカルロ法と呼ぶ。
（モンテカルロ市はモナコの首都、カジノのメッカ）
物理、化学：核・化学反応シミュレーション
生物科学：たんぱく質折りたたみシミュレーション
金融工学：株価変動、金融商品の価格付
遊び、芸術：ゲーム、漫画の模様（トーン）、コンピュータグラフィックス

5乱数を、どうやって生成するか…それが問題だ。
“決定的な動作しかしない計算機では、乱数は生成できない”物理乱数発生器:• 物理雑音から拾ってくる: もっとも素朴な方法• 例：サイコロ、熱雑音、株価の変動• 問題点: コスト、スピード、再現不能性.

6

再現不能性とは、同じ乱数列を再現するのに、それらをすべて記録しておかないとならないこと。
再現性が必要となる場合：

• 追試
• 最適化
参考：核シミュレーションでは乱数は何兆個も消費される⇒ それらをすべて記録しておくのは効率が悪い。

7

擬似乱数発生法:
漸化式を用いて、乱数のように見える数列を生成する方法
例: 線形合同法 Linear Congruential Generator(LCG, Lehmer ’60)

• ある整数 x1 を初期シードとして選ぶ。
• 次の例のような漸化式により x2, x3, . . . を次々に生成：
xn+1 = 1103515245xn + 12345 mod 232.
例 x1 = 3 ならば
3 × 1103515245 + 12345 =3310558080 (mod232) → 3310558080 = x23310558080 × 1103515245 + 12345 = 3653251310737941945 (mod232) → 465823161 = x3465823161 × 1103515245 + 12345 = 514042959637601790 (mod232) → 679304702 = x4679304702 × 1103515245 + 12345 = 749623094657194335 (mod232) → 2692258143 = x58

擬似乱数のメリット:• 漸化式と初期シードを記録しておけば、誰でも同じ数列を再現できる• 高速で低コスト
問題点：「乱数と呼んでいいのか」...
擬似乱数の創始者 von Neumann「漸化式で乱数を作るのはある種の罪」“Anyone who considers arithmetical methods ofproducing random digits is, of course, in a state of sin”

9

たとえば:• 先の線形合同法は、70 年代から 80 年代にかけて
ANSI-C などの標準擬似乱数 rand() であった。
いまでも教科書にのっていて、広く使われている。

• この数列の周期は、初期シードの選び方によらず 232。
• 現代のパソコンは数分で 232 個の乱数を使ってしまう• 生成される数列はかなり乱数っぽく見えるが、数千万個の出力を使うと、非乱数性が現れてくる

10この線形合同法 rand による
232 個の 3 次元ランダム点プロット（70 倍拡大図）
62 = 8.5 × 10−3.
赤領域体積 ≃ 7253
45 = 6.2 × 10−3.
青領域体積 ≃ 7253
(真の値: 8.333 · · · × 10−3)11メルセンヌ・ツイスター擬似乱数(MT, 松本-西村 ’98) による
3 次元ランダム点プロット
61 = 8.4 × 10−3.
赤領域 ≃ 7248
64 = 8.8 × 10−3.
青領域 ≃ 7248
(真の値: 8.333 · · · × 10−3)12擬似乱数の落とし穴
落とし穴その１：偶数と奇数が交互に出るものがある
先の生成法 rand()3 × 1103515245 + 12345 =3310558080 (mod232) → 3310558080 = x2
3310558080 × 1103515245 + 12345 = 3653251310737941945 (mod232) → 465823161 = x3
465823161 × 1103515245 + 12345 = 514042959637601790 (mod232) → 679304702 = x4
679304702 × 1103515245 + 12345 = 749623094657194335 (mod232) → 2692258143 = x5
2692258143 × 1103515245 + 12345 = 2970947904275902380 (mod232) → 3498995628 = x63498995628 × 1103515245 + 12345 = 3861195017686361205 (mod232) → 1035215989 = x71035215989 × 1103515245 + 12345 = 1142376625729264650 (mod232) → 273505290 = x8273505290 × 1103515245 + 12345 = 301817257103158395 (mod232) → 1018653819 = x9

• 奇数と偶数が交互に出る。
• この現象は、線形合同法を「 mod 偶数」で使った場合必ず起きる。

• そして、知っていなければ気づきにくい。
13(slashdot.jp より引用：)• 2006 年 12 月 06 日
「カルドセプトサーガ」にダイス目が偶数と奇数を繰り返すバグ
株式会社バンダイナムコゲームス発売の Xbox360 向けゲー
ムソフト「カルドセプトサーガ」に「次のダイス目が偶数か奇数か推測できる」という致命的バグが見つかりました。

• 2006 年 12 月 16 日
ダイス目が偶数と奇数を繰り返すバグが取り上げられた
Xbox360 向けゲーム「カルドセプトサーガ」ですが、CNETの記事によると店頭在庫の回収を行うことが決定したようです。
14

モンテカルロ法用擬似乱数への要請• 高速性
素粒子シミュレーションなどでは、全体の計算時間の４０％を擬似乱数生成が占めることもある(MT に変えたら全体が 17%高速化したという報告あり)• 乱数性
擬似乱数は、所詮漸化式による数列「乱数性」の実用的な定義がない⇒「乱数である」ことを証明できない
仕方がないので代わりに：– 周期, 分布
周期や高次元分布に着目して「良いもの」を使う
15線形合同法の問題点：

xn+1 = axn + c mod M
なる数列の周期は高々M（∵ xn が取りうる値は 0, 1, . . . , M − 1 の M 個しかない）
メリット：周期や分布を求めるアルゴリズムがある
限界：周期を大きくするには M を大きくするしかなく、掛け算や割り算が必要で生成が遅くなる

16擬似乱数の落とし穴その２：周期が短い
rand(), TurboC 1.5, Visual C++, Borland C++, エクセルのVisual Basic の標準乱数は線形合同法（2000 年ごろの和田維作氏による調査、現在は改善か？）
周期は：232 = 4294967296 以下、エクセルは 224 以下。
→プログラム開発初期段階ではうまく行くが、本格的に大量の乱数を使うとおかしくなるというやっかいな状態

17例：４次元球の体積のモンテカルロ積分誤差
左：rand(周期 231) 右：メルセンヌ・ツイスター (周期 > 106000)
提供:日本原子力研究開発機構 核融合研究開発部門 清水勝宏氏・慶応大学 理工学部 物理学科早川奈伊紀氏 (2014/8/27 松本宛電子メール「長周期乱数コードに感謝」より)
cpu

err

cpuerr

10-1

104

10-1

104

4th-Dimension

4th-Dimension

randomrandMT

10-2

103

10-2

103

102
limit

10-4

10-5
104

101

105

106

107

108

109

100
1010

10-3

102

10-4

101

10-5
104

Nsamp

105

106

107
Nsamp

18

108

109

100
1010

CPU (s)10-3

relative Err

1/sqrt(N)CPU (s)relative Err

1/sqrt(N)メルセンヌ・ツイスター法 (松本-西村拓士 ’98):• 周期：219937 − 1 ≒ 4.3 × 106001• 1 周期で 623 次元空間に均等分布することが証明（32 ビット精度で）• 生成速度は、近年の線形合同法（ mod 248）よりも高速• 多くの計算機言語で標準擬似乱数として採用 (Python, Ruby,R, PHP, MATLAB, C++（C++11 から）など)
他、広く用いられている（多くのソフト、ポケモンゲーム、任天堂 Wii など）• MT の Wikipedia も見てください19速度比較 32 ビット整数一個を生成するのにかかるCPU サイクル数
cycle

cycles per generation

random: ラグ付き
フィボナッチ周期 ∼ 263
SFMT SIMD Fast MT
MT: Mersenne Twister
周期 219937 − 1rand: LCG 周期 231
rand48: LCG 周期 248
MRG: L’Ecuyer
周期 ∼ 2186

604020random

SFMT

MT

rand

rand4820MRG

メルセンヌ・ツイスターの動作原理:1 + 1 = 0 の数学：
二元体 F2F2 := {0, 1} とおく。

0,1 の掛け算は普通に定義して、F2 からはみ出ない。
足し算は 1 + 1 = 2 だけが F2 からはみ出してしまうので、1+1=0と定義する（2 で割ったあまりを見ている）。

21

F2 の世界で漸化式
xn+3 = xn+2 + xn (n ≥ 1),

x1 = 0, x2 = 0, x3 = 1

で x1x2x3x4 · · · を計算すると周期 7 で循環

123456789
0011
00111
0011100011101
00111010
001110100
0011101001
00111010011
001110100111
22

同様の漸化式

xn+607 = xn+273 + xnで数列を作る（初期値として x1, x2, . . . , x607 に 607 個の 0 ま
たは 1 を与える）と周期が 2607 − 1 となることが示せる。
注：この数は約 5.3 × 10182 である。
１ペタ Flops の計算機は、１秒に 1015 回の計算しかできない。
計算機で一周期求めることはできない。
宇宙の素粒子の数でも 1080 くらいである。
実験ではなく、数学的証明によってのみ周期を示しえる。

23このような F2 上の高階線形漸化式を用いて 0-1 列を生成、擬似乱数として用いる方法を
Tausworthe 法あるいはLinear Feedbacked Shift Register 法（LFSR 法）と言う (Tausworthe, 1965)。
メリット：• 周期を長くしても、生成速度が遅くならない「二か所見て、足して、一か所に書く」という動作だから

xn+607 = xn+273 + xn
• 周期が極大なので、「全て 0」以外のビットパターンが(xn, . . . , xn+607−1) (n = 1, 2, . . . , 2607 − 1)のなかにちょうど一度だけ現れる（607 次元均等分布性）
24

ベクトル化GFSR (Lewis-Payne ’73) 計算機ワード長の F2 ベクトル列を⃗xn+p := ⃗xn+q + ⃗xnで生成 (+ は F2 ベクトルとしての和)例 4 ビット計算機なら 4 ビット同時に作る

0101 ⃗x10110⃗x21111 ⃗x31010 ⃗x4 ⃗x41100⃗x50011
1001
0101
25

整定数 p, q をうまく選ぶと周期 2p − 1 にできる

• p として 31〜607 が良く使われていた• 各桁は Tausworthe 法で生成される数列に一致• 高速だが、各桁の間に情報のやり取りがない• 乱数性に問題あり（特にランダムウォークで）思いつき=⇒ 桁の間に情報の攪拌 (Twist) を行ったら？

結果:非常に効果があった。

26

Twisted GFSR (松本-栗田良春 ’92, ’94):
Twister と呼ぶ F2 係数正方行列 A を導入する：⃗xn+p = ⃗xn+q + ⃗xnA.
• A は桁の間の情報を混ぜる⇒ より長周期: 2wp − 1 が達成可能 (w:ワード長)• A は次のようなものを選ぶ:定数ベクトル ⃗a により{shiftright(⃗x)(⃗x の最下位ビットが 0 の場合)
⃗xA =shiftright(⃗x) + ⃗a (⃗x の最下位ビットが 1 の場合)
(0 1 0 0)⇒ 高速に計算可能, 十分一般: A = 00 00 10 01 .
a3 a2 a1 a0• 高次元均等分布性を改善するため ⃗xnT を出力 (’94)
27

開発秘話その 1：’87 学部 4 年生のときに、当時指導教員の米田信夫教授からGFSR 法とその問題点を聞いた。
GFSR 法を統計検定しているつくば計量研究所の栗田良春氏を訪問、いろいろ改良をこころみた。
（当時のスーパーコンピュータ・クレイが使い放題。
段ボール箱ふた箱に検定結果をプリントアウトした。）
根津の飲み屋「車屋」で米田研で飲んでいるとき、
松本「GFSR の 1 ワードを、rotate してはどうでしょう」
米田「それは僕もやった。けど、あまり良くなかった。」
松本「では、companion 行列を掛けたらどうでしょう。」
米田「... それはいい考えだ。早速、また栗田さんのところに行きなさい。」
28

開発秘話その２：Twisted GFSR に関する投稿論文は、投稿した学術雑誌が紛失した。
投稿後 2 年たっても「査読中」と言われた。3 年目に、雑誌の編集長が変わったときに「紛失した」と言われた。
そのうえ、reject された。ので、ぐれてほったらかした。

’90 に擬似乱数の大家 Pierre L’Eecuyer 教授が来日したとき、直接 Twsited GFSR を説明する機会があり、日の目を見ることになった。
(この「他力」がなければ、MT はなかったかもしれない。)

29

Mersenne Twister: (松本-西村拓士 ’98):
TGFSR の周期を 2wp−r − 1 の形にするため、状態空間 wp ビット中 r 次元を「使わない」工夫：⃗xn+p = ⃗xn+q + ⃗xn+1B + ⃗xnC
⇒ 2wp−r − 1 を素数にするように選ぶことができる。
素数に対する周期判定は容易なため、超長周期が実現可能: 2624×32−31 − 1 = 219937 − 1.
注：2p − 1 の形の素数をメルセンヌ素数、そのときの p をメルセンヌ指数という。
現在 48 個知られており、19937 は 24 番目のメルセンヌ指数。
2013 年 1 月時点で知られている最大の指数は 257885161
30

開発秘話その３メモリの一部を使わないことによって一般性を増し、周期をメルセンヌ素数にする考えは’94 ごろ思いついた。
が、自分でパラメータ探索プログラムを書いてもうまく動かなかった。
当時すでに研究を純粋数学分野にシフトしていたため、このアイデアをほったらかした。

’95 慶応大学理工学部数理科学科講師として赴任。
榎本彦衛教授の修士学生であった西村拓士氏にMT のアイデアを話し、実験を始めた。
理論と実験が合わなかったとき：
僕の証明の間違い＝西村氏のプログラムが正しい確率 9 割。
⇐ 彼の力が無ければ MT は実現できなかった可能性大。
31

普及への道96 年国際会議で MT を口頭発表。

96 年ザルツブルグ大学の乱数研究グループホームページでニュースとして取り上げられる。
97 年 10 月 27 日朝日新聞夕刊記事（内村直之記者）

32

33

どうして新聞に載ったか？
⇒ 96 年夏ごろ朝日新聞に「秋の夜の数学」というシリーズ企画を内村さんが書いていた。
この企画で、慶応理工の教授が内村さんに取材を受けた。
その教授が、MT について内村さんに売り込んだ。
内村さんが取材に来た。大きな記事になってびっくりした。
（東大工学部の伏見正則教授に、その重要性を確かめてくれていた。）⇒ 慶応大学に電話での問い合わせが多数あった。
電話で回答するのには無理があった。
（プログラムを下さい、という依頼。）

34

普及への道（つづき）
回答できないのでホームページを立ち上げた。
ついでに英文のホームページも立ち上げて、無料でダウンロード可能とした。
ライセンスフリー、商用利用も許可なく可能とした。
⇒ 多数のボランティアが無償で多くの言語に移植し、高速バージョンを作って送ってくれた。
それをまた、ホームページに載せた。

35• なまけものなので、新聞記事にならなければホームページもつくらず忘れ去られた可能性大。

• インターネットがなければ、忘れ去られた可能性大。
（乱数の大御所の Knuth, Compagner, Marsaglia らが、「整数の方が 2 元体より乱数生成に適している」と主張しつづけており、MT の発表後も、古典的なアルゴリズムを推奨しつづけた。）• 有償にしていたら、普及しなかった可能性大。

36

普及への道（つづきその２）
ホームページで発表後、MT は世界中で使われ始めた。
多くの言語、ゲーム、商用プログラムで利用が広がった。

’98 年 JIS 規格の「ランダムサンプリング」の改訂にあたり、伏見正則教授が改訂委員長となり、「日本発の世界規格を目指して、MT を JIS 規格乱数の目玉としよう」という動きになり、JIS 規格に入った。
’07 年、ISO 規格で「ランダムサンプリング」を改訂するにあたり、統計数理研究所椿広計教授らが ISO 規格委員に入り、JIS 規格をもとにして MT を含めた標準疑似乱数の提案を行い、採択された。
(ただし、規格になる前にすでに広く普及していた。)
37

数学の効用周期と高次元分布性を保証するのに、ガロアに遡る近代数学が用いられている。
有限体、線形代数、メルセンヌ素数、Berlekamp-Massey 法、
格子縮約アルゴリズム、MacWilliams 恒等式などなど。
研究された当時は応用など思いもつかなかった 1 + 1 = 0 の数学が、今、実用されている。
終わり38

終わりに：今でも危険な疑似乱数が使われている
random: ’90-現在 UNIX 系 C 言語での標準的擬似乱数の最下位ビット 0-1 を見る。(原本博史-M ’07)
横軸：過去 31 回中の 1 の個数
縦軸：その条件下で、次の 8 回が全て 0 の確率：
0.014
0.012
0.01
0.008
0.006
0.004
0.002
1012141639

1820初期値を系統的に選んだときの、非乱数性
ranlux(カオス理論に基づく擬似乱数, Lüscher ’94 ) で、20, 25, 41 番目の出力（縦軸）を、初期シードを 1, 2, 3, . . . , 100（横軸）と動かしてプロット「20 番目の打者がいつもヒット」
0xFFFFFF
n=20
n=25
n=41

0xF00000
0xE00000
0xD00000
0xC00000
ranlux generated value

0xB00000
0xA00000
0x900000
0x800000
0x700000
0x600000
0x500000
0x400000
0x300000
0x200000
0x100000
0x000000
11020304050
seed40607080

90100

GNU Scientific Library に入っている58 個の擬似乱数発生法のうち、最新のものも含めて45 個にこのような問題が観測された。
Mersenne Twister では観測されなかった。
（M, 和田維作、倉本愛、芦原評 2007）⇒ モンテカルロ法には MT を使いましょう。
(注：MT は暗号用にはそのままでは使えません。)

41

蛇足

• 2008 年 5 月 16 日 Debian の OpenSSL に脆弱性:
SSL ／ SSH 暗号鍵の基になる乱数が推測可能⇐ 乱数発生アルゴリズムが悪かった。
• 2010 年 8 月 13 日 乱数生成器の問題によって、Android 版 BitCoin アプリで生成したウォレットは脆弱な秘密鍵を利用している可能性がある⇐ 初期値を十分ランダムにとっていなかった。

42

• 2010 年 10 月 02 日ポケットモンスター BW「乱数調整で強力なモンスターをゲット」⇐ MT を使っているが、初期値を MAC アドレスと時刻のみから決めていた。
時刻を選んで (MAC アドレスと MT を使って) 乱数を再現することで、望みのパラメータをもつモンスターの出現時期がわかる。
目的にあった疑似乱数を、正しく使いましょう。
御清聴ありがとうございました。

43

最近の研究• SIMD-oriented Fast Mersenne Twister(SFMT, 斎藤睦夫-M, 2006）
最近の CPU は 128 ビット演算命令(Single Instruction Multiple Data, SIMD)を持っている。

SIMD 命令やパイプライン処理を最大限に生かした漸化式を考案。

SFMT は、SIMD を使わなくても MT よりも高速。
SIMD を使うと 4 倍程度高速。
周期は MT と同じ、高次元均等分布性は MT より良い。
（ホームページで公開中、ダウンロード数万件）
44

• WELL (Panneton, L’Ecuyer, M ’06)
高次元均等分布性について、理論的上限を達成し、MT に近い高速性を持つ• CryptMT: 暗号用 MT, MT の出力 32 ビット整数を奇数化して整数積算 ( mod 232) し、最上位 8 ビットだけ使う(斎藤睦夫、西村拓士、萩田真理子-M, 2007)• Graphic Processor (高並列性・低機能マルチプロセッサ）
用の Mersenne Twister の開発、（斎藤睦夫 2010）• Dynamic Creator:多数のマシンで別々の疑似乱数を使うため、多数のパラメータを生成。(M-西村 98, 斎藤 2010)
特に原瀬らの「最短基底計算の高速化」が有用。

• WAFOM：準モンテカルロ法の新指標 (M-斎藤-Matoba2013)
45

総まとめ
擬似乱数は極めて大量に用いられる。微細な統計的偏りや、初期値へのわずかな依存性が、計算機の高速化・大規模化に伴いシミュレーションを狂わせる可能性がある。
この際、使用者は擬似乱数が原因だと中々気づけない。

⇒

精密なデタラメさを高速生成する必要性
それに答えるのが、「1+1=0 の数学」

講演者は、
「擬似乱数を MT 系に変えたらうまく動いた」というメールをたくさんもらっている。

46