Title: Why Floats Are Never Equal

Author: Martin Moene

Date: 08 September 2016 16:14:05 +01:00 or Thu, 08 September 2016 16:14:05 +01:00

Summary: Silas S. Brown tries his hand at optimising floating point equality comparisons.

Body: 

Optimiser 1, Silas Brown 0

You’re probably aware that comparing two floats or doubles for equality is a bad idea due to rounding error. If two numbers which you expect to be equal were calculated in two different ways, the two different calculations will likely have different rounding errors and you don’t quite end up getting an equal result, so you have to do something like:

  fabs(b-a) < .0001

where .0001 is a tolerance chosen appropriately for your application (there are ways of working it out in the general case, but if you know your application then you probably have a good idea of what sort of numbers you’re dealing with). Do use fabs rather than coding it yourself, as many compilers will convert fabs to a single instruction.

I was doing some voluntary programming for a cancer research team [1], and it involved a thermodynamics calculation on fragments of DNA at different alignments. Once all possible alignments had been checked and we knew which one gave the lowest result, I wanted to display information about it. I had three options:

1. Calculate display information about ALL the items, store it, and then display only the one we found to be minimum;
2. Store the loop counter of the minimum item so far so we can go back to it;
3. Just store the minimum value, then go back through the whole loop and find it again.

Option 1 would run slower and moreover would be a memory management hassle to code. Normally I would have picked option 2, but in this case each loop iteration was making incremental changes to quite a few variables along the way, and on the other hand there weren’t very many loop iterations, so I decided on option 3.

Now surely, I thought, this could be the one occasion where it is all right to save a couple of CPU cycles by comparing two floating-point values for equality. After all, the second version of the loop is exactly retracing the steps of the first, with exactly the same formulae and only the addition of some extra display code in between. The rounding error can’t possibly be any different this time, can it?

Wrong.

The extra display code was nothing to do with the floating-point calculation, but it did affect the optimiser’s decisions of which values to keep in registers and which ones to write back to memory. And it turned out the x86 CPU had 80-bit floating-point registers, which were being rounded down to 64-bit (or 32-bit) when writing to memory. So if the optimiser found enough room to keep an intermediate result in a register, it would be kept at 80-bit precision, but if it decided to spill that register to memory to make room for something else, precision would be lost. And optimisers these days have ways of using floating-point registers as holding bays for non-floating point data, so just because the extra code wasn’t doing any floating point itself, didn’t mean it wouldn’t be compiled to something that ‘wanted’ those registers. The presence of this display code was causing the floating-point registers to be saved out to memory, and precision to be lost, making the equality comparison still wrong even though I was looking at the same C formulae in both cases.

I didn’t spot this bug at first because, when I was developing, I happened to be using a Mac, which is based on BSD. BSD defaults to putting the x86 CPU into a mode where it always rounds its intermediate results, so there is no difference between a result held in a register and one written to memory. So I was getting away with it. But then we tried to run the same code on GNU/Linux, which doesn’t tell the CPU to round, because it prefers accuracy to predictability. And it got stuck in an infinite loop looking for a minimum that didn’t exist.

Of course there are ways to dodge the issue on x86 CPUs. You could find out how to change the CPU’s rounding mode yourself (it involves assembly). Or you could make sure to always use 80-bit variables (usually called ‘long double’ by x86 compilers), although if you don’t need that much precision then the extra memory operations would probably slow you down far more than a ‘nearly equal’ comparison would (and the same goes for GCC’s -ffloat-store option, which prevents any intermediate floating-point values from being held in registers at all). But what if someone wants to compile your code on some other kind of CPU? Can you be sure that there won’t be a similar problem?

Floats are never equal. Well they are sometimes, but you don’t know when, even if you think you do. Not unless you’re programming in assembly and you know exactly which intermediate results are held in which registers at which times, and the chances are your application doesn’t make it worth your time to go there. Floats are never equal.

Finally, if you really want to save 2 CPU instructions and are looking for a known minimum, then instead of doing fabs(x-minimum)<.0001, you can first add .0001 to minimum outside the loop, then simply test for x < minimum. Do check this tolerance is appropriate to your application and can be represented by your chosen floating-point precision.

GCC has a warning when you try to compare two floats for equality, but you have to switch it on with the -Wfloat-equal option, which is sadly not included in -Wall or -Wextra.

Reference

[1] http://people.ds.cam.ac.uk/ssb22/pooler

Notes: 

More fields may be available via dynamicdata ..