重新考虑浮点比较

Question

这个主题在StackOverflow上出现了很多次,但我相信这是一个新的看法.是的,我已经阅读了布鲁斯道森的文章和每个计算机科学家应该知道的关于浮点算术的内容和这个很好的答案.

据我了解,在一个典型的系统上,比较浮点数是否相等有四个基本问题:

1. 浮点计算不准确
2. 无论a-b是"小"取决于规模a和b
3. 是否a-b为"小"取决于类型a和b(例如浮动,双,长双)
4. 浮点通常具有+ -infinity,NaN和非规范化表示,其中任何一个都可能干扰天真的公式

这个答案 - 又名."Google方法" - 似乎很受欢迎.它确实处理了所有棘手的案件.并且它确实非常精确地缩放比较,检查两个值是否在彼此的固定数量的ULP内.因此,例如,非常大的数字将"几乎相等"与无穷大相比较.

然而:

• 在我看来,它非常混乱.
• 它不是特别便携,严重依赖于内部表示,使用union来从float中读取位等.
• 它只处理单精度和双精度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单精度和双精度的Google实现相同,(c)是完全标准的C++.

这些说法中的一个或多个几乎肯定是错误的.我将接受任何证明这样的答案,最好是修复.一个好的答案应该包括以下一个或多个:

• 特定输入的差异超过ulps最后一个位置的单位,但此函数返回true(差异越大越好)
• 特定输入相差ulps最后位置的单位,但此函数返回false(差异越小越好)
• 任何我错过的案例
• 此代码依赖于未定义的行为或中断的任何方式,具体取决于实现定义的行为.(如果可能,请引用相关规范.)
• 修复您识别的任何问题
• 任何方式来简化代码而不会破坏它

我打算在这个问题上提出一个非平凡的赏金.

Answer 1

"几乎等于"不是一个好的功能

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.由于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不会纠正这个.考虑将1和-0x1p-149f与设置为1/limits :: epsilon的ulps进行比较.除非您的有效位数有149位,否则减法结果将舍入为1,小于或等于1/limits :: epsilon ULP.然而,数学差异明显超过1.

次要说明

该表达式factor * limits::epsilon / 2将因子转换为浮点类型,这会导致不可精确表示的较大因子值的舍入误差.可能,例程并不打算用于如此大的值(浮点数百万个ULP),因此应该将其指定为例程的限制而不是bug.