Comparing floating point numbers

Bruce Dawson

This article is obsolete. Its replacement - which will fix some errors and better explain the relevant issues - is being crafted as a multi-part series here. Please update your links.

Ultimately this article will go away, once the series of articles is complete.

I mean it. Some of the problems with this code include aliasing problems, integer overflow, and an attempt to extend the ULPs based technique further than really makes sense. The series of articles listed above covers the whole topic, but the key article that demonstrates good techniques for floating-point comparisons can be found here. This article also includes a cool demonstration, using sin(double(pi)), of why the ULPs technique and other relative error techniques breaks down around zero.

In short, stop reading. Click this link.

Okay, you've been warned. The remainder of this article exists purely for historical reasons.

Comparing for equality

Floating point math is not exact. Simple values like 0.2 cannot be precisely represented using binary floating point numbers, and the limited precision of floating point numbers means that slight changes in the order of operations can change the result. Different compilers and CPU architectures store temporary results at different precisions, so results will differ depending on the details of your environment. If you do a calculation and then compare the results against some expected value it is highly unlikely that you will get exactly the result you intended.

In other words, if you do a calculation and then do this comparison:

if (result == expectedResult)

then it is unlikely that the comparison will be true. If the comparison is true then it is probably unstable – tiny changes in the input values, compiler, or CPU may change the result and make the comparison be false.

Comparing with epsilon – absolute error

Since floating point calculations involve a bit of uncertainty we can try to allow for this by seeing if two numbers are ‘close’ to each other. If you decide – based on error analysis, testing, or a wild guess – that the result should always be within 0.00001 of the expected result then you can change your comparison to this:

if (fabs(result - expectedResult) < 0.00001)

The maximum error value is typically called epsilon.

Absolute error calculations have their place, but they aren’t what is most often used. When talking about experimental error it is more common to specify the error as a percentage. Absolute error is used less often because if you know, say, that the error is 1.0 that tells you very little. If the result is one million then an error of 1.0 is great. If the result is 0.1 then an error of 1.0 is terrible.

With the fixed precision of floating point numbers in computers there are additional considerations with absolute error. If the absolute error is too small for the numbers being compared then the epsilon comparison may have no effect, because the finite precision of the floats may not be able to represent such small differences.

Let's say you do a calculation that has an expected answer of about 10,000. Because floating point math is imperfect you may not get an answer of exactly 10,000 - you may be off by one or two in the least significant bits of your result. If you're using 4-byte floats and you're off by one in the least significant bit of your result then instead of 10,000 you'll get +10000.000977. So we have:

float expectedResult = 10000;

float result = +10000.000977; // The closest 4-byte float to 10,000 without being 10,000

float diff = fabs(result - expectedResult);

diff is equal to 0.000977, which is 97.7 times larger than our epsilon. So, our comparison tells us that result and expectedResult are not nearly equal, even though they are adjacent floats! Using an epsilon value 0.00001 for float calculations in this range is meaningless – it’s the same as doing a direct comparison, just more expensive.

Absolute error comparisons have value. If the range of the expectedResult is known then checking for absolute error is simple and effective. Just make sure that your absolute error value is larger than the minimum representable difference for the range and type of float you’re dealing with.

Comparing with epsilon – relative error

An error of 0.00001 is appropriate for numbers around one, too big for numbers around 0.00001, and too small for numbers around 10,000. A more generic way of comparing two numbers – that works regardless of their range, is to check the relative error. Relative error is measured by comparing the error to the expected result. One way of calculating it would be like this:

relativeError = fabs((result - expectedResult) / expectedResult);

If result is 99.5, and expectedResult is 100, then the relative error is 0.005.

Sometimes we don’t have an ‘expected’ result, we just have two numbers that we want to compare to see if they are almost equal. We might write a function like this:

// Non-optimal AlmostEqual function - not recommended.

bool AlmostEqualRelative(float A, float B, float maxRelativeError)

