from cmath import *
g = 7
lanczos_coef = [ \
0.99999999999980993,
676.5203681218851,
-1259.1392167224028,
771.32342877765313,
-176.61502916214059,
12.507343278686905,
-0.13857109526572012,
9.9843695780195716e-6,
1.5056327351493116e-7]
def gamma(z):
z = complex(z)
if z.real < 0.5:
return pi / (sin(pi*z)*gamma(1-z))
else:
z -= 1
x = lanczos_coef[0]
for i in range(1, g+2):
x += lanczos_coef[i]/(z+i)
t = z + g + 0.5
return sqrt(2*pi) * t**(z+0.5) * exp(-t) * x
from pylab import *
y = []
x = arange(-4.5, 4.5, 0.1)
for point in x:
y.append(gamma(point))
plot(x, y, 'r-')
axis([-5, 5, -20, 25])
title('gamma plot using lanczos approximation')
grid(True)
savefig('gamma.png')
