Xorshift

Xorshiftは疑似乱数列生成法の1つである。George Marsagliaが2003年に提案した。演算が排他的論理和とビットシフトのみであるため高速であるなどの特徴がある。

実装例

XorshiftアルゴリズムのCによる実装例:

uint32_t xor(void) {
  static uint32_t y = 2463534242;
  y = y ^ (y << 13); y = y ^ (y >> 17);
  return y = y ^ (y << 15);
}

uint32_t xor64(void) {
  static uint64_t x = 88172645463325252ULL;
  x = x ^ (x << 13); x = x ^ (x >> 7);
  return x = x ^ (x << 17);
}

uint32_t xor96(void) {
  static uint32_t x = 123456789;
  static uint32_t y = 362436069;
  static uint32_t z = 521288629;
  uint32_t t;

  t = (x ^ (x << 3)) ^ (y ^ (y >> 19)) ^ (z ^ (z << 6));
  x = y; y = z;
  return z = t;
}

uint32_t xor128(void) { 
  static uint32_t x = 123456789;
  static uint32_t y = 362436069;
  static uint32_t z = 521288629;
  static uint32_t w = 88675123; 
  uint32_t t;
 
  t = x ^ (x << 11);
  x = y; y = z; z = w;
  return w = (w ^ (w >> 19)) ^ (t ^ (t >> 8)); 
}

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

64ビット整数を効率よく扱える計算機では、xor,shift操作の組を3つから2つにして計算負荷を減らしても、周期は[math]2^{64}-1[/math]に保たれる。[1]

uint64_t xor64(void) {
  static uint64_t x = 88172645463325252ULL;
  x = x ^ (x << 7);
  return x = x ^ (x >> 9);
}

周期の特定

シード値を128bitの二元横ベクトル[math]\beta[/math]、現在の状態から次状態への二元遷移行列を[math]T[/math]と置くと、Xorshiftアルゴリズムにより生成される乱数列は[math]\beta, \beta T, \beta T^{2}, \dots[/math]と表される。詳しい証明は省略するが、行列[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]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]に対する操作が

x ^= x << a;
x ^= x >> b;
x ^= x << c;
return x;

である場合、[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]を調べれば十分である。

別の例としては

t = x ^ (x << 11);
x = y; y = z; z = w;
return w = (w ^ (w >> 19)) ^ (t ^ (t >> 8));

に対しては[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]は以下の表のようになる。

[math]n[/math]	[math]32[/math]	[math]64[/math]	[math]128[/math]	[math]256[/math]
[math]d[/math]	0x0000ffff	0x00000280fffffd7f	0x0000000000042f00fffffffffffbd0ff	0x000000000000000000d3eafc3af14600ffffffffffffffffff2c1503c50eb9ff
	0x00ff00ff	0x0000ffff0000ffff	0x00000280fffffd7f00000280fffffd7f	0x000000000000013540775b48cc32ba00fffffffffffffecabf88a4b733cd45ff
	0x0f0f0f0f	0x00663d80ff99c27f	0x00003d30f19cd100ffffc2cf0e632eff	0x0000000000042f00fffffffffffbd0ff0000000000042f00fffffffffffbd0ff
	0x33333333	0x00ff00ff00ff00ff	0x0000ffff0000ffff0000ffff0000ffff	0x00000280fffffd7f00000280fffffd7f00000280fffffd7f00000280fffffd7f
	0x55555555	0x0f0f0f0f0f0f0f0f	0x00663d80ff99c27f00663d80ff99c27f	0x00003d30f19cd100ffffc2cf0e632eff00003d30f19cd100ffffc2cf0e632eff
		0x3333333333333333	0x00ff00ff00ff00ff00ff00ff00ff00ff	0x0000ffff0000ffff0000ffff0000ffff0000ffff0000ffff0000ffff0000ffff
		0x5555555555555555	0x0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f	0x00663d80ff99c27f00663d80ff99c27f00663d80ff99c27f00663d80ff99c27f
			0x33333333333333333333333333333333	0x00ff00ff00ff00ff00ff00ff00ff00ff00ff00ff00ff00ff00ff00ff00ff00ff
			0x55555555555555555555555555555555	0x0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f0f
				0x3333333333333333333333333333333333333333333333333333333333333333
				0x5555555555555555555555555555555555555555555555555555555555555555

参考文献

• Marsaglia (July 2003). “Xorshift RNGs”. Journal of Statistical Software Vol. 8 (Issue  14). http://www.jstatsoft.org/v08/i14/paper.