Пересмотрено сравнение с плавающей точкой

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

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

1. Вычисления с плавающей точкой не точны
2. Будь то a-b «маленький» зависит от масштаба a а также b
3. Будь то a-b «маленький» зависит от типа a а также b (например, float, double, long double)
4. Плавающая точка обычно имеет + -infinity, NaN и денормализованные представления, любое из которых может мешать наивной формулировке

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

Тем не мение:

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

Я хочу что-то подобное, но с использованием стандартного C ++ и обработки длинных пар. Под «стандартом» я подразумеваю C ++ 03, если это возможно, и C ++ 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;
}

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

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

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

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