81

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

if y == 1 or y == x:
    print(1)

if y > x:
    print(0)        
else:
    a = math.factorial(x)
    b = math.factorial(y)
    div = a // (b*(x-y))
    print(div)  

This binomial coefficient program works but when I input two of the same number which is supposed to equal to 1 or when y is greater than x it is supposed to equal to 0.

9
  • 1
    What do you need help with? The formula you're using for binomial coefficients doesn't look quite right, is that it?
    – Joni
    Oct 25, 2014 at 8:53
  • why you are using while ? you can merely use if !!
    – Mazdak
    Oct 25, 2014 at 8:55
  • when I input a number greater than x it comes up with an error or if x and y are equal to each other
    – jcoke
    Oct 25, 2014 at 8:56
  • Enter a value for x: 1 Enter a value for y: 1 1 Traceback (most recent call last): File "D:\CE151 Computer Programming\ass1.py", line 122, in <module> elif len(line)==1 and "1"<=line<="8": exlist[int(line)]() File "D:\CE151 Computer Programming\ass1.py", line 83, in ex4 div = (a//(b*(x-y))) ZeroDivisionError: integer division or modulo by zero
    – jcoke
    Oct 25, 2014 at 8:57
  • 1
    You're in luck that python has unlimited length integers, because in every other language computing two enormous factorials and then dividing to give a fairly modest result will lead to overflow and other problems. In other words: maybe use a different way of computing. Nov 14, 2018 at 13:54

15 Answers 15

162

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

2
  • 8
    I'll just add a caution that scipy.special.binom returns a floating point approximation. This is good enough for most applications, but might not suffice for theoretical purposes. Jun 27, 2016 at 1:30
  • 5
    Scipy also offers the comb function, which can be used to compute exact values. Dec 7, 2017 at 1:39
128

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
32

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
7
  • 1
    Bah, who cares if the formula is correct. As long as my code doesn't throw any exceptions, it's fine, right? Right? :) Oct 25, 2014 at 9:50
  • @TimPietzcker Hey, it's not your fault if the customer gives you the wrong specs for the software they want you to write for them. :)
    – PM 2Ring
    Oct 25, 2014 at 9:52
  • 3
    Your (second) code seems to work with integer division as well, which is a great feature, since for n greater than (around) 60, floats start to give uncorrect results due to precision errors.
    – mmj
    Jun 5, 2016 at 14:57
  • 1
    @gciriani My 2nd code block is essentially the same as your suggestions except it doesn't call math.factorial, which didn't exist when I originally wrote that code. The code in the 1st block is probably faster for most values of n & r, since it runs at C speed (Python for loops are relatively slow), although it will slow down a bit when the values get larger than the machine integer size. For really large factorials or binomial coefficients you can use libraries like mpmath &/or gmpy.
    – PM 2Ring
    Jan 17, 2020 at 15:18
  • 1
    @PM2Ring, I had missed the 2nd code. It's a great formula, so simple and effective. Thanks!
    – gciriani
    Feb 1, 2020 at 14:27
12

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.

3

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

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

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)))
4
  • 1
    Why is avoiding math.factorial an advantage?
    – BartoszKP
    Oct 25, 2014 at 9:48
  • @BartoszKP: May be it is not, but @pm-2ring in his answer pointed out that math.factorial doesn't exist in old Pythons, so I desided to avoid it just for fun. Anyway, I've defined product. Oct 25, 2014 at 9:53
  • @BartoszKP and firegurafiku : math.factorial() is running at C speed so it's probably much faster than solutions that use Python loops. OTOH, factorial() grows very quickly: factorial(13) is too big to fit into an int, so the much slower long arithmetic must be used. firegurafiku's algorithm is better on that score than the simple factorial-based algorithm, but it still ends up working with large numbers. Continued in next comment...
    – PM 2Ring
    Oct 25, 2014 at 11:40
  • 1
    @BartoszKP and firegurafiku : My loop-based version is the usual way to do it in languages that don't have big integers because the division in each loop keeps the cumulative value as small as possible. Also, by using min(r, n - r) it does the minimum number of loops.
    – PM 2Ring
    Oct 25, 2014 at 11:40
2

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

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.

2

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

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

