这个话题在StackOverflow上已经出现过很多次了，但是我相信这是一个新想法。是的，我已阅读Bruce Dawson's articles和What Every Computer Scientist Should Know About Floating-Point Arithmetic和this nice answer。
据我了解，在典型系统上，比较浮点数是否相等时存在四个基本问题:
a-b是否“小”取决于a和b的大小a-b是否为“small”取决于a和b的类型(例如float，double，long double)This answer-又名“Google方法”-似乎很流行。它确实可以处理所有棘手的情况。而且它确实非常精确地缩放比较，检查两个值是否在彼此的固定数量ULPs之内。因此，例如，非常大的数字将“几乎等于”与无穷大进行比较。
然而:
我想要类似的东西，但是使用标准C++并处理长 double 。所谓“标准”，是指如果可能的话是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单精度和 double 的Google实现相同，并且(c)是完全标准的C++。
这些主张中的一项或多项几乎可以肯定是错误的。我将接受任何证明这一点的答案，最好是通过修复来解决。一个好的答案应包括以下一项或多项:
ulps单位，但为此功能返回true(差异越大，效果越好)ulps单位，但为此功能返回false(差异越小越好)我打算对这个问题给予不小的悬赏。
最佳答案
“几乎等于”不是一个好函数
4不是适当的值:您指出的答案指出“因此，4应该足够用于普通用途”，但不包含该主张的依据。实际上，在许多情况下，即使通过精确数学计算得出的数值相等，许多情况下ULP也会以不同的方式在浮点数中计算得出的数值有所不同。因此，公差应没有默认值。应要求每个用户提供自己的代码，希望基于对代码的全面分析。
作为为什么默认4 ULP不好的一个示例，请考虑1./49*49-1。数学上精确的结果为0，但计算结果(64位IEEE 754二进制)为2–53，误差超过精确结果的10307 ULP和计算结果的1016 ULP。
有时，没有合适的值:在某些情况下，公差不能相对于要比较的值，无论是数学上精确的相对公差还是量化的ULP公差。例如，FFT中几乎每个输出值都会受到几乎每个输入值的影响，任何一个元素中的误差都与其他元素的大小有关。必须为“几乎等于”例程提供附加上下文以及有关潜在错误的信息。
“几乎等于”具有较差的数学属性:这表明“几乎等于”的缺点之一:缩放会改变结果。下面的代码显示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.745d1816•2-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中进行舍入；如果差值小于此值，则减去是准确的，Sterbenz的引理。有关于使用
long double纠正此问题的建议。但是，long double无法纠正此问题。考虑将ulps设置为1/limits::epsilon来比较1和-0x1p-149f。除非有效位数为149位，否则减法结果将舍入为1，小于或等于1/limits::epsilon ULP。然而数学上的差异显然超过了1。次音符
表达式
factor * limits::epsilon / 2将factor转换为浮点类型，这会导致无法精确表示的大factor值舍入错误。可能，例程不打算与如此大的值(float中的数百万个ULP)一起使用，因此应将其指定为例程的限制而不是错误。