1.0 - допустимый вывод std :: generate_canonical?

124

Я всегда думал, что случайные числа лежат между нулем и единицей, без него1 , т.е. это числа из полуоткрытого интервала [0,1). Справки о на 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

и лязг 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

cschwan
источник
21
@ Дэвид Лайвли 1.f == 1.fво всех случаях (какие там все случаи? Я даже не видел никаких переменных 1.f == 1.f; здесь только один случай:, 1.f == 1.fи это всегда true). Пожалуйста, не распространяйте этот миф дальше. Сравнение с плавающей запятой всегда точное.
Р. Мартиньо Фернандес
15
@DavidLively: Нет, это не так. Сравнение всегда точное. Это ваши операнды могут быть неточными, если они вычисляются, а не литералы.
Гонки легкости на орбите
2
@Galik любое положительное число ниже 1.0 является допустимым результатом. 1.0 нет. Это так просто. Округление не имеет значения: код получает случайное число и не округляет его.
R. Martinho Fernandes
7
@DavidLively он говорит, что есть только одно значение, которое сравнивается с 1.0. Это значение составляет 1,0. Значения, близкие к 1,0, не сравниваются с 1,0. Неважно, что делает функция генерации: если она вернет 1.0, она будет равна 1.0. Если он не вернет 1.0, он не будет сравнивать равным 1.0. Ваш пример с использованием abs(random - 1.f) < numeric_limits<float>::epsilonпроверяет, близок ли результат к 1.0 , что совершенно неверно в данном контексте: здесь есть числа, близкие к 1.0, которые являются действительными результатами, а именно все те, которые меньше 1.0.
R. Martinho Fernandes
4
@Galik Да, с этим будут проблемы. Но эту проблему должен решать разработчик. Пользователь никогда не должен видеть 1.0, и пользователь всегда должен видеть равное распределение всех результатов.
Р. Мартиньо Фернандес

Ответы:

121

Проблема заключается в отображении из codomain of std::mt19937( std::uint_fast32_t) в float; алгоритм, описанный в стандарте, дает неправильные результаты (несовместимые с его описанием вывода алгоритма), когда происходит потеря точности, если текущий режим округления IEEE754 отличается от округления до отрицательной бесконечности (обратите внимание, что по умолчанию используется округление -в-ближайший).

7549723-й вывод mt19937 с вашим семенем - 4294967257 ( 0xffffffd9u), который при округлении до 32-битного числа с плавающей запятой дает 0x1p+32, что равно максимальному значению mt19937, 4294967295 ( 0xffffffffu), когда оно также округляется до 32-битного числа с плавающей запятой.

Стандарт мог бы гарантировать правильное поведение, если бы он указывал, что при преобразовании вывода URNG в RealTypeof generate_canonicalдолжно выполняться округление в сторону отрицательной бесконечности; в этом случае это даст правильный результат. Как QOI, было бы хорошо, если бы libstdc ++ внесла это изменение.

После этого изменения 1.0больше не будет создаваться; вместо этого граничные значения 0x1.fffffep-Nдля 0 < N <= 8будут генерироваться чаще (примерно 2^(8 - N - 32)на N, в зависимости от фактического распределения MT19937).

Я бы рекомендовал не использовать floatwith 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.

ecatmur
источник
2
@user качество реализации - все, что делает одну совместимую реализацию лучше другой, например производительность, поведение в крайних случаях, полезность сообщений об ошибках.
ecatmur
2
@supercat: Чтобы немного отвлечься, на самом деле есть веские причины попытаться сделать синусоидальные функции как можно точнее для малых углов, например, потому что небольшие ошибки в sin (x) могут превратиться в большие ошибки в sin (x) / x (которые встречается довольно часто в реальных расчетах), когда x близко к нулю. «Дополнительная точность», кратная π, обычно является лишь побочным эффектом этого.
Илмари Каронен
1
@IlmariKaronen: Для достаточно малых углов sin (x) - это просто x. Мой крик о синусоидальной функции Java связан с углами, близкими к числу пи. Я бы сказал, что в 99% случаев, когда код просит sin(x), то на самом деле он хочет синуса (π / Math.PI), умноженного на x. Люди, поддерживающие Java, настаивают на том, что лучше иметь медленную рутинную математическую процедуру, чтобы сообщать, что синус Math.PI является разницей между π и Math.PI, чем сообщать значение, которое немного меньше, несмотря на то, что в 99% приложений он лучше бы ...
supercat 04
3
@ecatmur Предложение; обновите этот пост, чтобы упомянуть, что std::uniform_real_distribution<float>вследствие этого возникла та же проблема. (Чтобы люди, которые ищут uniform_real_distribution, увидели этот вопрос / ответ).
MM
1
@ecatmur, я не уверен, почему вы хотите округлить до отрицательной бесконечности. Поскольку generate_canonicalдолжно генерировать число в диапазоне [0,1), а мы говорим об ошибке, когда оно иногда генерирует 1.0, не будет ли округление до нуля столь же эффективным?
Маршалл Клоу,
39

По стандарту 1.0не действует.

C ++ 11 §26.5.7.2 Шаблон функции generate_canonical

Каждая функция, созданная из шаблона, описанного в этом разделе 26.5.7.2, сопоставляет результат одного или нескольких вызовов предоставленного генератора однородных случайных чисел gодному члену указанного RealType, так что, если значения g i, созданные с помощью g, равномерно распределены, результаты создания экземпляров t j , 0 ≤ t j <1 , распределяются настолько равномерно, насколько это возможно, как указано ниже.

Ю Хао
источник
25
+1 Я не вижу никаких недостатков в программе OP, поэтому я называю это ошибкой libstdc ++ и libc ++ ... что само по себе кажется маловероятным, но вот и все.
Lightness Races in Orbit
-2

Я только что столкнулся с аналогичным вопросом 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 раз. Мои эксперименты с полным перебором с libc ++ подтверждают это ожидание.

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 раз. У меня нет терпения проверить это ожидание. :)

Насколько я понимаю, такое поведение нормально. Это может оскорбить наше чувство «полуоткрытого-rangeness» , что распределение требуя вернуть номера «меньше , чем 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- из-за округления. Вы просто не можете уйти от этого; Так зачем нам в этом единственном случае делать вид, что вы можете?

Quuxplusone
источник
2
Я категорически не согласен с этой точкой зрения. Если стандарт диктует значения из полуоткрытого интервала, а реализация нарушает это правило, реализация неверна. К сожалению, как правильно указал в своем ответе ecatmur, стандарт также диктует алгоритм, в котором есть ошибка. Это также официально признано здесь: open-std.org/jtc1/sc22/wg21/docs/lwg-active.html#2524
cschwan
@cschwan: Я считаю, что реализация не нарушает правила. Стандарт диктует значения из [0,1); реализация возвращает значения из [0,1); некоторые из этих значений округляются до IEEE, 1.0fно это просто неизбежно, когда вы приводите их к числам с плавающей запятой IEEE. Если вам нужны чисто математические результаты, используйте систему символьных вычислений; если вы пытаетесь использовать числа с плавающей запятой IEEE для представления чисел, находящихся в пределах eps1, вы находитесь в состоянии греха.
Quuxplusone
Гипотетический пример, который может быть нарушен этой ошибкой: разделить что-то на canonical - 1.0f. Для каждого представимого поплавка [0, 1.0), x-1.0fне равно нуль. Имея ровно 1.0f, вы можете получить деление на ноль вместо очень маленького делителя.
Питер Кордес