Integer Arithmetic Continued

I continue my thing about integer arithmetic. The plan is to add so much math to this blog that nobody will read it, although to be fair, the amount of math that takes is “any” and nobody was reading to begin with.

avg §

Taking the average of two numbers is important. Super important. But how do you do it?

It’s been over seven years since Peter Norvig posted “Nearly All Binary Searches and Mergesorts are Broken” and still, the top result for “binary search in C” has these maths for the indices:

first = 0;
last = n - 1;
middle = (first+last)/2;

Of course, the problem here is that if first+last exceeds 2³¹-1 or 2147483647, then it overflows into the negative, and using the now-negative middle as an index will cause “undefined behavior^*.”

Before jumping into fixes, let’s take a second to appreciate that page. I didn’t expect anything that bad to come up as the first result. Let’s take in the link to a compiled EXE (!), a screenshot of it running (?), and this brilliant programmers’ admonishment about using binary search (?!):

Binary search is faster than linear search but list should be sorted, hashing is faster than binary search and perform searches in constant time.

Alright then.

In this case, the fix is really very simple: since both first and last are known to be positive, their difference can’t overflow the range of integers and you can do this:

middle = first + (last - first) / 2;

Guaranteed to not overflow. However, what happens when you need to take an average of two numbers that may not be positive? If first is very large and positive and last is very large and negative, last - first could overflow.

One naïve solution would be to use both formulae, and switch between them based on the inputs. After all, if (a+b)/2 works when a and b have different signs, and a+(b-a)/2 works when they’re both positive, then we just write a routine that applies them at the right times, right?

Except… I kind of ignored how the expression last - first is guaranteed to be positive. But if the order of our avg function’s two input values is not known, then that value can be negative. Why is this a problem?

As you might guess, it’s rounding. See, division in C++ is defined to truncate towards zero, and in C it’s a frequent default to do the same. So, negative numbers round up and positive numbers round down. You see the problem when you have an order-dependent operation like subtraction in the expression (b-a)/2?

Let’s try an example; let’s take the average of 2 and 7.

Yeah, we ended up with an avg function where avg(a, b) != avg(b, a). Not good. Sure, we could just correct for this with another conditional to check if (b-a) < 0 and subtract one from the result if one input is odd and the other even^†. But you’d agree that’s pretty complicated, right?

// note: this is to prove a point. It's not been tested.
int32_t avg(int32_t a, int32_t b) {
    if ( (a < 0) != (b < 0) ) {
        return (a + b) / 2;
    } else {
        return a + (b - a) / 2 + ( (b - a < 0) & ( (a & 1) != (b & 1) ) );
    }
}

Let’s try another approach. Maybe we should divide each input value first, then add the halves together.

int32_t avg(int32_t a, int32_t b) {
    return (a / 2) + (b / 2);
}

The error here should be obvious now. I’ll show you with examples.

Wait. So our avg function works as if avg(a, b) = avg(a + 1, b + 1)? That’s… kind of wrong. But at least avg(a, b) = avg(b, a)!

Let’s see here. So in order to build a function that works similarly to (a+b)/2, the gold standard of midpointing, we must satisfy the following:

1. avg(a, b) = avg(b, a)
2. avg(a + 1, b + 1) = avg(a, b) + 1
3. avg(-a, -b) = -avg(a, b)

That last one actually specifies a wholly different and yet unexplained requirement: symmetry about zero. I just feel that if you’re gonna be averaging signed numbers together, the function ought to work similarly on either side of the number line and not bias results in either direction. So it only makes sense to round the results toward zero.

Anyways, so how do with address our problem with the “sum of halves” solution? Can we fix it without adding a bunch of branches like with the former “conditional half of difference” solution?

Well, let’s think about why it happens. The error manifests when both halves have a half aftering halving: that is, when both a / 2 and b / 2 are something.5 (read that as “something point five”). This results in both parts being rounded down before getting summed. But, since both parts were decreased by .5, this causes a total decrease of 1 in the sum!

Of course, x / 2 needs only be rounded if x is odd, so the only time the above error manifests is then if both inputs are odd. So it’s a simple matter of adding one to the result when both are odd, right?

int32_t avg(int32_t a, int32_t b) {
    return (a / 2) + (b / 2) + (a & b & 1);
    // can also be           + (a % 2 && b % 2)
}

Kind of. See, this does take care of most of the errors, but there’s still a subtle bug that causes it to not satisfy requirement 2. avg(a + 1, b + 1) = avg(a, b) + 1. Can you spot it?

I’ll give it to ya straight: when a and b have opposite signs, they are rounded in opposite directions. Negatives are rounded up, and positives are rounded down. Not a big deal if both inputs are odd or both are even. In fact, the error only manifests when one input is odd and the other is even, as well as another condition that I haven’t even yet been able to divine. Suffice to say that it’s not really the right path to go down.

So what can we do? The easiest thing I figured out is to simply use a division that consistently rounds in the same direction, then round the result as needed. To do that, we’ll use an arithmetic right shift by 1 in order to perform a divide-by-2.

It’s a common misconception that the C expressions x / 2 and x >> 1 are equivalent. In fact, for unsigned integers, they are. However, the major distinction is that the latter rounds towards negative infinity on most platforms^‡! This is much more convenient for us to use, since corrections happen consistently across the range of integers.

int32_t avg(int32_t a, int32_t b) {
    return (a >> 1) + (b >> 1) + (a & b & 1);
    // can also be             + (a % 2 && b % 2)
}

Already, this is a correct midpoint function. However, it’s still just a bit off my mark: it doesn’t meet our third requirement, symmetry on the number line about zero. That’s a simple fix though. Toss in a 1 when parities don’t match but the sum of halves is negative, and voilà, we have a winner! Rounds towards zero and everything!

/* SPDX-License-Identifier: MIT OR Apache-2.0 */
/**
 * Overflow-safe average of two signed integers. The naive average function is
 * erroneous when the sum of the inputs overflows integer limits; this average
 * works by summing the halves of the input values and then correcting the sum
 * for rounding.
 *
 * @param a first value, as a 32-bit signed integer
 * @param b second value, as a 32-bit signed integer
 * @return signed average of the two values, rounded towards zero if their
 * average is not an integer
 */
static inline int32_t avg(int32_t a, int32_t b) {
    // shifts divide by two, rounded towards negative infinity
    const int32_t sumHalves = (a >> 1) + (b >> 1);
    // this has error of magnitude one if both are odd
    const uint32_t bothOdd = (a & b) & 1;
    // round toward zero; add one if one input is odd and sum is negative
    const uint32_t roundToZero = (sumHalves > 0) & (a ^ b);

    // result is sum of halves corrected for rounding
    return sumHalves + bothOdd + roundToZero;
}

• Real Talk: Integer Arithmetic
• Big Data EVT

• C++,
• Math