内部状態としてカウンターのみを持った疑似乱数生成器

TestU01のBigCrushを通すということ

疑似乱数の統計的検定プログラムとしてはTestU01というのが有名です（Wikipedia）。SmallCrush、Crush、BigCrushというテストの詰め合わせが用意されていて、順に実行時間が長くなっていき全てのテストに通るのも難しくなっていきます。
自作の疑似乱数生成器を試してみるのもとても簡単で、How to Test with TestU01にやりかたがまとめられています。呼ばれる度にuint32_tの疑似乱数を返す関数のポインタを渡してやるだけなので、最低限のC言語の知識さえあれば十分です。
このBigCrushに通れば安心して使える乱数なのかというとそこまで話は単純ではないのですが¹、10行程度のコードを書いては一喜一憂する材料としてはもってこいでしょう。

出力関数

先程のHow to Test with TestU01ですが、これはpermuted congruential generator (PCG) （Wikipedia）という生成器の作者のブログ記事です。2014年と比較的最近提唱された手法なので評価が定まっているのかはちょっとわからないのですが²、かいつまんでいうと伝統的な線形合同法（LCG, linear congruential generator）に「出力関数」をかませてやればBigCrushは通るよ、といったものです。
ここで出力関数とはどのようなものか見ていきます。

まず線形合同法では

x_{n+1} = A x_n + C\quad\mod M

という漸化式で数列$\{x_n\}$を作っているわけですが、$x_n$は「内部状態」でかつ「出力」になっています（64-bitの線形合同法で上32-bitのみを出力したりするのも一般的ですが、このことは一旦忘れて話を進めさせてください）。線形合同法では生の$x_n$を出力していて特に下位ビットの品質が悪かったので、それなら

y_n = f(x_n)

のように関数をかました$y_n$を使えば「疑似乱数の品質」も改善するだろう、というのが出力関数の考え方です。$f(x)$としてどんな関数を使うかというと、例えばuint64_t -> uint64_tな関数だとすると「可逆な（全単射になっている）」関数を使います。そもそも64-bitの線形合同法の作る数列は

\{0, 1, 2,...,2^{64}-1 \}

を並び替えてたものになっており、同じ値が複数回出てきたり出てこない値があったりしません³。関数$f$が可逆なら、数列$\{y_n\}$もこれを再び並び替えたものになっています。$f$の具体例は：

• 奇数定数での掛け算

uint64_t x;
x *= 0xdeadbeef;

オーバーフローのときが気にはなりますが、$2^{64}$を法とした整数では奇数には逆元がみつかります。いっぽう偶数を掛けてしまうと駄目なのはすぐにわかりますね。

• 定数での加算

 x += 0xdeadc0de;

これは同じ値を引けば元に戻ります。

• Shift & xor

x ^= x>>3;
x ^= x<<5;

Xorで汚染されなかったビットを再びshiftしてxorしていくことで元に戻すことができます。ビットベクトルに正則な行列を掛けていることにもなっています。
代数的構造がちゃんぽんにできるのがプログラミング言語の楽しいところですね。8-bitや16-bit整数のベクトルとみなしてシフトして掛けて足すなどの操作も、MMX/SSE2/AVX2で実装するなら面白いと思います。

Counter Basedな疑似乱数生成器

出力関数で線形合同法の品質を改善できるなら、もっともっと駄目な状態からでもBigCrushを通せないものかというのが人情です。例えば、

uint64_t stat(){
    static uint64_t x = 0;
    return x++;
}

は数列$\{0, 1, 2,...,2^{64}-1 \}$になりますが、どのみち内部状態が64-bitの疑似乱数はこれの並び替えなわけで⁴、出力関数で十分に並び替えを繰り返せば検定にも通りそうです。試してみようかと思っていたら先行研究が存在することがわかったので、実際に試しながら記事を書くことにしました。
このような方法はCounter-based random number generator (CBRNG)(en.wikipedia)といって2011年の論文ですからPCGよりも前ですね。"Parallel random numbers: as easy as 1, 2, 3"というおしゃれなタイトルになっています。筆頭著者はD. E. Shaw Researchに移ったJohn Salmonで重力多体問題の並列実装で有名な人です⁵。

書いてみたもの

とりあえずSmallCrushとCrushが通ってからこの記事を書き始めました。BigCrushは通っても通らなくても結果を晒すつもりでいましたが、書き終わる前に無事に通ってくれました。

uint32_t counter64(){
    static uint64_t status = 0;

    const uint64_t a = 6364136223846793005;
    const uint64_t c = 1442695040888963407;

    uint64_t x = status++;

    x = a*x + c;
    x ^= x>>32;
    x ^= x>>16;

    x = a*x + c;
    x ^= x>>8;
    x ^= x>>4;

    x = a*x + c;
    x ^= x>>2;
    x ^= x>>1;

    return x >> 32;
}

ヤマ感で適当に書いています。係数$a$と$c$は線形合同法から借りてきました。考える前にテストにかけてみることができるのもTestU01の作者のおかげです。

========= Summary results of BigCrush =========

 Version:          TestU01 1.2.3
 Generator:        Counter64
 Number of statistics:  160
 Total CPU time:   03:10:44.45

 All tests were passed

まとめ

線形合同法は疑似乱数生成法の中でもjumpが得意で、jumpの生成に状態更新の$\log n$倍程度、ジャンプの適用は状態更新と同じ計算でできてしまいます。しかしCBGはもうほぼほぼjumpという概念すら消え去ってしまっています。その代わりに出力関数は重めになりますが、並列計算機でスループットを出すということなのでしょう。

---

1. シングルコアで数時間流す程度のテストなので、周期が$2^{32}$では通らないものの$2^{36}$あれば通るもののようです。しかしこの程度の周期はシングルコアでも1分程度で食べ尽くせてしまいます。 ↩

2. 個人的な感想ですが「BigCrushを通すのを目標に作ればそりゃ通せるでしょ」という時代なのかなという感があります。 ↩

3. これは「乱数」としては非常に不自然なものですが、周期のある擬似乱数としては避けられません。その代りに周期を全部なめるのが現実的でないぐらいに大きく取ります。 ↩

4. xorshift等$0$は状態として含まず周期が$2^{64}-1$のものもあります。 ↩

5. 重力多体問題の人が分子動力学の研究所に行くというのはよくあることのように思われますが、そこで疑似乱数で論文を書いてしまうというのは凄いですね。 ↩