C++ floating points comparison why epsilon() * std::fabs(x+y)?

Question

This code sample from https://en.cppreference.com/w/cpp/types/numeric_limits/epsilon

typename std::enable_if<!std::numeric_limits<T>::is_integer, bool>::type
    almost_equal(T x, T y, int ulp)
{
    // the machine epsilon has to be scaled to the magnitude of the values used
    // and multiplied by the desired precision in ULPs (units in the last place)
    return std::fabs(x-y) <= std::numeric_limits<T>::epsilon() * std::fabs(x+y) * ulp
        // unless the result is subnormal
        || std::fabs(x-y) < std::numeric_limits<T>::min();
}

Can someone explain why epsilon is scaled to fabs(x+y) instead of std::fmax(std::fabs(x), std::fabs(y))?

If `x` and `y` are negative then `std::fmax(x, y)` is negative. Do you mean the absolute maximum: `std::fmax(std::fabs(x), std::fabs(y))`? `std::fabs(x + y)` might just be a faster approximation. — Goswin von Brederlow, Jun 11 '22 at 17:07
Related, the Knuth's floating point comparisons: https://stackoverflow.com/a/253874/4641116 — Eljay, Jun 12 '22 at 02:09
Note that there might be disastrous cancellation in `x+y` if x is approximately -y and both are normal. But in that case, `x-y` is not small anyway. You can't have disastrous cancellation in both `x-y` and `x+y`. — MSalters, Jun 13 '22 at 13:56

score 1 · Answer 1 · answered Jun 13 '22 at 14:21

We can distinguish two main cases: either x and y are almost_equal, or they aren't.

When x and y are almost equal (down to a few bits), x+y is almost equal to 2*x. It's likely that x and y have the same exponent, but it's also possible that their exponent differs by exactly one (e.g. 0.49999 and 0.500001).

When x and y are not even close, e.g. when they differ by more than a factor of 2, then almost_equal(x,y) will return true anyway - the difference between the current and proposed implementations doesn't really matter.

So, what we're really worried about is the edge case - where x and y differ in about ulp bits. So getting back to the definition of epsilon: it differs from 1.0 in exactly 1 Unit in the Last Place (1 ulp). Differing in N ULP's means differing by ldexp(epsilon, N-1).

But here we need to look carefully at the direction of epsilon. It's defined by the next bigger number: nextafter(1.0, 2.0)-1.0. But that's twice the value of 1.0-nextafter(1.0, 0.0) - the previous smaller number.

This is because the difference between adjacent floating-point numbers around x scales with approximately log(x), but it's a step function. At every power of 2, the step size changes also by a step of 2.

So, your proposed lower bound of approximately x indeed works for the best case (near powers of 2) but can be off by a factor of 2. And x+y, which is almost equal to 2*x when it matters, also works for the worst-case inputs.

So you could just multiply your bound by 2, but it's already a bit slower due to the hidden conditional.

C++ floating points comparison why epsilon() * std::fabs(x+y)?

1 Answers1