C++11の乱数ジェネレーターを自作してみた

この記事は、C++ (fork) Advent Calendar 2013 19日目の記事になります。
近々C++14が標準化されようという時期に、今更感が拭えませんが
C++11から、C言語由来のrand()の他に新しい乱数ライブラリが導入されたので。
自前の乱数ジェネレータを作ってみました。

この記事は、江添さんのブログを参考に書いています。
詳しい実装方法に付いては、江添さんのブログを読んでください、僕には説明できません。
本の虫

使い方自体は簡単で
乱数そのものを生成するengineクラスを、
乱数の範囲を調節するdistributionクラスに引数として渡すだけです

ソースコード

#include <random>

#include <QDebug>

int main(int argc, char *argv[])

{

    //乱数生成クラス

    std::mt19937 engine;

    //範囲調節クラス 0から100までの乱数を作る

    std::uniform_int_distribution<int> distribution(0,100) ;

    qDebug() << distribution(engine);

    return 0;

}

今回はこの std::mt19937 engine の代わり使えるクラスを実装します。

まずヘッダ定義

ヘッダ定義

#ifndef RDTSC_RANDOM_GENERATOR_H

#define RDTSC_RANDOM_GENERATOR_H

#include <climits>

#include <chrono>

#include <QDebug>

#include <random>

//乱数エンジン

class rdtsc_random_generator

{

public:

    rdtsc_random_generator();

    //乱数を格納する型

    using result_type = unsigned int;

    //乱数の最小値

    static constexpr result_type min(){ return 0 ;}

    //乱数の最大値

    static constexpr result_type max(){ return UINT_MAX;}

    //関数オブジェクトの呼び出し

    result_type operator () () {return randam();}   

    //タイマー読み出し

    unsigned long long timer();

    //指定した回数乱数を空読みする

    void discard(unsigned long long ppp);

private:

    //乱数生成関数

    result_type randam();

    //初期化

    inline void Init();

    //xorshit

    result_type x,y,z,w;

    //再初期化カウント

    int count;

};

#endif // RDTSC_RANDOM_GENERATOR_H

次に定義

関数定義

#include "rdtsc_random_generator.h"

rdtsc_random_generator::rdtsc_random_generator() :

    //初期化

    x(875456789), y(62436069), z(21288629), w(88675123)

{

    //タイマーで初期化

    Init();

    //かたよらないように 数回空読み

    discard(3);

    count=0;

}

//タイマーを使って初期化

inline void rdtsc_random_generator::Init()

{

    count=0;

    const unsigned int time=static_cast<unsigned int>(timer());

    x=x^time;

    //ただXORするのも芸がないので 適当に数値をいじる。

    unsigned int tmp1,tmp2;

    //28bit右にシフト

    tmp2=time>>28;

    //4bit左にシフト

    tmp1=time<<4;

    //合成

    tmp1=tmp1+tmp2;

    y=y^tmp1;

    //28bit右にシフト

    tmp2=tmp1>>28;

    //4bit左にシフト

    tmp1=tmp1<<4;

    //合成

    tmp1=tmp1+tmp2;

    z=z^tmp1;

    //28bit右にシフト

    tmp2=tmp1>>28;

    //4bit左にシフト

    tmp1=tmp1<<4;

    //合成

    tmp1=tmp1+tmp2;

    w=w^tmp2;

}

//タイマー読み出し

unsigned long long rdtsc_random_generator::timer()

{

    /*

    unsigned int a,b;

    asm("rdtsc":"=a"(a),"=d"(b));

    unsigned long long out=b;

    out=out<<32;

    out+=a;

    return out;

    */

    auto clock=std::chrono::high_resolution_clock::now();

    auto time=clock.time_since_epoch();

    auto count = std::chrono::duration_cast<std::chrono::nanoseconds>(time).count();

    return count;

}

//乱数発生 

rdtsc_random_generator::result_type rdtsc_random_generator::randam()

{

    //10回呼び出されたら タイマーで再初期化する

    if(count>10){

        Init();

    }else{

        count++;

    }

    //xorshift

    result_type t;

    t=x^(x<<11);

    x=y;

    y=z;

    z=w;

    w=(w^(w>>19))^(t^(t>>8));

    return w;

}

//空読み

void rdtsc_random_generator::discard(unsigned long long ppp)

{

    for(unsigned long long i=0;i<ppp;i++){

        randam();

    }

}

ソースコード

#include <QDebug>

#include "rdtsc_random_generator.h"

int main(int argc, char *argv[])

{

    //乱数生成クラス

    rdtsc_random_generator engine;

    //範囲調節クラス 0から100までの乱数を作る

    std::uniform_int_distribution<int> distribution(0,100) ;

    qDebug() << distribution(engine);

    return 0;

}

C++１１の<chrono>で取得したナノ秒単位の時刻を乱数の種に、
Xorshiftで擬似乱数を発生させるクラスになります。
コンストラクタで乱数の種を設定するので、手動で初期化する手間が省けます 。

高速が売りのXorshiftなので、ベンチマークを取ってみる。

ということで、標準の乱数ジェネレータと速度比較します。
乱数を10000回生成するのにかかった時間を比較するのですが、
単にループで回すだけだと、当然最適化で全部消されます。
最適化無しでの速度比較は、あんまり実体に即してない気がするのでやりません、
なので生成した乱数を配列に書き込み、
処理時間の出力時にランダム一つ出力することで、過度の最適化を回避します。
時間の単位はナノ秒です。

ソースコード

#include <random>

#include <QDebug>

#include "rdtsc_random_generator.h"

