Permuted congruential generatorのソースを表示

'''Permuted congruential generator''' ('''PCG''') は[[2014年]]に発表された[[擬似乱数]]生成器。[[線形合同法]]の出力を加工することで統計的性質を改善したものである。高速、省メモリかつ小さいコードサイズで非常に良い統計的性質を示す<ref>{{cite web
 |title=Testing non-cryptographic random number generators: my results
 |first=Daniel |last=Lemire
 |date=22 August 2017  |accessdate=2017-10-03
 |url=https://lemire.me/blog/2017/08/22/testing-non-cryptographic-random-number-generators-my-results/
}}</ref><ref>{{cite web
 |title=Testing the PCG random number generator
 |first=John D. |last=Cook
 |date=7 July 2017  |accessdate=2017-10-03
 |url=https://www.johndcook.com/blog/2017/07/07/testing-the-pcg-random-number-generator/
}}</ref><ref>{{cite web
 |title=Testing RNGs with PractRand
 |first=John D. |last=Cook
 |date=14 August 2017  |accessdate=2017-10-03
 |url=https://www.johndcook.com/blog/2017/08/14/testing-rngs-with-practrand/
}}</ref><ref name=pcg-passes>{{cite web
 |title=PCG Passes PractRand
 |first=M.E. |last=O'Neill
 |date=29 July 2017 |accessdate=2017-11-03
 |url=http://www.pcg-random.org/posts/pcg-passes-practrand.html
}}</ref><ref name=ONeill>{{cite techreport
 |title=PCG: A Family of Simple Fast Space-Efficient Statistically Good Algorithms for Random Number Generation
 |first=Melissa E. |last=O'Neill
 |publisher=[[Harvey Mudd College]]
 |id=HMC-CS-2014-0905
 |date=5 September 2014
 |url=http://www.pcg-random.org/pdf/hmc-cs-2014-0905.pdf
}}</ref>。

PCGは古典的な線形合同法と3つの点で異なる。
* 線形合同法で使われる状態が大きい（通常、出力の2倍）。
* [[除法|法]]として[[2の冪]]を使う（これにより偏りがなく最長の周期を効果的に実装できる）。
* 状態は直接出力されるのではなく、最も周期の長いビットによってビット回転やビットシフトの量が決まり、出力に影響する。
下位ビットの短い周期という線形合同法の弱点は状態に依存するビット回転量により解決する{{r|ONeill|p=31–34}}。

== PCGの種類 ==
PCGにはいくつかの種類がある。内部で使用する[[線形合同法]]は8ビットから128ビットのものが使用される。実用的な用途では64ビットと128ビットのものだけが推奨される。それより小さいものは統計的性質のテストに使用される。

線形合同法で加算される定数は異なる[[乱数列]]を生成するために異なる値を使用できる<ref name=critiquing>{{cite web
 |title=Critiquing PCG's streams (and SplitMix's too)
 |first=M.E. |last=O'Neill
 |date=10 August 2017 |accessdate=2017-11-03
 |url=http://www.pcg-random.org/posts/critiquing-pcg-streams.html
}}</ref>。定数は任意の奇数で、明示的に指定する必要はない。最下位のビットを1にすれば状態を保持する変数のアドレスを使ってもよい。

状態を出力に変換する方法がいくつかある。どの方法もうまく働くが、それぞれ余裕の幅が異なる{{r|ONeill|p=39}}。変換方法は以下の要素によって構築される。
* RR: ランダムな（入力に依存する）ビット回転で、入力の半分の長さの出力になる。<math>2^{b}</math>ビットの入力の最上位の<math>b-1</math>ビットが回転量を決めるのに使われ、そのすぐ下の<math>2^{b-1}</math>ビットが右に回転され出力され、残りの下位<math>2^{b-1}-(b-1)</math>ビットは捨てられる。
* RS: ランダムな（入力に依存する）ビットシフト。ビット回転の計算時間が長いときに使用する。これも入力の半分の長さの出力になる。<math>2^{b}</math>ビットの入力の最上位の<math>b-3</math>ビットによりシフト量が決定され、すぐ下の<math>2^{b-1}+2^{b-3}-1</math>ビットに適用され、その中の下位<math>2^{b-1}</math>ビットが出力される。残りの下位<math>2^{b-1}-2^{b-3}-(b-4)</math>ビットは捨てられる。
* XSH: [[Xorshift]]の操作 (<code>x ^= x >> (定数)</code>)。定数は次の操作で捨てられないビットの半分（端数切り捨て）となるようにする。
* XSL: XSHの (定数) の部分を全体の長さの半分にしたもの。
* RXS: XSHの (定数) の部分をランダムな（入力に依存する）量にしたもの。
* M: 定数の乗算。

