LoginSignup
28
23

More than 1 year has passed since last update.

@umireon(Kaito Udagawa)inリックソフト株式会社

10年前の自分に伝えたい、たった一つの擬似乱数生成器

Last updated at Posted at 2021-09-19

TL;DR

Xorshift (10年前)

George Marsagliaが考案しました。当時はC#で実装しました。Equidistributionの性質が不明なので、現代では言うまでもなく使用すべきではありません。

uint32_t xor128(void) { 
  static uint32_t x = 123456789;
  static uint32_t y = 362436069;
  static uint32_t z = 521288629;
  static uint32_t w = 88675123; 
  uint32_t t;

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

From https://ja.wikipedia.org/wiki/Xorshift

Mersenne TwisterとSFMT (7年前)

松本眞氏と西村拓士氏が考案しました。C++から使っていました。
Equidistributionが簡単に大きな値にできると言うメリットがあり、非暗号用途なら何にでも応用できるというものです。一方で、Equidistributionが大きいということは、当然巨大なステートをメモリ上に保持しなければならないということを意味しているので、CPUキャッシュヒットさせて高速なシミュレーションを行う上ではメルセンヌツイスタは障害となります。

void sfmt_gen_rand_all(sfmt_t * sfmt) {
    int i;
    __m256i r;
    w128_t * pstate = sfmt->state;

    r = LD2(&pstate[SFMT_N - 2]); 
    for (i = 0; i < SFMT_N - SFMT_POS1; i+=2) {
        r = mm256_recursion(LD2(&pstate[i]), LD2(&pstate[i + SFMT_POS1]), r);
        ST2(&pstate[i], r);
    }
    for (; i < SFMT_N; i+=2) {
        r = mm256_recursion(LD2(&pstate[i]), LD2(&pstate[i + SFMT_POS1 - SFMT_N]), r);
        ST2(&pstate[i], r);
    }
}

inline static __m256i mm256_recursion(__m256i a, __m256i b, __m256i c)
{
    __m256i x, y, z, w;
    __m256i mask = _mm256_broadcastsi128_si256(sse2_param_mask.si);

    x = _mm256_slli_si256(a, SFMT_SL2);
    y = _mm256_srli_epi32(b, SFMT_SR1);

    x = _mm256_ternarylogic_epi32(x, y, mask, 0x78); // A ^ (B & C)

    z = _mm256_srli_si256(c, SFMT_SR2);

    x = _mm256_ternarylogic_epi32(x, a, z, 0x96); // XOR3(A,B,C)

    w = _mm256_permute2f128_si256(c, x, 0x21);  
    w = _mm256_slli_epi32(w, SFMT_SL1); 
    x = _mm256_xor_si256(x, w);
    return x;
}

From https://github.com/MersenneTwister-Lab/SFMT/blob/master/SFMT-avx256.h

Xorshift*, Xorshift+ (5年前)

Sebastiano Vignaが考案しました。Cのリファレンス実装をC++で利用していました。
最小のステートで最大のEquidistributionを保証する擬似乱数生成器の走りだったと記憶しています（これ以前にも存在していたかもしれませんが私は認識していません）。
Chromeに採用されて話題になったのがこの擬似乱数生成器のファミリーです。欠陥があったらしくVignaによって新しいアルゴリズムに置き換えられました。
余談ですが、私はこのソースコードからCC0と言うライセンスを知りました。

uint64_t s[2];

uint64_t next(void) {
    uint64_t s1 = s[0];
    const uint64_t s0 = s[1];
    const uint64_t result = s0 + s1;
    s[0] = s0;
    s1 ^= s1 << 23; // a
    s[1] = s1 ^ s0 ^ (s1 >> 18) ^ (s0 >> 5); // b, c
    return result; 
}

From https://prng.di.unimi.it/xorshift128plus.c

Xoroshiro+ (5年前)

Sebastiano Vignaが考案しました。自力でC++の標準互換の実装を作成しました。
さらに高速なアルゴリズムが開発されたので、Vignaによって新しいアルゴリズムに置き換えられました。
私は最終的にはこの擬似乱数生成器を自分の研究で多用していました。

/**
 * @file xoroshift128plus.hpp
 * xoroshiro128plus-cpp11 - xoroshift128plus implementation for C++11
 *
 * Written in 2016 by Kaito Udagawa <umireon@gmail.com>
 *
 * To the extent possible under law, the author(s) have dedicated all copyright
 * and related and neighboring rights to this software to the public domain
 * worldwide. This software is distributed without any warranty.
 *
 * You should have received a copy of the CC0 Public Domain Dedication along
 * with this software. If not, see
 * <http://creativecommons.org/publicdomain/zero/1.0/>.
 */

#pragma once
#include "splitmix64.hpp"

#include <array>
#include <cstdint>
#include <vector>

/**
 * @class xoroshiro128plus
 * @brief A C++11 PRNG implementation of the xoroshiro128+ algorithm
 * This is a C++11 implementation of the xoroshiro128+ algorithm, which is
 * a pseudo-random number generator proposed by Sebastiano Vigna in 2016.
 *
 * This code is based on the reference implementation written by Sebastiano
 * Vigna and David Blackman, which is available at
 * <http://xoroshiro.di.unimi.it/xoroshiro128plus.c>.
 */
class xoroshiro128plus {
public:
  // concept Uniform Random Number Generator
  typedef std::uint64_t result_type;

  constexpr static result_type default_seed = UINT64_C(1);

  xoroshiro128plus(const xoroshiro128plus &) = default;
  xoroshiro128plus(xoroshiro128plus &&) noexcept = default;
  xoroshiro128plus &operator=(const xoroshiro128plus &) = default;
  xoroshiro128plus &operator=(xoroshiro128plus &&) noexcept = default;

  result_type operator()() noexcept {
    const result_type s0 = state[0];
    result_type s1 = state[1];
    const result_type result = s0 + s1;

    s1 ^= s0;
    state[0] = rotl(s0, 55) ^ s1 ^ (s1 << 14); // a, b
    state[1] = rotl(s1, 36);                   // c

    return result;
  }

  constexpr static result_type min() noexcept { return UINT64_C(1); }

  constexpr static result_type max() noexcept {
    return UINT64_C(0xFFFFFFFFFFFFFFFF);
  }

  // concept Pseudo-Random Number Generator
  explicit xoroshiro128plus(result_type i = default_seed) : state() { seed(i); }

  template <class SeedSeq> explicit xoroshiro128plus(SeedSeq &s) : state() {
    seed(s);
  }

  void seed(result_type i = default_seed) {
    splitmix64 rng(i);
    seed(rng);
  }

  template <class SeedSeq> void seed(SeedSeq &s) {
    auto *p = reinterpret_cast<std::uint32_t *>(&state[0]);
    s.generate(p, p + 4);
  }

  void discard(unsigned long long j) noexcept {
    for (auto i = j; i < j; i++)
      this->operator()();
  }

  void jump() noexcept {
    constexpr std::array<result_type, 2> JUMP{
        {0xbeac0467eba5facb, 0xd86b048b86aa9922}};

    result_type s0 = 0;
    result_type s1 = 0;
    for (std::size_t i = 0; i < JUMP.size(); i++) {
      for (int b = 0; b < 64; b++) {
        if (JUMP[i] & UINT64_C(1) << b) {
          s0 ^= state[0];
          s1 ^= state[1];
        }
        this->operator()();
      }
    }

    state[0] = s0;
    state[1] = s1;
  }

private:
  std::array<result_type, 2> state;

  constexpr static inline result_type rotl(const result_type x,
                                           int k) noexcept {
    return (x << k) | (x >> (64 - k));
  }
};

これは公開していません。

なぜコピーコンストラクタとムーブコンストラクタに明示的にデフォルトコンストラクタを指定しているのかは、今となってはわからないので有識者の方誰か教えてください。

xoshiro256** (2年前)

Sebastiano Vignaが考案しました。Vignaが考案した最新のアルゴリズムです。
私はこのファミリーの擬似乱数生成器をVimのrand()として実装してVimのコントリビューターになりました。

static inline uint32_t rotl(const uint32_t x, int k) {
    return (x << k) | (x >> (32 - k));
}


static uint32_t s[4];

uint32_t next(void) {
    const uint32_t result = rotl(s[1] * 5, 7) * 9;

    const uint32_t 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] = rotl(s[3], 11);

    return result;
}

