CODEWITHSUNDEEP
numericalnumerical-analysis
Waizung Taam
Waizung Taam

Reputation: 139

What's the best midpoint formula for floating point numbers?

The first formula

m = (a + b) / 2

is simple, but has a great risk of overflow. Besides, Numerical Analysis, 9th Edition by Burden and Faires points out that

when b - a is near the maximum precision of the machine, it is possible for (a + b) / 2 to return a midpoint the is not even in the interval [a, b].

though no further explanation is provided.

The second one

m = a + (b - a) / 2

is also correct, with a smaller chance of overflow. But for floating numbers, nearly equal values of a and b may lead to a loss of significance.

So, which formula is better in practice? Also, an explanation for the quote statement will be appreciated.

Upvotes: 8

Views: 980

Answers (3)

Hans Brende
Hans Brende

Reputation: 8567

Here's a JavaScript TL;DR of the very informative paper @vinc17 linked to:

function midpoint(a, b) {
    const mid = 0.5*(a + b);
    if (Number.isFinite(mid) || a !== a || b !== b) {
        // Handle the common case first, along with NaNs
        return mid;
    }
    if (Number.isFinite(a) !== Number.isFinite(b)) {
        // handle the case where exactly one operand is infinite
        return Math.sign(mid)*Number.MAX_VALUE;
    }
    // handle 2 infinities or 2 overflowing reals
    // NaN is falsy in javascript so the "|| 0" covers
    // the case where we have oppositely signed infinities.
    return (0.5*a + 0.5*b) || 0;
}

The result of this function is guaranteed to be correctly rounded, and to lie in the interval [min(a, b), max(a, b)] (a slight deviation from the paper in that a is not required to be less than b).

Python version:

def midpoint(a, b):
    mid = 0.5*(a + b)
    if math.isfinite(mid) or a != a or b != b:
        # common case and NaNs
        return mid
    if mid != mid:
        # opposite infinities
        return 0.0
    if math.isinf(a) is not math.isinf(b):
        # closest floating point number to infinity
        # (same as sys.float_info.max with the sign of mid)
        return math.nextafter(mid, 0)
    return 0.5*a + 0.5*b  # overflow

Note: if you don't care about the handling of infinite inputs, these functions can be simplified down to:

// JavaScript
function midpoint(a, b) {
    const mid = 0.5*(a + b);
    return Number.isFinite(mid) ? mid : 0.5*a + 0.5*b;
}

# Python
def midpoint(a, b):
    mid = 0.5*(a + b)
    return mid if math.isfinite(mid) else 0.5*a + 0.5*b

These simplified versions will yield NaN for the midpoint of -inf and +inf instead of 0, and the midpoint of a finite number with an infinite one will be infinite rather than finite (more mathematically "correct" but not as useful for binary search).

Upvotes: 0

vinc17
vinc17

Reputation: 3476

There is no best midpoint formula. This depends on the context and on the properties you would like to preserve. For binary floating-point arithmetic (base 2), I recommend Frédéric Goualard's article How do you compute the midpoint of an interval? A part of the abstract: "We review several implementations from prominent C/C++ interval arithmetic packages and analyze their potential failure to deliver correct results. We then highlight two implementations that avoid common pitfalls."

Upvotes: 4

Davis Herring
Davis Herring

Reputation: 39818

The simple (a+b)/2 is perhaps not as overflow-prone as you think—for IEEE 754 double precision, at least one of the operands has to be at least 8.988e307 (half the maximum finite value of 1.788e308) for a+b to overflow. Moreover, if it does not overflow it is correctly rounded (again, for 754) because at most one operation rounds (the division (potentially) rounds only for numbers smaller than 4.450e-308 (down to the absolute minimum of 5e-324), and no addition whose result is that close to 0 ever rounds). Since it is correctly rounded, it of course cannot be outside [a,b] since at least one of those would be closer to the true value.

If you might overflow, at least one of your values is very large, so you can just use a/2+b/2, which is then also correctly rounded (since each division is exact or else irrelevant). This is of course one more floating-point operation.

There is a caveat that the rounding mode can produce unexpected overflow or underflow with these formulae, but that’s not a common concern.

As for a+(b-a)/2, it’s just as bad for overflow if a and b might have different signs. It doesn’t, however, have “loss of significance” concerns: while the relative error in a small difference of large approximate values can of course be very large, such an operation is always exact in terms of the exact floating-point input values and thus does not contribute any numerical problems beyond those inherent in any such calculation.

Upvotes: 4

Related Questions