はじめに
駄目な疑似乱数として知られている線形合同法は次の式で値を更新します。
x_{n+1} = a x_n + c \mod m
$x_n$が疑似乱数列、$a$と$c$と$m$がその計算方法のパラメタ（定数）となります。
$m$に関しては計算機の語長に合わせる（つまり何も考えずに単に掛けたり足したりする）と$2^{32}$や$2^{64}$が自然な選択となります。メルセンヌ素数の$2^{31}-1$や$2^{61}-1$のこともあるようです。とりあえず一覧は英語版Wikipediaを見てください。
一度に複数回ジャンプ
以降$\mod m$を書くのは面倒なの全部剰余環$\mathbb Z / m \mathbb Z$上の計算ということにします1。
x_{n+2} = a (a x_n + c) + c = a^2 x_n + (a + 1)c
なので
\begin{align}
a^{(2)} =& a^2 \\
c^{(2)} =& (a+1)c
\end{align}
の係数を覚えておけば全く同じアルゴリズム$x_{n+2}=a^{(2)} x_n + c^{(2)}$で2つ飛び数列を作れます。再帰的に組み合わせたら3つでも4つでも作れそうですね2。
話を円滑に進めるために状態の更新操作をベクトル行列積のかたちで書いてみます。
\pmatrix{x_{n+1} \\ 1} = \pmatrix{a & c \\ 0 & 1} \pmatrix{x_n \\ 1}
縦ベクトルに突然$1$が現れることには驚くかもしれませんが、回転行列を拡大して並進も加えるときなどにも使われるテクニックです（参考：アフィン変換/平行移動だって変換行列で）。$2 \times 2$の上三角行列は$U$とでも呼ぶことにしましょう。
- 乗算：
\pmatrix{A & C \\ 0 & 1}
\pmatrix{a & c \\ 0 & 1}
=
\pmatrix{Aa & Ac + C \\ 0 & 1}
可換ではありませんが、行列$U$に対して$U^n U^k = U^k U^n = U^{n+k}$ぐらいは成り立っているので心配することはないでしょう。
- 冪乗： 乗算ができればとてもエレガントなアルゴリズムで冪乗を実装できます。
U^n = \begin{cases}
I & (n=0) \\
\left(U^{n/2}\right)^2 & (n \text{ is even}) \\
U U^{n-1} & (n \text{ is odd})
\end{cases}
ただし$I$は$a=1, c=0$の単位行列。
- 逆元： 状態を一歩手前に戻す操作はどうやったら求まるでしょうか。$m$が２冪で係数が適切なとき周期も$m$になります（適切でなくても周期は$m$の約数になるので$m$乗したら単位元に戻ります）。
U^{-1} = U^{m-1}
実装なのですが、this->pow(-1)だけで済んでしまいます。$-1$がuint_64では勝手に$2^{64}-1$になってくれるからで3、見た目が直感的でいいですね。
実装
というわけこんなコードになりました。係数の乗算則、冪乗、逆元、デフォルト値を実装してあります。
template<typename UINT>
struct RandCoeff{
UINT a;
UINT c;
static RandCoeff default_value();
RandCoeff operator*(const RandCoeff &rhs) const {
return {a * rhs.a, a * rhs.c + c};
}
RandCoeff pow(const UINT n) const {
if (0 == n){
return RandCoeff{1, 0};
} else if (0 == n%2){
return ((*this) * (*this)).pow(n/2);
} else {
return (*this) * this->pow(n-1);
}
}
RandCoeff inv() const {
return this->pow(-1);
}
};
template<>
RandCoeff<uint32_t> RandCoeff<uint32_t>::default_value(){
return {134775813, 1};
}
template<>
RandCoeff<uint64_t> RandCoeff<uint64_t>::default_value(){
return {6364136223846793005, 1};
}
inv()の結果を検算：
$ bc < "obase=16; 6364136223846793005 * 13877824140714322085"
4271E21C27CCEC3E0000000000000001
逆元がちゃんと求まってます。
状態の方はこんなコードです：
template<typename UINT>
struct RandStat{
UINT val;
RandStat next(const RandCoeff<UINT> &ac) const {
return {ac.a * val + ac.c};
}
void update(const RandCoeff<UINT> &ac){
(*this) = this->next(ac);
}
auto operator()() const {
return val;
}
void init(const UINT seed, const RandCoeff<UINT> &ac){
val = seed;
this->update(ac);
}
double to_double() const {
return double(val) * (1.0 / (1.0 + std::numeric_limits<UINT>::max())) ;
}
operator double() const {
return this->to_double();
}
};
template<>
double RandStat<uint64_t>::to_double() const{
return val * 5.4210108624275216e-20;
}
グローバル変数に状態を置いてAPIが呼ばれると状態の値を返しながらグローバルの状態をこっそり更新するという設計にはしていません。殆どの疑似乱数APIってここが駄目だと思うんです。いや「疑似」であることを忘れたいならrand()を呼べば乱数が返ってくるでいいのかもしれませんが。ひとつのAPI呼び出しに対して2つの動作があって、しかもその片方は副作用もあるとかやっぱり駄目でしょう4。
メイン関数の例は
int main(){
enum { seed = 42 };
RandStat<uint64_t> s;
RandCoeff<uint64_t> ac = ac.default_value();
s.init(seed, ac);
for(int i=0; i<51; i++){
auto u = s();
double d(s);
printf("%2d : %016llx : %20.16f %a\n", i, u, d, d);
s.update(ac);
}
puts("");
s.init(seed, ac);
auto ac5 = ac.pow(5);
for(int i=0; i<51; i+=5){
auto u = s();
double d(s);
printf("%2d : %016llx : %20.16f %a\n", i, u, d, d);
s.update(ac5);
}
}
$10 \times 5$の行列を乱数で埋めて、1列目だけ縦読みで読み出すような操作が無駄なくできるという例になっています。
何故このようなものを書いたのか
ここまで書いたのは、疑似乱数生成器がオプションとして提供できる機能として「ジャンプ」というものがあること、あと線形合同法ではとても簡単にこれが実装できるということです。線形合同法に限らず、漸化式が有限体上の行列ベクトル積で書けるものになっていれば行列行列積をひとつ実装すればジャンプは実現できそうですね。
それで、ジャンプが有り難くなるのはマルチスレッドやマルチプロセスで並列化するときです。1000万の要素を疑似乱数で埋める作業を10人で分担するとしましょう。0番目から9番目までの人が、100万$\times$自分の番号だけジャンプしてから100万個書き込めば並列化前とまったく同じ結果を保証できます。スレッド並列でこれをやる場合、係数や状態がグローバルに格納されているのは論外です。
また、ここまで厳密に並列化前との一致が必要でない場合、とりあえず全員SEEDだけ別の値を使っておけばという雑な考え方もあるかもしれません（MPIのようなプロセス並列の場合グローバル変数の状態もプロセス数分あるわけですし）。この方法の少し怖いところは、数列$x_n$の中であるプロセスのいる位置と別のプロセスのいる位置が十分離れている保証がなく、そもそも距離をどう測ればいいのかもわかりません5。（$2^{64}$は周期としてちょっと短く思えてきますが）ジャンプを使ってプロセスごとの距離を$2^{48}$ぐらい空けておくといった方法で衝突を回避できます。
応用
16進数ダンプの下の桁を見ると周期性がありそうです。
0 : 7d720f6e9086dd63 : 0.4900216717656137 0x1.f5c83dba421b6p-2
1 : 7800fec18e280768 : 0.4687651846453776 0x1.e003fb0638a01p-2
2 : 8f6ada53083de549 : 0.5602241947856252 0x1.1ed5b4a6107bcp-1
3 : 06619f9c471d84d6 : 0.0249271160638412 0x1.9867e711c7612p-6
4 : 12ec0efbfda4839f : 0.0739144673549550 0x1.2ec0efbfda483p-4
5 : a792cbe702c203f4 : 0.6545836867324012 0x1.4f2597ce0583fp-1
6 : 124ee7f67914bde5 : 0.0715165116062790 0x1.24ee7f67914bdp-4
7 : 667153181522fc42 : 0.4001666959567379 0x1.99c54c60548bep-2
8 : 4de67d927cb5159b : 0.3042982561716665 0x1.3799f649f2d44p-2
9 : 3f8f9eac58c3b140 : 0.2482852144974135 0x1.fc7cf562c61d8p-3
10 : 25374f19d994e841 : 0.1453751981840302 0x1.29ba78cecca73p-3
11 : b2b65a87983a126e : 0.6980949955979910 0x1.656cb50f30741p-1
12 : 93a81af5f35fcf57 : 0.5767838335481194 0x1.275035ebe6bf9p-1
13 : ab513794d1569b4c : 0.6692080248188523 0x1.56a26f29a2ad2p-1
14 : 8ae6e2e1b680005d : 0.5425855446869717 0x1.15cdc5c36cfffp-1
15 : 2ef6c82b66cf335a : 0.1834530931435523 0x1.77b6415b36799p-3
16 : d14ee04b7847acd3 : 0.8176098045589969 0x1.a29dc096f08f5p-1
17 : 8bdaa40cee250e18 : 0.5463049441879272 0x1.17b54819dc4a1p-1
18 : 4295366292796239 : 0.2600893011717509 0x1.0a54d98a49e58p-2
19 : c20e57d20e3d8b06 : 0.7580313575583570 0x1.841cafa41c7bp-1
20 : e9cc5c22bb476a0f : 0.9132745346826403 0x1.d398b845768ecp-1
21 : 1dd491b387e615a4 : 0.1165247977550268 0x1.dd491b387e615p-4
22 : 3572322550a229d5 : 0.2087737408734333 0x1.ab91912a85114p-3
23 : 8060c4ce343b0572 : 0.5014765742709919 0x1.00c1899c6876p-1
24 : aecd6d2a796d830b : 0.6828220585534178 0x1.5d9ada54f2dafp-1
25 : c0aaf74835a97df0 : 0.7526087333004962 0x1.8155ee906b52fp-1
26 : 5f73ba8a10f53331 : 0.3728596293521675 0x1.7dceea2843d4cp-2
27 : 7aaca634f7044e9e : 0.4791969184216366 0x1.eab298d3dc113p-2
28 : 63a4f47938b833c7 : 0.3892357631636691 0x1.8e93d1e4e2e0cp-2
29 : f45d950290e3d2fc : 0.9545529490723218 0x1.e8bb2a0521c79p-1
30 : ddc1b4dd1a631a4d : 0.8662369766046931 0x1.bb8369ba34c62p-1
31 : df2645faf848d28a : 0.8716777551383387 0x1.be4c8bf5f0919p-1
32 : 8e5c796f47917843 : 0.5560985466650601 0x1.1cb8f2de8f22ep-1
3896f452b01e7cda...
実際に$U^{(16)}$の値をみてみましょう。
ac.pow(16).dump();
すみません、dump()メソッドはさっき追加したので上の方のコードには含まれていません。
結果は16進数で：
8d5e2ddc895abe41, fd8341fcddebfcb0
一番下の桁を見ると$a^{(16)}=1, c^{(16)}=0$なわけですから、16回のジャンプが恒等写像になっています。
$2^{4n}$回の係数は：
0 : 1, 0
4 : 8d5e2ddc895abe41, fd8341fcddebfcb0
8 : 2078e0dd6db6401, 41f7c92a64676b00
12 : 469c6146fd364001, 5eeb75973616b000
16 : 902da3ff53640001, 39f376e3016b0000
20 : 7cc9cf7536400001, e61cddd016b00000
24 : bc2c775364000001, c73d7d016b000000
28 : 5247753640000001, e377d016b0000000
32 : a477536400000001, d77d016b00000000
36 : 4775364000000001, 77d016b000000000
40 : 7753640000000001, 7d016b0000000000
44 : 7536400000000001, d016b00000000000
48 : 5364000000000001, 16b000000000000
52 : 3640000000000001, 16b0000000000000
56 : 6400000000000001, 6b00000000000000
60 : 4000000000000001, b000000000000000
64 : 1, 0
綺麗ですね。つまり、上位1-bitをマスクするたびにその下は半分の周期を持っていることがわかります。法でもあり周期でもある$m$が半分になるので納得の行く挙動です（注意して使いましょう）。
筆算
2進数で奇数の二乗を筆算してみます。($x, y, z \in {0, 1}$)
\begin{array}{cccc}
& & x & 1 \\
& \times & x & 1 \\
\hline
& & x & 1 \\
& x^2 & x \\
\hline
y & z & 0 & 1
\end{array}
\qquad
\begin{array}{ccccc}
& & x & 0 & 1 \\
& \times & x & 0 & 1 \\
\hline
& & x & 0 & 1 \\
x^2 & 0 & x & & \\
\hline
x^2 & y & 0 & 0 & 1 \\
\end{array}
最下位の$1$を保存しながら$0$が上へと伸びていくのがわかります。
追記
uint32_t lcg64(){
static uint64_t x = 42;
const uint64_t a = 6364136223846793005;
const uint64_t c = 1;
x = a*x + c;
uint32_t y = x >> 32;
// return y;
y ^= y>>1;
y ^= y>>2;
y ^= y>>4;
y ^= y>>8;
y ^= y>>16;
return y;
}
こんな糞コードでもBigCrush通りました、わ〜い。
========= Summary results of BigCrush =========
Version: TestU01 1.2.3
Generator: LCG64
Number of statistics: 160
Total CPU time: 03:19:00.79
All tests were passed
実際は、LCG64をxorでかき混ぜたものです。
-
$m$が素数のときのみ有限体$\mathbb F_m$になる（0以外のすべての要素に逆元がみつかる）のでした。 ↩
-
ごめんなさい、今更一般型の存在に気が付きました。$a^{(n)}=a^n$と$c^{(n)} = c \left(\sum_0^{n-1} a^n \right)$になります。 ↩
-
言語仕様で保証されているかはわかりかねますが普通の計算機ならまず大丈夫でしょう。 ↩
-
ただ関数呼び出しの引数に状態と係数が入ってしまうと際限なく面倒くさいです。グローバル変数に
thread_local指定子を付けて置いておくぐらいが折衷案でしょうか。 ↩ -
「有限体における冪は計算量的にコストが掛からないのに対し、冪の逆である離散対数は計算量的にコストが掛かる」だそうです ↩