It's OK to compare floating-points for equality

NB: The title of this post is an intentional clickbait. Even though I do stand for its statement, a more honest one would be something like: It's NOT OK to compare floating-points using epsilons.

You've probably heard the mantra that you must never compare floating-point numbers for exact equality, and you absolutely must use some kind of epsilon-comparison instead, like

bool approxEqual(float x, float y)
{
return abs(x - y) < 1e-4f;
}

Over the 15+ years that I've been writing code, – which often deals with geometry, graphics, physics, simulations, etc, and thus has to work with floating-points on a daily basis, – I've encountered only one or two cases where such epsilon-comparison is actually a good solution. Pretty much always there is a better solution that either involves rewriting the code in some way, or simply compares floating-points just like x == y. And pretty much always the epsilon solution was actually one of the worst possible options.

I'll show a bunch of examples where adding some kind of epsilon might be your first instinct, but actually a much better – and often much simpler – solution exists. But first, let's talk about floating-point numbers.

The whole idea of epsilon-comparison seems to come from the general perception of floating-point numbers as some kind of random black-box machine that sometimes produces inexact results because the gods of computing force it to. In reality it is a pretty deterministic (modulo compiler options, CPU flags, etc) and highly standardized system.

Floating-point numbers are necessarily inexact in that they cannot represent all possible real numbers. In fact, no finite amount of memory can, because that's how maths works – there are just way too many real numbers (or even just rational numbers, fwiw). Given that we probably only want to allocate a fixed (and not just a finite) amount of bits per such a number, we're forced to accept that only a finite set of numbers will be representable (specifically, at most (2^{bits}) of them), and for all others we'll have to deal with approximations.

I don't want this to turn into a lecture on how floating-point numbers are represented, though; I think wikipedia does a good enough job. For a more in-depth source, see this classic by David Goldberg, or this more recent one that follows the same idea. What's more important is that this "inexactness" of floating-point numbers doesn't mean some "uncertainty" or "randomness" in its behavior!

For example, any single arithmetic operation (e.g. addition, multiplication, etc) on two floating-point numbers is required to produce a floating-point number which is the closest to the actual true answer if we treat the inputs as being exact (there are some rounding rules in case of ties, i.e. when two representable numbers are equally close to the true answer). Notice that, even though the result is approximate, it is still guaranteed to be as close to the truth as possible, and is very much deterministic.

However, it does mean that our usual mathematical formulas don't always hold for floating-points. For example, even though addition (and multiplication) are guaranteed to be commutative ((a + b = b + a)), they aren't necessarily associative: it may happen that ((a+b)+c\neq a+(b+c)). It's fairly easy to find such examples; here's one that works for 32-bit floating-points:

// Outputs 0.89999998
std::cout << std::setprecision(8) << ((0.2f + 0.3f) + 0.4f) << '\n';
// Outputs 0.90000004
std::cout << std::setprecision(8) << (0.2f + (0.3f + 0.4f)) << '\n';

Notice that, even though the expressions aren't equal, they are very close (the difference is about 6e-8), and the floating-point standard has some guarantees that we can use to predict how large the difference can be.

So, what's the problem, then? We know that the result of some computation is only approximate, and we compare it with some expected result approximately, sounds about right.

Or does it?

I have 3 major problems with those epsilon-comparisons:

The second point is probably the toughest. One part of the program treats 2D points as equal if they are separated in Manhattan distance by no more than 1e-4, another part of the program treats points as equal if they are separated in L-inf distance (max of coordinates) by no more than 1e-6, the input points themselves are generated using some other algorithm, and now all the input-output invariants are messed up, and your line rendering is broken but only in this specific scenario, with this specific data, only at night and in full moon. Good luck debugging that.

A rare wrong case of line rendering isn't that much of an issue, but it can manifest in a ton of other ways, up to crashing the program, and can be really nasty when combined with a lot of other geometric code. I've encountered many such cases, and mismatched epsilon comparisons were an extremely common root cause.

The problem is that these epsilons are pretty much always simply guessed, and there is no correct way to choose one epsilon out of many. If you have several such comparisons, no combination of epsilons will ensure that everything works correctly.

Another problem with epsilons is that such comparison isn't transitive. This might sound like a technical nitpick, but in reality most algorithms assume that things like comparisons are in fact transitive, and these algorithms can simply break (produce nonsense or even crash) if you use a non-transitive comparison with them.

So, what should we do, then? We need to think about the problem! Shocking, I know. Why are we comparing floating-points in the first place? Maybe we don't trust our algorithms? Maybe we don't trust the data? Maybe we don't trust the CPU? There's no single answer, so let's look at some examples.

Say, you have a turn-based game where units move on a grid. A unit has some movement points and can do several moves per turn, but for UX sake you only allow executing the next move after the previous one finished.

Now, you're probably interpolating the unit's position somehow and not just teleporting it to the target grid cell, so you need to check when exactly the move is finished before allowing the player to select the next move target.

You could just wait for a certain amount of time (and that would be a perfectly good solution in many cases!), but different units have different animations and thus different timings, there are some accessibility settings to reduce animations, etc, so you decide that relying on animation time isn't a good idea.

Then you realize that the move finishes exactly when the position of the unit coincides with the target cell's center. You write something like

void update() {
if (selectedUnit.position != targetCell.center)
return;
// Do the frame update
}

and after the move was executed, nothing happens and the game effectively hangs, because the condition selectedUnit.position != targetCell.center is never true. With a typical linear interpolation with easing, the code

vec2f getPosition(float time) {
return lerp(start, end, easing(time));
}

will produce enough floating-point roundings that the result will never be equal to end when time == 1.f. Heck, in some interpolation schemes it's not even supposed to be ever true!

Damn, stupid floating-points! — you grumble as you add the holy grail of floating-points — an epsilon-comparison:

void update() {
if (distance(selectedUnit.position, targetCell.center) > 1e-4)
return;
// Do the frame update
}

This solves the issue, so why is this so bad? For a number of reasons:

1e-4 might break if you decide to choose a different interpolation scheme. So, how do we solve this? One option is to use some sort of acceptance radius: once the unit is within, say, 0.25 of the target cell's center, we stop the animation and allow the player to issue the next commands. Then, we find the actually good value by a lot of testing.