Как сравнить float и double, учитывая потерю точности?

Будьте предельно осторожны при использовании любых других предложений. Всё зависит от контекста.

Я потратил много времени на поиск ошибок в системе, которая предполагала, что a==b, если |a-b|<epsilon. Основные проблемы заключались в следующем:

1. Неприметное предположение в алгоритме, что если a==b и b==c, то a==c.
2. Использование одного и того же epsilon для значений, измеряемых в дюймах, и для значений, измеряемых в милs (0.001 дюймов). То есть a==b, но при этом 1000a!=1000b. (По этой причине AlmostEqual2sComplement требует указания epsilon или максимального ULPS).
3. Применение одного и того же epsilon как для косинусов углов, так и для длин линий!
4. Использование такой функции сравнения для сортировки элементов в коллекции. (В этом случае использование встроенного оператора == для double в C++ давало корректные результаты.)

Как я уже сказал: всё зависит от контекста и ожидаемого размера a и b.

Кстати, std::numeric_limits<double>::epsilon() — это "машинный эпсилон". Это разница между 1.0 и следующим значением, представимым в формате double. Я полагаю, его можно использовать в функции сравнения, но только если ожидаемые значения меньше 1. (Это в ответ на комментарий @cdv...)

Также, если у вас фактически имеется арифметика int в double (здесь мы используем double для хранения целых значений в определённых случаях), ваши вычисления будут корректными. Например, 4.0/2.0 будет равно 1.0+1.0. Это верно, пока вы не выполняете операции, которые приводят к дробным значениям (4.0/3.0), или не выходите за пределы размера int.

Портативный способ получения эпсилон в C++ выглядит следующим образом:

#include <limits>
std::numeric_limits<double>::epsilon()

Теперь функция сравнения может быть реализована так:

#include <cmath>
#include <limits>

bool AreSame(double a, double b) {
    return std::fabs(a - b) < std::numeric_limits<double>::epsilon();
}

Такой подход позволяет корректно сравнивать два числа с плавающей запятой, учитывая погрешности, которые могут возникать в результате операций с ними.

Я провел значительное время, изучая материалы в этой отличной теме. Я сомневаюсь, что всем хочется тратить столько времени, поэтому выделю краткое резюме того, что я узнал, и решение, которое я реализовал.

Краткое резюме

1. Является ли 1e-8 приблизительно равным 1e-16? Если вы анализируете шумные данные с датчиков, то, возможно, да. Но если вы занимаетесь молекулярным моделированием, то, возможно, нет! Главное: всегда нужно учитывать значение допусков в контексте конкретного вызова функции, а не просто использовать его как общую жестко закодированную константу.

2. Для общих функций библиотек приятно иметь параметр с умолчательным допуском. Обычно выбирается numeric_limits::epsilon(), который равен FLT_EPSILON в float.h. Однако это проблематично, поскольку эпсилон для сравнения значений, таких как 1.0, не равен эпсилону для значений типа 1E9. FLT_EPSILON определяется для 1.0.

3. Очевидная реализация для проверки, находится ли число в пределах допуска, — это fabs(a-b) <= epsilon, однако это не работает, потому что умолчательный эпсилон определяется для 1.0. Нам нужно масштабировать эпсилон в зависимости от a и b.

4. Существует два решения этой проблемы: либо вы устанавливаете эпсилон пропорционально max(a,b), либо получаете ближайшие представимые числа вокруг a и проверяете, попадает ли b в этот диапазон. Первое называется "относительным" методом, а второе — методом ULP.

5. Оба метода на самом деле все равно не работают при сравнении с 0. В этом случае приложение должно предоставить корректный допуск.

Реализация утилитарных функций (C++11)

// реализует относительный метод - не используйте для сравнения с нулем
// используйте это большинство времени, допуск должен иметь смысл в вашем контексте
template<typename TReal>
static bool isApproximatelyEqual(TReal a, TReal b, TReal tolerance = std::numeric_limits<TReal>::epsilon())
{
    TReal diff = std::fabs(a - b);
    if (diff <= tolerance)
        return true;

    if (diff < std::fmax(std::fabs(a), std::fabs(b)) * tolerance)
        return true;

    return false;
}

// предоставьте допуск, который имеет смысл в вашем контексте
// например, умолчательный допуск может не работать, если вы сравниваете double с float
template<typename TReal>
static bool isApproximatelyZero(TReal a, TReal tolerance = std::numeric_limits<TReal>::epsilon())
{
    if (std::fabs(a) <= tolerance)
        return true;
    return false;
}

