calculating binomial coefficients

Question

$\binom{n}{i}=\displaystyle\prod_{j=1}^{n-i}\frac{i+j}{j}$

$\binom{5000}{4}=26010428123750$

$\binom{100000}{60}=1.18069197996257*10^{218}$

$\binom{1000}{500}=2.702882409454366*10^{299}$

n = [5000, 100000, 1000,]
i = [4, 60, 500]
for k in range(3):
   sum1 = 1
   for j in range(n[k]-i[k]):
       sum1 = sum1 * (i[k]+j+1)/(j+1)
   print (sum1)

My calculations:

$26010428123749.992$

$1.180691979962546e+218$

$2.7028824094543666e+299$

Is there a reason to why i get $...49.992$?

Why do I have to use floating numbers?

is it now possible to get "overflow" if the binomial coefficient that I have to calculate is less than the biggest floating number which can be reperesented on my machine?

Answer 1

The inaccuracy in the results is due to the limited resolution of floating-point, only able to represent numbers of $52+1$ bits, about $16$ exact digits.

Note that such a high accuracy is more than enough for real life applications.

Answer 2

I've no experience with Python and its integer-restrictions; but I think the following product formula should allow you to compute at least the first one correctly in integers (just reformulate the loop from Pari/GP-notation)

 n=5000
 i=4

 prd = 1 
 for(k= 1,i, prd = prd* (n+1-k)/k ) 

 print(prd," by product formula")
 print(binomial(n,i)," by binomial function")

giving the result

 26010428123750  by product formula
 26010428123750  by binomial function