モンテカルロ法による円周率の推定（その4 PRNG）

「モンテカルロ法による円周率の推定」もひととおり終わったので，今回は擬似乱数生成器（pseudo random number generator）の話。

math/rand の擬似乱数生成器

Go 言語では math/rand パッケージで擬似乱数を取り扱えることは「その1」で紹介した通り。 math/rand パッケージに実装されている擬似乱数生成器はラグ付フィボナッチ法（Lagged Fibonacci Generator）のバリエーションらしい。

If I am not mistaken again, the generator is an ALFG (Additive Lagged Fibonacci Generator, thats what Wikipedia calls it). Knuth describes the algorithm in Volume 2 of The art of computer programming in section 3.2.2 (around equation 7). Both Wikipedia and Knuth state the parameter combination 607,273 as possible combination with a period length of 2^(e-1) * (2^607-1) where e is the length of the random number in bits.
I actually found a few references examining its properties and it seems to be a good rng so faar, but there is still seems to be a lack of mathematical background and it is fairly easy to get into trouble by not seeding properly.

via [Announce] A rand package for high quality 64bit random numbers (possibly go2)

ラグ付フィボナッチ法は線形合同法（後述）を改善することを目的としたものでフィボナッチ数の生成法を元にしている。

Lagged Fibonacci generator - Wikipedia

ラグ付フィボナッチ法は以下の式で表される。

\begin{alignat*}{2} S_{n} \equiv S_{n-j} * S_{n-k} \pmod{m}, & \; & 0 \lt j \lt k \end{alignat*}

ラグ付フィボナッチ法は $ * $ 演算子によってバリエーションがあるが math/rand パッケージの実装では加算を使うため “Additive Lagged Fibonacci Generator” ということらしい。ソースコードで言うとこの部分かな。

// Int63 returns a non-negative pseudo-random 63-bit integer as an int64.
func (rng *rngSource) Int63() int64 {
    rng.tap--
    if rng.tap < 0 {
        rng.tap += _LEN
    }

    rng.feed--
    if rng.feed < 0 {
        rng.feed += _LEN
    }

    x := (rng.vec[rng.feed] + rng.vec[rng.tap]) & _MASK
    rng.vec[rng.feed] = x
    return x
}

擬似乱数生成器のバリエーション

math/rand パッケージでは以下の interface を持つ別の擬似乱数生成器を使うことができる。

// A Source represents a source of uniformly-distributed
// pseudo-random int64 values in the range [0, 1<<63).
type Source interface {
    Int63() int64
    Seed(seed int64)
}

以下に2つほど紹介する。

線形合同法

線形合同法（Linear Congruential Generator）は昔の擬似乱数ライブラリでよく使われていたアルゴリズムで，以下の式で表される。

\[ X_{n+1} = (A \times X_{n} + B) \bmod M \]

定数 $A$ および $B$ の与え方により幾つかバリエーションがある。

線形合同法のメリットは実装サイズが小さく計算量も少ない点だろうか。一方デメリットとしては，多次元で疎に分布する性質があり，周期も小さいため乱数を大量に発生させる必要がある科学技術シミュレーションなどには向かないと言われている。このためメモリサイズが限られるマイクロ・コントローラのようなものでもない限り線形合同法が使われることはなくなった。

Go 言語でわざわざ線形合同法を実装しているパッケージは少ないのだが，たとえば以下のパッケージがある¹。

davidminor/gorand: Basic golang implementation of a permuted congruential generator for pseudorandom number generation

Mersenne Twister

Mersenne Twister とは松本眞・西村拓士両氏によって1996年に発表された擬似乱数生成アルゴリズムである。他の擬似乱数生成アルゴリズムと比べて以下の特徴があるそうだ。（「Mersenne Twister とは?」より）

