This question is old but as it comes up high on search results I will point out that scipy has two functions for computing the binomial coefficients:
scipy.special.binom()scipy.special.comb()import scipy.special # the two give the same results scipy.special.binom(10, 5) # 252.0 scipy.special.comb(10, 5) # 252.0 scipy.special.binom(300, 150) # 9.375970277281882e+88 scipy.special.comb(300, 150) # 9.375970277281882e+88 # ...but with `exact == True` scipy.special.comb(10, 5, exact=True) # 252 scipy.special.comb(300, 150, exact=True) # 393759702772827452793193754439064084879232655700081358920472352712975170021839591675861424
Note that scipy.special.comb(exact=True) uses Python integers, and therefore it can handle arbitrarily large results!
Speed-wise, the three versions give somewhat different results:
num = 300
%timeit [[scipy.special.binom(n, k) for k in range(n + 1)] for n in range(num)]
# 52.9 ms ± 107 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit [[scipy.special.comb(n, k) for k in range(n + 1)] for n in range(num)]
# 183 ms ± 814 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)each)
%timeit [[scipy.special.comb(n, k, exact=True) for k in range(n + 1)] for n in range(num)]
# 180 ms ± 649 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
(and for n = 300, the binomial coefficients are too large to be represented correctly using float64 numbers, as shown above).
This question is old but as it comes up high on search results I will point out that scipy has two functions for computing the binomial coefficients:
scipy.special.binom()scipy.special.comb()import scipy.special # the two give the same results scipy.special.binom(10, 5) # 252.0 scipy.special.comb(10, 5) # 252.0 scipy.special.binom(300, 150) # 9.375970277281882e+88 scipy.special.comb(300, 150) # 9.375970277281882e+88 # ...but with `exact == True` scipy.special.comb(10, 5, exact=True) # 252 scipy.special.comb(300, 150, exact=True) # 393759702772827452793193754439064084879232655700081358920472352712975170021839591675861424
Note that scipy.special.comb(exact=True) uses Python integers, and therefore it can handle arbitrarily large results!
Speed-wise, the three versions give somewhat different results:
num = 300
%timeit [[scipy.special.binom(n, k) for k in range(n + 1)] for n in range(num)]
# 52.9 ms ± 107 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit [[scipy.special.comb(n, k) for k in range(n + 1)] for n in range(num)]
# 183 ms ± 814 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)each)
%timeit [[scipy.special.comb(n, k, exact=True) for k in range(n + 1)] for n in range(num)]
# 180 ms ± 649 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
(and for n = 300, the binomial coefficients are too large to be represented correctly using float64 numbers, as shown above).
Note that starting Python 3.8, the standard library provides the math.comb function to compute the binomial coefficient:
math.comb(n, k)
which is the number of ways to choose k items from n items without repetition n! / (k! (n - k)!):
import math
math.comb(10, 5) # 252
math.comb(10, 10) # 1
Videos
André Nicolas correctly identifies that we should maybe use: $$\binom{n}{k}=\frac{n}{k}\binom{n-1}{k-1}$$
But we should also use:
$$\binom{n}{k}=\binom{n}{n-k}$$
And this also seems important:
$$\binom{n}{0} = \frac{n!}{(n-0)!} = \frac{n!}{n!} = 1$$
Great, we have everything we need to implement the function recursively!
I like to implement n choose k in a functional language like so:
n `choose` k
| k > n = undefined
| k == 0 = 1
| k > (n `div` 2) = n `choose` (n-k)
| otherwise = n * ((n-1) `choose` (k-1)) `div` k
Or in an imperative language like so:
f(n, k) {
if(k > n)
throw some_exception;
if(k == 0)
return 1;
if(k > n/2)
return f(n,n-k);
return n * f(n-1,k-1) / k;
}
It is pretty easy to see that both implementations run in $O(k)$ time and avoids overflow errors of fixed precision numbers much better then calculating n!/(n-k)!.
Maybe use $$\binom{n}{k}=\frac{n}{k}\binom{n-1}{k-1}.$$