1

I was comparing the result of my following python calculation with Mathematica: https://www.wolframalpha.com/input?i=sum+%28500+choose+r+%29%28-1%29%5Er+%2F%28r%21%29+%2C+r%3D0+to+500

import numpy as np
from decimal import *
import scipy.special
from scipy.special import factorial
getcontext().prec = 30

i = 500
sum(np.array([scipy.special.comb(Decimal(i), (r), exact=True)*pow(-1, r)/Decimal(factorial(r, exact=False)) for r in range(i+1)]))

Both calculations are giving almost same value until i = 400 but failing to converge after that even though I set an arbitrary precision in my python code via the decimal module. Calculation with Mathematica seems correct. May I know, how can we get same result in python as that of Mathematica for larger i?

R S John
  • 507
  • 2
  • 9
  • 16
  • 1
    Note that the factorial is recomputed from scratch in the loop which is not very efficient. you can compute it yourself or even precompute it by building a list before the operation. This should be about 10x to 20x times faster in this case. Why are you using `np.array` here? `sum` of Python does the job already well. Using Numpy array will not make that faster (in fact, it is [often the opposite in such case](https://stackoverflow.com/questions/69584027/why-is-np-sumrangen-very-slow/69587888#69587888)) and is not more convenient either here. – Jérôme Richard May 13 '22 at 11:09

1 Answers1

5

The problem is that you have exact=False in the factorial.

>>> import numpy as np
>>> from decimal import *
>>> import scipy.special
>>> from scipy.special import factorial
>>> getcontext().prec = 30
>>> 
>>> i = 500
>>> sum(np.array([scipy.special.comb(Decimal(i), (r), exact=True)*pow(-1, r)/Decimal(factorial(r, exact=False)) for r in range(i+1)]))
Decimal('-7.13859062099388393889008217957')
>>> sum(np.array([scipy.special.comb(Decimal(i), (r), exact=True)*pow(-1, r)/Decimal(factorial(r, exact=True)) for r in range(i+1)]))
Decimal('0.196589352363439561009074161963')
VRehnberg
  • 431
  • 3
  • 9