初めに
背景
元々、C++でのメルセンヌツイスタ実装としてはboostのものがありましたが、乱数機能が貧弱であった標準ライブラリでは、C++11 規格で、このboostの実装を参考にしメルセンヌツイスタが導入されたようです。
※参考:本の虫：C++0xの新しい乱数ライブラリ、random
という話から、両者の実装は同じだろうと考えていたのですが、どうも違うようだ…ということで、調べてみました。
結論
ということなのですが、先に結論です。
性能に関しては、あくまで手元の環境で、の話です。
- boost::random::mt19937とg++標準のstd::mt19937の実装は少し違う
- boostの方が性能チューニングされていて速い
- boostは過去に実装の変更がされていて、この過去の実装と標準は共にメルセンヌツイスタのリファレンス実装に近い。
- とは言えboostの過去の実装よりも、標準の実装の方が速い ( 少し実装が違う )
- g++の実装は ( おそらく性能に関係ない範囲で ) 一度変更されているが、それ以外に変更はない
まあ、全て見切れている訳ではないので、あくまで観測できた範囲内で、ですが。
環境
- Windows10/WSL ( Ubuntu 16.04.2 LTS相当 )
- gcc/g++ 5.4.0 ( パッケージ版 )
- boost 1.58 ( パッケージ版 ), 1.46.1 ( 手動インストール版 )
なお、メルセンヌツイスタは調整可能なパラメータがありますが、ここではMT19937に絞って比較をしています。
メルセンヌツイスタ概要
リファレンス実装
メルセンヌツイスタの詳細についてはここでは割愛しますが、特徴としては次の通りです
- GFSR ( bit-xorによる n項間隣接漸化式により生成される数列を乱数とするもの ) に一部ひねりを加えて乱数列としての特性を改善したもの
- そのためGFSRと同様、漸化式を回すために以前の項を保存するバッファが必要。( MT19937で624ワード )
- 乱数を得る度に漸化式を回すのは効率が悪いためか、バッファの個数分の乱数を溜めて置き、なくなったところで一気に漸化式を回してバッファを刷新する
- バッファに保存された値を更に「調律」( tempering ) した値が最終的な乱数 ( ワードの取り得る整数値に対する均一分布乱数 ) となる
なお、具体的な実装については http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/CODES/mt19937ar.c のgenrand_int32()に見ることができます。
C++での実装概要
それでは続いてC++(boost/g++標準)での現時点での実装についてです
※過去の実装では異なる点もある
- メルセンヌツイスタは、リファレンス実装をほぼほぼ踏襲した、可変パラメータによるテンプレートクラス
mersenne_twister_engineとして実装されている - そのテンプレートクラスにMT19937相当のパラメータを適用したクラスに
mt19937という名前がつけられている - 乱数生成ルーチンの実体は ( テンプレートなのでまあ妥当と言うか ) ヘッダファイル内にある。コンパイル済みの外部ライブラリ ( .a や .so ファイル ) がリンクされるわけではない。
※裏を返せば、inline化等のコンパイラ時の最適化の影響が性能に現れる - 個々の乱数ワードの取り出しは、
operator()というメンバ関数が担当する - この
operator()はバッファからの取り出しとtemperingのみを担当する - バッファの刷新は、別のメンバ関数 ( boostの場合は
twist()、g++標準の場合は_M_gen_rand()) に分離されている - このバッファの刷新の実装が性能に大きく寄与する部分である
なお、乱数生成ルーチンの実態はboostならばboost/random/mersenne_twister.hppに、g++標準ならbits/random.tccにあります。
で、operator()自体の実装には差がありません。以下にboostの方の実装を挙げます
template<class UIntType,
std::size_t w, std::size_t n, std::size_t m, std::size_t r,
UIntType a, std::size_t u, UIntType d, std::size_t s,
UIntType b, std::size_t t,
UIntType c, std::size_t l, UIntType f>
inline typename
mersenne_twister_engine<UIntType,w,n,m,r,a,u,d,s,b,t,c,l,f>::result_type
mersenne_twister_engine<UIntType,w,n,m,r,a,u,d,s,b,t,c,l,f>::operator()()
{
if(i == n)
twist();
// Step 4
UIntType z = x[i];
++i;
z ^= ((z >> u) & d);
z ^= ((z << s) & b);
z ^= ((z << t) & c);
z ^= (z >> l);
return z;
}
バッファ内の乱数を幾つ消費したかをiで管理していて、使い切っていればtwist()で刷新、取得したバッファ内の1要素に色々bit演算をかますことで乱数値をtemperingしていることが分かります。
性能/実装の差異
性能測定コード
では、boostとg++標準の差異について、まずは性能を見てみます。
比較対象は、boost1.58, boost1.46.1, g++標準です。後述しますが、boost1.46.1はメルセンヌツイスタ実装が切り替わる直前のバージョンです。
しょぼい環境なのであまり真面目には測ってなくて、乱数生成1億回のコードを、-O3 -march=nativeコンパイル、timeコマンドでCPU時間測定、位です。
コマンドとしては次の通り。
$ g++ -std=c++11 -O3 -march=native -o mtbench-std mtbench.cpp
$ g++ -std=c++11 -O3 -march=native -DUSE_BOOST -o mtbench-boost mtbench.cpp
$ g++ -std=c++11 -O3 -march=native -DUSE_BOOST -I./boost_1_46_1 -o mtbench-boost1461 mtbench.cpp
$ time ./mtbench-std
$ time ./mtbench-boost
$ time ./mtbench-boost1461
ということで測定した結果です…が。なんと、今のboostの方が5倍近い速度が出ています。ホンマかいなという感じです。逆に昔のboostは標準よりもかなり遅いです。
| ケース | CPU時間(USER) |
|---|---|
| 標準 | 5.359s |
| boost(1.58) | 1.172s |
| boost(1.46.1) | 9.656s |
ちなみに、実行したコードは次の通りです。
何と言うか、クラス名もインターフェースも同じなので、マクロUSE_BOOSTでスイッチし、実質#includeとusingの内容を変えるだけで、あとは全く同じコードとして組むことができます。
なお、乱数のseedをデフォルトにしているので、毎回同じ乱数列が生成されます。一応ダミーとしてサイコロよろしく1～6の乱数にして、その1億回分の合計を適当な数で割った余りを出力していますが、どの実装でも毎回49994605が出力されます。
#include <iostream>
#include <cstdint>
#include <random>
#ifndef USE_BOOST
using std::mt19937;
#else
# include <boost/version.hpp>
# include <boost/random.hpp>
# if BOOST_VERSION < 104700
using boost::mt19937;
# else
using boost::random::mt19937;
# endif
#endif
int main() {
constexpr int iter=100000000;
mt19937 mt;
std::uniform_int_distribution<int> rand(1,6);
int64_t sum=0;
for ( int i=0; i<iter; i++ )
sum += rand(mt);
std::cout << sum%100000007 << std::endl;
}
実装の差異
さて、では上で見た性能差はどこから出てくるか。違いがあるのは乱数バッファの刷新の部分だけです。そこで、それぞれの実装を抜粋して見てみます。
-
リファレンス実装 の
genrand_int32()のバッファ刷新処理
for (kk=0;kk<N-M;kk++) {
y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL];
}
for (;kk<N-1;kk++) {
y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
}
y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
- g++標準の
_M_gen_rand()メンバ関数
void
mersenne_twister_engine<_UIntType, __w, __n, __m, __r, __a, __u, __d,
__s, __b, __t, __c, __l, __f>::
_M_gen_rand(void)
{
const _UIntType __upper_mask = (~_UIntType()) << __r;
const _UIntType __lower_mask = ~__upper_mask;
for (size_t __k = 0; __k < (__n - __m); ++__k)
{
_UIntType __y = ((_M_x[__k] & __upper_mask)
| (_M_x[__k + 1] & __lower_mask));
_M_x[__k] = (_M_x[__k + __m] ^ (__y >> 1)
^ ((__y & 0x01) ? __a : 0));
}
for (size_t __k = (__n - __m); __k < (__n - 1); ++__k)
{
_UIntType __y = ((_M_x[__k] & __upper_mask)
| (_M_x[__k + 1] & __lower_mask));
_M_x[__k] = (_M_x[__k + (__m - __n)] ^ (__y >> 1)
^ ((__y & 0x01) ? __a : 0));
}
_UIntType __y = ((_M_x[__n - 1] & __upper_mask)
| (_M_x[0] & __lower_mask));
_M_x[__n - 1] = (_M_x[__m - 1] ^ (__y >> 1)
^ ((__y & 0x01) ? __a : 0));
_M_p = 0;
}
なお、リファレンス実装でのmag01というのは、0 もしくは MATRIX_A の2要素からなる配列であって、g++標準では配列参照の代わりに条件演算子による値の選択に変わっています。__aというのがMATRIX_Aに相当します。( これはMTの線型漸化式における行列Aに相当する定数です )
両実装は、この部分しか根本的な違いはありません。
続いてboostの実装です。
- boostの
twist()メンバ関数
void
mersenne_twister_engine<UIntType,w,n,m,r,a,u,d,s,b,t,c,l,f>::twist()
{
const UIntType upper_mask = (~static_cast<UIntType>(0)) << r;
const UIntType lower_mask = ~upper_mask;
const std::size_t unroll_factor = 6;
const std::size_t unroll_extra1 = (n-m) % unroll_factor;
const std::size_t unroll_extra2 = (m-1) % unroll_factor;
// split loop to avoid costly modulo operations
{ // extra scope for MSVC brokenness w.r.t. for scope
for(std::size_t j = 0; j < n-m-unroll_extra1; j++) {
UIntType y = (x[j] & upper_mask) | (x[j+1] & lower_mask);
x[j] = x[j+m] ^ (y >> 1) ^ ((x[j+1]&1) * a);
}
}
{
for(std::size_t j = n-m-unroll_extra1; j < n-m; j++) {
UIntType y = (x[j] & upper_mask) | (x[j+1] & lower_mask);
x[j] = x[j+m] ^ (y >> 1) ^ ((x[j+1]&1) * a);
}
}
{
for(std::size_t j = n-m; j < n-1-unroll_extra2; j++) {
UIntType y = (x[j] & upper_mask) | (x[j+1] & lower_mask);
x[j] = x[j-(n-m)] ^ (y >> 1) ^ ((x[j+1]&1) * a);
}
}
{
for(std::size_t j = n-1-unroll_extra2; j < n-1; j++) {
UIntType y = (x[j] & upper_mask) | (x[j+1] & lower_mask);
x[j] = x[j-(n-m)] ^ (y >> 1) ^ ((x[j+1]&1) * a);
}
}
// last iteration
UIntType y = (x[n-1] & upper_mask) | (x[0] & lower_mask);
x[n-1] = x[m-1] ^ (y >> 1) ^ ((x[0]&1) * a);
i = 0;
}
リファレンス実装・g++標準と比べて、やっていることはほとんど同じなのですが、
- 0か
aかの選択を、配列参照 ( リファレンス実装 )、条件演算子 ( g++標準 ) でもなく、乗算でやっている
※つまり、0→0, 1→aなら×aの計算でも同じ、ということ - 2か所あるforループの繰り返し回数を、
unroll_factor=6の倍数回の主実行部と、残りの端数部分に分割している
という違いがあるようで、ここが恐らく性能に効いているのでしょう。
…ってあれ? g++ の-O3ってわざわざループを分けても、そもそもアンローリングはしてくれなかった気が、というのと、アンローリングなんてコンパイラに自動でやらせるべきでは…? という疑問が。実際、アンローリング関係の最適化オプション-funroll-loops,-fno-unroll-loopsでも性能に違いは見られなかったので、実質条件演算子か乗算か、だけでそれだけの性能差が出ていると考えざるを得ません。
…何それ。それくらいコンパイラでなんとかしてくれないの?
最後に、一番遅かったboost旧版の実装を見てみます。
- boost1.46.1の
twist()メンバ関数
template<class UIntType, int w, int n, int m, int r, UIntType a, int u,
int s, UIntType b, int t, UIntType c, int l, UIntType val>
void mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::twist(int block)
{
const UIntType upper_mask = (~0u) << r;
const UIntType lower_mask = ~upper_mask;
if(block == 0) {
for(int j = n; j < 2*n; j++) {
UIntType y = (x[j-n] & upper_mask) | (x[j-(n-1)] & lower_mask);
x[j] = x[j-(n-m)] ^ (y >> 1) ^ (y&1 ? a : 0);
}
} else if (block == 1) {
// split loop to avoid costly modulo operations
{ // extra scope for MSVC brokenness w.r.t. for scope
for(int j = 0; j < n-m; j++) {
UIntType y = (x[j+n] & upper_mask) | (x[j+n+1] & lower_mask);
x[j] = x[j+n+m] ^ (y >> 1) ^ (y&1 ? a : 0);
}
}
for(int j = n-m; j < n-1; j++) {
UIntType y = (x[j+n] & upper_mask) | (x[j+n+1] & lower_mask);
x[j] = x[j-(n-m)] ^ (y >> 1) ^ (y&1 ? a : 0);
}
// last iteration
UIntType y = (x[2*n-1] & upper_mask) | (x[0] & lower_mask);
x[n-1] = x[m-1] ^ (y >> 1) ^ (y&1 ? a : 0);
i = 0;
}
}
今のtwist()と違い、blockという引数があります。つまり、operator()の方で、状況に応じてblockの値を変えつつtwist()を呼んでいるということです。
実は、この実装では乱数バッファの要素数が倍用意されていて、前半が埋まっていて後半を刷新する時 ( block==0 ) はforループを1つにまとめて、後半が埋まっていて前半を刷新する場合は、標準実装のようにforループを2つに分けているのです。
乗算でなくて条件演算子を使っているのは標準実装と同じです。
…なんというか、余計な小細工で性能を落としているという感じでしょうか。
終わりに
まとめ
- boostとg++標準のメルセンヌツイスタ実装は少し違う
- 違いは些細ながらboostの方が大分速い
- 昔のboostはむしろ遅い
…というか、機能は同じなんだから、実装も同期した方が良いのでは、という気もしますが。
現時点では、乱数生成がボトルネックになるようなプログラムでは、boostを積極的に使うようにした方が良いかもしれません。
もっとも、インターフェースが同じですから、別記事標準ヘッダのオーバーライドのような手法を使えば、標準の<random>を#includeしつつその実boostの方の実装を使うようにできる気もしますが。
※コードに手を入れずに外から何とかするズボラ的解決法
参考：各実装の変遷
- g++標準版
- boost旧版