疑似乱数生成

2020年2月2日 #PRNG #MersenneTwister #xorshift

概要

乱数 (random number) はランダムに選択された数のこと。特定の出現パターンを持たず、選択される値が予測できないという性質を持っており、ゲームや科学技術シミュレーション、暗号セキュリティの分野で重要な役割を持っている。

ランダムなビット列は裏表に 0 と 1 のラベルが付けられた "公正な" コインによるコイントスの結果と考えることができる。次のビットが 0 になるか 1 になるかの確率は正しく $1/2$ であり、各コイントスは完全に独立していて前の結果が後の結果に影響することはない。しかし乱数生成のためにコイントスを使用することは明らかに非効率である。コンピュータの分野ではこのような真の乱数に近い理想的な乱数を演算によって生成する疑似乱数生成器 (pseudorandom number generator) が用いられている。

多くの暗号アプリケーションには乱数が必要である。例えば一般的な暗号システムでは乱数を使用して暗号用の鍵を生成しているし、署名や認証プロトコルのチャレンジなど、様々な局面で乱数を必要とする。

概要
シフトレジスタジェネレータ
1. Mersenne Twister
2. xorshift
線形合同法
その他の疑似乱数生成アルゴリズム
パフォーマンス
参照

コンピュータの分野での乱数列の生成には、(真の) 乱数生成器 (RNG; random number generator) と疑似乱数生成器 (pRNG; pseudorandom number generator) という 2 種類の基本的な生成器が使用される。これらの乱数生成器はいずれも 0 と 1 のビット列を生成し、それをバイト列や整数列に分割して乱数として扱う。

特に暗号用途では前方予測不可能性 (forward unpredictability) と後方予測不可能性 (backward unpredictability) の両方が満たされていなければならない。つまり、シードやある乱数から次に生成される乱数が予測不可能出なければならず、同様にある乱数からシードや前の乱数が予測不可能でなければならない。

真の乱数

真の乱数 (true random number) は非決定論的で予測や再現が完全に不可能な性質を持つ乱数である。真性乱数生成器 (TRNG; hardware true random number generator) と呼ばれる専用のデバイスを使ってランダムであると予想される物理的現象、例えば抵抗やダイオードのノイズ、熱雑音などを計測し、場合によって測定プロセスで発生する可能性のあるバイアスを補正して真の乱数とする。

このような真の乱数は Unix 系 OS であれば /dev/random から得ることができる。多くの OS で真の乱数を取得する手段を提供しているが、対象とする物理現象に十分なエントロピーが得られるまで処理をブロックすることから生成速度が遅く、一度に大量の乱数列を生成するのには向いていない。

真の乱数は研究目的のいくつかの Web サイトでも入手することができる。例えば原子崩壊に基づく HotBit や、大気ノイズに基づくrandom.org などがある。

疑似乱数

疑似乱数 (pseudo-random number) は疑似乱数生成器 (PRNG; pseudo-random number generator) と呼ばれる数理に基づいた計算アルゴリズムを使って生成されるランダム性の高い数値シーケンス。真の乱数ではないため疑似乱数と呼ばれる。この記事では疑似乱数生成アルゴリズムについて扱う。

疑似乱数生成アルゴリズムは大量の乱数を高速に発生させることができることから、コンピュータ上で行う統計分析やシミュレーション、テスト、ディスク領域をランダムビットで埋めるようなタスクのような、ある程度のランダム性だけを必要としているケースでは、決定論的で再現可能な性質を持つ疑似乱数が適しているケースが多い。

暗号論的疑似乱数

疑似乱数生成はさらに暗号や署名などのセキュリティ用途で利用することを意図した暗号論的疑似乱数生成器 (cryptographically secure PRNGs; CSPRNGs) かどうかに分類される。暗号論的疑似乱数はランダム性に加えて次の要件を満たしている必要がある。

アルゴリズムの出力シーケンスの最初の $k$ ビットが与えられたとき、次に出現するビットを 50% を大幅に超える確率で予測できる多項式時間アルゴリズムが存在しない。

