PCG(permuted congruential generator)とは
PCG自体の元ネタは PCG: A Family of Simple Fast Space-Efficient Statistically Good Algorithms for Random Number Generation にあります。(この論文は疑似乱数の統計的性質の比較とかいろいろやってるデカいやつなので、単にアルゴリズムだけ知りたいなら後半まで読み飛ばしてよい)
ざっくりいうと PCG は LCG と順列関数を組み合わせて出力品質を向上させたものです。
LCG(線形合同法)については前書いたのでこちらをどうぞ
↑でも触れているけど、LCGにはmを2べきにすると下位bitの周期が短くなる既知の性質がありますがこのへんも改善されてるよう。
順列関数にいくつか種類がありますが、今回は NumPy などにも採用されている DXSM といわれるものを採用したパターンについて整理します。
DXSM とは
いきなりなんですが、興味あって調べようとしたんだけど厳密な仕様が見つけられなかったです。始まる前から負けている。現時点だと実装か NumPy ドキュメントくらいしかまとまった情報ないんじゃないか...？
(最近 Go1.22 で PCG-DXSM の実装入ったんですが、そこのコードからも実装とかドキュメントへのリンクしかなかったから多分そう)
以下わかる範囲での説明を試みます。
DXSM は "double xorshift multiply" の略っぽい です。2019に考案されたようなので比較的あたらしめ。
頑張って挙動に関連するドキュメント探すと、NumPyのドキュメントに
https://numpy.org/devdocs/reference/random/bit_generators/pcg64dxsm.html#numpy.random.PCG64DXSM
The PCG64DXSM state vector consists of 2 unsigned 128-bit values, which are represented externally as Python ints. One is the state of the PRNG, which is advanced by a linear congruential generator (LCG). The second is a fixed odd increment used in the LCG.
とあるので、内部状態 128bit *2 をつかって64bit疑似乱数を吐けるやーつぽい。内部状態は LCG で作るやつがひとつと、もう一つが LCG のための固定 increment 値らしい。雰囲気LCGの状態遷移時にも多少気にすることがあるのかな。
やってる操作を見ると shift なりなんなりを掛け合わせて値をこねてる様子。操作もそんなに難しくないが、さっきのリンクのコメントから読み取れるノリとしては十分掻き回せていい感じっぽい。
かき回す箇所だけピックアップすると以下
https://github.com/imneme/pcg-cpp/blob/428802d1a5/include/pcg_random.hpp#L1043-L1048
lo |= 1; hi ^= hi >> (xtypebits/2); hi *= xtype(cheap_multiplier<itype>::multiplier()); hi ^= hi >> (3*(xtypebits/4)); hi *= lo; return hi;
みたかんじ
- hi 下位半分でXOR
- 定数で mul
- hi 上位1/4でXOR
- hi と下位1bitをたてた lo で mul
の操作をしている。なるほどそういう感じで置換しとるのね。
Go 1.22 の実装をみてみる
せっかくなので実装を眺めて雰囲気をつかんでいきます！(というか詳細な仕様を発見できなかったので実装に学ぶしかなかった)
他の実装と一貫してるかとかは特に確認しないのであしからず
構造
LCG ですすむ内部状態 128bitを持っているよう
https://cs.opensource.google/go/go/+/refs/tags/go1.22.1:src/math/rand/v2/pcg.go;l=17-20;bpv=0;bpt=0
type PCG struct { hi uint64 lo uint64 }
多分これが NumPy ドキュメントにあった One is the state of the PRNG, which is advanced by a linear congruential generator (LCG). かな
状態遷移
きになる状態遷移関数 next() はこんなかんじ
https://cs.opensource.google/go/go/+/refs/tags/go1.22.1:src/math/rand/v2/pcg.go;l=74-98;bpv=0;bpt=0
func (p *PCG) next() (hi, lo uint64) { // https://github.com/imneme/pcg-cpp/blob/428802d1a5/include/pcg_random.hpp#L161 // // Numpy's PCG multiplies by the 64-bit value cheapMul // instead of the 128-bit value used here and in the official PCG code. // This does not seem worthwhile, at least for Go: not having any high // bits in the multiplier reduces the effect of low bits on the highest bits, // and it only saves 1 multiply out of 3. // (On 32-bit systems, it saves 1 out of 6, since Mul64 is doing 4.) const ( mulHi = 2549297995355413924 mulLo = 4865540595714422341 incHi = 6364136223846793005 incLo = 1442695040888963407 ) // state = state * mul + inc hi, lo = bits.Mul64(p.lo, mulLo) hi += p.hi*mulLo + p.lo*mulHi lo, c := bits.Add64(lo, incLo, 0) hi, _ = bits.Add64(hi, incHi, c) p.lo = lo p.hi = hi return hi, lo }
cpp 実装のリンクが貼られてるが、定数はそっちとあわせてるよう。この incHi, incLo の二つの定数が NumPy ドキュメントで The second is a fixed odd increment used in the LCG. と言わせていたふたつめの内部状態かな。
コメントよるとなんか NumPy だと参照する mul 定数が64bitらしいんだけど、多倍長整数演算なのでどっちみち複数回やらないとあかんくてサボれるの1回くらいだから気にせんでよくね、だから Go では普通に128bitでやってるよ、みたいなことがかかれている(多分)
まあやってることとしては大したことではなく、普通に mul とって inc して LCG を進めてるだけかな。LCG の段階でも DXSM のための increment 定数を掛け合わせたりしてるので、そのへんは多少考慮ありといった感じだろうか。
面白ポイントとしては hi x hi がないこと。hi x hi は確実に hi の桁を超えるため算出する必要がない。かしこ〜い
多倍長整数演算チョットワカルなので計算の話とか、戻り値にcarry含まれてるのとか前置きなく出てきても無事読めたんですが、あまり前提知識ない方は過去に書いた記事もどうぞ。実装を簡素にした上で基礎的なことをまとめています。
DXSM
// Uint64 return a uniformly-distributed random uint64 value. func (p *PCG) Uint64() uint64 { hi, lo := p.next() // XSL-RR would be // hi, lo := p.next() // return bits.RotateLeft64(lo^hi, -int(hi>>58)) // but Numpy uses DXSM and O'Neill suggests doing the same. // See https://github.com/golang/go/issues/21835#issuecomment-739065688 // and following comments. // DXSM "double xorshift multiply" // https://github.com/imneme/pcg-cpp/blob/428802d1a5/include/pcg_random.hpp#L1015 // https://github.com/imneme/pcg-cpp/blob/428802d1a5/include/pcg_random.hpp#L176 const cheapMul = 0xda942042e4dd58b5 hi ^= hi >> 32 hi *= cheapMul hi ^= hi >> 48 hi *= (lo | 1) return hi }
まあそうだねといった趣。ふつうに cpp にある実装まんまという感じですね。
コメントに XSL-RR だと別の書き方になるけど、作者からもコメントあったし DXSM にするで(意訳)みたいなのもあってよかった。
めちゃめちゃ実装へのリンクつけられてるからやはり詳細な仕様はまだこの世にないのか...？
感想
気になって調べた結果PCG自体は概ね理解した。DXSM については詳細なドキュメントが発見できなかったが実装でどのような操作が行われているかは理解した。
あとPCGの元ネタ論文は疑似乱数を比較する上でいくつかの観点を整理していたりしてけっこう面白そうだった。めちゃめちゃ量あったので元気な時に比較まわりもちゃんと読んでみようかな...
