Never been to TextSnippets before?

Snippets is a public source code repository. Easily build up your personal collection of code snippets, categorize them with tags / keywords, and share them with the world (or not, you can keep them private!)

About this user

« Newer Snippets
Older Snippets »
1 total  XML / RSS feed 

Comparing floating point numbers

AlmostEqual2sComplement is an effective way of handling floating point comparisons. 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.

// 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;
}


// Support functions and conditional compilation directives for the
// master AlmostEqual function.
#define INFINITYCHECK
#define NANCHECK
#define SIGNCHECK

inline bool IsInfinite(float A)
{
    const kInfAsInt = 0x7F800000;

    // An infinity has an exponent of 255 (shift left 23 positions) and
    // a zero mantissa. There are two infinities - positive and negative.
    if ((*(int*)&A & 0x7FFFFFFF) == kInfAsInt)
        return true;
    return false;
}

inline bool IsNan(float A)
{
    // A NAN has an exponent of 255 (shifted left 23 positions) and
    // a non-zero mantissa.
    int exp = *(int*)&A & 0x7F800000;
    int mantissa = *(int*)&A & 0x007FFFFF;
    if (exp == 0x7F800000 && mantissa != 0)
        return true;
    return false;
}

inline int Sign(float A)
{
    // The sign bit of a number is the high bit.
    return (*(int*)&A) & 0x80000000;
}

// This is the 'final' version of the AlmostEqualUlps function.
// The optional checks are included for completeness, but in many
// cases they are not necessary, or even not desirable.
bool AlmostEqualUlpsFinal(float A, float B, int maxUlps)
{
    // There are several optional checks that you can do, depending
    // on what behavior you want from your floating point comparisons.
    // These checks should not be necessary and they are included
    // mainly for completeness.

#ifdef  INFINITYCHECK
    // If A or B are infinity (positive or negative) then
    // only return true if they are exactly equal to each other -
    // that is, if they are both infinities of the same sign.
    // This check is only needed if you will be generating
    // infinities and you don't want them 'close' to numbers
    // near FLT_MAX.
    if (IsInfinite(A) || IsInfinite(B))
        return A == B;
#endif

#ifdef  NANCHECK
    // If A or B are a NAN, return false. NANs are equal to nothing,
    // not even themselves.
    // This check is only needed if you will be generating NANs
    // and you use a maxUlps greater than 4 million or you want to
    // ensure that a NAN does not equal itself.
    if (IsNan(A) || IsNan(B))
        return false;
#endif

#ifdef  SIGNCHECK
    // After adjusting floats so their representations are lexicographically
    // ordered as twos-complement integers a very small positive number
    // will compare as 'close' to a very small negative number. If this is
    // not desireable, and if you are on a platform that supports
    // subnormals (which is the only place the problem can show up) then
    // you need this check.
    // The check for A == B is because zero and negative zero have different
    // signs but are equal to each other.
    if (Sign(A) != Sign(B))
        return A == B;
#endif

    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;

    // Now we can compare aInt and bInt to find out how far apart A and B
    // are.
    int intDiff = abs(aInt - bInt);
    if (intDiff <= maxUlps)
        return true;
    return false;
}

« Newer Snippets
Older Snippets »
1 total  XML / RSS feed