競プロ 乱数 速度調査

疑似乱数の生成速度について軽く調べた。まず結果から紹介する。それぞれ $10^8$ 回乱数を生成してかかった時間をAtCoderのコードテストで計測した。実験コードはここ

Name	Output bit	Time
LCG32	32bit	114ms
mt19937	32bit	325ms
mt19937_64	64bit	389ms
xorshift32	32bit	192ms
xorshift64	64bit	193ms
xoshiro128++	32bit	249ms
xoshiro128**	32bit	238ms
xoshiro256++	64bit	249ms
xoshiro256**	64bit	238ms
pcg32	32bit	134ms
pcg32_fast	32bit	105ms
pcg64	64bit	202ms
pcg64_fast	64bit	157ms
Mwc128XXA32	32bit	92ms
Mwc256XXA64	64bit	108ms

Output bitで分類すると次のようになる。

Name	Output bit	Time
LCG32	32bit	114ms
mt19937	32bit	325ms
xorshift32	32bit	192ms
xoshiro128++	32bit	249ms
xoshiro128**	32bit	238ms
pcg32	32bit	134ms
pcg32_fast	32bit	105ms
Mwc128XXA32	32bit	92ms

Name	Output bit	Time
mt19937_64	64bit	389ms
xorshift64	64bit	193ms
xoshiro256++	64bit	249ms
xoshiro256**	64bit	238ms
pcg64	64bit	202ms
pcg64_fast	64bit	157ms
Mwc256XXA64	64bit	108ms

LCG (線形合同法)

$X = (A \times X + B) \bmod M$ という形で表されるやつ。パラメーターはwikipediaから(あまりよくないパラメーターらしい)

// Reference:
// https://ja.wikipedia.org/wiki/%E7%B7%9A%E5%BD%A2%E5%90%88%E5%90%8C%E6%B3%95
namespace lcg32 {

inline static u32 state = 12345;

u32 next() { return state = 1'103'515'245 * state + 12'345; }

}  // namespace lcg32

mt19937 (メルセンヌ・ツイスター)

C++のSTLのmt19937 / mt19937_64をそのまま

// Reference:
// https://cpprefjp.github.io/reference/random/mt19937.html
namespace mt19937 {

inline static std::mt19937 engine(12345);

u32 next() { return (u32)engine(); }

}  // namespace mt19937

// Reference:
// https://cpprefjp.github.io/reference/random/mt19937_64.html
namespace mt19937_64 {

inline static std::mt19937_64 engine(12345);

u64 next() { return (u64)engine(); }

}  // namespace mt19937_64

xorshift

主観だと競プロで一番使われていそうなやつ。実装はwikipediaから殆どそのまま

// Reference:
// https://ja.wikipedia.org/wiki/Xorshift
namespace xorshift32 {

inline static u32 a = 12345;

u32 next() {
    u32 x = a;
    x ^= x << 13;
    x ^= x >> 17;
    x ^= x << 5;
    return a = x;
}

}  // namespace xorshift32

// Reference:
// https://ja.wikipedia.org/wiki/Xorshift
namespace xorshift64 {

inline static u64 a = 12345;

u64 next() {
    u64 x = a;
    x ^= x << 13;
    x ^= x >> 7;
    x ^= x << 17;
    return a = x;
}

}  // namespace xorshift64

xoshiro++ / xoshiro**

https://prng.di.unimi.it/ に詳しい紹介がある。速くて質もいいらしい。実装はこのサイトから殆どそのまま

// Reference:
// https://prng.di.unimi.it/
namespace xoshiro128plusplus {

inline static u32 s[4] = {123, 234, 345, 567};

u32 next() {
    const u32 result_starstar = std::rotl(s[0] + s[3], 7) + s[0];

    const u32 t = s[1] << 9;

    s[2] ^= s[0];
    s[3] ^= s[1];
    s[1] ^= s[2];
    s[0] ^= s[3];

    s[2] ^= t;

    s[3] = std::rotl(s[3], 11);

    return result_starstar;
}

}  // namespace xoshiro128plusplus

// Reference:
// https://prng.di.unimi.it/
namespace xoshiro128starstar {

inline static u32 s[4] = {123, 234, 345, 567};

u32 next() {
    const u32 result_starstar = std::rotl(s[1] * 5, 7) * 9;

    const u32 t = s[1] << 9;

    s[2] ^= s[0];
    s[3] ^= s[1];
    s[1] ^= s[2];
    s[0] ^= s[3];

    s[2] ^= t;

    s[3] = std::rotl(s[3], 11);

    return result_starstar;
}

}  // namespace xoshiro128starstar

// Reference:
// https://prng.di.unimi.it/
namespace xoshiro256plusplus {

inline static u64 s[4] = {123, 234, 345, 567};

u64 next() {
    const u64 result_starstar = std::rotl(s[0] + s[3], 23) + s[0];

    const u64 t = s[1] << 17;

    s[2] ^= s[0];
    s[3] ^= s[1];
    s[1] ^= s[2];
    s[0] ^= s[3];

    s[2] ^= t;

    s[3] = std::rotl(s[3], 45);

    return result_starstar;
}

}  // namespace xoshiro256plusplus