1
  • Floating point will lose precision at larger numbers.
    – qwr
    Oct 28 at 20:12
1

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.

EDIT: I decided to add an explanation to this old question. As mentioned in other answers there are ready-to-use library methods like math.comb(n, k), scipy.special.comb(), or scipy.special.binom(). However, the original question was about the actual implementation.

The formula:

comb(n,k)=n*(n-1)*...*(n-k+1) / (1*2*...*k) = (n!/(n-k)!) / k!

While it is possible to use formulas blindly, calculated factorials likely become very large and require big integer arithmetic very quickly. Instead, it makes more sense to rearrange the operands:

comb(n,k)=n/1 * (n-1)/2 * (n-2)/3 ... *(n-k+1)/k

Because the multiplication of k consecutive positive integers is always divisible by k!, the partial calculation result is always an integer never exceeds the actual value of comb(n,k).

So, the solution by @assafsha that uses the multiplicative formula, unnecessarily converts components to float. Solutions, using factorials and the direct formula from a schoolbook unnecessarily use big integer arithmetic. The solution with recursion unnecessarily use (potentially large amount of) stack and the user may end up with the out-of-memory error. This is because Python doesn't have tail recursion optimization.

1
  • For those curious about why this works, it's because comb(n,k+1) = comb(n,k) * (n-k)/(k+1). (source) Jan 19 at 10:51
0

For those of you interested, here is the math.comb CPython implementation in mathmodule.c. The basic recursive formula (no memoization) is

j = k // 2
C(n, k) = C(n, j) * C(n-j, k-j) // C(k, j)

This is a repeated use of the identity C(n, k) = (n/k) * C(n-1, k-1). The rationale given is that Karatsuba multiplication, a divide-and-conquer algorithm, runs faster multiplying numbers about the same size (therefore the choice of j).

/* Calculate P(n, k) or C(n, k) using recursive formulas.
 * It is more efficient than sequential multiplication thanks to
 * Karatsuba multiplication.
 */
static PyObject *
perm_comb(PyObject *n, unsigned long long k, int iscomb)
{
    if (k == 0) {
        return PyLong_FromLong(1);
    }
    if (k == 1) {
        return Py_NewRef(n);
    }

    /* P(n, k) = P(n, j) * P(n-j, k-j) */
    /* C(n, k) = C(n, j) * C(n-j, k-j) // C(k, j) */
    unsigned long long j = k / 2;
    PyObject *a, *b;
    a = perm_comb(n, j, iscomb);
    if (a == NULL) {
        return NULL;
    }
    PyObject *t = PyLong_FromUnsignedLongLong(j);
    if (t == NULL) {
        goto error;
    }
    n = PyNumber_Subtract(n, t);
    Py_DECREF(t);
    if (n == NULL) {
        goto error;
    }
    b = perm_comb(n, k - j, iscomb);
    Py_DECREF(n);
    if (b == NULL) {
        goto error;
    }
    Py_SETREF(a, PyNumber_Multiply(a, b));
    Py_DECREF(b);
    if (iscomb && a != NULL) {
        b = perm_comb_small(k, j, 1);
        if (b == NULL) {
            goto error;
        }
        Py_SETREF(a, PyNumber_FloorDivide(a, b));
        Py_DECREF(b);
    }
    return a;

error:
    Py_DECREF(a);
    return NULL;
}

There is also a specialization for small n, using some pre-computed values. Pre-computed factorials in general are useful if computing many binomial coefficients, including mod p.


/* Number of permutations and combinations.
 * P(n, k) = n! / (n-k)!
 * C(n, k) = P(n, k) / k!
 */