従来の様々な生成法の欠点を考慮して設計されています。
従来にない長周期，高次元均等分布を持ちます。（周期が $2^{19937}-1$ で²、623次元超立方体の中に均等に分布することが証明されています。）
生成速度がかなり速い。（処理系にもよりますが、パイプライン処理やキャッシュメモリのあるシステムでは、Cの標準ライブラリの rand() より高速なこともあります。なお、開発当時には cokus 版は rand() より4倍程度高速でしたが、その後 ANSI-C の rand() が LCG 法から lagged-fibonacci に変更されたこともあり、2002年現在 rand と MT の速度差はあまりありません。）
メモリ効率が良い。（32ビット以上のマシン用に設計された mt19937.c は、 624ワードのワーキングメモリを消費するだけです。 1ワードは32ビット長とします。

ちなみに Mersenne Twister のオリジナル・コードは BSD ライセンスで提供されている。

公式のリポジトリには C/C++ による実装のみのようだが， Go 言語で実装している人もいるようである。

seehuhn/mt19937: An implementation of Takuji Nishimura’s and Makoto Matsumoto’s Mersenne Twister pseudo random number generator in Go.
nutterts/randgen: Pseudo Random Number Generators implementing the Go(lang) math/rand.Source Interface
farces/mt19937_64: Mersenne Twister (int64) for Go
cuixin/goalg: golang algorithms
cpmech/gosl: Go scientific library : SFMT や TinyMT に対応。オリジナルのコードを cgo で結合しているのでクロス環境では注意

擬似乱数生成器を組み込む

では，先ほど紹介した擬似乱数生成器を今回のコードに組み込んでみることにしよう。こんな感じ。

package gencmplx

import (
    "math/rand"

    "github.com/davidminor/gorand/lcg"
    "github.com/seehuhn/mt19937"
)

//RNGs is kind of RNG
type RNGs int

const (
    NULL RNGs = iota
    GO
    LCG
    MT
)

//NewRndSource returns Source of random numbers
func NewRndSource(rng RNGs, seed int64) rand.Source {
    switch rng {
    case LCG:
        return lcg.NewLcg64(uint64(seed))
    case MT:
        mt := mt19937.New()
        mt.Seed(seed)
        return mt
    default:
        return rand.NewSource(seed)
    }
}

gencmplx.NewRndSource() 関数で rand.Source オブジェクトを生成する³。これを「その1」で作った gencmplx.New() 関数に渡せばよい。

CLI は以下のように調整してみた。

$ go run main.go estmt --help
Estimate of Pi with Monte Carlo method.

Usage:
  pi estmt [flags]

Flags:
  -e, --ecount int   Count of estimate (default 100)

Global Flags:
  -p, --pcount int       Count of points (default 10000)
  -r, --rsource string   Source of RNG (GO/LCG/MT) (default "GO")

で，それぞれの擬似乱数生成器で評価を行ってみようと思ったのだが，今回のケースに限ってはあまり違いが出ないようである。

まずは線形合同法の場合。

$ go run main.go estmt -e 10000 -p 100000 -r LCG > estmt100k-lcg.dat
random number generator: LCG
minimum value: 3.12204
maximum value: 3.16224
average value: 3.14164
standard deviation: 0.00524 (68.3%)

次は Mersenne Twister の場合。

$ go run main.go estmt -e 10000 -p 100000 -r MT > estmt100k-mt.dat
random number generator: MT
minimum value: 3.12380
maximum value: 3.16140
average value: 3.14165
standard deviation: 0.00517 (67.8%)

もっと多次元だったりすると変わってくるのかなぁ。

暗号技術用途の乱数生成器

Go 言語では暗号技術用途の乱数として crypto/rand パッケージが用意されている。これは math/rand とは互換性がない。

具体的には，UNIX 系のプラットフォームでは乱数生成に /dev/urandom デバイスを参照している⁴。また Windows プラットフォームでは CryptoAPI 2.0 の CryptGenRandom 関数を使っている⁵。

そもそも暗号技術用途の乱数生成器は科学技術シミュレーションやゲームで使う擬似乱数生成器とは要件が異なる。