暗号理論の文脈で多項式時間アルゴリズムとは、(ビット幅のような) 設定パラメータ $x$ に対して高々その定数乗 $O(x^a)$ 程度の計算で解けるアルゴリズムを意味する。

乱数の品質と検定

真にランダムな性質を持つビット列であれば列を観測して得た 0 の数は全体の $\frac{1}{2}$ に収束するはずである。また有限の計算資源と記憶領域に依存する疑似乱数生成アルゴリズムはある一定の期間で同じパターンが繰り返し出現する周期性を持っている。このようなビットの出現頻度や周期の長さは乱数アルゴリズムの品質指標の一つといえる。この品質は、例えば The Art of Computer Programming [7] では:

$\chi^2$ 検定 … 観測した事象がある多項分布に従っているかを評価する一般的な統計的仮説検定。例えば $[0,1)$ 範囲の乱数列を $0.1$ 刻みの区間にビニングして各ビンの頻度を算出して $\chi^2$ 検定を行うことができる。
Kolmogorov-Smirnov 検定 … 観測した事象がある確率分布にしたがっているかを、確率分布の累積密度関数と比較して評価する。無限の数列や $[0,1)$ 範囲に連続して分布するような乱数列に向いている。

といった統計的仮説検定の手法を基礎として、0 と 1 が均等に分布している等分布検定 (頻度検定) や、0 または 1 が期待値と同じ程度連続しているかの連検定といった 12 の観点での検定 (+ 理論的検定からスペクトル検定) が提案されている。NIST [5, 6] や Diehard、TestU01 といった検定スイートは、それらの個別の検定をまとめて (疑似) 乱数のランダム性の品質を実験的に検定するプログラムである。それぞれの検定がどのように機能するかの概要は文献 [8] を参照。

ただし、これらの検定スイートに含まれるアルゴリズムや多くのパラメータはヒューリスティックに選択されたものであり、また結果には統計的仮説検定に由来する蓋然性が含まている。したがって、これらに含まれる全ての検定に合格したから完全という保証はなく、また合格しなかった検定があるから不完全と明確に断言できるものではないが、それでもある程度の基準を満たすランダム性を持っていることを乱数の利用者が確証するためには有用である。

これらの検定スイートは暗号アプリケーションで使用するような強いランダム性を検定することを目的としており、自然科学や社会科学で観測される事象のランダム性を評価するのには向いていない。そのような用途でランダムさの一様性を簡易的に確かめるだけであればビニングした $\chi^2$ 検定だけでも有用である。また RMT 検定はそのような過度に強いランダム性を持たない事象向けと言われている (が、実験科学向けには必要な標本数が多い)。

シフトレジスタジェネレータ

過度に疎な多項式を使う線形合同法と異なり、シフトレジスタジェネレータ (shift-register generator) は線形帰還シフトレジスタの考え方に基づいて計算効率を向上させた疑似乱数ジェネレータである。

Mersenne Twister

Mersenne Twister [4] は $2^{19937}-1$ という非常に長い周期と 623 次元で均等分布をもつことで知られる疑似乱数生成アルゴリズム。現時点の Python, Ruby, R, MATLAB といったいくつかのプラットフォームでの標準 PRNG として利用されている。

現代的な CPU に実装されている SSE や AVX といった SIMD 拡張命令を使用することで速度を向上した SFMT、さらに GPU を使用して並列計算を行う MTGP、より小さなメモリ空間で動作する TinyMT などの変種が存在する。

xorshift

ビット演算のみを使用して $2^{32}-1$ から $2^{192}-1$ 周期の乱数を生成するアルゴリズム。線形合同法の結果に統計的補正を加える PCG としても使用されている。

論文翻訳: Xorshift RNGs

ビット演算のみを使用した非常に高速でコンパクトな Xorshift 擬似乱数生成アルゴリズムに関する 2003 年の論文。著者はキャリー付き乗算の論文にも携わっている。Abstract にあるようにこの論文自体はアイディアの説明であり、良い乱数・悪い乱数で何点かの間違いが指摘されている。…

