Python Binomial Coefficient

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.


Solution 1:

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

Solution 2:

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

Solution 3:

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