Сравнение плавающих точек

Эта тема много раз возникала в StackOverflow, но я считаю, что это новый подход. Да, я прочитал статьи Брюса Доусона и Что каждый компьютерный ученый должен Знайте о арифметике с плавающей точкой и этот хороший ответ.

Как я понимаю, в типичной системе есть четыре основные проблемы при сравнении чисел с плавающей запятой для равенства:

  • Расчеты с плавающей точкой не являются точными.
  • Независимо от того, a-b "маленький" зависит от шкалы a и b
  • Независимо от того, a-b "small" зависит от типа a и b (например, float, double, long double)
  • Плавающая точка обычно имеет + -инфекцию, NaN и денормализованные представления, любые из которых могут мешать разработке na & iuml; ve

Этот ответ - ака. "подход Google" - кажется, популярен. Он справляется со всеми сложными делами. И это очень точно масштабирует сравнение, проверяя, находятся ли два значения в пределах фиксированного числа ULPs друг друга. Так, например, очень большое число сравнивает "почти равное" с бесконечностью.

Однако:

  • Это очень грязно, на мой взгляд.
  • Он не особенно переносим, ​​сильно полагаясь на внутренние представления, используя объединение для чтения бит из поплавка и т.д.
  • Он обрабатывает только IEEE 754 с одноточечной и двойной точностью (в частности, не длинный двойной x86)

Мне нужно что-то подобное, но с использованием стандартного С++ и обработки длинных удвоений. Под "стандартным", я имею в виду С++ 03, если это возможно, и С++ 11, если это необходимо.

Вот моя попытка.

#include <cmath>
#include <limits>
#include <algorithm>

namespace {
// Local version of frexp() that handles infinities specially.
template<typename T>
T my_frexp(const T num, int *exp)
{
    typedef std::numeric_limits<T> limits;

    // Treat +-infinity as +-(2^max_exponent).
    if (std::abs(num) > limits::max())
    {
        *exp = limits::max_exponent + 1;
        return std::copysign(0.5, num);
    }
    else return std::frexp(num, exp);
}
}

template<typename T>
bool almostEqual(const T a, const T b, const unsigned ulps=4)
{
    // Handle NaN.
    if (std::isnan(a) || std::isnan(b))
        return false;

    typedef std::numeric_limits<T> limits;

    // Handle very small and exactly equal values.
    if (std::abs(a-b) <= ulps * limits::denorm_min())
        return true;

    // frexp() does the wrong thing for zero.  But if we get this far
    // and either number is zero, then the other is too big, so just
    // handle that now.
    if (a == 0 || b == 0)
        return false;

    // Break the numbers into significand and exponent, sorting them by
    // exponent.
    int min_exp, max_exp;
    T min_frac = my_frexp(a, &min_exp);
    T max_frac = my_frexp(b, &max_exp);
    if (min_exp > max_exp)
    {
        std::swap(min_frac, max_frac);
        std::swap(min_exp, max_exp);
    }

    // Convert the smaller to the scale of the larger by adjusting its
    // significand.
    const T scaled_min_frac = std::ldexp(min_frac, min_exp-max_exp);

    // Since the significands are now in the same scale, and the larger
    // is in the range [0.5, 1), 1 ulp is just epsilon/2.
    return std::abs(max_frac-scaled_min_frac) <= ulps * limits::epsilon() / 2;
}

Я утверждаю, что этот код (а) обрабатывает все соответствующие случаи, (б) выполняет то же самое, что и реализация Google для одно- и двухточечной обработки IEEE-754, и (в) является совершенно стандартным С++.

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

  • Конкретные входы, отличающиеся более чем на ulps Единицы в последнем месте, но для которых эта функция возвращает true (чем больше разница, тем лучше)
  • Конкретные входы, отличающиеся на ulps Единицы в последнем месте, но для которых эта функция возвращает false (чем меньше разница, тем лучше)
  • Любые случаи, когда я пропустил
  • Любой способ, которым этот код опирается на поведение undefined или ломается в зависимости от поведения, определенного реализацией. (Если возможно, процитируйте соответствующую спецификацию.)
  • Исправления для любых проблем, которые вы идентифицируете.
  • Любой способ упростить код, не нарушая его.

Я намерен разместить нетривиальную награду по этому вопросу.

Ответы

Ответ 1

"Почти равно" не является хорошей функцией

4 не является подходящим значением: Ответ, который вы указываете на состояния "Следовательно, 4 должно быть достаточно для обычного использования", но не содержит оснований для этого требования. Фактически, существуют обычные ситуации, когда числа, рассчитанные в плавающей запятой различными способами, могут различаться многими ULP, даже если они будут равны, если рассчитывать по точной математике. Следовательно, для допуска не должно быть значения по умолчанию; каждый пользователь должен будет предоставить свои собственные, мы надеемся, основываясь на тщательном анализе их кода.

