# How to evaluate binomial coefficients efficiently and as correctly as possible?

This question is more precisely about evaluation with a computer, of a binomial coefficient of the form $\binom{x}{m}$ where $x$ is a real number and $m$ a rational integer.

The reason why I ask is that I found out recently that sage is using the naive definition with the $\Gamma$ function, which means that it gets as a result NaN (not-a-number) with quite small parameters, for which the real result is pretty reasonable and should have been given (see the bug report).

I have proposed to change the implementation by returning zero in more cases than it already does, to reduce to a situation $\binom{x}{m}$ with $x\geq m\geq 0$, so we can write $x=m+k+u$ with $k$ a natural integer and $u\in[0;1[$, then computing the quotient $\Gamma(x+1)/\Gamma(m+1)$ with a Pochhammer symbol times the quotient $\Gamma(m+1+u)/\Gamma(m+1)$. For that last quotient, I was proposing a direct computation for small $m$ and a polynomial expansion in $u$ for big $m$.

There are two problems with this approach:

• I don't really know how big the error is, which for a numerical computation is a pretty big issue ;
• I used the naive code as a starting point, and added naive ideas to the mix : there may exist better approaches (it's also because of this point that I didn't try to evaluate the error more precisely).

It would be surprising if there existed no algorithm for this kind of computations, given how important those coefficients are in various situations...

• This question might be better suited for scicomp.stackexchange.com . – Emil Jeřábek Feb 8 '12 at 15:10
• Well, discussing a specific implementation might be better suited elsewhere, but a general discussion on algorithms seemed appropriate here. – Julien Puydt Feb 8 '12 at 15:36
• Does sage do better if you ask it to compute the beta function of the appropriate arguments? – Igor Rivin Feb 9 '12 at 15:04
• As I said, I'm more interested in discussing better algorithms than a specific implementation : I mostly wanted to explain why I came up with the question. And I'm not sure it has beta... and it's difficult to tell with such a name :-/ – Julien Puydt Feb 9 '12 at 18:20
• How about using the log-gamma function, and use subtraction instead of division? – Luis Mendo Sep 30 '13 at 11:44

## 1 Answer

Since you're using sage, I'd write a quick C or Cython function. Use the Arb library (I think this is already included in the latest release of sage) by computing the function in the naive way using the gamma function.

The benefit to doing so, is that the library uses "Ball" arithmetic and automatically finds rigorous error bounds for the computation.

Basically, the way I'd do it is by given a required precision...compute naively the result. If the relative accuracy in terms of bits isn't high enough, increase the working precision and compute again. Continue the process until the result is as accurate as required. Then, convert the arb type result to an mpfr variable by making an appropriate arb function call.