Является ли 1.0 действительным выходом из std:: generate_canonical?
Я всегда думал, что случайные числа будут лежать между нулем и одним, без 1
, т.е. это числа из полуоткрытого интервала [0,1). documento на cppreference.com std::generate_canonical
подтверждает это.
Однако, когда я запускаю следующую программу:
#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: g++ из git, похоже, страдает от той же проблемы. Я нахожусь на
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
Ответы
Ответ 1
Проблема заключается в отображении из кодомена в std::mt19937
(std::uint_fast32_t
) на float
; алгоритм, описанный в стандарте, дает неверные результаты (несовместимые с описанием выхода алгоритма), когда происходит потеря точности, если текущий режим округления IEEE754 представляет собой нечто вроде округлой-отрицательной-бесконечности (обратите внимание, что значение по умолчанию равно -в-ближайший).
Результат 7549723rd mt19937 с вашим семенем - 4294967257 (0xffffffd9u
), который при округлении до 32-битного поплавка дает 0x1p+32
, который равен максимальному значению mt19937, 4294967295 (0xffffffffu
), когда это также округляется до 32-битного поплавка.
Стандарт мог бы обеспечить правильное поведение, если бы он указывал, что при преобразовании с выхода URNG в RealType
of generate_canonical
округление должно выполняться к отрицательной бесконечности; это даст правильный результат в этом случае. Как QOI, было бы хорошо для libstdС++ внести это изменение.
При этом изменении 1.0
больше не будет сгенерировано; вместо этого граничные значения 0x1.fffffep-N
для 0 < N <= 8
будут генерироваться чаще (приблизительно 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
.
Ответ 2
Согласно стандарту, 1.0
недействителен.
С++ 11 §26.5.7.2 Шаблон функции generate_canonical
Каждая функция, созданная из шаблона, описанного в этом разделе 26.5.7.2, отображает результат одного или нескольких вызовов созданного равномерного генератора случайных чисел g
одному члену указанного RealType таким образом, что если значения g i, создаваемые g
, равномерно распределены, результаты создания экземпляров t j, 0 j < 1, распределяются как можно более равномерно, как указано ниже.
Ответ 3
Я столкнулся с аналогичным вопросом с 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
около 1 в 2 ^ 25 раз. Мое экспериментирование с libС++ подтвердило это ожидание.
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
около 1 в 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
- из-за округления. Вы просто не можете уйти от него; так почему бы нам притвориться в этом единственном экземпляре, что вы можете?