[解決済み] なぜ新しい random ライブラリは std::rand() よりも優れているのですか？

質問

という講演を見たのですが rand() は有害であると考えられている という講演があり、乱数生成のエンジン分布のパラダイムを使うことを提唱し、単純な std::rand() プラス係数パラダイムを使用することを提唱しました。

しかしながら、私は std::rand() の失敗を直接見てみたかったので、簡単な実験をしてみました。

1. 基本的に、私は2つの関数を書きました getRandNum_Old() と getRandNum_New() を使って0から5までの乱数を生成しています。 std::rand() と std::mt19937 + std::uniform_int_distribution をそれぞれ使用します。
2. 次に、quot;old" を使って 96万（6で割り切れる）個の乱数を生成し、0～5の数字の頻度を記録しました。そして、これらの頻度の標準偏差を計算しました。私が求めているのは、分布が本当に一様であればそうなるように、標準偏差をできるだけ小さくすることです。
3. そのシミュレーションを 1000 回実行し、各シミュレーションの標準偏差を記録しました。また、かかった時間もミリ秒単位で記録しました。
4. その後、まったく同じことをもう一度行いましたが、今度は "新しい" 方法で乱数を生成しました。
5. 最後に、古い方法と新しい方法の両方について、標準偏差のリストの平均と標準偏差を計算し、古い方法と新しい方法の両方について、かかった時間のリストの平均と標準偏差を計算しました。

以下はその結果です。

[OLD WAY]
Spread
       mean:  346.554406
    std dev:  110.318361
Time Taken (ms)
       mean:  6.662910
    std dev:  0.366301

[NEW WAY]
Spread
       mean:  350.346792
    std dev:  110.449190
Time Taken (ms)
       mean:  28.053907
    std dev:  0.654964

意外なことに、ロールの総延長はどちらの方法でも同じだった。すなわち std::mt19937 + std::uniform_int_distribution は、単純な std::rand() + % . もうひとつ観察したのは、新しい方法は古い方法よりも4倍ほど遅いということです。全体として、品質がほとんど向上しないのに、速度に大きなコストを払っているように思えました。

私の実験には何らかの欠陥があるのでしょうか。それとも std::rand() は本当にそんなに悪くないし、もしかしたらもっと良いのかもしれませんね？

参考までに、私が使用したコードをそのまま掲載します。

#include <cstdio>
#include <random>
#include <algorithm>
#include <chrono>

int getRandNum_Old() {
    static bool init = false;
    if (!init) {
        std::srand(time(nullptr)); // Seed std::rand
        init = true;
    }

    return std::rand() % 6;
}

int getRandNum_New() {
    static bool init = false;
    static std::random_device rd;
    static std::mt19937 eng;
    static std::uniform_int_distribution<int> dist(0,5);
    if (!init) {
        eng.seed(rd()); // Seed random engine
        init = true;
    }

    return dist(eng);
}

template <typename T>
double mean(T* data, int n) {
    double m = 0;
    std::for_each(data, data+n, [&](T x){ m += x; });
    m /= n;
    return m;
}

template <typename T>
double stdDev(T* data, int n) {
    double m = mean(data, n);
    double sd = 0.0;
    std::for_each(data, data+n, [&](T x){ sd += ((x-m) * (x-m)); });
    sd /= n;
    sd = sqrt(sd);
    return sd;
}

int main() {
    const int N = 960000; // Number of trials
    const int M = 1000;   // Number of simulations
    const int D = 6;      // Num sides on die

    /* Do the things the "old" way (blech) */

    int freqList_Old[D];
    double stdDevList_Old[M];
    double timeTakenList_Old[M];

    for (int j = 0; j < M; j++) {
        auto start = std::chrono::high_resolution_clock::now();
        std::fill_n(freqList_Old, D, 0);
        for (int i = 0; i < N; i++) {
            int roll = getRandNum_Old();
            freqList_Old[roll] += 1;
        }
        stdDevList_Old[j] = stdDev(freqList_Old, D);
        auto end = std::chrono::high_resolution_clock::now();
        auto dur = std::chrono::duration_cast<std::chrono::microseconds>(end-start);
        double timeTaken = dur.count() / 1000.0;
        timeTakenList_Old[j] = timeTaken;
    }

    /* Do the things the cool new way! */

    int freqList_New[D];
    double stdDevList_New[M];
    double timeTakenList_New[M];

    for (int j = 0; j < M; j++) {
        auto start = std::chrono::high_resolution_clock::now();
        std::fill_n(freqList_New, D, 0);
        for (int i = 0; i < N; i++) {
            int roll = getRandNum_New();
            freqList_New[roll] += 1;
        }
        stdDevList_New[j] = stdDev(freqList_New, D);
        auto end = std::chrono::high_resolution_clock::now();
        auto dur = std::chrono::duration_cast<std::chrono::microseconds>(end-start);
        double timeTaken = dur.count() / 1000.0;
        timeTakenList_New[j] = timeTaken;
    }

    /* Display Results */

    printf("[OLD WAY]\n");
    printf("Spread\n");
    printf("       mean:  %.6f\n", mean(stdDevList_Old, M));
    printf("    std dev:  %.6f\n", stdDev(stdDevList_Old, M));
    printf("Time Taken (ms)\n");
    printf("       mean:  %.6f\n", mean(timeTakenList_Old, M));
    printf("    std dev:  %.6f\n", stdDev(timeTakenList_Old, M));
    printf("\n");
    printf("[NEW WAY]\n");
    printf("Spread\n");
    printf("       mean:  %.6f\n", mean(stdDevList_New, M));
    printf("    std dev:  %.6f\n", stdDev(stdDevList_New, M));
    printf("Time Taken (ms)\n");
    printf("       mean:  %.6f\n", mean(timeTakenList_New, M));
    printf("    std dev:  %.6f\n", stdDev(timeTakenList_New, M));
}

