Handling large numbers to calculate the number of combinations

Question

I wrote the following code to calculate a probability. I know the formula I used is correct. The code gives reasonable results for small values of m and n. However, for the large values of m and n, the result, sumC, is outside the interval [0,1] which is not expected. The variable iter takes large values which may be the reason. How can I handle this problem?

import math
from scipy.special import comb
n = 50
m = 100
k = 15
P=[]
sumC = 0
for j in range(k, m+1):
    if not (n-j < 0):
        iter = (-1)**(j+k) * comb(j, k, exact=True) * comb(m, j, exact=True) * math.factorial(n) * (m-j)**(n-j) /( (m)**(n) * math.factorial(n-j))
        sumC = sumC + iter
print(sumC )

Answer 1

Using the Python mpmath arbitrary precision library with 50 digits of precision produces a value in [0, 1] range.

from scipy.special import comb
from mpmath import fac, mp

mp.dps = 50; mp.pretty = True

def compute_prob(k, m, n):
  """ Original summation, but using factorial ('fac') from mpmath """
  sumC = 0
  for j in range(k, m+1):
    if not (n-j < 0):
        iter_ = (-1)**(j+k) * comb(j, k, exact=True) * comb(m, j, exact=True) * fac(n) * (m-j)**(n-j) /( (m)**(n) * fac(n-j))
        sumC = sumC + iter_

  return sumC

n = 50
m = 100
k = 15
print(compute_prob(k, m, n))

Output

0.000054845306977312907595945622368606216050228369266162