{

if (A == B)

return true;

float relativeError = fabs((A - B) / B);

if (relativeError <= maxRelativeError)

return true;

return false;

}

The maxRelativeError parameter specifies what relative error we are willing to tolerate. If we want 99.999% accuracy then we should pass a maxRelativeError of 0.00001.

The initial comparison for A == B may seem odd – if A == B then won’t relativeError be zero? There is one case where this will not be true. If A and B are both equal to zero then the relativeError calculation will calculate 0.0 / 0.0. Zero divided by zero is undefined, and gives a NAN result. A NAN will never return true on a <= comparison, so this function will return false if A and B are both zero (on some platforms where NAN comparisons are not handled properly this function might return true for zero, but it will then return true for all NAN inputs as well, which makes this poor behavior to count on).

The trouble with this function is that AlmostEqualRelative(x1, x2, epsilon) may not give the result as AlmostEqualRelative(x2, x1, epsilon), because the second parameter is always used as the divisor. An improved version of AlmostEqualRelative would always divide by the larger number. This function might look like this;

// Slightly better AlmostEqual function – still not recommended

bool AlmostEqualRelative2(float A, float B, float maxRelativeError)

{

if (A == B)

return true;

float relativeError;

if (fabs(B) > fabs(A))

relativeError = fabs((A - B) / B);

else

relativeError = fabs((A - B) / A);

if (relativeError <= maxRelativeError)

return true;

return false;

}

Even now our function isn’t perfect. In general this function will behave poorly for numbers around zero. The positive number closest to zero and the negative number closest to zero are extremely close to each other, yet this function will correctly calculate that they have a huge relative error of 2.0. If you want to count numbers near zero but of opposite sign as being equal then you need to add a maxAbsoluteError check also. The function would then return true if either the absoluteError or the relativeError were smaller than the maximums passed in. A typical value for this backup maxAbsoluteError would be very small – FLT_MAX or less, depending on whether the platform supports subnormals.

// Slightly better AlmostEqual function – still not recommended

bool AlmostEqualRelativeOrAbsolute(float A, float B,

float maxRelativeError, float maxAbsoluteError)

{

if (fabs(A - B) < maxAbsoluteError)

return true;

float relativeError;

if (fabs(B) > fabs(A))

relativeError = fabs((A - B) / B);

else

relativeError = fabs((A - B) / A);

if (relativeError <= maxRelativeError)

return true;

return false;

}

Comparing using integers

There is an alternate technique for checking whether two floating point numbers are close to each other. Recall that the problem with absolute error checks is that they don’t take into consideration whether there are any values in the range being checked. That is, with an allowable absolute error of 0.00001 and an expectedResult of 10,000 we are saying that we will accept any number in the range 9,999.99999 to 10,000.00001, without realizing that when using 4-byte floats there is only one representable float in that range – 10,000. Wouldn’t it be handy if we could easily specify our error range in terms of how many floats we want in that range? That is, wouldn’t it be convenient if we could say “I think the answer is 10,000 but since floating point math is imperfect I’ll accept the 5 floats above and the 5 floats below that value.”

It turns out there is an easy way to do this.

The IEEE float and double formats were designed so that the numbers are “lexicographically ordered”, which – in the words of IEEE architect William Kahan means “if two floating-point numbers in the same format are ordered ( say x < y ), then they are ordered the same way when their bits are reinterpreted as Sign-Magnitude integers.”

This means that if we take two floats in memory, interpret their bit pattern as integers, and compare them, we can tell which is larger, without doing a floating point comparison. In the C/C++ language this comparison looks like this:

if (*(int*)&f1 < *(int*)&f2)

This charming syntax means take the address of f1, treat it as an integer pointer, and dereference it. All those pointer operations look expensive, but they basically all cancel out and just mean ‘treat f1 as an integer’. Since we apply the same syntax to f2 the whole line means ‘compare f1 and f2, using their in-memory representations interpreted as integers instead of floats’.

