这个话题在StackOverflow上已经被提及多次,但我相信这是一个新的观点。是的,我阅读了Bruce Dawson的文章和计算机科学家应该知道的浮点运算知识以及这个不错的答案。
据我理解,在典型系统上,比较浮点数是否相等时存在四个基本问题:
- 浮点计算并不精确
a-b
是否“小”取决于a
和b
的规模a-b
是否“小”取决于a
和b
的类型(例如:float、double、long double)- 浮点通常具有 +-infinity、NaN 和非规格化表示,其中任何一种都可能干扰天真的公式
这个答案--也被称为“谷歌方法”--似乎很受欢迎。它可以处理所有棘手的情况。它确实非常精确地比较两个值是否相差ULPs个单位。因此,例如,一个非常大的数与无穷大“几乎相等”。
然而:
- 在我看来,它非常混乱。
- 它不是特别可移植,严重依赖于内部表示,使用联合从浮点数中读取位等。
- 它只处理单精度和双精度IEEE 754(特别是没有x86长双精度)
我想要类似的东西,但使用标准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) 对于IEEE-754单精度和双精度与谷歌实现的功能相同,(c) 完全符合C++标准。
这些声明中至少有一个是错误的。我将接受任何证明这一点的答案,最好附带修复方案。一个好的答案应该包括以下内容之一:
- 具体输入差异超过ulps个单位,但此函数返回true(差异越大越好) - 具体输入差异不超过ulps个单位,但此函数返回false(差异越小越好) - 我所错过的任何情况 - 代码依赖未定义行为或因实现定义的行为而破坏的任何方式。(如果可能,请引用相关规范。) - 您确定的任何问题的修复程序 - 在不破坏代码的情况下简化代码的任何方法
我打算在这个问题上放置一个非微不足道的悬赏。