2020年2月2日(Sun) 2003年の論文 #Xorshift #PRNG

線形合同法

線形合同法 (linear congruential generator; LCG) は歴史の長さと実装の単純さから最もよく知られている疑似乱数生成アルゴリズム。ただし現在ではランダム性の品質も計算効率も悪いと見なされており、シミュレーションなど乱数を大量に使用する用途ではより良い選択肢を考慮する必要がある。

線形合同法での乱数列の生成アルゴリズムは前の値 $X_n$ に乗算、加算、余剰演算を行って次の乱数 $X_{n+1}$ を得る。\[ X_{n+1} = (\alpha \times X_n + \beta) \bmod M \] ここで $\alpha$, $\beta$, $M$ はそれぞれ定数であり、$0 \lt \alpha \lt M$, $0 \le \beta$ である必要がある。線形合同法の詳しい解説や法 $M$、乗数 $\alpha$ の適切な選び方は [7] で論議されている。

線形合同法はいくつかのプログラミング言語において標準の疑似乱数生成アルゴリズムとして使用されてきた。例えば Java 標準の Randomクラスは以下のような実装を行っている。

protected synchronized int next(int bits) {
  seed = (seed * 0x5DEECE66DL + 0xBL) & ((1L << 48) - 1);
  return (int) (seed >>> (48 - bits));
}

ただし Java 8 から標準で SplitMax に基づく SplittableRandomを利用することができるため、過去の互換性以外の目的で線形合同法に基づいた Random を使用する必要はない。

その他の疑似乱数生成アルゴリズム

Salsa20 / ChaCha20: ChaCha20 は RFC 7539 で仕様策定が行われているストリーム暗号アルゴリズム。

パフォーマンス

Table 1. Rust `rand` crate performance on
Algorithm	Footprint [bits]	Period	nsec/op	🔐	TestU01
ChaCha8			---	✔
ChaCha12			---	✔
ChaCha20	1,024	$2^{128}$	---	✔
HC-128			---	✔
ISAAC 32			---	✔
ISAAC 64			---	✔
OS provided PRNG			---	?
MT19937-64	20,032	$2^{19937}-1$	---
SFMT	20,032	$2^{19937}-1$	---
TinyMT32	224	$2^{128}-1$	---
TinyMT64	256	$2^{128}-1$	---
PCG32			---
PCG32			---
PCG32			---
SplitMix64	64	$2^{64}$	---		✔
xorshift64		$2^{64}-1$	---		✔
xoshiro128+	128	$2^{128}-1$	---		✔
xoshiro128++	128	$2^{128}-1$	---		✔
xoshiro128**	128	$2^{128}-1$	---		✔
xoshiro256+	256	$2^{256}-1$	---		✔
xoshiro256++	256	$2^{256}-1$	---		✔
xoshiro256**	256	$2^{256}-1$	---		✔
xoshiro512+	512	$2^{512}-1$	---		✔
xoshiro512++	512	$2^{512}-1$	---		✔
xoshiro512**	512	$2^{512}-1$	---		✔
xoroshiro64*	64	$2^{64}-1$	---		✔
xoroshiro64**	64	$2^{64}-1$	---		✔
xoroshiro128+	128	$2^{128}-1$	---		✔
xoroshiro128++	128	$2^{128}-1$	---		✔
xoroshiro128**	128	$2^{128}-1$	---		✔

Benchmark Source (click to open)

$ cat main.rs

// rustup install nightly
// cargo +nightly bench
#![feature(test)]

extern crate rand;
extern crate test;
extern crate tinymt;

use rand::{Rng, SeedableRng};
use rand_chacha::*;
use rand_core::OsRng;
use rand_hc::Hc128Rng;
use rand_isaac::{Isaac64Rng, IsaacRng};
use rand_pcg::*;
use rand_xorshift::XorShiftRng;
use rand_xoshiro::*;
use test::Bencher;