/* Calculate C(n, k) for n in the 63-bit range. */
static PyObject *
perm_comb_small(unsigned long long n, unsigned long long k, int iscomb)
{
    assert(k != 0);

    /* For small enough n and k the result fits in the 64-bit range and can
     * be calculated without allocating intermediate PyLong objects. */
    if (iscomb) {
        /* Maps k to the maximal n so that 2*k-1 <= n <= 127 and C(n, k)
         * fits into a uint64_t.  Exclude k = 1, because the second fast
         * path is faster for this case.*/
        static const unsigned char fast_comb_limits1[] = {
            0, 0, 127, 127, 127, 127, 127, 127,  // 0-7
            127, 127, 127, 127, 127, 127, 127, 127,  // 8-15
            116, 105, 97, 91, 86, 82, 78, 76,  // 16-23
            74, 72, 71, 70, 69, 68, 68, 67,  // 24-31
            67, 67, 67,  // 32-34
        };
        if (k < Py_ARRAY_LENGTH(fast_comb_limits1) && n <= fast_comb_limits1[k]) {
            /*
                comb(n, k) fits into a uint64_t. We compute it as

                    comb_odd_part << shift

                where 2**shift is the largest power of two dividing comb(n, k)
                and comb_odd_part is comb(n, k) >> shift. comb_odd_part can be
                calculated efficiently via arithmetic modulo 2**64, using three
                lookups and two uint64_t multiplications.
            */
            uint64_t comb_odd_part = reduced_factorial_odd_part[n]
                                   * inverted_factorial_odd_part[k]
                                   * inverted_factorial_odd_part[n - k];
            int shift = factorial_trailing_zeros[n]
                      - factorial_trailing_zeros[k]
                      - factorial_trailing_zeros[n - k];
            return PyLong_FromUnsignedLongLong(comb_odd_part << shift);
        }

        /* Maps k to the maximal n so that 2*k-1 <= n <= 127 and C(n, k)*k
         * fits into a long long (which is at least 64 bit).  Only contains
         * items larger than in fast_comb_limits1. */
        static const unsigned long long fast_comb_limits2[] = {
            0, ULLONG_MAX, 4294967296ULL, 3329022, 102570, 13467, 3612, 1449,  // 0-7
            746, 453, 308, 227, 178, 147,  // 8-13
        };
        if (k < Py_ARRAY_LENGTH(fast_comb_limits2) && n <= fast_comb_limits2[k]) {
            /* C(n, k) = C(n, k-1) * (n-k+1) / k */
            unsigned long long result = n;
            for (unsigned long long i = 1; i < k;) {
                result *= --n;
                result /= ++i;
            }
            return PyLong_FromUnsignedLongLong(result);
        }
    }
    else {
        /* Maps k to the maximal n so that k <= n and P(n, k)
         * fits into a long long (which is at least 64 bit). */
        static const unsigned long long fast_perm_limits[] = {
            0, ULLONG_MAX, 4294967296ULL, 2642246, 65537, 7133, 1627, 568,  // 0-7
            259, 142, 88, 61, 45, 36, 30, 26,  // 8-15
            24, 22, 21, 20, 20,  // 16-20
        };
        if (k < Py_ARRAY_LENGTH(fast_perm_limits) && n <= fast_perm_limits[k]) {
            if (n <= 127) {
                /* P(n, k) fits into a uint64_t. */
                uint64_t perm_odd_part = reduced_factorial_odd_part[n]
                                       * inverted_factorial_odd_part[n - k];
                int shift = factorial_trailing_zeros[n]
                          - factorial_trailing_zeros[n - k];
                return PyLong_FromUnsignedLongLong(perm_odd_part << shift);
            }

            /* P(n, k) = P(n, k-1) * (n-k+1) */
            unsigned long long result = n;
            for (unsigned long long i = 1; i < k;) {
                result *= --n;
                ++i;
            }
            return PyLong_FromUnsignedLongLong(result);
        }
    }

    /* For larger n use recursive formulas:
     *
     *   P(n, k) = P(n, j) * P(n-j, k-j)
     *   C(n, k) = C(n, j) * C(n-j, k-j) // C(k, j)
     */
    unsigned long long j = k / 2;
    PyObject *a, *b;
    a = perm_comb_small(n, j, iscomb);
    if (a == NULL) {
        return NULL;
    }
    b = perm_comb_small(n - j, k - j, iscomb);
    if (b == NULL) {
        goto error;
    }
    Py_SETREF(a, PyNumber_Multiply(a, b));
    Py_DECREF(b);
    if (iscomb && a != NULL) {
        b = perm_comb_small(k, j, 1);
        if (b == NULL) {
            goto error;
        }
        Py_SETREF(a, PyNumber_FloorDivide(a, b));
        Py_DECREF(b);
    }
    return a;

error:
    Py_DECREF(a);
    return NULL;
}
-1

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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.