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

Posted 2023-02-25

【中文标题】为啥新的随机库比 std::rand() 更好？

【英文标题】：Why is the new random library better than std::rand()?为什么新的随机库比 std::rand() 更好？

【发布时间】：2019-04-02 02:46:42

【问题描述】：

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

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

基本上，我编写了 2 个函数 getRandNum_Old() 和 getRandNum_New()，它们分别使用 std::rand() 和 std::mt19937+std::uniform_int_distribution 生成一个介于 0 和 5 之间的随机数。 然后我使用“旧”方式生成了 960,000 个（可被 6 整除）随机数，并记录了数字 0-5 的频率。然后我计算了这些频率的标准偏差。我正在寻找的是尽可能低的标准差，因为如果分布真正均匀，就会发生这种情况。 我运行该模拟 1000 次并记录每个模拟的标准偏差。我还记录了它所花费的时间，以毫秒为单位。 之后，我再次执行完全相同的操作，但这次以“新”方式生成随机数。 最后，我计算了旧方法和新方法的标准差列表的均值和标准差，以及旧方法和新方法所用时间列表的均值和标准差。

结果如下：

[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));

【问题讨论】：

这就是这个建议存在的原因。如果您不知道如何测试 RNG 是否有足够的熵或者它是否对您的程序很重要，那么您应该假设 std::rand() 不够好。 en.wikipedia.org/wiki/Entropy_(computing)

rand() 是否足够好的底线很大程度上取决于您使用随机数集合的目的。如果您需要特定类型的随机分布，那么库的实现当然会更好。如果您只需要随机数并且不关心“随机性”或产生什么类型的分布，那么rand() 就可以了。将合适的工具与手头的工作相匹配。

可能被骗：***.com/questions/52869166/… 我只是不想锤这个，所以我避免实际投票。

for (i=0; i<k*n; i++) a[i]=i%n; 产生与最好的 RNG 相同的精确平均值和标准差。如果这对您的应用程序来说足够好，只需使用此序列即可。

“尽可能低的标准偏差” - 不。那是错误的。您期望频率会有所不同 - 关于 sqrt(frequency) 是关于您期望的标准偏差。 n.m.的“递增计数器”。生产的标准差会低得多。 （并且是一个非常糟糕的rng）。

【参考方案1】：

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

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

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

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

编辑： @R. 正确指出 rand/srand 接口受到 srand 采用 @987654335 的事实的限制@，因此任何实现可能放在它们后面的生成器本质上都限于UINT_MAX 可能的起始种子（以及因此生成的序列）。确实如此，尽管 API 可以简单地扩展以使 srand 采用 unsigned long long，或添加单独的 srand(unsigned char *, size_t) 重载。

---

确实，rand() 的实际问题并不在于实现原则上，但是：

向后兼容性；许多当前的实现使用次优生成器，通常使用错误选择的参数；一个臭名昭著的例子是 Visual C++，它的 RAND_MAX 仅为 32767。但是，这不能轻易更改，因为它会破坏与过去的兼容性 - 使用带有固定种子的 srand 进行可重复模拟的人不会太高兴了（事实上，IIRC 上述实现可以追溯到 80 年代中期的 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 engstd::random_device();
    static std::uniform_int_distribution<int> dist0, 5;
    return dist(eng);

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

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

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

---

顺便说一句，我使用手动（但几乎符合标准库接口）XorShift64* 进行了相同的测试，它比rand() 快 2.3 倍（3.68 毫秒对 8.61 毫秒）；鉴于此，与 Mersenne Twister 和各种提供的 LCG 不同，它passes the current randomness test suites with flying colors 并且它的速度非常快，这让你想知道为什么它还没有包含在标准库中。

【讨论】：

这正是srand 和未指定算法的组合，导致std::rand 陷入困境。另见my answer to another question。

rand 在 API 级别基本上受到限制，因为种子（以及因此可以产生的可能序列的数量）受UINT_MAX+1 的限制。

请注意：minstd 是一个糟糕的 PRNG，mt19973 更好但不是很多：pcg-random.org/…（在该图表中，minstd==LCG32/64）。很遗憾，C++ 没有提供任何高质量、快速的 PRNG，例如 PCG 或 xoroshiro128+。