Kahan says that we can compare them if we interpret them as sign-magnitude integers. That’s unfortunate because most processors these days use twos-complement integers. Effectively this means that the comparison only works if one or more of the floats is positive. If both floats are negative then the sense of the comparison is reversed – the result will be the opposite of the equivalent float comparison. Later we will see that there is a handy technique for dealing with this inconvenience.

Because the floats are lexicographically ordered that means that if we increment the representation of a float as an integer then we move to the next float. In other words, this line of code:

(*(int*)&f1) += 1;

will increment the underlying representation of a float and, subject to certain restrictions, will give us the next float. For a positive number this means the next larger float, for a negative number this means the next smaller float. In both cases it gives us the next float farther away from zero.

We can apply this logic in reverse also. If we subtract the integer representations of two floats then that will tell us how close they are. If the difference is zero, they are identical. If the difference is one, they are adjacent floats. In general, if the difference is n then there are n-1 floats between them.

The chart below shows some floating point numbers and the integer stored in memory that represents them. It can be seen in this chart that the five numbers near 2.0 are represented by adjacent integers. This demonstrates the meaning of subtracting integer representations, and also shows that there are no floats between 1.99999988 and 2.0.

With this knowledge of the floating point format we can write this revised AlmostEqual implementation:

// Initial AlmostEqualULPs version - fast and simple, but

// some limitations.

bool AlmostEqualUlps(float A, float B, int maxUlps)

{

assert(sizeof(float) == sizeof(int));

if (A == B)

return true;

int intDiff = abs(*(int*)&A - *(int*)&B);

if (intDiff <= maxUlps)

return true;

return false;

}

It’s certainly a lot simpler, especially when you look at all the divides and calls to fabs() that it’s not doing!

The last parameter to this function is different from the previous AlmostEqual. Instead of passing in maxRelativeError as a ratio we pass in the maximum error in terms of Units in the Last Place. This specifies how big an error we are willing to accept in terms of the value of the least significant digit of the floating point number’s representation. maxUlps can also be interpreted in terms of how many representable floats we are willing to accept between A and B. This function will allow maxUlps-1 floats between A and B.

If two numbers are identical except for a one-bit difference in the last digit of their mantissa then this function will calculate intDiff as one.

If one number is the maximum number for a particular exponent – perhaps 1.99999988 – and the other number is the smallest number for the next exponent – 2.0 – then this function will again calculate intDiff as one.

In both cases the two numbers are the closest floats there are.

There is not a completely direct translation between maxRelativeError and maxUlps. For a normal float number a maxUlps of 1 is equivalent to a maxRelativeError of between 1/8,000,000 and 1/16,000,000. The variance is because the accuracy of a float varies slightly depending on whether it is near the top or bottom of its current exponent’s range. This can be seen in the chart of numbers near 2.0 – the gap between numbers just above 2.0 is twice as big as the gap between numbers just below 2.0.

Our AlmostEqualUlps function starts by checking whether A and B are equal – just like AlmostEqualRelative did, but for a different reason that will be discussed below.

Compiler issues

In our last version of AlmostEqualUlps we use pointers and casting to tell the compiler to treat the in-memory representation of a float as an int. There are a couple of things that can go wrong with this. One risk is that int and float might not be the same size. A float should be 32 bits, but an int can be almost any size. This is certainly something to be aware of, but every modern compiler that I am aware of has 32-bit ints. If your compiler has ints of a different size, find a 32-bit integral type and use it instead.

Another complication comes from aliasing optimizations. Strictly speaking the C/C++ standard says that the compiler can assume that different types do not overlap in memory (with a few exceptions such as char pointers). For instance, it is allowed to assume that a pointer to an int and a pointer to a float do not point to overlapping memory. This opens up lots of worthwhile optimizations, but for code that violates this rule—which is quite common—it leads to undefined results. In particular, some versions of g++ default to very strict aliasing rules, and don’t like the techniques used in AlmostEqualUlps.

Luckily g++ knows that there will be a problem, and it gives this warning:

warning: dereferencing type-punned pointer will break strict-aliasing rules

There are two possible solutions if you encounter this problem. Turn off the strict aliasing option using the -fno-strict-aliasing switch, or use a union between a float and an int to implement the reinterpretation of a float as an int. The documentation for -fstrict-aliasing gives more information.

Complications

Floating point math is never simple. AlmostEqualUlps doesn’t properly deal with all the peculiar types of floating point numbers. Whether it deals with them well enough depends on how you want to use it, but an improved version will often be needed.

IEEE floating point numbers fall into a few categories:

• Zeroes
• Subnormals
• Normal numbers
• Infinities
• NANs

Zeroes

AlmostEqual is designed to deal with normal numbers, and it is there that it behaves its best. Its first problem is when dealing with zeroes. IEEE floats can have both positive and negative zeroes. If you compare them as floats they are equal, but their integer representations are quite different – positive 0.0 is an integer zero, but negative zero is 0x80000000! (in decimal this is -2147483648). The chart below shows the positive and negative floats closest to zero, together with their integer representations.

In AlmostEqualUlps I deal with this by checking for A and B being exactly equal, thus handling the case where one input is positive zero and the other is negative zero. However this still isn’t perfect. With this implementation positive zero and the smallest positive subnormal will be calculated as being one ulp apart, and therefore will generally count as being equal. However negative zero and the smallest positive subnormal will be counted as being far apart, thus destroying the idea that positive and negative zero are identical.

A more general way of handling negative numbers is to adjust them so that they are lexicographically ordered as twos-complement integers instead of as signed magnitude integers. This is done by detecting negative numbers and subtracting them from 0x80000000. This maps negative zero to an integer zero representation – making it identical to positive zero – and it makes it so that the smallest negative number is represented by negative one, and downwards from there. With this change the representations of our numbers around zero look much more rational.

Once we have made this adjustment we can no longer treat our numbers as IEEE floats – the values of the negative numbers will be dramatically altered – but we can compare them as ints more easily, in our new and convenient representation.

This technique has the additional advantage that now the distance between numbers can be measured across the zero boundary. That is, the smallest subnormal positive number and the smallest subnormal negative number will now compare as being very close – just a few ulps away. This is probably a good thing – it’s equivalent to adding an absolute error check to the relative error check. Code to implement this technique looks like this:

// Usable AlmostEqual function

bool AlmostEqual2sComplement(float A, float B, int maxUlps)

{

// Make sure maxUlps is non-negative and small enough that the

// default NAN won't compare as equal to anything.

assert(maxUlps > 0 && maxUlps < 4 * 1024 * 1024);

int aInt = *(int*)&A;

// Make aInt lexicographically ordered as a twos-complement int

if (aInt < 0)

aInt = 0x80000000 - aInt;

// Make bInt lexicographically ordered as a twos-complement int

int bInt = *(int*)&B;

if (bInt < 0)

bInt = 0x80000000 - bInt;

int intDiff = abs(aInt - bInt);

if (intDiff <= maxUlps)

return true;

return false;

}

Subnormals

The next potential issue is subnormals, also known as denormals. Subnormals are numbers that are so small that they cannot be normalized. This lack of normalization means that they have less precision – the closer they get to zero, the less precision they have. This means that when comparing two subnormals, an error of one ulp can imply a significant relative error – as great as 100%. However the interpretation of the ulps error as a measure of the number of representable floats between the numbers remains. Thus, this variation in the relativeError interpretation is probably a good thing – yet another advantage to this technique of comparing floating point numbers.

Infinities

IEEE floating point numbers have a special representation for infinities. These are used for overflows and for the result of divide by zeroes. The representation for infinities is adjacent to the representation for the largest normal number. Thus, the AlmostEqualUlps routine will say that FLT_MAX and infinity are almost the same. This is reasonable in some sense – after all, there are no representable floats between them – but horribly inaccurate in another sense – after all, no finite number is ‘close’ to infinity.

