Python Binomial Coefficient

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

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

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:
        return fac(x) // fac(y) // fac(x - y)
    except ValueError:
        return 0


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

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

Your program will continue with the second if statement in the case of y == x, causing a ZeroDivisionError. You need to make the statements mutually exclusive; the way to do that is to use elif ("else if") instead of if:

import math
x = int(input("Enter a value for x: "))
y = int(input("Enter a value for y: "))
if y == x:
    print(1)
elif y == 1:         # see georg's comment
    print(x)
elif y > x:          # will be executed only if y != 1 and y != x
    print(0)
else:                # will be executed only if y != 1 and y != x and x <= y
    a = math.factorial(x)
    b = math.factorial(y)
    c = math.factorial(x-y)  # that appears to be useful to get the correct result
    div = a // (b * c)
    print(div)  

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

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)

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.

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

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)

A bit shortened multiplicative variant given by PM 2Ring and alisianoi. Works with python 3 and doesn't require any packages.

def comb(n, k):
  # Remove the next two lines if out-of-range check is not needed
  if k < 0 or k > n:
    return None
  x = 1
  for i in range(min(k, n - k)):
    x = x*(n - i)//(i + 1)
  return x

Or

from functools import reduce
def comb(n, k):
  return (None if k < 0 or k > n else
    reduce(lambda x, i: x*(n - i)//(i + 1), range(min(k, n - k)), 1))

The division is done right after multiplication not to accumulate high numbers.

For everyone looking for the log of the binomial coefficient (Theano calls this binomln), this answer has it:

from numpy import log
from scipy.special import betaln

def binomln(n, k):
    "Log of scipy.special.binom calculated entirely in the log domain"
    return -betaln(1 + n - k, 1 + k) - log(n + 1)

(And if your language/library lacks betaln but has gammaln, like Go, have no fear, since betaln(a, b) is just gammaln(a) + gammaln(b) - gammaln(a + b), per MathWorld.)

import math

def binomial_coefficients(n,k):
      product = 1
      for i in range(k):
           product = math.floor(((product * (n - i))/(i + 1)) 
 return product 

in the calculation of the binomial coefficients, you should not calculate the finite product n(n-1) ... (n - k +1) for (n, k) and k!. This could cause an overflow error. Therefore, using a bit of number theory we can assume that the inputs will always be in integer form (since the combination of (n, k) only accepts integers)) we can see that for an integer 'i' in a product of consecutive integers, any term u in the product will always be divisible by i.

NOTE: you can do this without importing the math module. math.floor(a/b) is equivalent to a // b