どのように解決するのですか？

古い実装でも大丈夫です。 rand() を使用します。 LCG これらは一般に最高のジェネレータではありませんが、通常、このような基本的なテストで失敗することはありません。

bad"のよくある失敗例 - しかし十分よくあることです。 rand() の実装はそうです。

• 低次ビットの低ランダム性。
• 短周期。
• 低い RAND_MAX ;
• 連続した抽出の間に何らかの相関がある (一般に、LCG は限られた数の超平面上にある数を生成しますが、これは何らかの方法で緩和することができます)。

の API に固有のものではありません。 rand() . 特定の実装では、xorshift-family ジェネレータを srand / rand というように、アルゴリズム的にはインターフェイスの変更なしに最新の PRNG が得られるので、あなたが行ったようなテストでは、出力に弱点があることはわかりません。

編集してください。 @R. は、正しく rand / srand インタフェースが制限されるのは srand は unsigned int を取るので、実装がその後ろに置くジェネレータは本質的に UINT_MAX に制限されます。これは確かにそうなのですが、APIを拡張して srand を取ることができます。 unsigned long long を追加したり、別の srand(unsigned char *, size_t) をオーバーロードします。

確かに、実際に問題になっているのは rand() の実際の問題は、実装の多くではなく 原則的に が、しかし

• 後方互換性。現在の多くの実装では、最適とは言えないジェネレータを使用しており、通常、パラメータの選択が不適切です。 RAND_MAX をわずか 32767 に設定しています。しかし、これは過去との互換性を壊すことになるため、簡単には変更できません。 srand は、再現可能なシミュレーションのための固定シードで、あまり幸せではないでしょう (実際、前述の実装は 80 年代半ばの Microsoft C 初期バージョン、あるいは Lattice C にまでさかのぼります)。

• 単純化されたインターフェイス。 rand() は、プログラム全体のグローバルな状態を持つ単一のジェネレータを提供します。これは多くの単純な使用例では全く問題ありませんが (そして実際に非常に便利です)、問題を提起します。

  • マルチスレッド コード: これを修正するには、グローバルなミューテックスが必要です - これは無意味にすべての速度を低下させます。 と またはスレッドローカルの状態が必要です。この最後のものは、いくつかの実装 (特に Visual C++) で採用されています。
  • グローバルな状態に影響を与えない、プログラムの特定のモジュールへの、quot;private" で再現可能なシーケンスが必要な場合に使用します。

最後に rand の状態になります。

• は実際の実装を指定しないので（C 標準は実装のサンプルを提供するだけ）、異なるコンパイラ間で再現可能な出力（またはある既知の品質の PRNG を期待する）を生成することを意図するプログラムは、独自の生成器をロールバックしなければなりません。
• は適切な種を得るためのクロスプラットフォームな方法を提供しません ( time(NULL) は十分な粒度がなく、しばしば - RTC のない組み込みデバイスを考えると - 十分なランダム性さえないためです)。

そのため、新しい <random> ヘッダは、この混乱を修正しようとするもので、アルゴリズムを提供するものです。

• 完全に指定された (したがって、クロスコンパイラで再現可能な出力と保証された特性 (たとえば、ジェネレータの範囲) を持つことができます)。
• 一般的に最先端の品質である ( ライブラリが設計された時点の 下記参照)。
• クラスにカプセル化されている (そのため、グローバルな状態が強制されることはなく、スレッディングや非局所性の問題を完全に回避できる)。

... そして、デフォルトの random_device もシードするようにしました。

さて、もし私に言わせれば、私は も のために、この上に構築されたシンプルな API が欲しいです (Python がどのように "complicated" API を提供するかと同様に、些細な random.randint ビンゴカードの数字を抽出するたびにランダムなデバイス/エンジン/アダプタ/その他に溺れたくない、複雑ではない私たちのために、グローバルで事前にシードされたPRNGを使用します）、しかし現在の設備の上に自分で簡単に構築できることは事実です（一方、単純なものの上にquot; full" APIを構築することは不可能でしょう）。

最後に、パフォーマンスの比較に戻ります。他の人が指定したように、あなたは高速な LCG とより遅い (しかし一般的にはより高品質と考えられる) Mersenne Twister を比較しています。LCG の品質に問題がない場合は std::minstd_rand の代わりに std::mt19937 .

確かに、この関数に手を加えて std::minstd_rand を使うようにし、初期化に無駄な静的変数を使わないようにすると

int getRandNum_New() {
    static std::minstd_rand eng{std::random_device{}()};
    static std::uniform_int_distribution<int> dist{0, 5};
    return dist(eng);
}

9ms (古い) vs 21ms (新しい)です。 dist を取り除き (これは古典的なモジュロ演算子と比較して、入力範囲の倍数でない出力範囲の分布の歪みを処理します)、あなたが行っていることに戻ると getRandNum_Old()

int getRandNum_New() {
    static std::minstd_rand eng{std::random_device{}()};
    return eng() % 6;
}

6ms (つまり30%速い) になりました。 rand() , std::minstd_rand の方がインライン化しやすい。

ちなみに、同じテストを手打ちの（でも標準ライブラリのインターフェースにかなり準拠した） XorShift64* と比べて2.3倍速くなりました。 rand() (3.68 ms vs 8.61 ms) です。Mersenne Twister や提供されている様々な LCG とは異なり は現在のランダム性テストスイートを見事にクリアしています。 と であり、なぜまだ標準ライブラリに含まれていないのか不思議なくらいに高速です。