以下に挙げるものは推奨される組み合わせである。擬似コードも示す。
* XSH-RR: xorshiftで上位のビットを下位のビットに混ぜ、ビット 63-59 で回転量を決定しビット 27-58 に適用する。
*:（64→32ビット）<code>count = (int)(x >> 59); x ^= x >> 18; return rotr32((uint32_t)(x >> 27), count);</code>
* XSH-RS: 上の方法よりシフト量を決めるのに使われるビットが少ない。
*:（64→32ビット）<code>count = (int)(x >> 61); x ^= x >> 22; return (uint32_t)(x >> (22 + count));</code>
* XSL-RR: XSH-RRを単純化したもので、これは[[64ビット]]マシンで128ビットの状態を持つ実装に最適化されている。
*:（128→64ビット）<code>count = (int)(x >> 122); x64 = (uint64_t)(x ^ (x >> 64)); return rotr64(x64, count);</code>
* RXS-M-XS: 状態の半分を出力する場合、この方法は最も遅いが最も強い方法である。以下に示すように状態の全体を出力する場合は最も弱い方法になる。これを使う場合は状態のサイズが32ビットまたは64ビットに制限される。
*:（32→32ビット）<code>count = (int)(x >> 28); x ^= x >> (4 + count); x *= 277803737u; return x ^ (x >> 22);</code>
*:（64→64ビット）<code>count = (int)(x >> 59); x ^= x >> (5 + count); x *= 12605985483714917081u; return x ^ (x >> 43);</code>
* XSL-RR-RR: RXS-M-XSと同様に128ビットの状態から128ビットの状態に変換する。
*:（128→128ビット）<code>count = (int)(x >> 122); low64 = rotr64((uint64_t)(x ^ (x >> 64)), count); high64 = rotr((uint64_t)(x >> 64), low64 & 63); return (uint128_t)high64 << 64 | low64;</code>

それぞれの変換は元に戻すことができる（つまり[[全単射]]）か[[切り捨て]]である。したがってそれらを[[写像の合成|合成]]したものはある出力が得られる入力の個数は常に同じである。

もし<math>2^{128}</math>より長い周期が必要になった場合は複数の生成器を使って拡張することができる。例えば生成器を2つ使う場合、生成器1の出力が0になるごとに生成器2の出力を1つ進めて、最終的な出力は生成器1の出力と生成器2の出力の和とすることでそれぞれの周期の積の周期になる。

== コード例 ==
ほとんどの場合に推奨される生成器は64ビットの状態を持ち32ビットを出力するPCG-XSH-RRである{{r|ONeill|p=43}}。具体的な実装の一例を示す。

<syntaxhighlight lang="C">
#include <stdint.h>
static uint64_t       state      = 0x4d595df4d0f33173;		// 何らかの初期状態
static uint64_t const multiplier = 6364136223846793005u;
static uint64_t const increment  = 1442695040888963407u;	// 任意の奇数

static uint32_t rotr32(uint32_t x, unsigned r)
{
	return x >> r | x << (-r & 31);
}

uint32_t pcg32(void)
{
	uint64_t x = state;
	unsigned count = (unsigned)(x >> 59);		// 59 = 64 - 5

	state = x * multiplier + increment;
	x ^= x >> 18;								// 18 = (64 - 27)/2
	return rotr32((uint32_t)(x >> 27), count);	// 27 = 32 - 5
}

