Python (Numpy and SciPy) calculation with big integers

Question

I have been trying to use Python to calculate probability using the 'comb' function in the SciPy package. The code is as follows:

p = 0
for k in np.arange(8) + 1:
    p  += comb(8,k, exact = True)*k*8**(2*(2-k))*k**(2*k)
    print(p) ## used to check the value of p in each iteration
p /= 2**59

However, I got the result as nan. The following are the values of p for each iteration: result：

8.881784197e-16
2.44249065418e-15
5.76210086589e-15
1.35336620383e-14
3.16282379379e-14
-5.16631424912e-15
-5.20085316827e-15
nan

To make it work, I tried to do the division in each iteration:

p = 0
for k in np.arange(8) + 1:
    p  += comb(8,k, exact = True)*k*8**(2*(2-k))*k**(2*k)/ 2**59
    print(p)

But this didn't help, and I still got nan in the end:

8.881784197e-16
2.44249065418e-15
5.76210086589e-15
1.35336620383e-14
3.16282379379e-14
-5.16631424912e-15
-5.20085316827e-15
nan

I'm assuming this is caused by different representations of integers: the negative value and nan are probably caused by interpreting Long as float. Indeed, when I do the calculation step by step, it seems to work:

In [155]:
k = 8
comb(8,k, exact = True)*k*8**(2*(8-k))*k**(2*k) / 2**59
Out[155]:
0.00390625

In [156]:
k = 7
comb(8,k, exact = True)*k*8**(2*(8-k))*k**(2*k) / 2**59 + _
Out[156]:
0.008122931679330314

In [158]:
k = 6
comb(8,k, exact = True)*k*8**(2*(8-k))*k**(2*k) / 2**59 + _
Out[158]:
0.010721382431305493
...

Any idea how to solve this problem? Any suggestion is greatly appreciated!

Thanks!

Answer 1

Just change np.arange(8) + 1 to range(1,9).

Added explanation: The problem you are running into is that np.arange is producing int32 numbers, with a maximum value just over 2 billion. When you compute 8**np.int_(-11) this is computed as 1. / (8**np.int_(11)) and the expression in parentheses is substantially more than 2 billion, so you get an overflow.

hth.