1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
|
from cmath import exp, pi, sin, sqrt
p = [
676.5203681218851,
-1259.1392167224028,
771.32342877765313,
-176.61502916214059,
12.507343278686905,
-0.13857109526572012,
9.9843695780195716e-6,
1.5056327351493116e-7,
]
EPSILON = 1e-07
def drop_imag(z):
if abs(z.imag) <= EPSILON:
z = z.real
return z
def gamma(z):
z = complex(z)
if z.real < 0.5:
y = pi / (sin(pi * z) * gamma(1 - z)) # Reflection formula
else:
z -= 1
x = 0.99999999999980993
for (i, pval) in enumerate(p):
x += pval / (z + i + 1)
t = z + len(p) - 0.5
y = sqrt(2 * pi) * t ** (z + 0.5) * exp(-t) * x
return drop_imag(y)
"""
The above use of the reflection (thus the if-else structure) is necessary, even though
it may look strange, as it allows to extend the approximation to values of z where
Re(z) < 0.5, where the Lanczos method is not valid.
"""
print(gamma(1))
print(gamma(5))
print(gamma(0.5))
print(gamma(0.5* (1 + 1j)))
print(gamma(-0.5))
|