If treating infinity as being ‘close’ to FLT_MAX is undesirable then an extra check is needed.

NANs

IEEE floating point numbers have a series of representations for NANs. These representations – sharing an exponent with infinity but marked by their non-zero mantissa – are numerically adjacent to the infinities when compared as ints. Therefore it is possible for an infinite result, or a FLT_MAX result, to compare as being very close to a NAN. If your code produces NAN results then this could be very bad. However, two things protect against this problem. One is that most floating point code is designed to not produce NANs, and in fact most floating point code should treat NANs as an error by enabling floating point divide by zero and floating point illegal operation exceptions. The other reason this should not be an issue is that usually only one NAN value is generated. On x87 compatible processors this value is 0xFFC00000, which has a value separated by four million from the nearest infinity. This value is particularly well placed because another risk with NAN comparisons is that they could wrap around. A NAN with a value of 0xFFFFFFFF could compare as being very close to zero. The translation of negative numbers used by AlmostEqual2sComplement avoids this by moving the NANs where they can only wrap around to each other, but the NAN value 0xFFC00000 also avoids this problem since it keeps the NAN value four million ulps away from wrapping around.

One other complication is that comparisons involving NANs are always supposed to return false, but AlmostEqual2sComplement will say that two NANs are equal to each other if they have the same integer representation. If you rely on correct NAN comparisons you have to add extra checks.

Limitations

maxUlps cannot be arbitrarily large. If maxUlps is four million or greater then there is a risk of finding large negative floats equal to NANs. If maxUlps is sixteen million or greater then the largest positive floats will compare as equal to the largest negative floats.

As a practical matter such large maxUlps values should not be needed. A maxUlps of sixteen million means that numbers 100% larger and 50% smaller should count as equal. A maxUlps of four million means that numbers 25% larger and 12.5% smaller should count as equal. If these large maxUlps values are needed then separate checking for wrap-around above infinity to NANs or numbers of the opposite sign will be needed. To prevent accidental usage of huge maxUlps values the comparison routines assert that maxUlps is in a safe range.

AlmostEqual2sComplement is very reliant on the IEEE floating point math format, and assumes twos-complement integers of the same size. These limitations are the norm on the majority of machines, especially consumer machines, but there are machines out there that use different formats. For this reason, and because the techniques used are tricky and non-obvious, it is important to encapsulate the behavior in a function where appropriate documentation, asserts, and conditional checks can be placed.

Summary

AlmostEqual2sComplement is an effective way of handling floating point comparisons. Its behavior does not map perfectly to AlmostEqualRelative, but in many ways its behavior is arguably superior. To summarize, AlmostEqual2sComplement has these characteristics:

• Measures whether two floats are ‘close’ to each other, where close is defined by ulps, also interpreted as how many floats there are in-between the numbers
• Treats infinity as being close to FLT_MAX
• Treats NANs as being four million ulps away from everything (assuming the default NAN value for x87), except other NANs
• Accepts greater relative error as numbers gradually underflow to subnormals
• Treats tiny negative numbers as being close to tiny positive numbers

If the special characteristics of AlmostEqual2sComplement are not desirable then they can selectively be checked for. A version with the necessary checks, #ifdefed for easy control of the behavior, is available here.

AlmostEqual2sComplement works best on machines that can transfer values quickly between the floating point and integer units. This often requires going through memory and can be quite slow. This function can be implemented efficiently on machines with vector units that can do integer or floating point operations on the same registers. This allows reinterpreting the values without going through memory.

The same techniques can be applied to doubles, mapping them to 64-bit integers. Because doubles have a 53-bit mantissa a one ulp error implies a relative error of between 1/4,000,000,000,000,000 and 1/8,000,000,000,000,000.

References

IEEE Standard 754 Floating Point Numbers by Steve Hollasch

Lecture Notes on the Status of IEEE Standard 754 for Binary Floating-Point Arithmetic by William Kahan

Source code for compare functions and tests

Other papers...