// используйте это, когда хотите быть на безопасной стороне
// например, не запускайте ровера, пока сигнал не будет выше 1
template<typename TReal>
static bool isDefinitelyLessThan(TReal a, TReal b, TReal tolerance = std::numeric_limits<TReal>::epsilon())
{
    TReal diff = a - b;
    if (diff < tolerance)
        return true;

    if (diff < std::fmax(std::fabs(a), std::fabs(b)) * tolerance)
        return true;

    return false;
}

template<typename TReal>
static bool isDefinitelyGreaterThan(TReal a, TReal b, TReal tolerance = std::numeric_limits<TReal>::epsilon())
{
    TReal diff = a - b;
    if (diff > tolerance)
        return true;

    if (diff > std::fmax(std::fabs(a), std::fabs(b)) * tolerance)
        return true;

    return false;
}

// реализует метод ULP
// используйте это, когда вас беспокоит только проблема точности чисел с плавающей точкой
// например, если вы хотите проверить, является ли a равным 1.0, проверяя, находится ли оно в
// пределах 10 ближайших представимых чисел с плавающей точкой вокруг 1.0.
template<typename TReal>
static bool isWithinPrecisionInterval(TReal a, TReal b, unsigned int interval_size = 1)
{
    TReal min_a = a - (a - std::nextafter(a, std::numeric_limits<TReal>::lowest())) * interval_size;
    TReal max_a = a + (std::nextafter(a, std::numeric_limits<TReal>::max()) - a) * interval_size;

    return min_a <= b && max_a >= b;
}

Надеюсь, это поможет! Если у вас есть дополнительные вопросы, не стесняйтесь спрашивать.

Код, который вы написали, содержит ошибку:

return (diff < EPSILON) && (-diff > EPSILON);

Правильный код будет выглядеть так:

return (diff < EPSILON) && (diff > -EPSILON);

(...да, это действительно другое условие)

Что касается использования fabs, стоит задуматься, не приведет ли это к потере ленивой оценки в некоторых случаях. Я бы сказал, это зависит от компилятора. Рекомендуется попробовать оба варианта. Если они в среднем эквивалентны, то вы можете выбрать реализацию с использованием fabs.

Если у вас есть информация о том, какое из двух чисел с плавающей запятой более вероятно будет больше, вы можете изменить порядок сравнения, чтобы лучше использовать ленивую оценку.

Наконец, вы могли бы получить лучшие результаты, если инлайнить эту функцию, но вряд ли это даст значительное улучшение...

Правка: О, спасибо за исправление вашего кода. Я удалил свой комментарий соответственно.

Вы правы, ваш код работает, если:

• порядок величины ваших входных данных не изменяется значительно;
• очень маленькие числа противоположных знаков можно считать равными.

Однако в противном случае это может привести к проблемам. Числа с плавающей запятой двойной точности имеют разрешение примерно 16 десятичных знаков. Если два числа, которые вы сравниваете, больше по величине, чем EPSILON * 1.0E16, то фактически вы можете просто писать:

return a == b;

Я рассмотрю другой подход, предполагая, что вам нужно учитывать первую проблему, и предполагая, что вторая не вызывает беспокойства в вашей задаче. Решение могло бы выглядеть так:

#define VERYSMALL  (1.0E-150)
#define EPSILON    (1.0E-8)
bool AreSame(double a, double b)
{
    double absDiff = fabs(a - b);
    if (absDiff < VERYSMALL)
    {
        return true;
    }

    double maxAbs  = fmax(fabs(a), fabs(b));
    return (absDiff / maxAbs) < EPSILON;
}

Такой код может быть дорогостоящим с точки зрения вычислений, но иногда это необходимо. Именно так обстоят дела в нашей компании, поскольку мы работаем с инженерной библиотекой, и входные данные могут варьироваться на несколько порядков величины.

В любом случае, суть в следующем (и это относится практически к каждой программной задаче): оцените, какие у вас требования, а затем найдите решение, соответствующее вашим нуждам — не предполагая, что простое решение будет достаточным. Если после вашей оценки выясняется, что fabs(a-b) < EPSILON подходит, отлично — используйте его! Но будьте готовы к его недостаткам и другим возможным решениям.