Integer Arithmetic Continued

I con­tin­ue my thing about inte­ger arith­metic. The plan is to make this blog com­plete­ly unread­able by the time I get to Sil­i­con Val­ley, because screw pop­u­lar­i­ty. I’m on vaca­tion.

avg

Tak­ing the aver­age of two num­bers is impor­tant. Super impor­tant. But every­one does it wrong.

avg

No, you can’t just add them up and divide by 2.

It’s been over sev­en years since Peter Norvig (Peter frick­in’ Norvig!) post­ed “Near­ly All Bina­ry Search­es and Merge­sorts are Bro­ken” and still, the top result for “bina­ry search in C” has the­se maths for the indices:

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

Of course, the prob­lem here is that if first+last exceeds 231-1 or 2147483647, then it over­flows into the neg­a­tive, and using the now-neg­a­tive middle as an index will cause “unde­fined behav­ior1” through mem­o­ry cor­rup­tion.

I’ll get to how to fix that in a sec­ond, but I feel the need to WTF at that page for a sec­ond. I didn’t expect any­thing that bad to come up as the first result. What shocked me, besides the oth­er gra­tu­itous prob­lems with the code, are the link to a com­piled EXE (!), a screen­shot of it run­ning in a Win­dows Com­mand Prompt (?), and this bril­liant pro­gram­mers’ admon­ish­ment about using bina­ry search (?!):

Bina­ry search is faster than lin­ear search but list should be sort­ed, hash­ing is faster than bina­ry search and per­form search­es in con­stant time.

?!

Alright then. May­be I woke up on the wrong side of key­board today, but I just don’t feel like explain­ing how that string of words do not make sense2. Any­ways.

In this case, the fix is real­ly very sim­ple: since both first and last are known to be pos­i­tive, their dif­fer­ence can’t over­flow the range of inte­gers and you can just do this:

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

Guar­an­teed to not over­flow. How­ev­er, what hap­pens when you need to take an aver­age of two num­bers that may not be pos­i­tive? If first is very large and pos­i­tive and last is very large and neg­a­tive, last - first could over­flow.

One naïve solu­tion would be to use both for­mu­lae, and switch between them based on the inputs. After all, if (a+b)/2 works when a and b have dif­fer­ent signs, and a+(b-a)/2 works when they’re both pos­i­tive, then we just write a rou­tine that applies them at the right times, right?

Except… I kind of ignored how the expres­sion last - first is guar­an­teed to be pos­i­tive. But if the order of our avg function’s two input val­ues is not known, then that val­ue can be neg­a­tive. Why is this a prob­lem?

As you might guess, it’s round­ing. See, divi­sion in C++ is defined to trun­cate towards zero, and in C it’s a fre­quent default to do the same. So, neg­a­tive num­bers round up and pos­i­tive num­bers round down. You see the prob­lem when you have an order-depen­dent oper­a­tion like sub­trac­tion in the expres­sion (b-a)/2?

Let’s try an exam­ple; let’s take the aver­age of 2 and 7.

\begin{aligned} \text{avg}(2, 7) &= 2 + \frac{7 - 2}{2} = 2 + \lfloor 2.5\rfloor = 2 + 2 = 4 \\ \text{avg}(7, 2) &= 7 + \frac{2 - 7}{2} = 7 + \lceil -2.5\rceil = 7 + -2 = 5 \end{aligned}

Yeah, we end­ed up with an avg func­tion where avg(a, b) != avg(b, a). Not good. Sure, we could just cor­rect for this with anoth­er con­di­tion­al to check if (b-a) < 0 and sub­tract one from the result if one input is odd and the oth­er even3. But you’d agree that’s pret­ty com­pli­cat­ed, 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) ) );
    }
}

WATMAN

What the hell, Bat­man?

Let’s try anoth­er approach. May­be we should divide each input val­ue first, then add the halves togeth­er.

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

The error here should be obvi­ous now. I’ll show you with exam­ples.