macro_rules! rng_bench {
  ($identifier:ident, $generator:expr) => {
    #[bench]
    fn $identifier(b: &mut Bencher) {
      let mut rng = $generator;
      b.iter(|| rng.gen_range(0.0, 1.0));
    }
  };
}


// rand::random() is simply a shortcut for `thread_rng().gen()`.
//#[bench]
//fn bench_prng_random(b: &mut Bencher) {
//  b.iter(|| {
//    random::<f64>()
//  });
//}

// rng_bench!(bench_prng_thread_rng, thread_rng());   // thread-local version of StdRng
rng_bench!(bench_prng_xorshift, XorShiftRng::from_entropy());
rng_bench!(bench_prng_xoshiro128p, Xoshiro128Plus::from_entropy());
rng_bench!(bench_prng_xoshiro128pp, Xoshiro128PlusPlus::from_entropy());
rng_bench!(bench_prng_xoshiro128ss, Xoshiro128StarStar::from_entropy());
rng_bench!(bench_prng_xoshiro256p, Xoshiro256Plus::from_entropy());
rng_bench!(bench_prng_xoshiro256pp, Xoshiro256PlusPlus::from_entropy());
rng_bench!(bench_prng_xoshiro256ss, Xoshiro256StarStar::from_entropy());
rng_bench!(bench_prng_xoshiro512p, Xoshiro512Plus::from_entropy());
rng_bench!(bench_prng_xoshiro512pp, Xoshiro512PlusPlus::from_entropy());
rng_bench!(bench_prng_xoshiro512ss, Xoshiro512StarStar::from_entropy());
rng_bench!(bench_prng_xoroshiro64s, Xoroshiro64Star::from_entropy());
rng_bench!(bench_prng_xoroshiro64ss, Xoroshiro64StarStar::from_entropy());
rng_bench!(bench_prng_xoroshiro128p, Xoroshiro128Plus::from_entropy());
rng_bench!(bench_prng_xoroshiro128pp, Xoroshiro128PlusPlus::from_entropy());
rng_bench!(bench_prng_xoroshiro128ss, Xoroshiro128StarStar::from_entropy());
rng_bench!(bench_prng_isaac, IsaacRng::from_entropy());
rng_bench!(bench_prng_isaac64, Isaac64Rng::from_entropy());
rng_bench!(bench_prng_os, OsRng::default());
// rng_bench!(bench_prng_std, StdRng::from_entropy());    // 実行環境で最も適切な PRNG が選ばれる (ChaCha20 など)
rng_bench!(bench_prng_chacha8, ChaCha8Rng::from_entropy());
rng_bench!(bench_prng_chacha12, ChaCha12Rng::from_entropy());
rng_bench!(bench_prng_chacha20, ChaCha20Rng::from_entropy());
rng_bench!(bench_prng_hc128, Hc128Rng::from_entropy());
rng_bench!(bench_prng_splitmix64, SplitMix64::from_entropy());
rng_bench!(bench_prng_pcg32, Lcg64Xsh32::from_entropy());
rng_bench!(bench_prng_pcg64, Lcg128Xsl64::from_entropy());
rng_bench!(bench_prng_pcg64fast, Mcg128Xsl64::from_entropy());
rng_bench!(bench_prng_sfmt, sfmt::SFMT::from_entropy());
rng_bench!(bench_prng_tinymt32, tinymt::TinyMT32::from_entropy());
rng_bench!(bench_prng_tinymt64, tinymt::TinyMT64::from_entropy());

$ cat Cargo.toml

[package]
name = "prng-benchmark"
version = "0.1.0"
authors = ["TAKAMI Torao <koiroha@gmail.com>"]
edition = "2018"

# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html

[dependencies]
rand = "0.7"
rand_xorshift = "0.2"
rand_chacha = "0.2"
rand_isaac = "0.2"
rand_core = "0.5"
rand_xoshiro = "0.4"
rand_hc = "0.2"
rand_pcg = "0.2"
sfmt = "0.6"
tinymt = "1.0"