Xorshift 乱数生成アルゴリズムを書いてみた。

前回 BarrageLL の乱数について触れたので、今回もその続き。


計算機の応用において、擬似乱数というのは大切です。
なので、世の中には、いろいろな擬似乱数の生成アルゴリズムが存在するわけですが、
今回はその中で比較的新しいアルゴリズム、 xorshift 法の実装を書いてみました。


このアルゴリズムは、詳しくは検索すれば論文が出てくるのでそれを参照してもらうとして、大雑把には

  • 内部状態が小さい。他の有名な乱数生成アルゴリズムである mt19937 が 623*32 bit の内部状態を持つのに対し、 xorshift アルゴリズムの代表的な例である xor128 の内部状態は 128bit のみである。
  • 計算が速い。このアルゴリズムはビットシフトと xor 演算のみで構成されており、かつ条件分岐を含まないため、殆どのアーキテクチャで高速に実行できる。
  • 上記の性質を持つにも関わらず、周期が長い。 xor128 の周期は 2^128 - 1 であり、これは「全ての bit が 0 でない限り、乱数を生成し続ければ、特定の内部状態へといつかは遷移する」ということ。
  • 一方で乱数の分布の質は mt19937 に比べると見劣りする、とはいえ実用上は十分。

といった性質を持ちます。
特に大事なのは生成速度と周期の長さで、学術的に分布が一様であることが求められているならともかく、
それ以外の用途では、速度と質を兼ね備えた、十二分に優秀な乱数生成アルゴリズムなのです。


で、そんな xorshift アルゴリズムを、試しに C++ で実装してみました。
参考はこちら。 Mr.Exception 乱数生成器XorshiftでRandomクラスを実装して速度を測ってみた

// N-bit rotate するコード
#include <limits> // for std::numeric_limits
#include <boost/integer/integer_mask.hpp>

template<std::size_t N, class Int>
inline Int rotate_left( Int x )
{
  typedef std::numeric_limits<Int> limits;
  static const std::size_t Int_bits = limits::digits + ( limits::is_signed ? 1 : 0 );
  static const Int mask = static_cast<Int>( boost::low_bits_mask_t<N>::sig_bits );
  return ( x << N ) | ( mask & ( x >> ( Int_bits - N ) ) );
}

#include <boost/cstdint.hpp>
#include <cassert>
struct Xor128
{
  typedef boost::uint32_t result_type;
  typedef boost::uint32_t int_t;
  
  // デフォルトの初期状態で乱数生成器を構築 
  // これらの数値は論文に書いてあった値。全部 0 でないならなんでもいい。
  Xor128()
    : x(123456789), y(362436069), z(521288629), w(88675123) {}
  
  // 指定された内部状態を持つ乱数生成器を構築
  Xor128( int_t x_, int_t y_, int_t z_, int_t w_ )
    : x(x_), y(y_), z(z_), w(w_)
  {
    // 全部 0 だと困るのでそれだけチェック
    assert( x != 0 || y != 0 || z != 0 || w != 0 );
  }
  
  // 指定された種から乱数生成器を構築する
  explicit Xor128( int_t seed )
    : x(123456789), y(362436069), z(521288629), w(88675123)
  {
    // 種を埋め込む
    x ^= seed;
    y ^= rotate_left< 8>( seed );
    z ^= rotate_left<24>( seed );
    w ^= rotate_left<16>( seed );
  }
  
  // 乱数を実際に生成
  result_type operator()()
  {
    int_t const t = x ^ ( x << 11 );
    x = y; y = z; z = w;
    w = ( w ^ ( w >> 19 ) ) ^ ( t ^ ( t >> 8 ) );
    return w;
  }
  
 private:
  int_t x, y, z, w;
};

見たとおり、実装は特に難しいものではないのですが、一点、工夫した点があります。
それは乱数の種を与えて初期化するときの処理で、種を 8bit ずつ rotate したものを用いています。
これは、主に乱数の種として time 関数の戻り値が使われた場合を想定したもので、
その場合、種の下位 bit の方は多様性がありますが、上位 bit の方は殆ど一定と考えられるので、
単純に与えられた種から設定するのではなく、多様性のある部分が少しでもバラけるよう配慮してみました。
またその際、単純に 8bit ずつ rotate したものを順に埋め込んでいくのではなく、
乱数生成する時と同じ順にならないように一箇所 捻ることで、少しでも多様性を出そうとしています。


これらの多様性を出そうとする工夫は、乱数の周期が長い以上、一見して要らないようですが、
高々 32bit 程度の多様性しか持たない種から、 128bit 分の多様性を持った乱数列を生み出す以上、
ある程度は工夫しないと、最初のうちに似たような乱数列が生み出される結果となることは予想できるので、
工夫しておくに越したことはない…筈です。定量的に考察したわけではないですが。
誰か数学に詳しい人、定量的に考察してくれませんかね?(丸投げ

// まぁ普通は自前の乱数アルゴリズムではなく Boost や C++0x のライブラリを使うべきですが
// というか初期化された後にしばらく乱数列を読み飛ばせば偏りは簡単に回避できるわけですが

初期化について試行錯誤した結果、このようになりました。
http://d.hatena.ne.jp/gintenlabo/20100930/1285859540