Argthtjtr

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:

1. scipy.special.binom()


2. 
   

   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).

Here's a version that actually uses the correct formula . :)

#! /usr/bin/env python



''' Calculate binomial coefficient xCy = x! / (y! (x-y)!)

'''



from math import factorial as fac





def binomial(x, y):

    try:

        binom = fac(x) // fac(y) // fac(x - y)

    except ValueError:

        binom = 0

    return binom





#Print Pascal's triangle to test binomial()

def pascal(m):

    for x in range(m + 1):

        print([binomial(x, y) for y in range(x + 1)])





def main():

    #input = raw_input

    x = int(input("Enter a value for x: "))

    y = int(input("Enter a value for y: "))

    print(binomial(x, y))





if __name__ == '__main__':

    #pascal(8)

    main()

...

Here's an alternate version of binomial() I wrote several years ago that doesn't use math.factorial(), which didn't exist in old versions of Python. However, it returns 1 if r is not in range(0, n+1).

def binomial(n, r):

    ''' Binomial coefficient, nCr, aka the "choose" function 

        n! / (r! * (n - r)!)

    '''

    p = 1    

    for i in range(1, min(r, n - r) + 1):

        p *= n

        p //= i

        n -= 1

    return p

So, this question comes up first if you search for "Implement binomial coefficients in Python". Only this answer in its second part contains an efficient implementation which relies on the multiplicative formula. This formula performs the bare minimum number of multiplications. The function below does not depend on any built-ins or imports:

def fcomb0(n, k):

    '''

    Compute the number of ways to choose $k$ elements out of a pile of $n.$



    Use an iterative approach with the multiplicative formula:

    $$frac{n!}{k!(n - k)!} =

    frac{n(n - 1)dots(n - k + 1)}{k(k-1)dots(1)} =

    prod_{i = 1}^{k}frac{n + 1 - i}{i}$$



    Also rely on the symmetry: $C_n^k = C_n^{n - k},$ so the product can

    be calculated up to $min(k, n - k).$



    :param n: the size of the pile of elements

    :param k: the number of elements to take from the pile

    :return: the number of ways to choose k elements out of a pile of n

    '''



    # When k out of sensible range, should probably throw an exception.

    # For compatibility with scipy.special.{comb, binom} returns 0 instead.

    if k < 0 or k > n:

        return 0



    if k == 0 or k == n:

        return 1



    total_ways = 1

    for i in range(min(k, n - k)):

        total_ways = total_ways * (n - i) // (i + 1)



    return total_ways

Finally, if you need even larger values and do not mind trading some accuracy, Stirling's approximation is probably the way to go.

For Python 3, scipy has the function scipy.special.comb, which may produce floating point as well as exact integer results

import scipy.special



res = scipy.special.comb(x, y, exact=True)

See the documentation for scipy.special.comb.

For Python 2, the function is located in scipy.misc, and it works the same way:

import scipy.misc



res = scipy.misc.comb(x, y, exact=True)

What about this one? :) It uses correct formula, avoids math.factorial and takes less multiplication operations:

import math

import operator

product = lambda m,n: reduce(operator.mul, xrange(m, n+1), 1)

x = max(0, int(input("Enter a value for x: ")))

y = max(0, int(input("Enter a value for y: ")))

print product(y+1, x) / product(1, x-y)

Also, in order to avoid big-integer arithmetics you may use floating point numbers, convert
product(a[i])/product(b[i]) to product(a[i]/b[i]) and rewrite the above program as:

import math

import operator

product = lambda iterable: reduce(operator.mul, iterable, 1)

x = max(0, int(input("Enter a value for x: ")))

y = max(0, int(input("Enter a value for y: ")))

print product(map(operator.truediv, xrange(y+1, x+1), xrange(1, x-y+1)))

Here is a function that recursively calculates the binomial coefficients using conditional expressions

def binomial(n,k):

    return 1 if k==0 else (0 if n==0 else binomial(n-1, k) + binomial(n-1, k-1))

I recommend using dynamic programming (DP) for computing binomial coefficients. In contrast to direct computation, it avoids multiplication and division of large numbers. In addition to recursive solution, it stores previously solved overlapping sub-problems in a table for fast look-up. The code below shows bottom-up (tabular) DP and top-down (memoized) DP implementations for computing binomial coefficients.

def binomial_coeffs1(n, k):

    #top down DP

    if (k == 0 or k == n):

        return 1

    if (memo[n][k] != -1):

        return memo[n][k]



    memo[n][k] = binomial_coeffs1(n-1, k-1) + binomial_coeffs1(n-1, k)

    return memo[n][k]



def binomial_coeffs2(n, k):

    #bottom up DP

    for i in range(n+1):

        for j in range(min(i,k)+1):

            if (j == 0 or j == i):

                memo[i][j] = 1

            else:

                memo[i][j] = memo[i-1][j-1] + memo[i-1][j]

            #end if

        #end for

    #end for

    return memo[n][k]



def print_array(memo):

    for i in range(len(memo)):

        print('t'.join([str(x) for x in memo[i]]))



#main

n = 5

k = 2



print("top down DP")

memo = [[-1 for i in range(6)] for j in range(6)]

nCk = binomial_coeffs1(n, k)

print_array(memo)

print("C(n={}, k={}) = {}".format(n,k,nCk))



print("bottom up DP")

memo = [[-1 for i in range(6)] for j in range(6)]

nCk = binomial_coeffs2(n, k)

print_array(memo)

print("C(n={}, k={}) = {}".format(n,k,nCk))

Note: the size of the memo table is set to a small value (6) for display purposes, it should be increased if you are computing binomial coefficients for large n and k.

It's a good idea to apply a recursive definition, as in Vadim Smolyakov's answer, combined with a DP (dynamic programming), but for the latter you may apply the lru_cache decorator from module functools:

import functools



@functools.lru_cache(maxsize = None)

def binom(n,k):

    if k == 0: return 1

    if n == k: return 1

    return binom(n-1,k-1)+binom(n-1,k)

The simplest way is using the Multiplicative formula. It works for (n,n) and (n,0) as expected.

def coefficient(n,k):

    c = 1.0

    for i in range(1, k+1):

        c *= float((n+1-i))/float(i)

    return c

Multiplicative formula