From https://prng.di.unimi.it/xoshiro128starstar.c

Vignaの問題

ここにはあまり深く書きませんが、Vignaはコミュニティに対して（特にPCGを支持するものに対して）FUDを仕掛けるのが得意なので、コミュニティからは嫌われています。
あちこちでやっているのでググれば出てきます。
何年もやっているので、もう直らないと思います。

PCG (今)

正直なところ、私はPCGに真剣に向き合ったことがないので技術的詳細は知りません。しかしnumpyのコミュニティで支持されていることと論文が提出されていてそれに対する反論は（Vigna以外からは）ほとんど行われていないことから、将来有力な擬似乱数生成器だと考えます。特にストリームの機能が使いやすいです。C++の公式実装もあるので苦しまずに済みます。

// *Really* minimal PCG32 code / (c) 2014 M.E. O'Neill / pcg-random.org
// Licensed under Apache License 2.0 (NO WARRANTY, etc. see website)

typedef struct { uint64_t state;  uint64_t inc; } pcg32_random_t;

uint32_t pcg32_random_r(pcg32_random_t* rng)
{
    uint64_t oldstate = rng->state;
    // Advance internal state
    rng->state = oldstate * 6364136223846793005ULL + (rng->inc|1);
    // Calculate output function (XSH RR), uses old state for max ILP
    uint32_t xorshifted = ((oldstate >> 18u) ^ oldstate) >> 27u;
    uint32_t rot = oldstate >> 59u;
    return (xorshifted >> rot) | (xorshifted << ((-rot) & 31));
}

From https://www.pcg-random.org/download.html

結論

とりあえずPCGを使っておけばいいです。numpyを使えば自動でPCGになるので、数値シミュレーションにはnumpyを使用するべきなのかもしれません。

28
23
2

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
28
23