従前の観点から統計的にテストされた乱雑性は、セキュリティ用途に要求される予測困難性と同等ではありません。
例えば、（CRC Standard Mathematical Tables からのランダムテーブルのような）広く利用可能な一定のシーケンスの利用は、攻撃者に対して非常に弱いです。これを学習したり、推測する攻撃者は、容易に（過去・未来を問わず）そのシーケンス [CRC] に基づいて、すべてのセキュリティを破ることができます。他の例として、AES を 1, 2, 3 ... のような連続した整数を暗号化する一定の鍵と共に使うことは、優れた統計的乱雑性をもつが予測可能な出力を作り出します。他方、6 面のサイコロを連続して転がして、その結果の値を ASCII にエンコードすることは、実質的に予測困難なコンポーネントをもちながらも「統計的に貧弱な出力」を作り出します。それゆえ、「統計的テストの合否は、『何かが予測不可能であるか否か、あるいは、予測可能であるか否か』を表さないこと」に注意してください。

RFC 4086 - セキュリティのための乱雑性についての要件より

暗号技術用途の乱数生成器は，暗号分野においては中核技術のひとつであるが，一度に大量の乱数を生成させる必要のある科学技術シミュレーションなどの用途には向かない。上手く使い分けてほしい。

ブックマーク

Go 言語に関するブックマーク集はこちら。

参考図書

プログラミング言語Go (ADDISON-WESLEY PROFESSIONAL COMPUTING SERIES): Alan A.A. Donovan (著), Brian W. Kernighan (著), 柴田芳樹 (翻訳); 丸善出版 2016-06-20; 単行本（ソフトカバー）; 4621300253 (ASIN), 9784621300251 (EAN), 4621300253 (ISBN), 9784621300251 (ISBN); 評価

著者のひとりは（あの「バイブル」とも呼ばれる）通称 “K&R” の K のほうである。この本は Go 言語の教科書と言ってもいいだろう。

reviewed by Spiegel on 2016-07-13 (powered by PA-APIv5)

暗号技術入門第3版　秘密の国のアリス: 結城浩 (著); SBクリエイティブ 2015-08-25 (Release 2015-09-17); Kindle版; B015643CPE (ASIN); 評価

SHA-3 や Bitcoin/Blockchain など新しい知見や技術要素を大幅追加。暗号技術を使うだけならこれ1冊でとりあえず無問題。

reviewed by Spiegel on 2015-09-20 (powered by PA-APIv5)

ただし github.com/davidminor/gorand/lcg には不具合があって Int63() 関数で負の値を出力する場合がある。とりあえず fork 版の github.com/spiegel-im-spiegel/gorand/lcg で修正している。 Pull request も出したけど，古いコードだし，もうメンテしてないかなぁ。 ↩︎
$2^{19937}-1$ はメルセンヌ素数で Mersenne Twister の名前の由来になっている。 Mersenne Twister では周期サイズごとに複数の実装があるが， $2^{19937}-1$ がポピュラーな実装として広く使われているようだ。 ↩︎
Go 言語におけるオブジェクトの多態性については「Go 言語における「オブジェクト」」を参考にどうぞ。 ↩︎
/dev/urandom はハードウェア・デバイスから十分なエントロピー源が得られない場合は内部で疑似乱数生成器を使用する。このため一時は /dev/urandom の脆弱性が疑われたが，現時点では事実上は問題ないとされている。一方で，スマートデバイスのような場合はハードウェア・デバイスからのエントロピー源だけでは外部から推測され易いため，性能のよい疑似乱数生成器を組み合わせるほうが有効になる場合もあるようだ。 ↩︎
CryptGenRandom 関数の内部実装は公開されていないが，やはりキーボードやマウス等のデバイスの挙動をエントロピー源とし， NIST の SP800-90 勧告に従った実装をしているようである。余談だが SP800-90 は乱数生成の一部のアルゴリズムで脆弱性が発見され（これがまた NSA 絡みだったものだから大騒ぎになった），現在は修正版の SP800-90A Revision 1が発行されている。（参考：擬似乱数生成アルゴリズム Dual_EC_DRBG について） ↩︎