Implementing banker's rounding for integer division

Question

What would be the simplest formula to implement int divround(int a, int b) {...}, in C, where output is a/b with banker's rounding (round half to even)?

For example, divround(3,2) and divround(5,2) both evaluate to 2.

I'm writing embedded code, so I cannot rely on libraries. I want the code to be generic for ARM and RISC-V, so no assembly as well. I'm trying to mimick the behavior of np.around(a/b) in NumPy (which performs round half to even), so I can exactly compare output test vectors from Python and my embedded application.

Answer 1

Assuming a and b are nonnegative, and b is less than INT_MAX / 2, the following is a simple implementation:

int divround(int a, int b)
{
    int q = a / b;
    int r = a % b;

    if (2*r > b) {
        return q + 1;
    } else if (2*r < b) {
        return q;
    } else if (q % 2 == 0) {
        return q;
    } else {
        return q + 1;
    }
}

For an explanation:

• If the remainder of division is more than half of b round up and if it's less round down.
• If the remainder is b/2 exactly, then check if the quotient is even, and if so round down otherwise round up.

If you also want to handle negative values, first determine the sign of the result and normalize the values to nonnegative, then perform the above while multiplying the sign to the result. Works assuming b is between INT_MIN / 2 and INT_MAX / 2.

int divround(int a, int b)
{
    int s_a = 1, s_b = 1;

    if (a < 0) {
        a = -a;
        s_a = -1;
    }
    if (b < 0) {
        b = -b;
        s_b = -1;
    }

    int q = a / b;
    int r = a % b;
    int sign = s_a * s_b;

    if (2*r > b) {
        return sign * (q + 1);
    } else if (2*r < b) {
        return sign * q;
    } else if (q % 2 == 0) {
        return sign * q;
    } else {
        return sign * (q + 1);
    }
}

And of course, don't divide by 0.

Answer 2

If the result is non-negative (a and b have same sign, i.e. both are non-negative or both are negative), then the normally rounded division result is

  x = (a + b/2) / b;

where an exact half is rounded towards +infinity.

If the result is negative, then the normally rounded division result is:

  x = (a - b/2) / b;

where an exact half is rounded towards -infinity.

This is equal to the result of "banker's rounding" except if x is odd and the rounding was from an exact half. An exact half is obtained if "a/b=x-½" or equivalently, "b * x = a + b/2" (in the "non-negative" case) or "b * x = a - b/2" (in the "negative" case) and if that happens, then x must be decremented/incremented by 1 to get the correct result.

Notice that if b is odd, then with a and b being integers, it is not possible for "a/b" to have a fractional part of exactly ½. So the intermediate result only may need correction if b is even and in that case, the conditions above involve only integer terms.

All in all, we have:

int divround(int a, int b)
{
    int x;
    if((a<0) == (b<0))
    {   // a and b have same sign
        x = (a + b/2) / b;
        if((b%2 == 0) && (x%2 != 0) && (b*x == a + b/2))  // b even, x odd and a/b=x-1/2
        {
            x--;
        }
    }
    else
    {   // a and b have different sign
        x = (a - b/2) / b;
        if((b%2 == 0) && (x%2 != 0) && (b*x == a - b/2))  // b even, x odd and a/b=x+1/2
        {
            x++;
        }
    }
    return x;
}

Of course, if a and b are assumed to be non-negative, then only the first case in the if-else-statement is needed.

A few test cases:

divround(3, 2) == 2
divround(-3, 2) == -2
divround(3, -2) == -2
divround(-3, -2) == 2
divround(5, 2) == 2
divround(-5, 2) == -2
divround(5, -2) == -2
divround(-5, -2) == 2
divround(5, 1) == 5
divround(5, -1) == -5
divround(5, 3) == 2
divround(5, -3) == -2
divround(0, 2) == 0
divround(0, -2) == 0

Note: this solution relies on a+b/2 and a-b/2 to not overflow or underflow, so do not walk too close to the edge.

Answer 3

For b>0 and any a, the following works:

return (a+(b/2) - (~(b|a/b) &1))/ b

Logic:

• Do normal rounding: (a+(b/2))/b
• If b is even ~b&1 AND a/b is even ~(a/b)&1, (put together: (~(b|a/b) &1) ), round down

Answer 4

Here's a branchless version of @dbush's answer (with test cases), which handles positive or negative a and b. If you only have non-negative a and positive b, remove the sign handling part. (Or make that part branchy if negative inputs are rare for most uses.)

This has the same limitations as dbush's version, specifically that INT_MIN/2 < b < INT_MAX/2 (or thereabouts) is required. This is fewer instructions on ARM with clang or GCC than @njuffa's fully general-case version, but still quite a few total to handle the signs as well as the rounding.

#include <stdio.h>