// Reference:
// https://prng.di.unimi.it/
namespace xoshiro256starstar {

inline static u64 s[4] = {123, 234, 345, 567};

u64 next() {
    const u64 result_starstar = std::rotl(s[1] * 5, 7) * 9;

    const u64 t = s[1] << 17;

    s[2] ^= s[0];
    s[3] ^= s[1];
    s[1] ^= s[2];
    s[0] ^= s[3];

    s[2] ^= t;

    s[3] = std::rotl(s[3], 45);

    return result_starstar;
}

}  // namespace xoshiro256starstar

PGC

https://www.pcg-random.org/index.html に詳しい紹介がある。速くて質もいいらしい。

標準のpcg32 / pcg64と、ストリーム機能(とちょっと短い周期)の代わりに速度に特化したpcg32_fast / pcg64_fastがある。

pcg32 / pgc32_fast は wikipedia の実装を参考に、 pcg64 / pcg64_fast はrustの rand_pcg crateの実装を参考にした。

// Reference:
// https://www.pcg-random.org/download.html
// https://en.wikipedia.org/wiki/Permuted_congruential_generator
namespace pcg32 {

const u64 MULT = 6364136223846793005;
const u64 INC = 1442695040888963407;

inline static u64 state = 123;

u32 next() {
    u64 x = state;
    state = state * MULT + INC;

    u32 count = x >> 59;
    x ^= x >> 18;

    return std::rotr(u32(x >> 27), count);
}

}  // namespace pcg32

// Reference:
// https://www.pcg-random.org/download.html
// https://en.wikipedia.org/wiki/Permuted_congruential_generator
namespace pcg32_fast {

const u64 MULT = 6364136223846793005;

inline static u64 state = 123;

u32 next() {
    u64 x = state;
    state = state * MULT;

    u32 count = x >> 61;
    x ^= x >> 22;

    return (u32)(x >> (22 + count));
}

}  // namespace pcg32_fast

// Reference:
// https://www.pcg-random.org/download.html
// https://crates.io/crates/rand_pcg
namespace pcg64 {

const u128 MULT = u128(2549297995355413924) << 64 | u128(4865540595714422341);
const u128 INC = u128(6364136223846793005) << 64 | u128(1442695040888963407);

inline static u128 state = 123;

u64 next() {
    u128 x = state;
    state = state * MULT + INC;

    u32 rot = x >> 122;
    return std::rotr((u64)(x >> 64) ^ (u64)x, rot);
}

}  // namespace pcg64

// Reference:
// https://www.pcg-random.org/download.html
// https://crates.io/crates/rand_pcg
namespace pcg64_fast {

const u128 MULT = u128(2549297995355413924) << 64 | u128(4865540595714422341);

inline static u128 state = 123;

u64 next() {
    u128 x = state;
    state = state * MULT;

    u32 rot = x >> 122;
    return std::rotr((u64)(x >> 64) ^ (u64)x, rot);
}

}  // namespace pcg64_fast

Mwc128xxa32 / Mwc256xxa64

このブログで紹介されている。質が良くてPCGより速いらしい。

実装は githubから殆どそのまま。変更点はx1とcをまとめている。

// Reference:
// https://tom-kaitchuck.medium.com/designing-a-new-prng-1c4ffd27124d
// https://github.com/tkaitchuck/Mwc256XXA64
namespace mwc128xxa32 {

const u32 MULT = 3487286589;

inline static u32 x2 = 12345;
inline static u32 x3 = 0xcafef00d;
inline static u64 c_x1 = u64(0xd15ea5e5) << 32 | 23456;

u32 next() {
    u64 x = (u64)(x3)*MULT;
    u32 result = (x3 ^ x2) + ((u32)(c_x1) ^ (u32)(x >> 32));
    x3 = x2;
    x2 = (u32)(c_x1);
    c_x1 = x + (c_x1 >> 32);
    return result;
}

}  // namespace mwc128xxa32

// Reference:
// https://tom-kaitchuck.medium.com/designing-a-new-prng-1c4ffd27124d
// https://github.com/tkaitchuck/Mwc256XXA64
namespace mwc256xxa64 {

const u64 MULT = 0xfeb3'4465'7c0a'f413;

inline static u64 x2 = 12345;
inline static u64 x3 = 0xcafef00dd15ea5e5;
inline static u128 c_x1 = u128(0x1405'7B7E'F767'814F) << 64 | 23456;

u64 next() {
    u128 x = (u128)(x3)*MULT;
    u64 result = (x3 ^ x2) + ((u64)(c_x1) ^ (u64)(x >> 64));
    x3 = x2;
    x2 = (u64)(c_x1);
    c_x1 = x + (c_x1 >> 64);
    return result;
}

}  // namespace mwc256xxa64