这个话题在 StackOverflow 上出现过很多次，但我相信这是一个新的尝试。是的，我已经阅读了Bruce Dawson 的文章和What Every Computer Scientist Should Know About Floating-Point Arithmetic以及这个不错的答案。
据我了解，在典型系统上比较浮点数是否相等时有四个基本问题：
- 浮点计算不准确
- 是否
a-b“小”取决于规模a和b - 是否
a-b“小”取决于aand的类型b（例如float、double、long double） - 浮点通常具有 +-infinity、NaN 和非规范化表示，其中任何一个都可能干扰朴素的公式
这个答案——又名。“谷歌方法”——似乎很流行。它确实处理了所有棘手的情况。它确实非常精确地扩展了比较，检查两个值是否在彼此固定数量的ULP内。因此，例如，一个非常大的数字将“几乎等于”与无穷大进行比较。
然而：
- 在我看来，这非常混乱。
- 它不是特别便携，严重依赖内部表示，使用联合从浮点数中读取位等。
- 它只处理单精度和双精度 IEEE 754（特别是没有 x86 long double）
我想要类似的东西，但使用标准 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最后位置的 Units，但此函数返回 false（差异越小越好） - 我错过的任何案例
- 此代码依赖于未定义行为或中断的任何方式取决于实现定义的行为。（如果可能，请引用相关规范。）
- 修复您发现的任何问题
- 任何在不破坏代码的情况下简化代码的方法
我打算在这个问题上放一个不平凡的赏金。