Counting combinations modulo power of 2

It is pretty simple to calculate the number of ways to choose k items out of n: There are n ways to choose the first item, n-1 ways to choose the second item, and so on, until there are (n-k+1) ways to choose the kth item. By the rule of product there are n (n-1)···(n-k+1) ways to choose the k items in total.

But hold on, this process considers picking first A and then B distinct from picking first B and then A. If you don’t care about in which order the items are picked you have to remove the duplicates by dividing by the number of ways there are to order k items: this is k! = k·(k-1)···3·2·1. So, the number of ways you can choose k items out of n is:

$\binom{n}{k} = \frac{n(n-1) \cdots (n-k+1)}{k(k-1)\cdots 3 \cdot 2 \cdot 1} = \frac{n!}{k!(n-k)!}$

Of course, implementing this formula on a computer as it is is not a very good idea. If you use integers, the multiplications rapidly overflow fixed sized integers, and then the division doesn’t return the correct value because integer division is not properly defined modulo a power of two. If you use floating point types you’ll have a different set of problems. So let’s devise an algorithm that calculates $$\binom{n}{k}$$ correctly whenever it fits into a computer word, or in general, modulo a power of two.

We’ll replace the division with multiplication by the inverse modulo 2m. Here it is important to keep in mind that even numbers don’t have a multiplicative inverse, so we remove any powers of 2 from the numerator and denominator and keep track of the exponent of 2 separately. Done this, we can safely accumulate the denominator and numerator through multiplication and calculate only one modular inverse, in the end.

Without further talk, here’s a C program that implements the algorithm for 32 bit ints.

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

unsigned int inv(unsigned int n);   /* multiplicative inverse mod 2^32 */
int ord2(int n);                    /* exponent of 2 in n */

unsigned int binom(int n, int k)
{
unsigned int num=1, den=1, aux, i;
int ord=0;
if (n-k < k) k = n-k;
for (i = 1; i <= k; i++) {
aux = ord2(i);
ord -= aux;
den = den*(i>>aux);

aux = ord2(n+1-i);
ord += aux;
num = num*((n+1-i) >> aux);
}
return num*inv(den) << ord;
}

int main(int argc, char *argv[])
{
int a=3, b=1;
if (argc > 2) {
a = atoi(argv);
b = atoi(argv);
printf("C(%d,%d) = %u (mod 2^32)\n", a, b, binom(a,b));
} else  {
printf("%s a b - calculate C(a,b) mod 2^32", argv);
}
return 0;
}

unsigned int inv(unsigned int n) {
unsigned int x = 1;
x = x*(2-x*n); /* Newton's method */
x = x*(2-x*n); /* 5 steps enough for 32 bits */
x = x*(2-x*n);
x = x*(2-x*n);
x = x*(2-x*n);
return x;
}

int ord2(unsigned int n) {
unsigned int i, e;
i = n & -n; /* isolate rightmost 1-bit */
for (e=0; i > 1; e++) i >>= 1;
return e;
}

1. j_random_hacker
Posted 2013/08/10 at 23:37 | Permalink

Beautiful! I had a vague notion that multiplication can sometimes be used to effect a division (from time to time I’ve seen “magic constants” for dividing by 3, 5 etc.), but I had no idea such a straightforward implementation of choose(n, k) was possible.

A couple of questions:

1. Suppose we have two odd numbers x and y. If x is *not* a multiple of y, what happens when we multiply it by y’s inverse? Does it round down, as the usual integer division would?

2. I’m particularly surprised to see that plain old Newton’s method works to find the inverse, given that we’re working not only with integers but with integers that *suddenly wrap around*! I recall that Newton’s method roughly doubles the number of accurate digits with each iteration; does that mean that 6 iterations will suffice for 64 bits? (Side question: are there other starting points that guarantee correctness with fewer iterations?)

Many thanks. Also with your permission I’d like to post this to Reddit.

2. Joni
Posted 2013/08/12 at 16:53 | Permalink

If x is not a multiple of y, you get an ugly value when multiplying by y’s inverse. Let’s denote the inverse of y by y’. If you write x=ay+b where a is the quotient and b is the remainder when dividing x by y, you get xy’=ayy’+by’=a+by’. That is, xy’ is the correct quotient, plus the remainder times y’… and that’s not a nice number.

About Newton’s method, you’re right: with 6 iterations you would get the inverse for 64 bit integers. You can start iterating from any odd number because what matters is that at least the last bit is correct. The second to last bit is correct in half of the cases. If you can detect this can could stop the iteration early, but since the iterations are so cheap it’s not really worth the effort.

You can post this to Reddit, it’s on the public Internet :)

3. j_random_hacker
Posted 2013/08/16 at 00:25 | Permalink

Thanks for responding, that all makes sense. I suppose a well-behaved “rounding down” behaviour was too much to hope for :-P

I’ve now posted this to Reddit: http://www.reddit.com/r/programming/comments/1kg4am/smart_but_simple_way_to_avoid_overflow_when/. It doesn’t seem to have appeared yet on the main “Programming” page, but I expect it will get there in time.

Thanks!