Xorshiftのソースを表示

'''Xorshift'''は[[疑似乱数]]列生成法の1つである。George Marsaglia（[[w:George Marsaglia]]）が2003年に提案した。演算が[[ビット演算#XOR|排他的論理和]]と[[ビット演算#ビットシフト|ビットシフト]]のみであるため高速である<ref>2003年の論文執筆時点で、1800MHzのPCで毎秒2.2億個の擬似乱数を生成できた。</ref>
<!--、かつ[[線形合同法]]よりも長周期--><!-- ← 内部状態の大きさから可能な最長周期となる生成法の作り方が確立している点は線形合同法と同じです。定量的ではなく（間違った）感覚的な話であるが、何か勘違いをしている -->
などの特徴がある。

== 実装例 ==
Xorshiftアルゴリズム<ref name="origin"> Marsaglia, 2003 </ref>の[[C言語|C]]による実装例<ref>Cでは "<code>^</code>" は[[ビット演算#XOR|排他的論理和]]を、"<code>&lt;&lt;</code>" と "<code>&gt;&gt;</code>" は[[ビット演算#ビットシフト|ビットシフト]]を表す。</ref>:

<!-- 論文中の "unsigned long" は、ここでは "uint32_t" に変更してある。
なぜなら、64 ビット環境では "unsigned long" は 64 ビットであることもあるからだ。
一方マーサグリアは論文中では "unsigned long long" を 64 ビット整数として使用している。
どちらにしろ、サイズをはっきりさせた方が分かりやすいだろう。 -->
<syntaxhighlight lang="c">
#include <stdint.h>

struct xorshift32_state {
  uint32_t a;
};

/* The state word must be initialized to non-zero */
uint32_t xorshift32(struct xorshift32_state *state)
{
	/* Algorithm "xor" from p. 4 of Marsaglia, "Xorshift RNGs" */
	uint32_t x = state->a;
	x ^= x << 13;
	x ^= x >> 17;
	x ^= x << 5;
	return state->a = x;
}

struct xorshift64_state {
  uint64_t a;
};

uint64_t xorshift64(struct xorshift64_state *state)
{
	uint64_t x = state->a;
	x ^= x << 13;
	x ^= x >> 7;
	x ^= x << 17;
	return state->a = x;
}

struct xorshift128_state {
  uint32_t x[4];
};

/* The state array must be initialized to not be all zero */
uint32_t xorshift128(struct xorshift128_state *state)
{
	/* Algorithm "xor128" from p. 5 of Marsaglia, "Xorshift RNGs" */
	uint32_t t  = state->x[3];
    
    uint32_t s  = state->x[0];  /* Perform a contrived 32-bit shift. */
	state->x[3] = state->x[2];
	state->x[2] = state->x[1];
	state->x[1] = s;

	t ^= t << 11;
	t ^= t >> 8;
	return state->x[0] = t ^ s ^ (s >> 19);
}

/* The Xorshift128 algorithm can be reworded to use a pair of uint64_t state,
   or for systems with such support, a uint128_t state. */
</syntaxhighlight>

このアルゴリズムの周期はそれぞれ<math>2^{32}-1, 2^{64}-1, 2^{96}-1, 2^{128}-1</math> で、{{仮リンク|Diehardテスト|en|Diehard tests}}をパスしている。

64ビット整数を効率よく扱える計算機では、xor,shift操作の組を3つから2つにして計算負荷を減らしても、周期は<math>2^{64}-1</math>に保たれる。<ref>{{Cite web|和書|author=和田維作 |url=http://isaku-wada.my.coocan.jp/rand/rand.html |title=良い乱数・悪い乱数 |accessdate=2023-08-28}} 但し、周期が保たれる組み合わせは(7,9)と(9,7)のみである。</ref>
<syntaxhighlight lang="C">
struct xorshift64_state {
  uint64_t a;
};

uint64_t xorshift64(struct xorshift64_state *state)
{
	uint64_t x = state->a;
	x ^= x << 7;
	x ^= x >> 9;
	return state->a = x;
}
</syntaxhighlight>

== 周期の特定 ==
シード値を128bitの二元横ベクトル<math>\beta</math>、現在の状態から次状態への二元遷移行列を<math>T</math>と置くと、Xorshiftアルゴリズムにより生成される乱数列は<math>\beta, \beta T, \beta T^{2}, \dots</math>と表される。詳しい証明は省略するが<ref name="origin" />、行列<math>T</math>のOrderが<math>k=2^n-1</math>で表される時、かつその時に限り、任意の0でない<math>\beta</math>に対し<math>\beta, \beta T, \beta T^2, \dots, \beta T^{k-1}</math>は全て異なる値となる。
<!--
原文では<math>\beta, \beta T, \beta T^2, \dots</math>が(n次二元ベクトルとして)可能な全ての値を取る、という表現となっている。
-->

予備的なテストとしては、今周期<math>k</math>について<math>k=2^n-1</math>となることを期待しているため、期待する<math>n</math>が<math>T^{2^n}=T</math>を満たす最小の<math>n</math>であるかを調べる、というものがある。これは<math>T</math>を<math>n</math>回二乗すれば良いため、乱数列を調べるのと比較すると十分に速く調べられる。

完全なテストとしては、期待する周期<math>k</math>の約数<math>d</math>について<math>T^d \neq I \iff k \neq d</math>を示せば良い。<!-- 原文には明示的な記載は無いが、自明である。 -->例えば、<math>n</math> bitのビット列<math>x</math>に対する操作が<syntaxhighlight lang="c">
x ^= x << a;
x ^= x >> b;
x ^= x << c;
return x;
</syntaxhighlight>である場合、<math>T = (I+L^a)(I+R^b)(I+L^c)</math>である。但し、<math>L^a_{i,j}=
\begin{cases}
1 & i = j + a \\
0 & \text{otherwise}
\end{cases}</math>及び<math>R^a_{i,j}=\begin{cases}
1 & i = j - a \\
0 & \text{otherwise}
\end{cases}</math>である。

<math>n=32</math>の場合、<math>T^d \neq I \Longleftarrow
d=\begin{cases}
65535 \\
16711935 \\
252645135 \\
858993459 \\
1431655765
\end{cases}</math>及び<math>T^{2^{32}} = T</math>を調べれば十分である。

<math>n=64</math>の場合、<math>2^{64}-1</math>が<math>641</math>の倍数であるということに注意すると、<math>T^d \neq I \Longleftarrow
d = \begin{cases}
2753074036095 \\
281470681808895 \\
28778071877862015 \\
71777214294589695 \\
1085102592571150095 \\
3689348814741910323 \\
6148914691236517205
\end{cases}</math>及び<math>T^{2^{64}}=T</math>を調べれば十分である。

別の例としては<syntaxhighlight lang="c">
t = x ^ (x << 11);
x = y; y = z; z = w;
return w = (w ^ (w >> 19)) ^ (t ^ (t >> 8)); 
</syntaxhighlight>に対しては<math>(
\begin{smallmatrix}
x_{n+1} & y_{n+1} & z_{n+1} & w_{n+1}
\end{smallmatrix}
)=
(
\begin{smallmatrix}
x_n & y_n & z_n & w_n
\end{smallmatrix}
) \cdot T</math>のように立式し、<math>T=\left(
\begin{smallmatrix}
O & O & O & (I+L^{11})(I+R^8) \\
I & O & O & O \\
O & I & O & O \\
O & O & I & (I+R^{19}) \\
\end{smallmatrix}
\right)</math>について同様に解く。各変数が32 bitだとすれば<math>n=128</math>, 64 bitだとすれば<math>n=256</math>であり、対応する<math>d</math>は以下の表のようになる。
{| class="wikitable"
!<math>n</math>
!<math>32</math>
!<math>64</math>
!<math>128</math>
!<math>256</math>
|-
| rowspan="11" |<math>d</math>
|<code>0x0000ffff</code>
|<code>0x00000280fffffd7f</code>
|<code>0x0000000000042f00fffffffffffbd0ff</code>
|<code>0x000000000000000000d3eafc3af14600ffffffffffffffffff2c1503c50eb9ff</code>
|-
|<code>0x00ff00ff</code>
|<code>0x0000ffff0000ffff</code>
|<code>0x00000280fffffd7f00000280fffffd7f</code>
|<code>0x000000000000013540775b48cc32ba00fffffffffffffecabf88a4b733cd45ff</code>
|-
|<code>0x0f0f0f0f</code>
|<code>0x00663d80ff99c27f</code>
|<code>0x00003d30f19cd100ffffc2cf0e632eff</code>
|<code>0x0000000000042f00fffffffffffbd0ff0000000000042f00fffffffffffbd0ff</code>
|-
|<code>0x33333333</code>
|<code>0x00ff00ff00ff00ff</code>
|<code>0x0000ffff0000ffff0000ffff0000ffff</code>
|<code>0x00000280fffffd7f00000280fffffd7f00000280fffffd7f00000280fffffd7f</code>
|-
|<code>0x55555555</code>
|<code>0x0f0f0f0f0f0f0f0f</code>
|<code>0x00663d80ff99c27f00663d80ff99c27f</code>
|<code>0x00003d30f19cd100ffffc2cf0e632eff00003d30f19cd100ffffc2cf0e632eff</code>
|-
|
|<code>0x3333333333333333</code>
|<code>0x00ff00ff00ff00ff00ff00ff00ff00ff</code>
|<code>0x0000ffff0000ffff0000ffff0000ffff0000ffff0000ffff0000ffff0000ffff</code>
|-
|
|<code>0x5555555555555555</code>
|<code>0x0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f</code>
|<code>0x00663d80ff99c27f00663d80ff99c27f00663d80ff99c27f00663d80ff99c27f</code>
|-
|
|
|<code>0x33333333333333333333333333333333</code>
|<code>0x00ff00ff00ff00ff00ff00ff00ff00ff00ff00ff00ff00ff00ff00ff00ff00ff</code>
|-
|
|
|<code>0x55555555555555555555555555555555</code>
|<code>0x0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f</code>
|-
|
|
|
|<code>0x3333333333333333333333333333333333333333333333333333333333333333</code>
|-
|
|
|
|<code>0x5555555555555555555555555555555555555555555555555555555555555555</code>
|}

== 脚注 ==
{{Reflist}}

== 参考文献 ==
*{{Cite journal
| last=Marsaglia
| title=Xorshift RNGs
| journal=[[Journal of Statistical Software]]
| volume=Vol.&nbsp;8
| issue=Issue&nbsp; 14
| month=July
| year=2003
| url=http://www.jstatsoft.org/v08/i14/paper
}}

{{Applied-math-stub}}

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