void pcg32_init(uint64_t seed)
{
	state = seed + increment;
	(void)pcg32();
}
</syntaxhighlight>

この実装ではスーパースカラーのプロセッサーで[[命令レベルの並列性]]を最大化するために変換を線形合同法の計算後に行っている{{r|ONeill|p=43}}。

加算を使用せず、XSH-RSを使うことで若干高速化した実装を示す。この実装では周期は<math>2^{62}</math>になり、XSH-RRより弱くなる。

<syntaxhighlight lang="C">
static uint64_t mcg_state = 0xcafef00dd15ea5e5u;	// 奇数でなければならない

uint32_t pcg32_fast(void)
{
	uint64_t x = mcg_state;
	unsigned count = (unsigned)(x >> 61);	// 61 = 64 - 3

	mcg_state = x * multiplier;
	x ^= x >> 22;
	return (uint32_t)(x >> (22 + count));	// 22 = 32 - 3 - 7
}

void pcg32_fast_init(uint64_t seed)
{
	mcg_state = 2*seed + 1;
	(void)pcg32_fast();
}
</syntaxhighlight>

時間のかかる64ビットの乗算が残っているため高速化の効果は小さく、極端な場合を除き最初に挙げた実装のほうが良い。ただし、この高速化した実装も統計的検定を満足する{{r|pcg-passes}}。

[[32ビット]]プロセッサーで実行する場合、64ビット同士の乗算は32ビット同士の乗算3回が必要になる。これを2回に減らすには、32ビットの範囲に収まる乗数（0xf13283adなど{{r|critiquing}}）を使用する。

== 他の生成器との比較 ==
PCGは[[TestU01]]のBigCrushに合格するのに必要な最小限の状態のビット数を決定する過程で開発された<ref>{{cite web
 |title=Visualizing the Heart of Some PRNGs
 |first=M.E. |last=O'Neill
 |date=20 August 2017 |accessdate=2017-11-03
 |url=http://www.pcg-random.org/posts/visualizing-the-heart-of-some-prngs.html
}}</ref>。BigCrushは周期2^35の列を検出するのに十分なデータをテストするため、理想的な乱数生成器でも36ビットが必要である。十分に大きい状態の場合、[[平方採中法]]のような非常に悪い生成器でもBigCrushに合格することがある<ref>{{cite web
 |title=Too Big to Fail
 |first=M.E. |last=O'Neill
 |date=20 August 2017 |accessdate=2017-11-03
 |url=http://www.pcg-random.org/posts/too-big-to-fail.html
}}</ref>。小さい状態で合格することはそのアルゴリズムの質を示し、余裕の幅がどれだけあるかも示す。

PCG-RXS-M-XSはBigCrushに36ビット（理論上最小）で合格し、PCG-XSH-RRは39ビット、PCG-XSH-RSは49ビット必要である。ちなみに、xorshift*は40ビット{{r|ONeill|p=19}}、そして[[メルセンヌ・ツイスタ]]は19937ビットもの状態を持つのにもかかわらずBigCrushに合格できない<ref>{{cite journal
 |first1=Pierre |last1=L'Ecuyer  |first2=Richard |last2=Simard
 |url=http://www.iro.umontreal.ca/~lecuyer/myftp/papers/testu01.pdf
 |title=TestU01: A C library for empirical testing of random number generators
 |journal=[[ACM Transactions on Mathematical Software]]
 |volume=33 |issue=4 |page=22-1–22-40 |date=August 2007
 |doi=10.1145/1268776.1268777 |citeseerx=10.1.1.499.2830
}}</ref>。

== 関連項目 ==
* [[線形合同法]]
* [[TestU01]]

== 参考文献 ==
{{Reflist}}

== 外部リンク ==
* [http://www.pcg-random.org/ PCG, A Family of Better Random Number Generators] 公式サイト
* [https://www.reddit.com/r/programming/comments/2momvr/pcg_a_family_of_better_random_number_generators/ PCG, A Family of Better Random Number Generators — inspired by /r/programming!] 開発者によるReddit上の議論

[[Category:乱数]]
[[Category:数学に関する記事]]