int divround( int a, int b ) {
    printf( "%3d / %2d", a, b );
    int sign = ((a<0)^(b<0))*-2+1; // if 'a' and 'b' opposite signs then -1, else 1

  // This can have signed overflow; making both negative also works, see below
    a = a<0 ? -a : a;    // abs(a)
    b = b<0 ? -b : b;    // abs(b)

    int q = a / b; // quotient (will be +'ve)
    int r = a % b * 2; // remainder ... doubled!!
    return sign * ( q + (r==b)*(q%2) + (r>b) ); // see comments below...
}

int main( void ) {
    typedef int P[2];
    P arr[] = { { 12, 4 }, { 13, 4 }, { 14, 4 }, { 15, 4 } };

    for( P *p = arr; p < arr + sizeof arr/sizeof arr[0]; p++ ) {
        printf( " = %d\n",   divround( +(*p)[0], +(*p)[1] ) );
        printf( " = %d\n",   divround( +(*p)[0], -(*p)[1] ) );
        printf( " = %d\n",   divround( -(*p)[0], +(*p)[1] ) );
        printf( " = %d\n\n", divround( -(*p)[0], -(*p)[1] ) );
    }

    return 0;
}

Output:

 12 /  4 = 3
 12 / -4 = -3
-12 /  4 = -3
-12 / -4 = 3

 13 /  4 = 3
 13 / -4 = -3
-13 /  4 = -3
-13 / -4 = 3

 14 /  4 = 4
 14 / -4 = -4
-14 /  4 = -4
-14 / -4 = 4

 15 /  4 = 4
 15 / -4 = -4
-15 /  4 = -4
-15 / -4 = 4

Concerning return sign * ( q + (r==b)*(q%2) + (r>b) );

