比较两个双精度浮点数或两个浮点数最有效的方法是什么?
简单地这样做是不正确的:
bool CompareDoubles1 (double A, double B)
{
return A == B;
}
比如:
bool CompareDoubles2 (double A, double B)
{
diff = A - B;
return (diff < EPSILON) && (-diff < EPSILON);
}
似乎是浪费加工。
有人知道更聪明的浮点比较器吗?
不幸的是,即使您的“浪费”代码也是不正确的。EPSILON是可以添加到1.0并更改其值的最小值。值1.0非常重要——更大的数字在添加到EPSILON时不会改变。现在,您可以将这个值缩放到您正在比较的数字,以判断它们是否不同。比较两个双精度对象的正确表达式是:
if (fabs(a - b) <= DBL_EPSILON * fmax(fabs(a), fabs(b)))
{
// ...
}
这是最小值。一般来说,你会想要在计算中考虑噪声,并忽略一些最不重要的位,所以更现实的比较应该是这样的:
if (fabs(a - b) <= 16 * DBL_EPSILON * fmax(fabs(a), fabs(b)))
{
// ...
}
如果比较性能对您非常重要,并且您知道值的范围,那么您应该使用定点数字。
'返回fabs(a - b) < EPSILON;
这是可以的,如果:
输入的数量级变化不大
极少数相反的符号可以被视为相等
否则就会给你带来麻烦。双精度数的分辨率约为小数点后16位。如果您正在比较的两个数字在量级上大于EPSILON*1.0E16,那么您可能会说:
return a==b;
我将研究一种不同的方法,假设您需要担心第一个问题,并假设第二个问题对您的应用程序很好。解决方案应该是这样的:
#define VERYSMALL (1.0E-150)
#define EPSILON (1.0E-8)
bool AreSame(double a, double b)
{
double absDiff = fabs(a - b);
if (absDiff < VERYSMALL)
{
return true;
}
double maxAbs = max(fabs(a) - fabs(b));
return (absDiff/maxAbs) < EPSILON;
}
这在计算上是昂贵的,但有时是需要的。这就是我们公司必须做的事情,因为我们要处理一个工程库,输入可能相差几十个数量级。
无论如何,关键在于(并且适用于几乎所有的编程问题):评估你的需求是什么,然后想出一个解决方案来满足你的需求——不要认为简单的答案就能满足你的需求。如果在您的评估后,您发现fabs(a-b) < EPSILON将足够,完美-使用它!但也要注意它的缺点和其他可能的解决方案。
我最终花了相当多的时间在这个伟大的线程通过材料。我怀疑每个人都想花这么多时间,所以我将强调我所学到的总结和我实施的解决方案。
快速的总结
Is 1e-8 approximately same as 1e-16? If you are looking at noisy sensor data then probably yes but if you are doing molecular simulation then may be not! Bottom line: You always need to think of tolerance value in context of specific function call and not just make it generic app-wide hard-coded constant.
For general library functions, it's still nice to have parameter with default tolerance. A typical choice is numeric_limits::epsilon() which is same as FLT_EPSILON in float.h. This is however problematic because epsilon for comparing values like 1.0 is not same as epsilon for values like 1E9. The FLT_EPSILON is defined for 1.0.
The obvious implementation to check if number is within tolerance is fabs(a-b) <= epsilon however this doesn't work because default epsilon is defined for 1.0. We need to scale epsilon up or down in terms of a and b.
There are two solution to this problem: either you set epsilon proportional to max(a,b) or you can get next representable numbers around a and then see if b falls into that range. The former is called "relative" method and later is called ULP method.
Both methods actually fails anyway when comparing with 0. In this case, application must supply correct tolerance.
实用函数实现(c++ 11)
//implements relative method - do not use for comparing with zero
//use this most of the time, tolerance needs to be meaningful in your context
template<typename TReal>
static bool isApproximatelyEqual(TReal a, TReal b, TReal tolerance = std::numeric_limits<TReal>::epsilon())
{
TReal diff = std::fabs(a - b);
if (diff <= tolerance)
return true;
if (diff < std::fmax(std::fabs(a), std::fabs(b)) * tolerance)
return true;
return false;
}
//supply tolerance that is meaningful in your context
//for example, default tolerance may not work if you are comparing double with float
template<typename TReal>
static bool isApproximatelyZero(TReal a, TReal tolerance = std::numeric_limits<TReal>::epsilon())
{
if (std::fabs(a) <= tolerance)
return true;
return false;
}
//use this when you want to be on safe side
//for example, don't start rover unless signal is above 1
template<typename TReal>
static bool isDefinitelyLessThan(TReal a, TReal b, TReal tolerance = std::numeric_limits<TReal>::epsilon())
{
TReal diff = a - b;
if (diff < tolerance)
return true;
if (diff < std::fmax(std::fabs(a), std::fabs(b)) * tolerance)
return true;
return false;
}
template<typename TReal>
static bool isDefinitelyGreaterThan(TReal a, TReal b, TReal tolerance = std::numeric_limits<TReal>::epsilon())
{
TReal diff = a - b;
if (diff > tolerance)
return true;
if (diff > std::fmax(std::fabs(a), std::fabs(b)) * tolerance)
return true;
return false;
}
//implements ULP method
//use this when you are only concerned about floating point precision issue
//for example, if you want to see if a is 1.0 by checking if its within
//10 closest representable floating point numbers around 1.0.
template<typename TReal>
static bool isWithinPrecisionInterval(TReal a, TReal b, unsigned int interval_size = 1)
{
TReal min_a = a - (a - std::nextafter(a, std::numeric_limits<TReal>::lowest())) * interval_size;
TReal max_a = a + (std::nextafter(a, std::numeric_limits<TReal>::max()) - a) * interval_size;
return min_a <= b && max_a >= b;
}
你写的代码有bug:
return (diff < EPSILON) && (-diff > EPSILON);
正确的代码应该是:
return (diff < EPSILON) && (diff > -EPSILON);
(…是的,这是不同的)
我想知道晶圆厂是否会让你在某些情况下失去懒惰的评价。我会说这取决于编译器。你可能想两种都试试。如果它们在平均水平上是相等的,则采用晶圆厂实现。
如果你有一些关于两个浮点数中哪一个比另一个更大的信息,你可以根据比较的顺序来更好地利用惰性求值。
最后,通过内联这个函数可能会得到更好的结果。不过不太可能有太大改善……
编辑:OJ,谢谢你纠正你的代码。我相应地删除了我的评论