为什么新的随机库比 std::rand() 更好？

Question

所以我看到一个名为 rand() 被认为有害 的演讲，它提倡使用随机数生成的引擎分布范式，而不是简单的 std::rand() 加模范式。

但是，我想亲眼看看 std::rand() 的失败之处，所以我做了一个快速实验：

1. Basically, I wrote 2 functions getRandNum_Old() and getRandNum_New() that generated a random number between 0 and 5 inclusive using std::rand() and std::mt19937 + std::uniform_int_distribution 分别。
2. 然后我使用“旧”方式生成了 960,000 个（可被 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));
}

原文由 rcplusplus 发布，翻译遵循 CC BY-SA 4.0 许可协议

Answer 1

几乎所有“旧”的实现 rand() 都使用 LCG ；虽然它们通常不是最好的生成器，但通常你不会看到它们在这样的基本测试中失败 - 即使是最差的 PRNG，平均偏差和标准偏差也通常是正确的。

“坏”的常见故障 - 但足够常见 - rand() 实现是：

• 低位随机性低；
• 短期内;
• 低 RAND_MAX ；
• 连续提取之间的一些相关性（通常，LCG 产生的数字位于有限数量的超平面上，尽管这可以通过某种方式减轻）。

尽管如此，这些都不是特定于 rand() 的 API。一个特定的实现可以在 srand / rand 后面放置一个 xorshift 系列生成器，并且从算法上讲，在不更改接口的情况下获得最先进的 PRNG，因此没有像这样的测试您确实会在输出中显示出任何弱点。

编辑：@R。 correctly notes that the rand / srand interface is limited by the fact that srand takes an unsigned int , so any generator an implementation may put它们后面本质上仅限于 UINT_MAX 可能的起始种子（以及因此生成的序列）。确实如此，尽管 API 可以简单地扩展以使 srand 采用 unsigned long long 或添加单独的 srand(unsigned char *, size_t) 过载。

---

实际上， rand() 的实际问题 原则上并没有太多的实现， 但是：

• 向后兼容性；许多当前的实现使用次优生成器，通常使用错误选择的参数；一个臭名昭著的例子是 Visual C++，它的 RAND_MAX 只有 32767。然而，这不能轻易改变，因为它会破坏与过去的兼容性——人们使用 srand 可重复模拟的种子不会太高兴（实际上，IIRC 上述实现可以追溯到 Microsoft C 早期版本 - 甚至是 Lattice C - 从八十年代中期开始）；

• 简单的界面； rand() 为整个程序提供具有全局状态的单个生成器。虽然这对于许多简单的用例来说非常好（实际上非常方便），但它带来了问题：

  • 使用多线程代码：要修复它，您要么需要一个全局互斥锁 - 这会无缘无故减慢一切 并 杀死任何可重复性的机会，因为调用序列本身变得随机 - 或线程本地状态；最后一个已被多个实现（尤其是 Visual C++）采用；
  • 如果您想要将“私有”、可重现的序列放入程序的特定模块中，并且不会影响全局状态。

最后， rand 状态：

• 没有指定实际的实现（C 标准仅提供示例实现），因此任何旨在跨不同编译器产生可重现输出（或期望具有某种已知质量的 PRNG）的程序都必须滚动其自己的生成器；
• 没有提供任何跨平台方法来获得一个像样的种子（ time(NULL) 不是，因为它不够精细，而且通常 - 认为没有 RTC 的嵌入式设备 - 甚至不够随机）。

因此，新的 <random> 标头试图修复这种混乱，提供以下算法：

• 完全指定（因此您可以获得交叉编译器可重现的输出和保证的特性 - 例如，生成器的范围）；
• 通常具有最先进的质量（ _从设计图书馆时开始_；见下文）；
• 封装在类中（所以没有全局状态被强加给你，这完全避免了线程和非局部性问题）；

…以及默认的 random_device 以及播种它们。

现在，如果你问我，我 也 希望在此之上构建一个简单的 API，用于“简单”、“猜数字”的情况（类似于 Python 提供“复杂”API 的方式，但也很简单 random.randint & Co. 为我们那些不想被随机设备/引擎/适配器/任何我们每次想为宾果卡提取数字时不被淹没在随机设备/引擎/适配器/任何东西中的简单的人使用全球预播 PRNG），但确实，您可以在当前设施上轻松地自己构建它（而在简单的设施上构建“完整”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);
}

我得到 9 毫秒（旧）和 21 毫秒（新）；最后，如果我摆脱 dist （与经典的模运算符相比，它处理输出范围的分布偏斜而不是输入范围的倍数）并回到你正在做的事情 getRandNum_Old()

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

我把它降低到 6 毫秒（所以，快 30%），可能是因为，不像调用 rand() ， std::minstd_rand 更容易内联。

---

顺便说一句，我使用手动滚动（但几乎符合标准库接口） XorShift64* 进行了相同的测试，它比 rand() 快 2.3 倍（3.68 ms vs 8.61 ms );鉴于这一点，与 Mersenne Twister 和 提供的各种 LCG 不同，它 通过了当前的随机测试套件， 而且速度非常快，这让你想知道为什么它还没有包含在标准库中。

原文由 Matteo Italia 发布，翻译遵循 CC BY-SA 4.0 许可协议

Answer 2

如果您以大于 5 的范围重复您的实验，那么您可能会看到不同的结果。当您的范围明显小于 RAND_MAX 对于大多数应用程序来说没有问题。

例如，如果我们的 RAND_MAX 为 25，那么 rand() % 5 将产生具有以下频率的数字：

 0: 6
1: 5
2: 5
3: 5
4: 5

由于 RAND_MAX 保证大于 32767，并且最不可能和最可能之间的频率差异仅为 1，因此对于大多数用例而言，对于小数字，分布足够接近随机。

原文由 Alan Birtles 发布，翻译遵循 CC BY-SA 4.0 许可协议