Fixed-precision (n choose k) and overflow
I recently found myself needing to compute (n choose k) with 64-bit integers. Recall that (n choose k) is equal to n!/[k!(n-k)!]. Mathematically, this is not a difficult computation, but when considered in the context of integer overflow, the problem becomes much harder.
To illustrate the problem, consider the computation of (9 choose 4) using 8-bit signed integers. We can start off by doing some straightforward cancellation, which leaves us with [9*8*7*6]/[4*3*2]. Where do we go from here though? If we multiply all of the terms in the numerator first, we get an intermediate result of /[4*3*2], which clearly does not fit in the [-128..127] range. The method that we are taught on paper is to cancel factors until none are left in the denominator, then multiply the remaining factors in the numerator to get the answer. We can write a program that effectively does the same thing, but do we really have to create vectors of terms and duplicate the hand method?
I searched high and low for information about how best to implement (n choose k) with fixed precision integers, without success. While considering the mechanics of coding the hand method, I realized that computing greatest common divisors (GCDs) would be a critical component. I then began to wonder if there might be an iterative algorithm that does not require manipulating vectors of integers. Here is what I came up with. (n choose k) is [n*(n-1)*...*(n-k+1)] / [k*(k-1)*...*1]. Let us call the vectors of terms in the numerator and denominator [C] and [D], respectively, so (n choose k) is [C]/[D].
- If (k > n/2), set k <-- n-k. This does not change the result, but it reduces the computational overhead for later steps.
- Initialize accumulators A and B for the numerator and denominator of the result to 1, so that A/B is 1/1. Note that upon completion, B will always be 1, thus leaving the result in A, but during computation, B may be greater than 1.
- While possible without overflowing A (and while [C] is non-empty), repeatedly merge the first term of [C] (call it c) into A and remove the term from [C]. This is achieved via the following steps:
- Divide g <-- GCD(c, B).
- B <-- B/g and c <-- c/g. This removes common factors.
- A <-- A*c.
- While possible without overflowing B (and while [D] is non-empty), repeatedly merge the first term of [D] into B and remove the term from [D]. This is achieved using the same algorithm as for step 3.
- If no progress was made in steps 3 or 4, fail due to unavoidable overflow.
- If [C] or [D] is non-empty, go back to step 3.
Since implementing the algorithm, I have been troubled by a seemingly simple question: does this algorithm ever fail even though the final result can be expressed without overflow? My intuition is that the algorithm always succeeds, but a proof has thus far eluded me. I have exhaustively tested the algorithm for 32-bit integers, and the algorithm never fails. Unfortunately, I really need to move on to other work, since the algorithm certainly works well enough for my needs.