How to numerically evaluate large binomial coeffcient terms?

Question

I am looking at a linear program (Eq (72) of https://arxiv.org/pdf/1807.05354.pdf). $r$ is an $n$-dimensional vector with elements $r_k$ and $d=2$. One of the constraints is

$$\sum_{k=0}^{n}\left(\begin{array}{l} n \\ k \end{array}\right)\left(\frac{1}{d}\right)^{k}\left(d-\frac{1}{d}\right)^{n-k} r_{k}=1$$

The authors are able to evaluate such terms for $n = 10,000$. However, when I choose even $n=100$, MATLAB is unable to store the coefficients accurately since the binomial coefficient is just too large and the solver starts to generate warnings.

Is there a trick to numerically evaluate these for a linear program solver?

Answer 1

As $k$ increases, $\binom nk$ may increase, but the other factors $(1/2)^k$ and $(3/2)^{n-k}$ decrease. So I suggest combining the calculations to keep the product at a moderate size: $$\sum_ka_kr_k=1$$ $$a_k=\binom nk\left(\frac12\right)^k\left(\frac32\right)^{n-k}$$ $$a_0=\left(\frac32\right)^n,\quad a_{k+1}=a_k\cdot\frac{n-k}{k+1}\cdot\frac13$$ I haven't tested this, though. And for $n=10000$ the leading term $a_0$ is still greater than $\sqrt2^{10000}=2^{5000}$, so standard 64-bit floating-point numbers are not sufficient.

Answer 2

Writing $$ {n \choose k} = \frac{n}1\ \frac{n-1}2\ \frac{n-2}3\ \ldots \ \frac{n-(k-1)}k $$ I would write a subroutine that starts by multiplying the first few terms above until it gets bigger than some predefined threshold. I would then turn to successively multiplying the above by $1/2$ until it gets below a predefined $\varepsilon >0$. In between one could also sprinkle some factors $3/2$ in order to balance the growth of the product.

Thus, unless the final result of ${n \choose k} (1/2)^k(3/2)^{n-k}$ is too big for your system to handle, you will never face any overflow.