Floating-point oblivious way to compute multiset numbers

I have to compute $$R = \left(\!\!{n + 1\choose k}\!\!\right)$$, which happens to be:

$$R = \left(\!\!{n+1\choose k }\!\!\right) = \binom{n+k}{k} = \frac{(n + k)!}{n!k!} = \frac{(n+1)(n+2)\cdots(n+k)}{1\cdot 2\ \cdots\ k}.$$

The problem is that, if you compute the numerator and denominator separately, you will run very soon into overflow territory, so a better way to calculate this division, since both parts of the fraction have same number of components, is by: $$R = \left(\frac{n+1}{1}\right)\left(\frac{n+2}{2}\right)\cdots\left(\frac{n+k}{k}\right) = \prod_{i=1}^k\frac{n+i}{i} = \prod_{i=1}^k\left(\frac{n}{i} + 1\right)$$

The problem of this approach is that you have to apply floating point arithmetic to do the divisions, which are way more inefficient than using simple integral arithmetic.

But, knowing that the final result would be integral, is there any way to split the big fraction into a sequence of products of fractions so that each fraction has a numerator that is a multiple of its denominator? The idea is to remove any floating point operation to calculate $$R$$.

Has this formula any property that guarantees that such sequence exists?

 Ok, maybe not "any" way would be valid. If the process to avoid the floating point arithmetic would cost more than the floating point arithmetic itself, it won't be worthy.

• @D.W. I think it's ok. The multiset number $\left(\!\!{n + 1\choose k}\!\!\right)$ translates as $\binom{(n + 1) + k - 1}{k} = \binom{n+k}{k}$ which in turns translates as $\frac{(n+k)!}{(n+k-k)!k!} = \frac{(n+k)!}{n!k!}$. Is there something I'm not seeying here? Sep 19 '20 at 8:18
• Sorry, my fault. I read too quickly.
– D.W.
Sep 19 '20 at 8:33
• @D.W. No problem. Sep 19 '20 at 8:44

One common way to work with binomial coefficients without overflow is to use prime factorisation.

Legendre discovered the formula for this, and it is worth proving it for yourself. The factorial of a number can be factorised as powers of primes:

$$n! = p_1^{q_1} p_2^{q_2} \cdots$$

where:

$$q_i = \sum_{k=1}^{\left\lfloor \log_{p_i} n \right\rfloor} \left\lfloor \frac{n}{p_i^k}\right\rfloor$$

Instead of multiplying and dividing, therefore, you can just add and subtract exponents of primes.

There's the trick that is typically used when FFT is used for multiplication of large integers: If your calculation produces a floating-point number x, and you know that the mathematically correct result is an integer, and you can prove through careful analysis that the total rounding error in your result is less than 0.5, then x rounded to the nearest integer is the correct result.

In your case, assuming you use IEEE754 double precision floating point which can represent integers <= $$2^{53}$$ exactly, I wouldn't evaluate all the values (n + i) / i, but calculate (n + i) (n + i + 1) (n + i + 2) ... / ((i (i+1) (i + 2)) where you multiply as many numbers as fit into 53 bit, then calculate the quotients with a single floating-point division giving a single rounding error, and multiply the products with a single rounding error.

Say n = 100, k = 50, n+k ≤ 150, then the product of 7 consecutive numbers is less than $$2^{53}$$. So you have 8 groups of 2 times up to 7 integer products, 8 floating point division, and 7 floating point multiplications.

• I think there is no guarantee that the expression you list is an integer, if you only multiple some of them. Am I misunderstanding something?
– D.W.
Sep 19 '20 at 3:52
• Yes, it's only the final product that is mathematically an integer. Instead of dividing (n + i) / i using k floating-point divisions, you can use fewer divisions, which reduces the amount of floating-point operations and therefore the rounding error and therefore an upper bound for the rounding error. Sep 19 '20 at 9:32
• Isn't there a risk that the partial rounding errors you are evicting accumulates so that the final result is not the correct one? Sep 19 '20 at 14:11