Xorshift

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

などの特徴がある。

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


```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. */
```


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

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

```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;
}
```


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

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

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

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

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

別の例としては
```c
t = x ^ (x << 11);
x = y; y = z; z = w;
return w = (w ^ (w >> 19)) ^ (t ^ (t >> 8));
```
に対しては(
\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のように立式し、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)について同様に解く。各変数が32 bitだとすればn=128, 64 bitだとすればn=256であり、対応するdは以下の表のようになる。

## 脚注

## 参考文献
*