擬似乱数生成
\( \def\vector#1{\boldsymbol{#1}} \) \( \newcommand{\Z}{{\mathcal{Z}}} \)

論文翻訳: Xorshift RNGs

Takami Torao 2003年の論文 #Xorshift #PRNG
  • このエントリーをはてなブックマークに追加

George Marsaglia∗
The Florida State University

概要

\(k=32,64,96,128,160,192\) に対して \(2^k-1\) の周期を持つ、単純で非常に高速な乱数ジェネレータ (RNGs) の一種を説明する。これらはランダム性テストに非常によく適合しているようである。

Table of Contents

  1. 概要
  2. 1 導入
  3. 2 理論
    1. 2.1 すべての非ゼロバイナリベクトルを生成する行列 \(T\)
  4. 3 Xorshift RNG への適用
    1. 3.1 \(n=96,128,160,\ldots\) 次元のバイナリベクトル空間
  5. 4 概要
  6. References
  7. 翻訳抄
  • ∗While the author is now Professor Emeritus, portions of the research for this article were done under by grants from The National Science Foundation.

1 導入

xorshift 乱数生成器 (xorshift RNG) は単純な計算構成である。コンピュータワードのシフト演算と排他的論理和 (xor) を繰り返し適用することによって \(2^{32}-1\) 個の整数シーケンス \(x\)、\(2^{64}-1\) 個のペアシーケンス \(x,y\)、\(2^{96}-1\) 個のトリプルシーケンス \(x,y,z\) などを生成する。C 言語での基本的な操作は、左シフトの場合 y^(y<<a)、右シフトの場合 y^(y>>a) である (Fortran では単一の形式 ieor(y,ishft(y,a)) だけで、右シフトの場合は負の a、左シフトの場合は非負で期待した結果が得られる)。

このような xorshift 操作を様々なシフト数や引数と組み合わせると、ランダム性テストで非常に良好となる非常に高速で単純な RNG を得ることができる。xorshift 操作の能力と有効性を理解するために、4 つの乱数シード x,y,z,w が与えられたとき、呼び出しごとに 3 つの xorshift 操作で \(2^{128}-1\) 個のランダムな 32 ビット整数を提供する C プロシジャの重要な部分を次に示す:

tmp = (x ˆ (x << 15));
x = y;
y = z;
z = w;
return w = (w ˆ (w >> 21)) ˆ (tmp ˆ (tmp >> 4));

このような手順は非常に高速で一般的に秒間 2 億回以上であり、結果のランダム整数はすべてのランダム性テスト、特に [3] の "タフ" テストと Diehard Battery [2] の新しいバージョンに合格している。

\(2^{160}\) までの期間を持つ様々な xorshift RNG を提供する選択肢を確立するための背景理論をここに示す。例えばキャリー付き乗算 RNG がもたらすようなようなより長い周期を使用することもできるが、それらは整数の乗算と最新値を保存する (場合によっては大きくなる) テーブルを保持する必要がある。より詳細な比較は最後の要約セクションに記載している。

2 理論

ほとんどの RNG における数学的モデルは次の形式で記述できる: \(m\)-タプル \((x_1,x_2,\ldots,x_m)\) で構成されるシード集合 \(\mathcal{Z}\) と、\(\mathcal{Z}\) 上の一対一の関数 \(f(\ )\) とする。一般的には \(\mathcal{Z}\) は単なる整数の集合だが、より良い RNG ではペア、トリプルなどの集合である。\(\mathcal{Z}\) から一様なランダムさで \(z\) を選択したとき、RNG の出力はシーケンス \(f(z),f^2(z),f^3(z),\ldots\) となる。ここで \(f^2(z)\) は \(f(f(z))\) を意味している。\(f\) は \(\mathcal{Z}\) 上で一対一であるため、乱数 \(f(z)\) 自身は \(\mathcal{Z}\) 上で均一であり、\(f^2(z)\) も同様である。実際、シーケンス \(f(z),f^2(z),\ldots\) のそれぞれの要素はシード集合 \(\mathcal{Z}\) 上で一様に分布するが、それらは独立ではない。

xorshift RNG では、シード集合 \(\mathcal{Z}\) はゼロベクトルを除く \(1\times n\) 個のバイナリベクトル \(\beta=(b_1,b_2, \ldots,b_n)\) の集合である。通常 \(n\) は 32, 64, 96 などであり 32 ビットのコンピュータワードに隣接させて要素を構成する。\(\Z\) 内のベクトル \(\beta\) の要素は体 \(\{0,1\}\) にあるため、構成要素の 32 ビットパートを xor することでバイナリベクトルの追加を実装できる。我々の xorshift RNG では \(\Z\) 上の可逆関数が必要である。そのために我々は非特異 (正則) \(n\times n\) バイナリ行列 \(T\) で特徴付けられるバイナリベクトル空間上の線形変換を使用する。もし \(\beta\) が \(\Z\) から一様にランダムに選択されている (つまりシード) なら、シーケンス \(\beta T,\beta T^2, \beta T^3, \ldots\) の各要素も \(\Z\) 上で一様に分布しているため、\(\Z\) からの一様な要素である ID (identically distributed) のシーケンスを持つが、それらはIID (independent identically distributed) ではない。しかし多くの RNG のケースと同様に、しばしば ID 要素の関数は IID シーケンスの要素の同じ関数の分布に非常に近い分布を持つことがある。これは過去 50 年間のコンピュータでの有用性を正当化する、シード集合 \(\Z\) 上の関数 \(f(\ )\) の特定の選択肢の注目すべき特性である。

2.1 すべての非ゼロバイナリベクトルを生成する行列 \(T\)

最初に [1] で与えられた主旨は以下の通り:

定理. 非特異 \(n \times n\) バイナリ行列 \(T\) が非ゼロ初期 \(1 \times n\) バイナリベクトル \(\beta\) に対してシーケンス \(\beta T, \beta T^2, \beta T^3, \ldots\) ですべての非ゼロ \(1 \times n\) バイナリベクトルを生成するには、非特異 \(n \times n\) バイナリ行列の群において \(T\) の次数が \(2^n-1\) であることが必要かつ十分である。

証明. まず必要条件: もし \(\beta T, \beta T^2, \ldots\) の周期が \(k=2^n-1\) であれば各 \(1 \times n\) バイナリベクトル \(\beta\) に対して \(\beta T^k = \beta\) であり、行列 \(T^k+I\) の null 空間は全空間であるため \(T^k+I\) はゼロ行列、\(T^k=I\) でなければならない。もしある \(j \lt k\) に対して \(T^j=I\) ならば \(\beta T, \beta T^2, \ldots\) は \(2^n-1\) 未満となるだろう。

次に十分条件: \(T\) の次数が \(k=2^n-1\) の場合、行列 \(T,T^2,T^3,\ldots,T^k\) は非特異で明瞭 (district) であり、\(T\) の固有多項式 (characteristic polynomial) とユークリッドアルゴリズムを通じて次数 \(\lt n\) の \(T\) の多項式として表現される。次数 \(\lt n\) の \(T\) には \(k=2^n-1\) の非 null 多項式があるため、それらは何らかの順序で明確な非特異行列 \(T,T^2,\ldots,T^k\) でなければならない。特に \(T\) の多項式が特異行列である場合、\(T\) の固有多項式を介してゼロ行列に縮退しなければならない。これは、ある非 null \(\beta\) と \(j \lt k\) に対して \(\beta T^j=\beta\) が \(T^j+I\) が特異であることを意味するため、\(\beta T, \beta T^2, \ldots\) の周期は \(k=2^n-1\) でなければならない。

3 Xorshift RNG への適用

バイナリベクトル \(y\) に対して、コンピュータでの演算 \(yT\) は行列 \(T\) が特別な形式でない限り高コストとなる可能性が高い。\(L\) がバイナリベクトル \(y\) 上のある位置の左シフトに影響する \(n \times n\) バイナリ行列である場合、つまり \(L\) が主要な複対角線 (principal subdiagonal)上が 1 であることを除いてすべて 0 であり \(T=I+L^a\) である場合、C 言語での xorshift 操作 y^(y<<a) は線形変換 \(yT\): 加算、mod 2、バイナリベクトル \(y\) 自身への左シフトバージョンの線形変換を生成する。同様に \(R\) が右シフト 1 行列 (\(L\) の転置) であれば xorshift 演算 y^(y>>b) は \(T=I+R^b\) で \(yT\) を形成できる。

\(n \times n\) 行列 \(I + L^a\) と \(I + R^b\) は非特異である。これは、例えば \(L^n=0\) であり、したがって有限級数 \(I+L^a+L^{2a}+L^{3a}+\ldots\) が \((I+L^a)\) の逆行列となるためである。したがって位数 \(2^n-1\) を持つ行列の明らかな候補は \(T \equiv (I+L^a)(I+R^b)\) または \(T \equiv (I+R^b)(I+L^a)\) である。残念ながら \(n\) が 32 または 64 の場合、\(a\) および \(b\) の選択肢はそのような \(T\) に必要な次数を提供しない。(候補の事前検査は \(T\) を \(n\) 回二乗することである。結果が \(T\) でない場合 \(T\) は次数 \(2n-1\) ではない。バイナリ行列の行をコンピュータワードとして表し、(左) 乗算バイナリベクトルが 1 であるすべtネオ行の xor 演算を行うことにより、非常に高速なバイナリ行列積を単純な C または Fotrran プログラムで表することができ、\(T=(I+L^a)(I+R^b)\) となる可能性のある全ての \(a\), \(b\) を確認するのに数分しかかからない。)

ただし \(N=32\) または \(n=64\) となるような次数 \(2^n-1\) を持つ行列 \(T=(I+L^a)(I+R^b)(I+L^c)\) の \(a\), \(b\), \(c\) には多くの選択肢がある。次数 \(2^{32}-1\) を持つ \(32 \times 32\) のバイナリ行列 \(T=(I+L^a)(I+R^b)(I+L^c)\)に対して 81 個のトリプル \((a,b,c), a \lt c\) が存在する。

| 1, 3,10| 1, 5,16| 1, 5,19| 1, 9,29| 1,11, 6| 1,11,16| 1,19, 3| 1,21,20| 1,27,27|
| 2, 5,15| 2, 5,21| 2, 7, 7| 2, 7, 9| 2, 7,25| 2, 9,15| 2,15,17| 2,15,25| 2,21, 9|
| 3, 1,14| 3, 3,26| 3, 3,28| 3, 3,29| 3, 5,20| 3, 5,22| 3, 5,25| 3, 7,29| 3,13, 7|
| 3,23,25| 3,25,24| 3,27,11| 4, 3,17| 4, 3,27| 4, 5,15| 5, 3,21| 5, 7,22| 5, 9,7 |
| 5, 9,28| 5, 9,31| 5,13, 6| 5,15,17| 5,17,13| 5,21,12| 5,27, 8| 5,27,21| 5,27,25|
| 5,27,28| 6, 1,11| 6, 3,17| 6,17, 9| 6,21, 7| 6,21,13| 7, 1, 9| 7, 1,18| 7, 1,25|
| 7,13,25| 7,17,21| 7,25,12| 7,25,20| 8, 7,23| 8,9,23 | 9, 5,1 | 9, 5,25| 9,11,19|
| 9,21,16|10, 9,21|10, 9,25|11, 7,12|11, 7,16|11,17,13|11,21,13|12, 9,23|13, 3,17|
|13, 3,27|13, 5,19|13,17,15|14, 1,15|14,13,15|15, 1,29|17,15,20|17,15,23|17,15,26|

\(a \lt c\) となる 81 個のトリプルのうちトリプル \((c,b,a)\) も全周期 (full period) \(T\) を提供し、合計で 162 となる。さらに、転置を行うことで 162 個の全周期 \(T\) が生じる: トリプル \((a,b,c)\) を使用して \(T=(I+R^c)(I+L^b)(I+R^c)\) を形成し合計 324 である。そして最後に \((I+L^a)(I+R^b)(I+L^c)\) または \((I+R^a)(I+L^b)(I+R^c)\) 形式の 324 個の \(T\) それぞれに対応する行列 \((I+L^a)(I+L^c)(I+R^b)\) および \((I+R^a)(I+R^c)(I+L^b)\) もまた \(2^{32}-1\) の周期を有することが分かるため合計 648 個の選択肢がある。

要約すると、上記の 81 個のトリプル \(a\), \(b\), \(c\) それぞれが \(a \lt c\) である場合、以下に C で示す 8 行のいずれか一つは周期 \(2^{32}-1\) を持つ 32 ビットの xorshift RNG の基礎として機能する。

yˆ=y<<a; yˆ=y>>b; yˆ=y<<c;
yˆ=y<<c; yˆ=y>>b; yˆ=y<<a;
yˆ=y>>a; yˆ=y<<b; yˆ=y>>c;
yˆ=y>>c; yˆ=y<<b; yˆ=y>>a;
yˆ=y<<a; yˆ=y<<c; yˆ=y>>b;
yˆ=y<<c; yˆ=y<<a; yˆ=y>>b;
yˆ=y>>a; yˆ=y>>c; yˆ=y<<b;
yˆ=y>>c; yˆ=y>>a; yˆ=y<<b;

64 ビット整数の場合、次数 \(2^{64}-1\) となる \(64 \times 64\) 行列 \(T=(I+L^a)(I+R^b)(I+L^c)\) は次の 275 トリプルによってもたらされる:

| 1, 1,54| 1, 1,55| 1, 3,45| 1, 7, 9| 1, 7,44| 1, 7,46| 1, 9,50| 1,11,35| 1,11,50|
| 1,13,45| 1,15, 4| 1,15,63| 1,19, 6| 1,19,16| 1,23,14| 1,23,29| 1,29,34| 1,35, 5|
| 1,35,11| 1,35,34| 1,45,37| 1,51,13| 1,53, 3| 1,59,14| 2,13,23| 2,31,51| 2,31,53|
| 2,43,27| 2,47,49| 3, 1,11| 3, 5,21| 3,13,59| 3,21,31| 3,25,20| 3,25,31| 3,25,56|
| 3,29,40| 3,29,47| 3,29,49| 3,35,14| 3,37,17| 3,43, 4| 3,43, 6| 3,43,11| 3,51,16|
| 3,53, 7| 3,61,17| 3,61,26| 4, 7,19| 4, 9,13| 4,15,51| 4,15,53| 4,29,45| 4,29,49|
| 4,31,33| 4,35,15| 4,35,21| 4,37,11| 4,37,21| 4,41,19| 4,41,45| 4,43,21| 4,43,31|
| 4,53, 7| 5, 9,23| 5,11,54| 5,15,27| 5,17,11| 5,23,36| 5,33,29| 5,41,20| 5,45,16|
| 5,47,23| 5,53,20| 5,59,33| 5,59,35| 5,59,63| 6, 1,17| 6, 3,49| 6,17,47| 6,23,27|
| 6,27, 7| 6,43,21| 6,49,29| 6,55,17| 7, 5,41| 7, 5,47| 7, 5,55| 7, 7,20| 7, 9,38|
| 7,11,10| 7,11,35| 7,13,58| 7,19,17| 7,19,54| 7,23, 8| 7,25,58| 7,27,59| 7,33, 8|
| 7,41,40| 7,43,28| 7,51,24| 7,57,12| 8, 5,59| 8, 9,25| 8,13,25| 8,13,61| 8,15,21|
| 8,25,59| 8,29,19| 8,31,17| 8,37,21| 8,51,21| 9, 1,27| 9, 5,36| 9, 5,43| 9, 7,18|
| 9,19,18| 9,21,11| 9,21,20| 9,21,40| 9,23,57| 9,27,10| 9,29,12| 9,29,37| 9,37,31|
| 9,41,45|10, 7,33|10,27,59|10,53,13|11, 5,32|11, 5,34|11, 5,43|11, 5,45|11, 9,14|
|11, 9,34|11,13,40|11,15,37|11,23,42|11,23,56|11,25,48|11,27,26|11,29,14|11,31,18|
|11,53,23|12, 1,31|12, 3,13|12, 3,49|12, 7,13|12,11,47|12,25,27|12,39,49|12,43,19|
|13, 3,40|13, 3,53|13, 7,17|13, 9,15|13, 9,50|13,13,19|13,17,43|13,19,28|13,19,47|
|13,21,18|13,21,49|13,29,35|13,35,30|13,35,38|13,47,23|13,51,21|14,13,17|14,15,19|
|14,23,33|14,31,45|14,47,15|15, 1,19|15, 5,37|15,13,28|15,13,52|15,17,27|15,19,63|
|15,21,46|15,23,23|15,45,17|15,47,16|15,49,26|16, 5,17|16, 7,39|16,11,19|16,11,27|
|16,13,55|16,21,35|16,25,43|16,27,53|16,47,17|17,15,58|17,23,29|17,23,51|17,23,52|
|17,27,22|17,45,22|17,47,28|17,47,29|17,47,54|18, 1,25|18, 3,43|18,19,19|18,25,21|
|18,41,23|19, 7,36|19, 7,55|19,13,37|19,15,46|19,21,52|19,25,20|19,41,21|19,43,27|
|20, 1,31|20, 5,29|21, 1,27|21, 9,29|21,13,52|21,15,28|21,15,29|21,17,24|21,17,30|
|21,17,48|21,21,32|21,21,34|21,21,37|21,21,38|21,21,40|21,21,41|21,21,43|21,41,23|
|22, 3,39|23, 9,38|23, 9,48|23, 9,57|23,13,38|23,13,58|23,13,61|23,17,25|23,17,54|
|23,17,56|23,17,62|23,41,34|23,41,51|24, 9,35|24,11,29|24,25,25|24,31,35|25, 7,46|
|25, 7,49|25, 9,39|25,11,57|25,13,29|25,13,39|25,13,62|25,15,47|25,21,44|25,27,27|
|25,27,53|25,33,36|25,39,54|28, 9,55|28,11,53|29,27,37|31, 1,51|31,25,37|31,27,35|
|33,31,43|33,31,55|43,21,46|49,15,61|55, 9,56|

32 ビットの場合と同様に、64 ビットのシーケンスの場合は 275 個の \(a\), \(b\), \(c\) のいずれかを選択し、上記の 8 行の C のコードのいずれかを選択すると、64 ビットワードに対して合計 \(8 \times 275 = 2200\) 個の選択肢に対して周期 \(2^{64}-1\) となる xorshift RNG を提供する。以下は 32 ビットのシード値 \(y\) を取る基本的な 32 ビット xorshift の C プロシジャである。

unsigned long xor(){
  static unsigned long y = 2463534242;
  y ˆ= (y << 13);
  y = (y >> 17);
  return (y ˆ= (y << 5));
}

個人的に気に入っている選択肢の一つである \([a,b,c]=[13,17,5]\) を使用して Diehard [2] のバイナリランクテストを除くほとんどすべてのランダム性テストに合格するだろう。(長周期の xorshift RNG は必然的に非特異行列変換を使用するため、連続する \(n\) 個のベクトルがすべて線形独立でなければならないが、真にランダムなバイナリベクトルは時間の 30% だけ線形独立である。) 我々はこれらのうちいくつかをテストしただけだが、上記の 648 個の選択肢のいずれかが非常に速くシンプルで高品質な RNG を提供する可能性がある。

64 ビット整数を持つ C コンパイラの場合、64 ビットシード \(x\) が与えられると次の式は周期 \(2^{64}-1\) の優れた RNG を提供するだろう:

unsigned long long xor64(){
  static unsigned long long x=88172645463325252LL;
  x ˆ= (x << 13);
  x ˆ= (x >> 7);
  return (x ˆ= (x << 17));
}

ただし上記の 2200 個の選択肢のいずれかでも同様である。

3.1 \(n=96,128,160,\ldots\) 次元のバイナリベクトル空間

32 次元のベクトル空間要素を表すために 32 ビットコンピュータワードを使用したり、64 ビット整数を許容するコンパイラで 64 次元を使用することは便利だが、より長い xorshift 周期を得るためにはより高い次元のベクトル空間要素を表す方法が必要である。これを行う良い方法は、例えば 32 ビット成分 \((x,y,z)\) で構成される \(1 \times 96\) ベクトルや 32 ビット成分 \((x,y,z,w)\) で構成される \(1 \times 128\) ベクトルなどを可能にすることである。

自然な選択肢は \(T\) をブロック形式のコンパニオン行列にすることである。つまり \(n=64,96,128\) に対して \[ T = \left( \begin{array}{cc} 0 & A \\ I & B \end{array} \right), \ \ \ T = \left( \begin{array}{ccc} 0 & 0 & A \\ I & 0 & C \\ 0 & I & B \end{array} \right), \ \ \ T = \left( \begin{array}{cccc} 0 & 0 & 0 & A \\ I & 0 & 0 & C \\ 0 & I & 0 & D \\ 0 & 0 & I & B \end{array} \right) \ \ \ \]

次に、例えば \((x,y,z)T = (y,z,xA+yC+zB)\) とすると、\(32 \times 32\) の行列 \(A,B,C\) を求めることができるため、32 ビットの演算 \(xA\), \(yC\), \(zB\) は容易であり、\(T\) は \(96 \times 96\) の非特異バイナリ行列の群で次数 \(2^{92}-1\) を持つ。(ブロックの最後の列を \(A,B\)、\(A,C,B\) や \(A,C,D,B\) としたのは、全周期の選択肢 \(2^{64}-1\), \(2^{96}-1\), \(2^{128}-1\) があることが判明したからです。\(A=(I+L^a)(I+R^b)\) や \(B=(I+R^c)\) を選択するだけで他のブロックはすべてゼロの行列をブロックする。) したがってトリプル \([a,b,c]\) の適切な選択のためにこれらの行列はすべてそれぞれ \(2^{64}-1\), \(2^{96}-1\), \(2^{128}-1\) という必要な次数を持つ: \[ \left( \begin{array}{cc} 0 & (I+L^a)(I+R^b) \\ I & (I+R^c) \end{array} \right), \ \ \ \left( \begin{array}{ccc} 0 & 0 & (I+L^a)(I+R^b) \\ I & 0 & 0 \\ 0 & I & (I+R^c) \end{array} \right), \ \ \ \left( \begin{array}{cccc} 0 & 0 & 0 & (I+L^a)(I+R^b) \\ I & 0 & 0 & 0 \\ 0 & I & 0 & 0 \\ 0 & 0 & I & (I+R^c) \end{array} \right) \] トリプル \([a,b,c]\) のいくつか (4つ) の選択肢は: \(n=64\) に対して [10, 13, 10], [8, 9, 22], [2, 7, 3], [23, 3, 24]; \(n=96\) に対して [10, 5, 26], [13, 19, 3], [1, 17, 2], [10, 1, 26]; \(n=128\) に対して [5, 14, 1], [15, 4, 21], [23, 24, 3], [5, 12, 29] である。いずれの場合もブロック行列 \(T\) の次数は \(2^n-1\) であり、必要な長さのシーケンスを生成するための C プロシジャのエッセンスは static unsigned long の x, y, z, w とテンポラリの t である:

// 2^64-1 周期
t = (x ˆ (x << a));
x = y;
return y = (y ˆ (y >> c)) ˆ (t ˆ (t >> b));
// 2^96−1 周期
t = (x ˆ (x << a));
x = y;
y = z;
return z = (z ˆ (z >> c)) ˆ (t ˆ (t >> b));
// 2^128−1 周期
t = (x ˆ (x << a));
x = y;
y = z;
z = w;
return w = (w ˆ (w >> c)) ˆ (t ˆ (t >> b));

これらの例はパラメータ [a,b,c] = [2,1,4], [7,13,6], [1,1,20] を選択して同じプロモーションスキームを使用して \(2^{160}-1\) 周期の xorshift RNG を提供するが (static unsigned long の) シード x, y, z, w, v を必要とする。

t = (x ˆ (x >> a));
x = y;
y = z;
z = w;
w = v;
return v = (v ˆ (v >> c)) ˆ (t ˆ (t >> b));

必要に応じてブロックコンパニオン行列の最後の列を \(I+L^a\), \(I+R^b\), \(I+L^c\), \(I+R^d\) に選択することによって長い周期を見つけることもできる。この C の生成コードは:

// 2^96−1 周期
t = (x ˆ (x << 3)) ˆ (y ˆ (y >> 19)) ˆ (z ˆ (z << 6));
x = y;
y = z;
return (z = t);
// 2^128−1 周期
t = (x ˆ (x << 20)) ˆ (y ˆ (y >> 11)) ˆ (z ˆ (z << 27)) ˆ (w ˆ (w >> 6));
x = y;
y = z;
z = w;
return (w = t);

上記すべてのプロシジャは 32 ビット整数を返すが、プロシジャは \(n=64,96,128\) 次元のバイナリベクトル空間上に最大周期 \(2^n-1\) のペア \((x,y)\)、トリプル \((x,y,z)\)、クアドラプル \((z,y,z,w)\) のシーケンスを生成する。上記パラメータのすべての選択に対する xorshift RNG の結果は Diehard のすべてのテストに合格し、非常に高速で全てが 4～6 ナノ秒、つまり毎秒約 2 億の乱数を取得できる。これらは非常に高速であるため支配的な部分はリンケージ、つまりサブルーチン呼び出しで使用されるレジスターの保存前と復元後である。

上の例では最後の値 \((x,y)\) の \(y\)、\((x,y,z)\) の \(z\)、\((x,y,z,w)\) の \(w\) を返値としているが必ずしもそうである必要はない事に注意。単に代入するだけであればペア、トリプル、クアドラプルの完全なシーケンスは引き続き生成されるが、それらの値の関数を返すことができる。これによって興味深い RNG に様々な可能性がもたらされる。例えば、単に 69069*x を返すだけで \(2^n-1\) 周期の合同 RNG のようなものが得られ、一方 \(2^{32}(2^n-1)\) 周期は正規の xorshift と、我々が Weyl シーケンスと呼ぶものを組み合わせる (+ または xor) ことによって得られる: d+=362437; (362437 を任意の奇数定数に置き換えた) \(2^{32}\) 周期。ここに例を示す。

unsigned long xorwow(){
  static unsigned long x = 123456789, y = 362436069, z = 521288629,
                       w =  88675123, v =   5783321, d =   6615241;
  unsigned long t;
  t = (x ˆ (x >> 2));
  x = y;
  y = z;
  z = w;
  w = v;
  v = (v ˆ (v << 4)) ˆ (t ˆ (t << 1));
  return (d += 362437) + v;
}

シンプルで非常に高速 (1億2500万/秒)、\(2^{192}\)～\(2^{32}\) サイクルの要素は Diehard のすべてのテストを容易に合格する。

4 概要

xorshift 演算を様々な方法で組み合わせることによって単純で非常に高速な様々な RNG を開発することができるが、その単純性にもかかわらず生成される乱数はランダム性のタスとで非常に優れた特性を示している。上記で示した数百の \([a,b,c]\) に対してもっとも単純なものは基本的な C の操作である x^=(x<<a);x^=(x>>b);x^=(x<<c); を使用して、32 ビットワードの場合は \(2^{32}-1\) 周期、64 ビットワードの場合は \(2^{64}-1\) 周期を生成する。それらは単独で、また他の方法と組み合わせて使用するのに適している。

\(2^32\) 周期は現代の基準では低いと考えられている。非常にわずかな計算時間の追加だけで、xorshift 演算を使用して \(2^{64}-1\) 周期のペア \((x,y)\)、\(2^{96}-1\) 周期のトリプル \((x,y,z)\)、\(2^{128}-1\) 周期のクアドラプル \((x,y,z,w)\) または \(2^{192}-1\) 周期のクインティプル \((x,y,z,w,v)\) のシーケンスを生成することができる。より高次な \(k\)-タプルへの拡張も可能だが、前の \(k\) 次の作用として新しい値を計算する一般的なプロシジャから x=y;y=z,z=w; などのプロモートを行うことは、前に生成された \(k\) 個の値の循環テーブルを保持する方法よりも効率が低下する。このようなテーブルではキャリー付き乗算 (multiply-with-carry; MWC) やキャリー付き相補乗算 (complimentary-multiply-with-carry) の方法は \(2^{131102}\) という大きな周期を提供することができる。私はニュースグループへの投稿を通じて様々なバージョンを説明してきており、一般的な説明は [4] にある。

ここで \(2^{128}-1\) 周期の xorshift RNG と比較可能な周期を持つキャリー付き乗算 RNG とを比較する。まず xorshift:

unsigned long xor128(){
  static unsigned long x = 123456789, y = 362436069, z = 521288629, w = 88675123;
  unsigned long t;
  t = (x ˆ (x << 11));
  x = y;
  y = z;
  z = w;
  return (w = (w ˆ (w >> 19)) ˆ (t ˆ (t >> 8)));
}

次にキャリー付き乗算 (MWC):

unsigned long mwc(){
  static unsigned long x = 123456789, y = 362436069, z = 77465321, c = 13579;
  unsigned long long t;
  t = 916905990LL * x + c;
  x = y;
  y = z;
  c = (t >> 32);
  return z = (t & 0xffffffff);
}

MWC RNG は \(a = 916905990\), \(b = 2^{32}\) で \(x_n = ax_{n-3} + {\rm carry} \bmod b\) シーケンスを生成する。3 つ前の値 \(x,y,z\) および現在のキャリー \(c\) を保持し \(t=ax+c\) を 64 ビットで生成し、\(x=y,y=z\) をプロモートする。新しい \(c\) (キャリー) は \(t\) の上位 32 ビットであり、新しい \(z\) は下位 32 ビットである。周期は \((ab^3-1)/2 \approx 2^{125}\) (素数 \(ab^3-1\) に対する \(b\) の次数) である。

どちらのルーチンも Diehard の一連のテスト [2] ですべて合格する。

どちらもほんの数個の C 命令を使うだけである。xorshift RNG は最新の 4 つの値を保持する必要がある。MWC は最新の 3 つとキャリー \(c\) を保持する必要がある。

xor128 のシード集合はいずれもゼロではない 4 つの 32 ビット整数 x,y,z,w だが、MWC のシード集合は 3 つの 32 ビット整数 x,y,z と初期 c<a である (x=y=z=c=0 と x=y=z=b-1 かつ c=a-1 の 2 つのケースを除く)。

しかし xor128() は mwc() よりはるかに高速である。1800MHz の PC では xor128() は 4.4 ナノ秒 (つまり秒間 2 億 2000 万回を超える) かかるが、mwc() は 21 ナノ秒 (秒間 4800 万回) を必要とする。ただし、当然のことながら乱数を必要とするシミュレーションの問題では秒間 4800 万回はボトルネックとは考えられない。

非常に長い周期を持つ RNG に関心がある場合は最大 \(2^{131102}\) の周期を持つ MWC または CMWC 法を検討できるが、数百または数千の \(k\) の最新値のテーブルを保持する必要がある。しかし、もしそうではなくて周期が \(2^{160}\), \(2^{128}\), \(2^{96}\) 以下で満足なのであれば、単に最新の \(x\) や \(x,y\) や \(x,y,z\) を単純な形で保持するだけで大部分のアプリケーションには十分な大きさと思われる周期を提供する xorshift RNG の一つを検討する価値はある。特にこの xorshift RNG は非常に単純で、非常に高速で、ランダム性テストで非常にうまく機能する。

References

  1. Marsaglia, George and Tsay, L. H., 1985, Matrices and the structure of random number sequences, Linear Algebra and its Applications, 67, 147–156.
  2. Marsaglia, George, 1995, The Marsaglia Random Number CDROM, with The Diehard Battery of Tests of Randomness, produced at Florida State University under a grant from The National Science Foundation. Access available at www.stat.fsu.edu/pub/diehard, and a revised version of the Diehard tests atwww.csis.hku.hk/˜diehard.
  3. Marsaglia, George and Tsang, Wai Wan, 2002, Some difficult-to-pass tests of randomness, Journal Statistical Software, 7, Issue 3.
  4. Marsaglia, George, 2003, Random number generators, Journal of Modern Applied Statistical Methods, 2 No. 2.

翻訳抄

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

  1. Marsaglia, 2003, Xorshift RNGs, Journal of Statistical Software Vol.8.
  2. Google Chromeが採用した、擬似乱数生成アルゴリズム「xorshift」の数理