When evaluating geometric predicates with floating-point arithmetic, roundoff error can cause programs to make incorrect or inconsistent branching decisions. Jonathan Shewchuk refined and applied arbitrary precision floating-point techniques developed by Douglas Preist and others to the most commonly used geometric predicates. His adaptively evaluated predicates guarantee that their result is consonant with an exact evaluation and allow for the robust implementation of algorithms designed for Real-RAM machines.

Shewchuk’s technique is based on a structure called an expansion. An expansion is a set of regular floating point numbers whose exact sum is equal to the number that the expansion represents. He requires expansions to be nonoverlapping and sorted by magnitude ( largest, smallest). Two floating-point values and are nonoverlapping if the least significant nonzero bit of is more significant than the most significant nonzero bit of , or vice versa. The number zero does not overlap any number. Formally, and are nonoverlapping if there exist integers and such that and , or and .

Testing the Nonoverlapping Property

Implementing expansion logic from scratch is both illuminating and rather complicated. To aid in this process, I wanted a function that would help me unit test expansion manipulation functions to verify that they maintained the nonoverlapping requirement. My implementation is almost definitely not the fastest way to perform this check, but since I’m just using it for unit testing, speed is not critical.

One tricky part of Shewchuk’s definition is the notion of most/least significant nonzero bit. In the context of integers, this notion is easy to comprehend, but floating-point numbers have a complex structure that requires careful parsing. First, let’s review this diagram from Wikipedia:

The sign bit is not relevant for whether or not two numbers overlap. The mantissa (labeled as fraction above) is clearly relevant. We have to account for two other pieces: the implicit leading bit and the exponent. For normal numbers, the MSNZB will always be the leading 1 at bit position 23. For subnormal numbers, the leading bit is zero, so the MSNZB will be found somewhere in the mantissa (we explicitly filter out zero at the top of the function). Simply adding the exponent to the LSNZB/MSNZB is sufficient to account for it.

bool XYAreNonoverlapping(const float x, const float y)
{
    if (x == 0.0f || y == 0.0f)
        return true;

    unsigned int unX = *((unsigned int*)&x);
    unsigned int unY = *((unsigned int*)&y);

    int nXExp = ((unX >> 23) & 255) - 127;
    int nYExp = ((unY >> 23) & 255) - 127;

    // Find the LSNZB in the unmodified mantissa. For normal numbers, the MSNZB will
    // always be the leading 1 at bit position 23. For subnormal numbers, the leading 
    // bit is zero, so the MSNZB will be found somewhere in the mantissa (we explicitly
    // filter out zero at the top of the function).
    // Examples:
    // 22 21 20 19 18 17 16 15 14 13 12 11 10 09 08 07 06 05 04 03 02 01 00
    //  0  1  0  0  0  0  0  1  1  0  0  0  0  0  0  1  0  0  0  0  0  0  0
    // LSNZB = 07
    //  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
    // LSNZB = 22
    int lsnzbX = 23, msnzbX = nXExp == -127 ? -1 : 23;
    int lsnzbY = 23, msnzbY = nYExp == -127 ? -1 : 23;
    for (int i = 0; i < 23; ++i)
    {
        if (unX & (1 << i))
        {
            lsnzbX = std::min(lsnzbX, i);
            msnzbX = std::max(msnzbX, i);
        }

        if (unY & (1 << i))
        {
            lsnzbY = std::min(lsnzbY, i);
            msnzbY = std::max(msnzbY, i);
        }
    }

    lsnzbX += nXExp;
    msnzbX += nXExp;
    lsnzbY += nYExp;
    msnzbY += nYExp;

    return lsnzbX > msnzbY || lsnzbY > msnzbX;
}

Test Cases

The following test cases exercise normal and subnormal numbers, using the Catch test framework.

TEST_CASE("XYAreNonoverlapping", "[expansion]")
{
    float x = 1.5f; // x = 3 * 2^(-1)
    float y = 2.0f; // x = 2^(1)

    // x and y are initially nonoverlapping because 1.5 < 2.0.
    REQUIRE(XYAreNonoverlapping(x, y));

    // The range (2.0f, 0.5f] overlaps 1.5f because x = 3 * 2^(-1) and nothing in the
    // range is < 2^(-1). Since the nonoverlapping property is symmetrical, we note
    // that y will take on values that are of the form r*2^s, but all s values will 
    // be <= -1, and 1.5 > 2^(-1).
    y = std::nextafterf(y, 0.0f);
    while (y >= 0.5f)
    {
        REQUIRE_FALSE(XYAreNonoverlapping(x, y));
        y = std::nextafterf(y, 0.0f);
    }

    // The range (0.5f, 0.125f] does not overlap 1.5f because x = 3 * 2^(-1) and 
    // nothing in the range is > 2^(-1).
    while (y >= 0.125f)
    {
        REQUIRE(XYAreNonoverlapping(x, y));
        y = std::nextafterf(y, 0.0f);
    }

    // The range (2^(-127), 0.0f) does not overlap 1.5f because x = 3 * 2^(-1) and 
    // nothing in the range is > 2^(-1).
    y = std::nextafterf(std::numeric_limits<float>::min(), 0.0f);
    while (y > 0.0f)
    {
        REQUIRE(XYAreNonoverlapping(x, y));
        y = std::nextafterf(y, 0.0f);
    }
}