Xorshift

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

Xorshiftアルゴリズム^[2] のCによる実装例^[3] :

#include<stdint.h>
structxorshift32_state{
uint32_ta;
};
/* The state word must be initialized to non-zero */
uint32_txorshift32(structxorshift32_state*state)
{
/* Algorithm "xor" from p. 4 of Marsaglia, "Xorshift RNGs" */
uint32_tx=state->a;
x^=x<<13;
x^=x>>17;
x^=x<<5;
returnstate->a=x;
}
structxorshift64_state{
uint64_ta;
};
uint64_txorshift64(structxorshift64_state*state)
{
uint64_tx=state->a;
x^=x<<13;
x^=x>>7;
x^=x<<17;
returnstate->a=x;
}
structxorshift128_state{
uint32_tx[4];
};
/* The state array must be initialized to not be all zero */
uint32_txorshift128(structxorshift128_state*state)
{
/* Algorithm "xor128" from p. 5 of Marsaglia, "Xorshift RNGs" */
uint32_tt=state->x[3];

uint32_ts=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;
returnstate->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. */

このアルゴリズムの周期はそれぞれ${\displaystyle 2^{32}-1,2^{64}-1,2^{96}-1,2^{128}-1}${\displaystyle 2^{32}-1,2^{64}-1,2^{96}-1,2^{128}-1} で、Diehardテスト (英語版)をパスしている。

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

structxorshift64_state{
uint64_ta;
};
uint64_txorshift64(structxorshift64_state*state)
{
uint64_tx=state->a;
x^=x<<7;
x^=x>>9;
returnstate->a=x;
}

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

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

完全なテストとしては、期待する周期${\displaystyle k}${\displaystyle k}の約数${\displaystyle d}${\displaystyle d}について${\displaystyle T^{d}\neq I\iff k\neq d}${\displaystyle T^{d}\neq I\iff k\neq d}を示せば良い。例えば、${\displaystyle n}${\displaystyle n} bitのビット列${\displaystyle x}${\displaystyle x}に対する操作が

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

である場合、${\displaystyle T=(I+L^{a})(I+R^{b})(I+L^{c})}${\displaystyle T=(I+L^{a})(I+R^{b})(I+L^{c})}である。但し、${\displaystyle L_{i,j}^{a}={\begin{cases}1&i=j+a\0円&{\text{otherwise}}\end{cases}}}${\displaystyle L_{i,j}^{a}={\begin{cases}1&i=j+a\0円&{\text{otherwise}}\end{cases}}}及び${\displaystyle R_{i,j}^{a}={\begin{cases}1&i=j-a\0円&{\text{otherwise}}\end{cases}}}${\displaystyle R_{i,j}^{a}={\begin{cases}1&i=j-a\0円&{\text{otherwise}}\end{cases}}}である。

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

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

別の例としては

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

に対しては${\displaystyle ({\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}${\displaystyle ({\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}のように立式し、${\displaystyle 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)}${\displaystyle 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だとすれば${\displaystyle n=128}${\displaystyle n=128}, 64 bitだとすれば${\displaystyle n=256}${\displaystyle n=256}であり、対応する${\displaystyle d}${\displaystyle d}は以下の表のようになる。

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

• 乱数
• 数学に関する記事

• すべてのスタブ記事
• 応用数学関連のスタブ項目