動的メモリ確保を行わない std::seed_seq を作ってみた。

C++0x の話。 基本的なこととかは 本の虫: C++0xの新しい乱数ライブラリ、random を参照してください。


さて、上記の記事で登場している std::seed_seq ですが、
N3126 の 26.5.7.1 Class seed_seq によると、このクラスは内部に std::vector を保持しています。


つまり、このクラスは動的メモリ管理を行っている、ということで、これは微妙に効率が悪いです。
といっても、通常は このことが直ちにパフォーマンスの悪化につながるわけではなく*1
むしろ普通の用途には動的メモリ確保したほうが色々と扱いやすいので、特に文句はないのですが、

それは別として、動的メモリ確保をしないようなバリエーションがあっても悪くないよね、ってことで、
規格における seed sequence の要件を満たしたクラスを作る練習も兼ね、自分で作ってみることにします。


どのような設計にするかですが、
まず動的メモリ確保をしないので、 std::array を参考に、

template<std::size_t N>
class seed_array;

的な感じのクラスでいいでしょう。

で、実際の生成処理は、規格に書いてある std::seed_seq の動作と同じものにします。
つまり、 std::distance( first, last ) == N の場合、

std::seed_seq sseq1( first, last );
seed_array<N> sseq2( first, last );

std::mt19937 engine1( sseq1 );
std::mt19937 engine2( sseq2 );

というコードによって初期化されたエンジンが engine1 == engine2 となるようにします*2

問題は std::distance( first, last ) が N 以外の場合ですが、
その場合は適当に 0 で補うなり xor で合成するなりすれば良いでしょう。
規格を読むと、 seed sequence のコンストラクタで与えられたデータは 適当に読み捨てていいようですしね。


…で、そんな感じで作った俺々 seed sequence が、こちらです:

#include <cstddef>
#include <cstdint>
#include <iterator>
#include <type_traits>
#include <algorithm>
#include <initializer_list>

namespace etude {
  
  // std::seed_seq の動的確保しない版
  template<std::size_t N>
  struct seed_array
  {
    static_assert( N > 0, "N > 0" );
    
    typedef std::uint_least32_t result_type;
    
    // 構築
    seed_array()
      : v() {}
    
    template<class T>
    seed_array( std::initializer_list<T> il ) {
      init_( il.begin(), il.end() );
    }
    template<class Iter>
    seed_array( Iter first, Iter last ) {
      init_( first, last );
    }
    
    // 引数の個数を返す
    static std::size_t size(){ return N; }
    // 与えられた引数を出力
    template<class Iterator>
    void param( Iterator dest ) const {
      std::copy( &v[0], &v[N], dest );
    }
    
    // 生成本体
    template<class Iter>
    void generate( Iter first, Iter last ) const
    {
      typedef typename std::iterator_traits<Iter>::value_type value_type;
      static_assert(
        std::is_unsigned<value_type>::value && sizeof(result_type) <= sizeof(value_type),
        "value_type shall denote an unsigned integer type"
        "capable of accommodating 32-bit quantities."
      );
      
      if( first == last ){ return; }
      
      std::size_t const s = N;
      std::size_t const n = last - first;
      value_type const mask = result_type(0xffffffffu);
      auto const T = []( value_type x ){ return x ^ ( x >> 27 ); };
      
      // 初期化
      std::fill( first, last, static_cast<value_type>(0x8b8b8b8bu) );
      
      // パラメータ設定
      std::size_t const t = ( n >= 623 ) ? 11 :
                            ( n >=  68 ) ?  7 :
                            ( n >=  39 ) ?  5 :
                            ( n >=   7 ) ?  3 :
                                ( n - 1 ) / 2 ;
      std::size_t const p = ( n - t ) / 2;
      std::size_t const q = p + t;
      std::size_t const m = (std::max)( s + 1, n );
      
      // transform the elements of the range
      for( std::size_t k = 0; k < m; ++k )
      {
        std::size_t const k_mod_n = k % n;
        std::size_t const k_plus_p = ( k + p ) % n;
        std::size_t const k_plus_q = ( k + q ) % n;
        std::size_t const k_minus_1 = ( k + n - 1 ) % n;
        
        value_type const arg = first[k_mod_n] ^ first[k_plus_p] ^ first[k_minus_1];
        // arg &= mask; // xor なので要らない
        
        value_type const r1 = ( 1664525u * T( arg ) ) & mask;
        
        value_type const r2 = ( r1 + (
          ( k == 0 ) ? s :
          ( k <= s ) ? k_mod_n + v[k-1] :
                       k_mod_n
        ) ) & mask;
        
        first[k_plus_p] = ( first[k_plus_p] + r1 ) & mask;
        first[k_plus_q] = ( first[k_plus_q] + r2 ) & mask;
        first[k_mod_n] = r2;
      }
      
      std::size_t const m_plus_n = m + n;
      for( std::size_t k = m; k < m_plus_n; ++k )
      {
        std::size_t const k_mod_n = k % n;
        std::size_t const k_plus_p = ( k + p ) % n;
        std::size_t const k_plus_q = ( k + q ) % n;
        std::size_t const k_minus_1 = ( k + n - 1 ) % n;
        
        value_type const arg =
          ( first[k_mod_n] + first[k_plus_p] + first[k_minus_1] ) & mask;
        
        value_type const r3 = ( 1566083941u * T( arg ) ) & mask;
        
        value_type const r4 = ( r3 - k_mod_n ) & mask;
        
        // xor なので mask は不要
        first[k_plus_p] ^= r4;
        first[k_plus_q] ^= r3;
        first[k_mod_n] = r4;
      }
    }
    
    // コピーはなし
    seed_array( seed_array const& ) = delete;
    void operator=( seed_array const& ) = delete;
    
   private:
    result_type v[N];
    
    template<class Iter>
    void init_( Iter first, Iter last ) {
      // まず v を埋める
      for( std::size_t i = 0; i < N; ++i, ++first ) {
        // 範囲が足りないようなら残りは 0 で埋める
        if( first == last ) {
          for( ; i < N; ++i ) {
            v[i] = 0;
          }
          return;
        }
        
        // コピー
        v[i] = *first;
      }
      
      // 余るようなら xor で合成していく
      for( std::size_t i = 0; first != last; ++first ) {
        v[i] ^= *first;
        if( ++i >= N ){ i -= N; }
      }
    }
    
  };
  
} // namespace etude

#include <random>
#include <cassert>

int main()
{
  // 細かいテストは日記に載せる都合上カット
  std::seed_seq sseq1 = { 1, 2, 3, 4 };
  etude::seed_array<4> sseq2 = { 1, 2, 3, 4 };
  
  std::mt19937 engine1( sseq1 );
  std::mt19937 engine2( sseq2 );
  assert( engine1 == engine2 );
}

このクラスは標準規格における seed sequence の要件を満たしている…はずです。

まぁ gcc-4.5.0 だと assert で失敗するんですけどね*3
原因は gcc-4.5.0 の実装では x ^ ( x >> 27 ) とすべき部分で x ^ ( x << 27 ) としている為。
こちらの実装を gcc に合わせて変更すれば、きちんと assert が通ることは確認済みです。

追記: 右シフトと左シフトを間違えてたので修正(現状が正しいです)。

*1:なんといっても乱数エンジンの初期化は最初に一回行えば済む処理なので

*2:実際には gcc4.5.0 の std::seed_seq の実装は、規格によって規定されるものと異なるのですが、今回は規格の方に合わせます

*3:余談ですが、失敗する assert は assert じゃないので、このコードはあまり良いとは言えません