Xorshiftを統計テストで落とそう！
はじめに
For execution in software, xorshift generators are among the fastest PRNGs, requiring very small code and state. However, they do not pass every statistical test without further refinement. This weakness is amended by combining them with a non-linear function, as described in the original paper. Because plain xorshift generators (without a non-linear step) fail some statistical tests, they have been accused of being unreliable.[3]: 360
素のXorshift[1]は非線化を挟まないと特定の統計テストに落ちるとよく言われます。
今回は実際に状態数32bitのXorshiftをMatrixRank[2]で落としてみましょう！[3]
Xorshiftの復習
C言語のXorshiftの実装コードを確認します。
#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;
}
ご存知の通り、このXorshift特有の
x ^= x << 13 /* なにか定数 */;
はビットベクトル上の行列積として表せます。
ということで
x ^= x << 13;
x ^= x >> 17;
x ^= x << 5;
全部合わせて一個の行列とビットベクトルの積で表せます。
じゃあなんなの？
本来理想的な乱数は独立しているべきですがXorshiftの出力はビットベクトルとして見ると他の出力にめちゃくちゃ依存しています。
Matrix Rank
Matrix RankはPRNGの出力をビット行列にして行列のランクを取ります。
その分布を理想的な分布と比較してカイ二乗検定[4]すればだめなやつがふるい落とされるというわけです。
以下32x32のビット行列として話をします
ねらい
一般的にベクトルを線形変換した場合前と後のベクトルは互いに独立している場合もそうでない場合があります。
ですが少なくとも理想的な乱数と同じようなランクの分布を再現するのは難しいはずです。
実装
ビットベクトルのランクは掃き出し法で簡単に求めることができます。
// Rust
// Matrix rank of bit matrix
pub fn matrix_rank(matrix: &mut [u32]) -> usize {
let n = matrix.len();
let mut rank = 0;
for bit in 0..u32::BITS {
let mut pivot = None;
for i in rank..n {
if (matrix[i] >> bit) & 1 == 1 {
pivot = Some(i);
break;
}
}
if let Some(pivot) = pivot {
matrix.swap(rank, pivot);
for i in 0..n {
if i != rank && (matrix[i] >> bit) & 1 == 1 {
matrix[i] ^= matrix[rank];
}
}
rank += 1;
}
}
rank
}
理想的な乱数からビット行列を作ったときのランクの分布はこんな感じになるはずです[5]
横軸がランクの数で縦軸が頻度です。
ちなみに/dev/randomのランクの分布はこんな感じです
10e6回分の分布です（以下全部1e6回）
よさそう
結果
Xorshiftのランクはこんな感じです
なんと最高ランクしか出ません！
緩和策
素のXorshiftはビットの世界で線形なのでMatrixRankがよくありませんでした。
ですがビットの世界で非線形な変換、例えばmod 2^32の世界で加算すれば…
fn xorshift(state: &mut u32) -> u32 {
*state ^= *state << 13;
*state ^= *state >> 17;
*state ^= *state << 5;
state.wrapping_add(114514)
}
おお
Discussion