\begin{aligned} \text{avg}(2, 4) &= \frac{2}{2} + \frac{4}{2} = 1 + 2 = 3 \\ \text{avg}(3, 5) &= \frac{3}{2} + \frac{5}{2} = \lfloor 1.5\rfloor + \lfloor 2.5\rfloor = 1 + 2 = 3 \end{aligned}

Wait. So our avg func­tion 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 func­tion that works sim­i­lar­ly to (a+b)/2, the gold stan­dard of mid­point­ing, we must sat­is­fy the fol­low­ing:

  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 actu­al­ly spec­i­fies a whol­ly dif­fer­ent and yet unex­plained require­ment: round­ing the result towards zero. I just feel that if you’re gonna be aver­ag­ing signed num­bers togeth­er, the func­tion ought to be at least sym­met­ri­cal about zero. It only makes sense.

Any­ways, so how do with address our prob­lem with the “sum of halves” solu­tion? Can we fix it with­out adding a bunch of branch­es like with the for­mer “con­di­tion­al half of dif­fer­ence” solu­tion?

Sum of Halves

Well, let’s think­ing about why it hap­pens. The error man­i­fests when both halves have a half after­ing halv­ing: that is, when both a / 2 and b / 2 are something.5 (read that as “some­thing point five”). This results in both parts being round­ed down before get­ting summed. But, since both parts were decreased by .5, this caus­es a total decrease of 1 in the sum!

Of course, x / 2 needs only be round­ed if x is odd, so the only time the above error man­i­fests is then if both inputs are odd. So it’s a sim­ple mat­ter 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 sub­tle bug that caus­es it to not sat­is­fy require­ment 2. avg(a + 1, b + 1) = avg(a, b) + 1. Can you spot it?

\begin{aligned} \text{avg}(0, 3) &= 0 + \frac{3}{2} = 0 + \lfloor1.5\rfloor = 1 \\ \text{avg}(-1, 2) &= \frac{-1}{2} + 1 = \lceil-.5\rceil + 1 = 0 + 1 = 1 \end{aligned}

I’ll give it to ya straight: when a and b have oppo­site signs, they are round­ed in oppo­site direc­tions. Neg­a­tives are round­ed up, and pos­i­tives are round­ed down. Not a big deal if both inputs are odd or both are even. In fact, the error only man­i­fests when one input is odd and the oth­er is even, as well as anoth­er con­di­tion that I haven’t even yet been able to divine. Suf­fice to say that it’s not real­ly the right path to go down.

What a mess.

So what can we do? The eas­i­est thing I fig­ured out is to sim­ply use a divi­sion that con­sis­tent­ly rounds in the same direc­tion, then round the result as need­ed. To do that, we’ll use an arith­metic right shift by 1 in order to per­form a divide-by-2.

It’s a com­mon mis­con­cep­tion that the C expres­sions x / 2 and x >> 1 are equiv­a­lent. In fact, for unsigned inte­gers, they are. How­ev­er, the major dis­tinc­tion is that the lat­ter rounds towards neg­a­tive infin­i­ty on most plat­forms4! This is much more con­ve­nient for us to use, since cor­rec­tions hap­pen con­sis­tent­ly across the range of inte­gers.

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 cor­rect mid­point func­tion. How­ev­er, it’s still just a bit off my mark: it doesn’t meet our third require­ment, sym­me­try on the num­ber line about zero. That’s a sim­ple fix though. Toss in a 1 when par­i­ties don’t match but the sum of halves is neg­a­tive, and voilà, we have a win­ner! Rounds towards zero and every­thing!

/**
 * 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;
}
  1. C pro­gram­mer speak for shit hit­ting spin­ning blades of air waf­tage. Java pro­gram­mers just get Index­Out­Of­Bound­sEx­cep­tion. []
  2. Pro­pellers are more roman­tic than pota­toes but can’t be French fries, fish­ing is more roman­tic than pro­pellers and pro­vides pro­tein in con­stant time. []
  3. (a&1) != (b&1) or (a%2) != (b%2) []
  4. Sad­ly, this isn’t real­ly defined by the C or C++ stan­dards. []