浮点比较重新审视
这个话题在StackOverflow上多次出现，但我相信这是一个新的观点。 是的，我已阅读Bruce Dawson的文章以及每位计算机科学家应该了解的浮点算术和这个很好的答案。
据我了解，在一个典型的系统中，当比较浮点数时，存在四个基本问题：
ab是否“小”取决于a和b ab是否“小”取决于a和b的类型（例如float，double，long double） 这个答案 - 又名。 “谷歌方法” - 似乎很受欢迎。 它处理所有棘手的情况。 它确实比较精确地比较比较结果，检查两个值是否在固定数量的ULP之内。 因此，例如，一个非常大的数字比较“几乎相等”到无穷大。
然而：
我想要类似的东西，但使用标准的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 in Last Place，但该函数返回false（差异越小越好） 我打算在这个问题上提出一个不平凡的慷慨解囊。
“几乎等于”不是一个好功能
4不是一个合适的价值：你指出的答案是“因此，4应该足够用于普通用途”，但不包含该要求的基础。 事实上，在通常的情况下，通过不同的方式在浮点数中计算的数字可能会因许多ULP而有所不同，尽管如果通过精确数学计算它们是相等的。 因此，公差应该没有默认值; 应该要求每个用户提供他们自己的，希望基于他们的代码的彻底分析。
作为4 ULP默认值不好的例子，考虑1./49*49-1 。 数学上精确的结果为0，但计算结果（64位IEEE 754二进制）为-0x1p-53，超过1e307 ULP的精确结果的误差以及近似1e16的计算结果的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（即0x10p-4）中减去这个值将产生0xe.8ba2e8p-4。 由于1的ULP为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发生舍入， max_frac-scaled_min_frac是差异（即使舍入）超过ulps ; 如果差异小于这个数字，Sterbenz的引理就可以得出相减的结果。
有一个关于使用long double来纠正这个问题的建议。 然而， long double不会纠正这一点。 考虑比较1和-0x1p-149f与ulps设置为1 / limits :: epsilon。 除非您的有效位数为149位，否则相减结果会舍入为1，小于或等于1 / limits :: epsilon ULP。 然而，数学差异显然超过1。
小注
表达式factor * limits::epsilon / 2将因子转换为浮点类型，这会导致大数值因子的舍入误差不完全可表示。 有可能，这个例程并不打算用于这么大的值（数百万个浮点数ULP），所以这应该被指定为对例程的限制而不是bug。
简化：通过首先丢弃非限定案例，您可以避免my_frexp：
if( ! std::isfinite(a) || ! std::isfinite(b) )
return a == b;
至少在C ++ 11中似乎是有限的
编辑但是，如果意图是有limits::infinity()在1 ulp limits::max()
那么上面的简化不成立，但不应该my_frexp（）返回*exp limits::max_exponent+1 ，而不是max_exponent + 2？
未来的证明 ：如果你想将这种比较延伸到将来的十进制浮点数http://en.wikipedia.org/wiki/Decimal64_floating-point_format，并假设ldexp（）和frexp（）将以正确的基数处理这种类型，然后严格来说，0.5 return std::copysign(0.5, num); 应该被替换为T(1)/limits::radix() - 或std::ldexp(T(1),-1)或其他...（我在std :: numeric_limits中找不到方便的常量）
编辑 Nemo评论说，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 ulps*limits::epsilon()而不是ulps * epsilon（）/ 2
这也需要一个标准的> = C ++ 11
链接地址: http://www.djcxy.com/p/14999.html