@MatteoItalia 我们没有意见分歧。这也是 Bjarne 的观点。我们真的希望在标准中使用<random>，但我们也想要一个“只给我一个我现在可以使用的体面的实现”选项。对于 PRNG 和其他东西。

几个注意事项： 1. 用eng() % 6 替换std::uniform_int_distribution<int> dist0, 5(eng); 重新引入了std::rand 代码所遭受的偏斜因素（在这种情况下，诚然是轻微偏斜，其中引擎有2**31 - 1 输出，并且您将它们分配到 6 个存储桶中）。 2. 关于“srand 采用unsigned int”的注释，它限制了可能的输出，如所写，仅使用std::random_device() 为您的引擎播种有同样的问题； you need a seed_seq to properly initialize most PRNGs.

【参考方案2】：

正确答案是：这取决于您所说的“更好”。

“新”<random> 引擎在 13 年前被引入 C++，因此它们并不是真正的新引擎。 C 库 rand() 是几十年前推出的，在当时对许多事情都非常有用。

C++ 标准库提供三类随机数生成器引擎：线性同余（rand() 是一个示例）、滞后斐波那契和 Mersenne Twister。每个类都有权衡，每个类在某些方面都是“最好的”。例如，LCG 的状态非常小，如果选择正确的参数，在现代桌面处理器上相当快。 LFG 具有更大的状态，仅使用内存获取和加法操作，因此在缺乏专门数学硬件的嵌入式系统和微控制器上速度非常快。 MTG的状态巨大，速度慢，但可以有非常大的非重复序列，具有优异的光谱特性。

如果提供的生成器都不足以满足您的特定用途，C++ 标准库还为硬件生成器或您自己的自定义引擎提供接口。没有一个生成器打算单独使用：它们的预期用途是通过一个分布对象，该对象提供具有特定概率分布函数的随机序列。

<random> 相对于rand() 的另一个优势是rand() 使用全局状态，不是可重入或线程安全的，并且允许每个进程有一个实例。如果您需要细粒度的控制或可预测性（即能够在给定 RNG 种子状态的情况下重现错误），那么rand() 毫无用处。 <random> 生成器是本地实例化的，并且具有可序列化（和可恢复）的状态。

【参考方案3】：

首先，令人惊讶的是，答案会根据您使用随机数的目的而变化。如果要驱动，比如说，一个随机的背景颜色转换器，使用 rand() 是非常好的。如果您使用随机数来创建随机扑克牌或加密安全密钥，那就不行了。

可预测性：序列 012345012345012345012345... 将提供样本中每个数字的均匀分布，但显然不是随机的。对于一个随机的序列，n+1 的值不能通过 n 的值（甚至通过 n、n-1、n-2、n-3 等的值）轻易预测。显然是重复序列相同数字是退化的情况，但任何线性同余生成器生成的序列都可以进行分析；如果您使用公共库中公共 LCG 的默认开箱即用设置，恶意人员可以毫不费力地“破坏序列”。过去，一些在线赌场（和一些实体赌场）因使用不良随机数生成器的机器而遭受损失。就连本该更了解的人也被追上了；来自多家制造商的 TPM 芯片已被证明比密钥的位长预测的更容易破解，因为密钥生成参数选择不当。

分布：正如视频中提到的，取模数 100（或任何不能整除序列长度的值）将保证某些结果的可能性至少比其他结果稍微高一些。在以 100 为模的 32767 个可能的起始值的宇宙中，数字 0 到 66 的出现频率将比值 67 到 99 多 328/327 (0.3%)；一个可能为攻击者提供优势的因素。

【讨论】：

“可预测性：序列 012345012345012345012345... 将通过您的“随机性”测试，因为您的样本中的每个数字都会均匀分布“实际上，不是真的；他正在测量的是运行之间的 stddev stddevs，即本质上各种运行直方图是如何展开的。使用 012345012345012345... 生成器，它将始终为零。

好点；恐怕我读过 OP 的代码有点太快了。编辑了我的答案以反映。

呵呵我知道，因为我也想进行那个测试，我注意到我得到了不同的结果?

【参考方案4】：

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

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

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

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

【讨论】：

这在 STL 的第二张幻灯片中有解释

好的，但是... STL 是谁？什么幻灯片？ （严肃的问题）

@kebs，Stephan Lavavej，请参阅问题中的 Youtube 参考资料。