## Feb 17, 2012

### on the algorithms of binomial coefficients

factorial based algorithms:
The most naive definition of binomial coefficients is: ${ n \choose m } = \frac{n!}{m! (n-m)!},$ where $n!$ is factorial of $n$. However this formula is not a practical algorithm. The reason is the value of factorial increases extremely fast. For example if we want to calculate ${13 \choose 6 }$. Using above algorithm, $13! = 6227020800 > 2^{32} = 4294967296$ overflows the typical unsigned 4-byte integer. On the other hand, ${13 \choose 6 }= 1716$. In addition, the time complexity is $T(n, m) = O( n )$.

A more practical variation of this formula is: ${ n \choose m } = \prod_{k=1}^m \frac{n-m+k}{k} = \frac{n}{m} \cdot \frac{n-1}{m-1} \cdots \frac{n-m+1}{1}.$ Also notice that ${ n \choose m } = { n \choose n - m }$. One can choose the smaller one to do computing: ${ n \choose m } = \prod_{k=1}^{\bar{m}} \frac{n-\bar{m}+k}{k}$ where $\bar{m} = \min\{ m, n-m \}$. It can be shown the time complexity is only $T(n, m) = O( \bar{ m } )$.

recursive algorithms:
A recursive algorithm based on above algorithm is, ${ n \choose m } = \frac{n}{m} \cdot { n - 1 \choose m - 1}$ with initial values ${k \choose 0} = { 1 \choose l } = 1 \quad k, l \in \mathbb{N}, l \le 1$. $T(n,m) = O(\bar{m})$. Another recursive algorithm, known as Pascal's triangle (Meru-prastaara, Khayyam triangle, Yang Hui's triangle, Tartaglia's triangle) is, ${n \choose m} = { n - 1 \choose m - 1} + { n - 1 \choose m }$ with the same initial values above. $T(n,m) = O( n^2 )$.

falling factorial and gamma function based algorithm for generalized binomial coefficients:
Another thing one often encounters is the generalized binomial coefficients. For instance, the generalized binomial theorem states: $(x + y)^\alpha = \sum_{k=0}^\infty {\alpha \choose k} x^k y^{\alpha - k} \qquad ( x, y \in \mathbb{C}, | x | < | y |, \alpha \in \mathbb{R} )$ where ${ \alpha \choose k } = \frac{\alpha \cdot (\alpha - 1) \cdots ( \alpha - k + 1)}{ k !} \equiv \frac{(\alpha)_k}{k!}$. Symbol $(\cdot )_k$ is termed falling factorial. The algorithm for computing falling factorial fashioned binomial coefficients is the same as the one mentioned above.

The most general definition of binomial coefficients is via Gamma function or equivalent limit of some Gamma function: ${ \alpha \choose \beta } \equiv \frac{\Gamma(\alpha + 1)}{\Gamma(\beta + 1) \Gamma(\alpha - \beta + 1)}$ with $\alpha, \beta \in \mathbb{C}$. There are several issues for using gamma function in practical computing. First of all, they are computationally complex. Therefore, only if $m$ in ${n \choose m}$ is not an integer (which is rare in practice), Gamma function definition could be used. In that case, one actually only need the initial values for $\Gamma(s) \quad (0 \le s \le 1)$. Calculation of gamma function can be obtained by Lanczos approximation:
$\Gamma(z) = \sqrt{2 \pi} (z + g + \frac{1}{2})^{z + \frac{1}{2}}e^{- (z + g + \frac{1}{2})} A_g(z)$ where, $g$ is an arbitrary constant subject to restriction $\text{Re}( z + g + \frac{1}{2}) > 0$. The reflection formula $\Gamma(z) \Gamma(1-z) = \frac{\pi}{\sin(\pi z)}$ can be used to achieve this. $A_g(z) = c_0 + \sum_{k=1}^\infty \frac{c_k}{z + k}$ where $c_k$ are pre-calculable constants. GNU Scientific Library choose $g = 7$, and the first 9 coefficients are given such that 15 correct digits are guaranteed:

$c_0 = 0.99999999999980993$,
$c_1 = 676.5203681218851$,
$c_2 = -1259.1392167224028$,
$c_3 = 771.32342877765313$,
$c_4 = -176.61502916214059$,
$c_5 = 12.507343278686905$,
$c_6 = -0.13857109526572012$,
$c_7 = 9.9843695780195716\text{E}-6$,
$c_8 = 1.5056327351493116\text{E}-7$.

Stirling's approximation and Spouge's approximation may also be used to calculate gamma function numerically, under different circumstance.

The second issue of gamma function definition, also reviewed in Lanczos approximation, is singularities at negative integers. However, binomial coefficients do not suffer from singularities. M. J. Kronenburg  (arXiv:1105.3689v1 [math.CO], 2011) extended the definition of binomial coefficients to arbitrary integer coefficients,  adopted by Mathematica 8. Let $n, m$ be integers, ${ n \choose m } = \left\{ \begin{array}{l l l} & \frac{n!}{m! (n-m)!} &\qquad (n \ge 0, 0\le m \le n ); \\ & (-1)^m { -n + m - 1 \choose m } &\qquad (n < 0, m \ge 0 ); \\ & (-1)^{n-m} { -m-1 \choose n - m} &\qquad (n < 0, m \le n ); \\ & 0 &\qquad \text{ others }. \\ \end{array} \right.$

logarithm method:
${ n \choose m }$ is bounded in, $\left(\frac{n}{k} \right) ^k \le { n \choose k } \le \left(\frac{n \text{e}}{k} \right)^k.$ This still may be a large number. In actuality, binomial coefficients emerging in intermediate processes may be very large, yet several binomial coefficients cancels out and result a relative small number at the end of the day. Thus it is inappropriate to evaluate binomial coefficients directly. The solution is to calculate $\log {n \choose m } = \sum_{k=1}^\bar{m} \log (n-\bar{m}+k) - \log(k)$ for each binomial coefficients, and add/subtract them depending on multiplication/division demand. At the end of the calculation, take exponential. Generalized binomial coefficients may not be positive. When using logarithm method, the absolute value and the sign should be considered separately.