Cのrand()よりmt19937の方が速いことがあるという話

おはようございます。2年ぶりの記事ですね。
もう1月程前になってしまいましたが、id:sleepy_yoshi:20130720　で id:sleepy_yoshi さんが高速な非復元抽出をやっておられ、その中で、Cのrand関数を使っておられました。僕は、普段、std::mt19937を使っていたので、ちょっと比較してみた、という記事です。

C++11では、大別して、2つの擬似乱数生成の方法があります。１つはC（cstdlib）のrand関数で、高速ですが乱数の質が低く、もう１つはrandomヘッダのmt19937（メルセンヌ・ツイスタ）で、低速ですが乱数の質が高い（科学実験に適する）と、一般には思われていると思います。この高速・低速ですが、mt19937を使うことがボトルネックになるほど遅いことは殆どない、というのが今までの実感でした。なので、僕は、非復元抽出のような処理では、特にボトルネックにならない限り、mt19937の方を用いてきました。擬似乱数の質に起因する数値的なバグが万が一出たりすると、デバッグが困難、というのがその理由です。

が、ちゃんと速度比較をしたことがなかったので、比較をしてみました。すると…定説に反して、Cのrand()よりmt19937の方が、使い方によっては速い、という結果になりました。

「使い方」というのはこういう事です。今、[0,1)の範囲の実数の乱数が得られるとします。これをそのまま使う事もありますが、よくあるのは、ランダムに要素を選ぶために、[0,N)のような範囲の整数の乱数が欲しい、という場合です。このような整数の乱数を、実数の乱数から得る時によく使われているテクニックとしては、N * [0,1)の範囲の実数、を計算した上、整数型に変換して、[0,N)の範囲の乱数を求める、という方法があります。一方、mt19937を使う場合は、C++11ではそのような変換は必要なく、std::uniform_int_distributionという専用の関数があります。
比較してみました。

ソースコード(github)

#include <iostream>
#include <random>
#include <chrono>
#include <cstdlib>

namespace{
  class mtrand_t{
    std::mt19937 mt_;
    public:
    inline std::string name() const{ return "mt19937"; }
    inline double      zeroonerand() const { std::uniform_real_distribution<double> gen(0.0, 1.0); return gen(mt_);} 
    inline size_t      nrand(size_t n_rand) const { std::uniform_int_distribution<size_t> gen(0, n_rand-1); return gen(mt_);} 
    inline void        seed(int seed) { mt_.seed(seed); } 
  };
  class crand_t{
    public:
    inline std::string name() const{ return "std::rand"; }
    inline double      zeroonerand() const { return (((double) std::rand()-1)/((double)RAND_MAX)); } 
    inline size_t      nrand(size_t n_rand) const{ return n_rand * zeroonerand();}
    inline void        seed(int seed) { std::srand(seed); }
  };

  class zeroone_caller_t{
    public:
    template<class P>
    inline void call(const P randgen) const{ randgen->zeroonerand( ); }
  };
  class nrand_caller_t{
    size_t n_rand_;
    public:
    template<class P>
    inline void call(const P randgen) const{ randgen->nrand( n_rand_); }
    inline size_t n_rand()const{ return n_rand_;}
    inline void set_n_rand(size_t n_rand){ n_rand_ = n_rand; }
    nrand_caller_t(size_t n_rand): n_rand_(n_rand) {}
  };
  template<class T, class CALLER>
  void measure_time(T* randgen, const CALLER& caller, size_t n_iter, int seed){
    randgen->seed(seed);
    auto start = std::chrono::system_clock::now();
    for(size_t i = 0 ; i < n_iter ; ++i){  caller.call(randgen);  }
    auto d = std::chrono::system_clock::now() - start;
    std::cout << randgen->name() << ":\t" << std::chrono::duration<double>(d).count() << std::endl;
  }
}// namespace

int main(int argc, char* argv[]){
  size_t n_iter = 1000*1000*1000;
  size_t n_rand = 1000;
  zeroone_caller_t zeroone_caller;
  crand_t crand;
  measure_time(&crand, zeroone_caller, n_iter,  7);
  mtrand_t mtrand;
  measure_time(&mtrand, zeroone_caller, n_iter, 7);

  nrand_caller_t nrand_caller(n_rand);
  measure_time(&crand, nrand_caller, n_iter,  7);
  measure_time(&mtrand, nrand_caller, n_iter, 7);
}

実験結果：単位（秒）

std::rand: 7.17306
mt19937: 8.87919
std::rand: 7.18293
mt19937: 4.70852

ソースコードでは、1セットにつき、1000^3=1G個=10億個の乱数を生成するのに必要な時間を測っています。std::randがCのrand関数、mt19937が、std::19937です。最初の2行が、[0,1)の範囲の実数の乱数生成の実行時間、次の2行が、[0,N)の範囲の整数の乱数生成の実行時間です。[0,N)範囲の整数の乱数生成の実行時間は、mt19937の方が約2.5秒も短い事が分かりました。

もちろん、実装・環境依存だと思います。が、いずれにしても数秒/1G回程度の速度差を理由に、mt19937の使用を回避しなければならない場合は殆どないはずです。従って、殆どの場合で、mt19937の使用が勧められるという結論になりました。気になる人は、ソースコードをコンパイルして実行して、自分の環境でどうなのか、試してみて下さい。乱数の質は、mt19937の方が明らかに高いです。また、最近では、SIMD命令を使った高速化版のSFMTも出ています。

なお、rand()は元々[0,RAND_MAX]の整数を返してくるのに、(((double) std::rand()-1)/*1;のように、わざわざ割り算を使ってdoubleに変換いるから遅くなるのだ、整数のままmodulo演算(%)を使って、rand()%Nとすればいいじゃないか、という人がいるかも知れません。が、その場合、modulo biasというのが出ます。仮に、RAND_MAX=4、N=3だとすると、RAND_MAXがNで割り切れないので、rand()%3では、0,1は2/5の確率で、2は1/5の確率で出る事になり、0,1,2が等確率で出なくなってしまいます。一方、同じrand()を使う方法でも、挙げたコードのように、割り算を使って一回doubleに変換すれば、そのような事にはなりません。

2013/08/19 1:14追記

rand()を使って[0,1)の範囲の乱数を作るのには、もっと深い話があるようです id:takeda25 さんが関連記事を書いてくださりました。
いつからその方法で偏りのない乱数が得られると錯覚していた？ - アスペ日記

また、id:yatt さん他何人かの方から指摘されましたが、そもそも、mt19937のページに、「Cの標準ライブラリのrand()より高速なこともあります。」と書かれているようです。

*1:double)RAND_MAX