1.0 是 std::generate_canonical 的有效输出吗？

Question

我一直认为随机数会介于零和一之间，没有1，即它们是来自半开区间 [0,1) 的数字。 std::generate_canonical 的 documention on cppreference.com 证实了这一点。

但是，当我运行以下程序时：

#include <iostream>
#include <limits>
#include <random>

int main()
{
    std::mt19937 rng;

    std::seed_seq sequence{0, 1, 2, 3, 4, 5, 6, 7, 8, 9};
    rng.seed(sequence);
    rng.discard(12 * 629143 + 6);

    float random = std::generate_canonical<float,
                   std::numeric_limits<float>::digits>(rng);

    if (random == 1.0f)
    {
        std::cout << "Bug!\n";
    }

    return 0;
}

它给了我以下输出：

Bug!

即它为我生成了一个完美的1，这会导致我的 MC 集成出现问题。这是有效的行为还是我这边有错误？这给出了与 G++ 4.7.3 相同的输出

g++ -std=c++11 test.c && ./a.out

和clang 3.3

clang++ -stdlib=libc++ -std=c++11 test.c && ./a.out

如果这是正确的行为，我该如何避免1？

编辑 1：来自 git 的 G++ 似乎也遇到了同样的问题。我来了

commit baf369d7a57fb4d0d5897b02549c3517bb8800fd
Date:   Mon Sep 1 08:26:51 2014 +0000

使用~/temp/prefix/bin/c++ -std=c++11 -Wl,-rpath,/home/cschwan/temp/prefix/lib64 test.c && ./a.out 编译得到相同的输出，ldd 产生

linux-vdso.so.1 (0x00007fff39d0d000)
libstdc++.so.6 => /home/cschwan/temp/prefix/lib64/libstdc++.so.6 (0x00007f123d785000)
libm.so.6 => /lib64/libm.so.6 (0x000000317ea00000)
libgcc_s.so.1 => /home/cschwan/temp/prefix/lib64/libgcc_s.so.1 (0x00007f123d54e000)
libc.so.6 => /lib64/libc.so.6 (0x000000317e600000)
/lib64/ld-linux-x86-64.so.2 (0x000000317e200000)

编辑 2：我在这里报告了该行为：https://gcc.gnu.org/bugzilla/show_bug.cgi?id=63176

编辑 3：clang 团队似乎意识到了这个问题：http://llvm.org/bugs/show_bug.cgi?id=18767

Answer 1

问题在于从std::mt19937 (std::uint_fast32_t) 的共域映射到float；如果当前 IEEE754 舍入模式不是舍入到负无穷大（注意默认值是舍入- 到最近）。

带有种子的 mt19937 的第 7549723 个输出是 4294967257 (0xffffffd9u)，当四舍五入到 32 位浮点数时，得到 0x1p+32，等于 mt19937 的最大值，4294967295 (0xffffffffu)也四舍五入为 32 位浮点数。

如果标准规定在从 URNG 的输出转换为 generate_canonical 的 RealType 时，向负无穷大进行舍入，则可以确保正确的行为；在这种情况下，这将给出正确的结果。作为 QOI，libstdc++ 做这个改变会很好。

进行此更改后，1.0 将不再生成；相反，0 < N <= 8 的边界值 0x1.fffffep-N 将更频繁地生成（大约为 2^(8 - N - 32) 每 N，取决于 MT19937 的实际分布）。

我建议不要直接将float 与std::generate_canonical 一起使用；而是在double 中生成数字，然后向负无穷大舍入：

    double rd = std::generate_canonical<double,
        std::numeric_limits<float>::digits>(rng);
    float rf = rd;
    if (rf > rd) {
      rf = std::nextafter(rf, -std::numeric_limits<float>::infinity());
    }

std::uniform_real_distribution<float> 也可能出现此问题；解决方案是相同的，在double 上专门分配分布并将结果向float 中的负无穷大舍入。

Answer 2

根据标准，1.0 无效。

C++11 §26.5.7.2 函数模板 generate_canonical

从本节 26.5.7.2 中描述的模板实例化的每个函数都将提供的统一随机数生成器g 的一次或多次调用的结果映射到指定 RealType 的一个成员，这样，如果值 g_{g 产生的 >i} 是均匀分布的，实例化的结果 t_j , 0 ≤ t_j 是尽可能均匀分布，如下所示。

Answer 3

我刚刚遇到了与uniform_real_distribution 类似的问题，以下是我对标准关于该主题的简洁措辞的解释：

标准总是根据 math 定义数学函数，从不根据 IEEE 浮点（因为标准仍然假装浮点 可能不 表示 IEEE浮点）。因此，每当您在标准中看到数学用语时，它都是在谈论真正的数学，而不是 IEEE。

标准规定uniform_real_distribution<T>(0,1)(g) 和generate_canonical<T,1000>(g) 都应返回半开范围[0,1) 内的值。但这些是数学值。当您在半开范围 [0,1) 中取一个实数并将其表示为 IEEE 浮点数时，很大一部分时间它将舍入到 T(1.0)。

当T 是float（24 个尾数位）时，我们预计会看到uniform_real_distribution<float>(0,1)(g) == 1.0f 大约2^25 次中的1。 My brute-force experimentation with libc++ confirms this expectation.

template<class F>
void test(long long N, const F& get_a_float) {
    int count = 0;
    for (long long i = 0; i < N; ++i) {
        float f = get_a_float();
        if (f == 1.0f) {
            ++count;
        }
    }
    printf("Expected %d '1.0' results; got %d in practice\n", (int)(N >> 25), count);
}

int main() {
    std::mt19937 g(std::random_device{}());
    auto N = (1uLL << 29);
    test(N, [&g]() { return std::uniform_real_distribution<float>(0,1)(g); });
    test(N, [&g]() { return std::generate_canonical<float, 32>(g); });
}

示例输出：

Expected 16 '1.0' results; got 19 in practice
Expected 16 '1.0' results; got 11 in practice

当T 是double（53 个尾数位）时，我们预计会看到uniform_real_distribution<double>(0,1)(g) == 1.0 大约2^54 次中的1 个。我没有耐心去检验这个期望。 :)

我的理解是这种行为很好。它可能会冒犯我们“半开放范围”的感觉，声称返回数字“小于 1.0”的分布实际上可以返回 等于 到 1.0 的数字；但这是“1.0”的两种不同含义，明白吗？第一个是数学 1.0；第二个是IEEE单精度浮点数1.0。几十年来，我们一直被教导不要比较浮点数的精确相等性。

无论您将随机数输入何种算法，都不会在意它是否有时会恰好得到1.0。除了数学运算之外，您无法执行浮点数，并且一旦您执行一些数学运算，您的代码将不得不处理舍入。即使您可以合理地假设generate_canonical<float,1000>(g) != 1.0f，您仍然不能假设generate_canonical<float,1000>(g) + 1.0f != 2.0f——因为四舍五入。你就是无法摆脱它；那么我们为什么要在这个单一的例子中假装你可以呢？