#include <QVector>

#include <chrono>

int main(int argc, char *argv[])

{

    //配列の長さを設定

    const int row=10000;

    //乱数生成

    std::uniform_int_distribution<int> distribution(0,1000);

    //最適化回避用

    std::uniform_int_distribution<int> distribution2(0,row-1);

    //時間計測開始

    auto start=std::chrono::high_resolution_clock::now();

    //標準乱数生成クラス

    std::mt19937 engine;

    QVector<int> a1(row),a2(row);

    for(int i=0;i<row;i++){

        a1[i]=distribution(engine);

    }

    //計測完了

    auto end=std::chrono::high_resolution_clock::now();

    auto time=end-start;

    qDebug() << "標準乱数生成クラス 処理時間" << time.count() << "最適化回避" << a1[distribution2(engine)];

    //時間計測開始

    auto start2=std::chrono::high_resolution_clock::now();

    //自作乱数生成クラス

    rdtsc_random_generator myengine;

    for(int i=0;i<row;i++){

        a2[i]=distribution(myengine);

    }

    //計測完了

    auto end2=std::chrono::high_resolution_clock::now();

    auto time2=end2-start2;

    qDebug() << "自作乱数生成クラス 処理時間" << time2.count() << "最適化回避"<< a2[distribution2(myengine)];

    return 0;

}

実行

出力

標準乱数生成クラス 処理時間 392787 最適化回避 965 

自作乱数生成クラス 処理時間 1037074 最適化回避 859 

ありゃりゃ
とまあ、3倍近い速度差が出てしまいました。
流石標準ライブラリは速いですね、
皆さんも、標準ライブラリにある機能は、自前で再実装せず標準ライブラリの機能を活用しましょう。














とまあ、ここで締めてもいいのですが、

カンの良い方なら気付いていると思いますが上の乱数ジェネレータが遅い理由は、
randam()が10回呼び出される度に、乱数の種の初期化を行うからです。
c++の<chrono>は地味に遅いのです。
なので速度が必要な場合、初期化の頻度を減らすか、消し去るかすれば大丈夫です。

関数定義

//乱数発生 

rdtsc_random_generator::result_type rdtsc_random_generator::randam()

{

    //10回ではなく100回にする

    if(count>100){

        Init();

    }else{

        count++;

    }

    //xorshift

    result_type t;

    t=x^(x<<11);

    x=y;

    y=z;

    z=w;

    w=(w^(w>>19))^(t^(t>>8));

    return w;

}

出力

標準乱数生成クラス 処理時間 396279 最適化回避 965 

自作乱数生成クラス 処理時間 245562 最適化回避 737

こんな感じで初期化の頻度を落とすと、標準ライブラリのstd::mt19937より速くなります。
std::mt19937の擬似乱数アルゴリズムメルセンヌ・ツイスタなので。
Xorshiftを使う自作の乱数ジェネレータの方が速いのは当然ではありますけどね、

さて、
乱数ジェネレータの名前で感づいた人も多いのでしょうが、
このジェネレータの名前はrdtsc_random_generatorです。
 chronoでもなく、high_resolution_clockでもなく、rdtscです。
つまりtscで初期化する環境依存しまくりのオレオレ乱数ジェネレータを記事のネタの為に書き換えたものです。
名前と実体が一致していませんが、今更クラスの名前を変えるのも面倒くさいので、
tsc版のソースも公開します。

関数定義

//rdtscを取得

unsigned long long rdtsc_random_generator::timer()

{

    unsigned int a,b;

    asm("rdtsc":"=a"(a),"=d"(b));

    unsigned long long out=b;

    out=out<<32;

    out+=a;

    return out;

}

出力

標準乱数生成クラス 処理時間 390203 最適化回避 965 

自作乱数生成クラス 処理時間 153232 最適化回避 94 

std::mt19937より、ずっとはやい!!

ですが、環境依存しまくりの外法です、ご利用は計画的に、

インラインアセンブリを使っているので、x86とx64限定かつコンパイラにも依存します。

gccとclangでは動きましたがMSVCでは試してません。

まあ、途中で初期化するとか中途半端なことやらなければいい話です
門外漢の浅知恵なので有効かどうかは知りません、
追記、x,y,z,w全てが0にならないかぎり限り、何で初期化してもいいみたいです

それとstd::chrono::high_resolution_clock::now()はtscを直に読むより二桁ぐらい遅いみたいです。
 追記、gccのchronoの実装を確認したところ、内部でシステムコールの
 clock_gettime(CLOCK_REALTIME, &hoge)とか
 syscall(SYS_clock_gettime, CLOCK_REALTIME, &hoge)あたりを呼び出しているらしい。
 プロセス切替えのオーバヘッドが大きいいらしく、ナノ秒精度の結果を返すが、読み出しに1マイクロ秒程度かかる。

ついでなので、c++11標準の乱数ジェネレータの速度比較もします。
条件は今までと同様に、乱数を10000回生成するのにかかった時間で比較します。
時間の単位はマイクロ秒。使用したコンパイラはｇｃｃです。

処理時間の比較

乱数ジェネレータ	処理時間 μs	注釈
xorshift	153	自作
mt19937	393	メルセンヌ・ツイスタ
minstd_rand0	355	線形合同法
minstd_rand	363	線形合同法の改良版
mt19937_64	384	メルセンヌ・ツイスタの64ビット版
ranlux24	1262	RANLUX法のレベル3
ranlux48	3907	RANLUX法のレベル4
knuth_b	401	KnuthのリオーダーアルゴリズムB

基本的にはstd::mt19937で充分だと思います。