1. Multiplying computed quotient by sign restoring possible negative.
2. 'q' is the calculated (integer division) quotient (that will be +'ve).
3. (r==b)*(q%2) If 2xRemainder r equals divisor, then this is exactly 0.5 remainder. Add 1 to q, but only if q is odd.
4. Otherwise, add 1 to q when 2xRemainder is GT the divisor.

---

Concerned about trying to change INT_MIN to a +'ve value? The following changes work around that problem. Instead of dividing +'ve numbers, divide -'ve numbers instead...

    a = a*(a<0) + -a*(a>=0); // if 'a' +'ve, then make -'ve
    b = b*(b<0) + -b*(b>=0); // same for 'b'
    /* ... */
    return sign * ( q + (r==b)*(q%2) + (r<b) ); // Now 'r<b'

---

An earlier version of this answer did the absolute-value part with the following code:

    a = -a*(a<0) + a*(a>=0); // change sign of 'a' if negative
    b = -b*(b<0) + b*(b>=0); // same for 'b'

But that doesn't compile as well with current GCC or clang for some ISAs. For example, GCC for x86-64 was using two CMOV instructions for each. And clang for ARM was using a mul and multiple other instructions, instead of a simple cmp / rsbmi reg, #0 to conditionally subtract from 0. See it on Godbolt.

Another way to avoid signed overflow UB is to use unsigned math in generating the abs result, like unsigned ua = a<0 ? 0U - a : a;. This relies on the well-defined range-reduction semantics for conversion of negative integer types to unsigned, so is well-defined everywhere, not relying on 2's complement int.

---

Compilers generally make branchless asm from ternary operators when the expressions in both halves are simple and have no side-effects, but technically only the chosen side of a ternary is evaluated in the abstract machine. (If that's a problem for more complex stuff, assign inputs to temporary variables and select between them with a ternary.)

Compilers will sometime make branchy asm from C source like return q + ( (r==b)*(q%2) + (r>b) );, as shown for the bankers_div_Fe2O3 function in the Godbolt link. (That version doesn't have sign handling; it assumes both are non-negative.)

---

Handling of sign could be optimized some if you're willing to assume 2's complement and that >> is an arithmetic right shift, as described in comments below.

Get 0 or -1 from (a^b) >> CHAR_BIT * sizeof(a) assuming a 2's complement system and arithmetic right shift.

Alternatives include a^b < 0 ? -1 : 0. Use it with a 2's complement identity bithack like you would for abs: (untested) return mask ^ (x + mask); which is either ~(x-1) or (x+0)

Compilers do see through ((a<0)^(b<0))*-2+1 and actually compile it to ((a^b)>>31) | 1 (for 32-bit int). But if you're assuming 2's complement, it simplifies nicely to just XORing and checking the sign bit of the result to see if they're different.

It saves one instruction to just get 0 or -1 (leave out the OR of 1), but it costs ADD/XOR at the end instead of imul. So same number of total instructions, but no multiply. This is a bigger gain on embedded systems with slow multipliers.

Answer 5

Below is an ISO-C99 implementation of divround() that (best I can determine) is free of spurious signed integer overflow in intermediate computations. It may contain code that relies on the two's complement representation of signed integers. On asker's intended target architectures int uses a two's complement representation, so this should not be a limitation in practice.

The division first computes a preliminary integer quotient q and a corresponding remainder r. Based on the magnitude of r relative to the divisor b a correction of magnitude 1 may need to be applied to q, where the sign of this correction is identical to that of the mathematical quotient. The rounding to nearest-or-even is implemented via three status bits that correspond directly to the conditions of the rounding mode, as well-known from floating-point computation: a round bit, a sticky bit, and the least significant bit of the unrounded result (the truncated quotient).

The round bit and the sticky bit take a bit of care to compute. Because r can be as large as b-1, one cannot use 2*r to check whether the fractional part of the quotient is equal to or greater than half. Instead, we need to compare the magnitude of the remainder r with b/2 while taking the discarded least significant bit of the divisor b into account separately.

My test scaffolding below uses a x86-specific reference function. By computing in extended precision with 64 significand bits, an accurate result is computed free of double-rounding issues: Samuel A. Figueroa. "When is Double Rounding Innocuous?" SIGNUM Newsl., 30(3):21–26, July 1995. Integer division typically fails and raises an exception when (1) the divisor is zero, or (2) the dividend is INT_MIN while the divisor is -1. My divround() does not check for these cases, and the test scaffolding avoids them.

The divround() implementation below should result in branchless code on any commonly-used compiler when compiled with full optimizations. My historical experience is that a /- operator and a %-operator sharing the same divisor do not necessarily result in a single idiv operation on architectures where an integer division returns both results, so where efficiency is required one would want to double check the generated assembly code.

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <limits.h>

int divround (int a, int b) 
{
    int q = a / b;
    int r = a % b;
    int abs_b_half = abs (b / 2);
    int quot_sign = ((a < 0) ^ (b < 0)) ? (-1) : 1;
    int abs_r = abs (r);
    int b_is_odd = b & 1;
    int q_is_odd = q & 1;
    int q_round  = abs_r >= (abs_b_half + b_is_odd);
    int q_sticky = b_is_odd | (abs_r != abs_b_half);
    q = q + quot_sign * (q_round & (q_sticky | q_is_odd));
    return q;
}

/* x86 specific! Computing in extended precision guarantees correct result */
int divround_ref (int a, int b) 
{
    uint16_t origcw;
    uint16_t newcw = 0x37f; // PC = extended precision, RC = to nearest or even
    int r;

    __asm fstcw  [origcw];
    __asm fldcw  [newcw];
    __asm fild   dword ptr [a];
    __asm fild   dword ptr [b];
    __asm fdivp  st(1), st;
    __asm frndint;
    __asm fistp  dword ptr [r];
    __asm fldcw  [origcw];
    
    return r;
}

// George Marsaglia's KISS PRNG, period 2**123. Newsgroup sci.math, 21 Jan 1999
// Bug fix: Greg Rose, "KISS: A Bit Too Simple" http://eprint.iacr.org/2011/007
static uint32_t kiss_z=362436069, kiss_w=521288629;
static uint32_t kiss_jsr=123456789, kiss_jcong=380116160;
#define znew (kiss_z=36969*(kiss_z&65535)+(kiss_z>>16))
#define wnew (kiss_w=18000*(kiss_w&65535)+(kiss_w>>16))
#define MWC  ((znew<<16)+wnew )
#define SHR3 (kiss_jsr^=(kiss_jsr<<13),kiss_jsr^=(kiss_jsr>>17), \
              kiss_jsr^=(kiss_jsr<<5))
#define CONG (kiss_jcong=69069*kiss_jcong+1234567)
#define KISS ((MWC^CONG)+SHR3)

int main (void)
{
    int a, b, res, ref;
    unsigned long long int c = 0;
    unsigned long long int m1 = 0x5555;
    unsigned long long int m2 = 0xaaaa;

    do {
        c++;

        a = (int)KISS;
        if ((c & m1) == m1) a = INT_MIN; // increase likelihood minimum dividend
        do {
            b = (int)KISS;
            if ((c & m2) == m2) b = INT_MIN; // incr. likelihood minimum divisor
        } while ((b == 0) ||                     // division by zero
                 ((a == INT_MIN) && (b == -1))); // integer division overflow

        res = divround (a, b);
        ref = divround_ref (a, b);
        if (res != ref) {
            printf ("mismatch: a=% 11d b=% 11d res=% 11d ref=% 11d %.16f\n", a, b, res, ref, (double)a/(double)b);
            return EXIT_FAILURE;
        }

        if ((c & 0xffffff) == 0) printf ("\r%llx ", c);
    } while (c);
    return EXIT_SUCCESS;
}