Python Problem with a gamma function/error message

[ChrisVer]

I wrote the following python script in order to try and look how $$\sum_{i=1}^\infty f(n) = \sum_{i=1}^\infty \frac{\Gamma(n+1/2)}{n^2 \Gamma(n)}$$ behaves. However it appears to have some problem with my function f(n) for n=172 and above throwing this error:

RuntimeWarning: invalid value encountered in double_scalars
return gamma(n+0.5)/(n**2 * gamma(n))

I thought that there was a problem with the function [as division with zero] and so I tried wolframalpha which returns me a non zero value for \Gamma(172) and also a non-infinite value for f(172).
I even tried a user-defined gamma function I had built in the past, and I got the same result for n=172... [it's highly unlike that my method's the same as scipy's]

from scipy.special import gamma

def f(n):
    return gamma(n+0.5)/(n**2 * gamma(n))

def sume(f,N):
    result=0.0
    for i in range(1,N+1): result+= f(i)
    return result

print f(172)
'''
for i in range(100,200):
    print i," : ",sume(f,i)
'''

 

[FactChecker]

gamma(172) = 171! is a huge number. It is probably just too large for the floating point double precision number format of the computer. The limit is about 1.8x10³⁰⁸. gamma(172) is about 1.2x10³⁰⁹. So that is where the number gets too large for the computer to handle.

 

[stevendaryl]

I wouldn't use the built-in gamma function for large arguments. Dividing two huge numbers is usually very inaccurate.

I would instead just calculate the expression inside the summation recursively

f(n) = \dfrac{\Gamma(n+\frac{1}{2})}{n^2 \Gamma(n)}

So f(n+1) = \dfrac{\Gamma(n+\frac{1}{2}+1)}{(n+1)^2 \Gamma(n+1)}

Then you use that \Gamma(x+1) = x \Gamma(x) to write this as follows:

So f(n+1) = \dfrac{(n+\frac{1}{2})\Gamma(n+\frac{1}{2})}{(n+1)^2 n \Gamma(n)} = \dfrac{n(n+\frac{1}{2})}{(n+1)^2} f(n)

 

[jim mcnamara]

Maybe consider an arbitrary precision approach - you declare your precision. http://www.sympy.org/en/index.html