疑似乱数xoshiro/xoroshiroについて(1) by @STDNUL

xoshiro/xoroshiroについては以前に鍋谷氏が書いているので前提は省いて
実装にはいる。

マルチインスタンス対応の実装

公式実装だと思う具合に並列化できないので、状態空間をグローバル変数ではなくローカル変数にして複数インスタンスを持てるように改造し、かつ、xoshiro256+/xoshiro256++/xoshiro256**のコードを一つに纏めたのが以下のコードであるが、そんな欲張り仕様でもこんなにコンパクトに収まるので有る。

#include <stdint.h>

struct XOSHIRO256 { 
    uint64_t s[4]; 
};
typedef struct XOSHIRO256 xoshiro256_t;

inline void next_common(xoshiro256_t *state) {
    const uint64_t t = state->s[1] << 17;

    state->s[2] ^= state->s[0];
    state->s[3] ^= state->s[1];
    state->s[1] ^= state->s[2];
    state->s[0] ^= state->s[3];

    state->s[2] ^= t;

    state->s[3] = rotl(state->s[3], 45);
}
inline uint64_t next_ss(xoshiro256_t *state) {
    const uint64_t result = rotl(state->s[1] * 5, 7) * 9;
    next_common(state);
    return result;
}
inline uint64_t next_p(xoshiro256_t *state) {
    const uint64_t result = state->s[0] + state->s[3];
    next_common(state);
    return result;
}
inline uint64_t next_pp(xoshiro256_t *state) {
	const uint64_t result = rotl(state->s[0] + state->s[3], 23) + state->s[0];
    next_common(state);
    return result;
}

static void jump_common(xoshiro256_t *state, const uint64_t jump_array[4]) {
    uint64_t s0, s1, s2, s3;
    s0 = s1 = s2 = s3 = 0;
    for(int i = 0; i < 4; i++) {
        for(int b = 0; b < 64; b++) {
            if (jump_array[i] & UINT64_C(1) << b) {
                s0 ^= state->s[0];
                s1 ^= state->s[1];
                s2 ^= state->s[2];
                s3 ^= state->s[3];
            }
            next_common(state);	
        }
    }
    state->s[0] = s0;
    state->s[1] = s1;
    state->s[2] = s2;
    state->s[3] = s3;
}

void jump(xoshiro256_t *state) {
    static const uint64_t JUMP[] = { 
        0x180ec6d33cfd0aba, 
        0xd5a61266f0c9392c, 
        0xa9582618e03fc9aa, 
        0x39abdc4529b1661c 
    };
    jump_common(state, JUMP);
}

void long_jump(xoshiro256_t *state) {
    static const uint64_t LONG_JUMP[] = { 
        0x76e15d3efefdcbbf,
        0xc5004e441c522fb3,
        0x77710069854ee241,
        0x39109bb02acbe635
    };
    jump_common(state, LONG_JUMP);
}
/* 
usage:
  xoshiro256_t state, state2;
  for (int i = 0; i < 4; i++) {
      state.s[i] = RDSEED64(); // example
  }

  memcpy(&state2, &state, sizeof(state));
  jump(&state2);

  uint64_t r1 = next_pp(&state);
  uint64_t r2 = next_pp(&state2);
*/

jump関数とマルチスレッド・SIMD化

xoshiro256の場合は256回のnext()関数の呼び出しとXORだけで2^128回もしくは2^192回先の状態に遷移することができる。
マルチスレッドやSIMDで複数のインスタンスを同時に扱いたい場合に有用である。
つまりjump関数はオリジナルの実装ではほとんど役に立たない。
複数のスレッドでSIMDを駆使して複数の状態ベクトルを並列に操作して大量に乱数を生成したいに、生成空間が被らないようにするために使うことができる。

SFMTにおけるjumpについて

疑似乱数におけるジャンプという概念を語るにあたり、その先鞭をつけた存在であろうSFMTは欠かせない存在で有ると思う。
SFMTは後付け的にjump関数を生成するためのツールが用意されたのだが、NTLというライブラリの導入の必要があったり、如何せん環境構築が面倒である。

その点においてxoshiro/xoroshiroのようにjump関数が最初から用意されていてすぐに使えるのは非常に利便性が良い。

もっとも、SFMTくらい状態空間が長いと、Intelの命令でRDSEEDなどのハードウェア乱数で初期ベクトルを初期化してやったらそうそう被らない気がするし、「MT/SFMTで乱数が被りやすい」という問題は、たった32ビット(2^32通り)のSEEDを線形合同法で伸ばしてベクトルを生成してそのまま使っているなど、系列を被らせないための使い方を使用者が認識していないなどの問題の方がよっぽど大きいと思っている。

並列実装こそxoshiro/xoroshiroの真骨頂

次回はxoshiroのSIMDを使った並列実装について語りたい。

続く