我一直认为随机数应该在0和1之间,没有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
Edit1: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
发布于 2014-09-04 23:34:48
问题是从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
(根据MT19937的实际分布情况,约为每个N
的2^(8 - N - 32)
)。
我建议不要在std::generate_canonical
中直接使用float
;而是在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
中将结果舍入到负无穷大。
发布于 2017-09-03 13:16:05
我刚刚在uniform_real_distribution
上遇到了一个类似的问题,下面是我如何解释标准在这个问题上的简洁措辞:
标准总是根据数学定义数学函数,而不是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次。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.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
--因为四舍五入。你就是无法摆脱它;那么为什么我们要在这个单一的实例中假装你可以呢?
https://stackoverflow.com/questions/25668600
复制相似问题