В качестве примера того, почему значение по умолчанию для 4 ULP является плохим, рассмотрим 1./49*49-1. Математически точный результат равен 0, но вычисленный результат (64-разрядный двоичный код IEEE 754) равен -0x1p-53, ошибка превышает 1e307 ULP точного результата и почти 1e16 ULP вычисленного результата.

Иногда значение не подходит: В некоторых случаях допуск не может быть относительным по отношению к сравниваемым значениям, ни к математически точной относительной толерантности, ни к квантованному допускам ULP. Например, почти каждое выходное значение в БПФ зависит почти от каждого входного значения, а ошибка в любом одном элементе связана с величиной других элементов. Подпрограмма "почти равных" должна быть снабжена дополнительным контекстом информацией о потенциальной ошибке.

"Почти равный" имеет плохие математические свойства:. Это показывает один из недостатков "почти равно": масштабирование изменяет результаты. Код ниже печатает 1 и 0.

double x0 = 1.1;
double x1 = 1.1 + 3*0x1p-52;
std::cout << almostEqual(x0, x1) << "\n";
x0 *= .8;
x1 *= .8;
std::cout << almostEqual(x0, x1) << "\n";

Другая неудача заключается в том, что она не транзитивна; almostEqual(a, b) и almostEqual(b, c) не означает almostEqual(a, c).

Ошибка в экстремальных случаях

almostEqual(1.f, 1.f/11, 0x745d17) неверно возвращает 1.

1.f/11 - 0x1.745d18p-4. Вычитая это из 1 (0x10p-4), получаем 0xe.8ba2e8p-4. Так как ULP 1 равно 0x1p-23, то есть 0xe.8ba2e8p19 ULP = 0xe8ba2e.8/2 ULP (сдвинутые 20 бит и деленные на 2, сетка 19 бит) = 0x745d17.4 ULP. Это превышает заданный допуск 0x745d17, поэтому правильный ответ будет равен 0.

Эта ошибка вызвана округлением в max_frac-scaled_min_frac.

Легко избавиться от этой проблемы - указать, что ulps должно быть меньше .5/limits::epsilon. Тогда округление происходит в max_frac-scaled_min_frac, только если разность (даже округленная) превышает ulps; если разность меньше, то вычитание является точным, по лемме Стербенца.

Было высказано предположение об использовании long double, чтобы исправить это. Однако long double не исправит это. Рассмотрим сравнение 1 и -0x1p-149f с установленными значениями ulps равными 1/пределам: epsilon. Если ваше значение не имеет 149 бит, результат вычитания округляется до 1, что меньше или равно 1/пределу:: epsilon ULP. Тем не менее математическая разница явно превышает 1.

Незначительное примечание

Выражение factor * limits::epsilon / 2 преобразует фактор в тип с плавающей точкой, что приводит к ошибкам округления для больших значений коэффициента, которые не являются точно представляемыми. Вероятно, подпрограмма не предназначена для использования с такими большими значениями (миллионы ULP в float), поэтому это должно быть указано как ограничение на рутину, а не на ошибку.

Ответ 2

Упрощение:. Вы можете избежать my_frexp, сначала отбросив не конечные случаи:

if( ! std::isfinite(a) || ! std::isfinite(b) )
    return a == b;

Кажется, что isfinite находится в С++ 11 по крайней мере

EDIT Однако, если намерение состоит в том, чтобы limits::infinity() в пределах 1 ulp limits::max()
то выше упрощения не выполняется, но не должно my_frexp() возвращать limits::max_exponent+1 в *exp, а не max_exponent + 2?

Ответ 3

БУДУЩЕЕ ПРОФИЛИРОВАНИЕ. Если вы когда-нибудь захотите расширить такое сравнение до десятичного поплавка http://en.wikipedia.org/wiki/Decimal64_floating-point_format в будущем, и предполагая, что ldexp() и frexp() будут обрабатывать такой тип с правильным основанием, тогда выражение striclty, 0.5 in return std::copysign(0.5, num); должно быть заменено на T(1)/limits::radix() - или std::ldexp(T(1),-1) или что-то в этом роде (я не смог найти удобная константа в std:: numeric_limits)

EDIT. Как заметил Немо, предположения, что ldexp и frexp будут использовать правильный FLOAT_RADIX, являются ложными, они придерживаются 2...

Итак, портативная версия Future Proof также должна использовать:

  • std::scalbn(x,n) вместо std::ldexp(x,n)

  • exp=std::ilogb(std::abs(x)),y=std::scalbn(x,-exp) вместо y=frexp(x,&exp)

  • теперь, когда выше y in является [1, FLOAT_RADIX) вместо [T (1)/Float_Radix, 1), верните copysign(T(1),num) вместо 0,5 для бесконечного случая my_frexp и проверьте вместо ulps*limits::epsilon() of ulps * epsilon()/2

Для этого также требуется стандарт >= С++ 11