1
+ """
2
+ Gamma function is a very useful tool in math and physics.
3
+ It helps calculating complex integral in a convenient way.
4
+ for more info: https://en.wikipedia.org/wiki/Gamma_function
5
+ In mathematics, the gamma function is one commonly
6
+ used extension of the factorial function to complex numbers.
7
+ The gamma function is defined for all complex numbers except
8
+ the non-positive integers
9
+ Python's Standard Library math.gamma() function overflows around gamma(171.624).
10
+ """
1
11
import math
2
12
3
13
from numpy import inf
4
14
from scipy .integrate import quad
5
15
6
16
7
- def gamma (num : float ) -> float :
17
+ def gamma_iterative (num : float ) -> float :
8
18
"""
9
- https://en.wikipedia.org/wiki/Gamma_function
10
- In mathematics, the gamma function is one commonly
11
- used extension of the factorial function to complex numbers.
12
- The gamma function is defined for all complex numbers except the non-positive
13
- integers
14
- >>> gamma(-1)
19
+ Calculates the value of Gamma function of num
20
+ where num is either an integer (1, 2, 3..) or a half-integer (0.5, 1.5, 2.5 ...).
21
+
22
+ >>> gamma_iterative(-1)
15
23
Traceback (most recent call last):
16
24
...
17
25
ValueError: math domain error
18
- >>> gamma (0)
26
+ >>> gamma_iterative (0)
19
27
Traceback (most recent call last):
20
28
...
21
29
ValueError: math domain error
22
- >>> gamma (9)
30
+ >>> gamma_iterative (9)
23
31
40320.0
24
32
>>> from math import gamma as math_gamma
25
- >>> all(.99999999 < gamma (i) / math_gamma(i) <= 1.000000001
33
+ >>> all(.99999999 < gamma_iterative (i) / math_gamma(i) <= 1.000000001
26
34
... for i in range(1, 50))
27
35
True
28
- >>> gamma (-1)/math_gamma(-1) <= 1.000000001
36
+ >>> gamma_iterative (-1)/math_gamma(-1) <= 1.000000001
29
37
Traceback (most recent call last):
30
38
...
31
39
ValueError: math domain error
32
- >>> gamma (3.3) - math_gamma(3.3) <= 0.00000001
40
+ >>> gamma_iterative (3.3) - math_gamma(3.3) <= 0.00000001
33
41
True
34
42
"""
35
43
if num <= 0 :
@@ -42,7 +50,66 @@ def integrand(x: float, z: float) -> float:
42
50
return math .pow (x , z - 1 ) * math .exp (- x )
43
51
44
52
53
+ def gamma_recursive (num : float ) -> float :
54
+ """
55
+ Calculates the value of Gamma function of num
56
+ where num is either an integer (1, 2, 3..) or a half-integer (0.5, 1.5, 2.5 ...).
57
+ Implemented using recursion
58
+ Examples:
59
+ >>> from math import isclose, gamma as math_gamma
60
+ >>> gamma_recursive(0.5)
61
+ 1.7724538509055159
62
+ >>> gamma_recursive(1)
63
+ 1.0
64
+ >>> gamma_recursive(2)
65
+ 1.0
66
+ >>> gamma_recursive(3.5)
67
+ 3.3233509704478426
68
+ >>> gamma_recursive(171.5)
69
+ 9.483367566824795e+307
70
+ >>> all(isclose(gamma_recursive(num), math_gamma(num))
71
+ ... for num in (0.5, 2, 3.5, 171.5))
72
+ True
73
+ >>> gamma_recursive(0)
74
+ Traceback (most recent call last):
75
+ ...
76
+ ValueError: math domain error
77
+ >>> gamma_recursive(-1.1)
78
+ Traceback (most recent call last):
79
+ ...
80
+ ValueError: math domain error
81
+ >>> gamma_recursive(-4)
82
+ Traceback (most recent call last):
83
+ ...
84
+ ValueError: math domain error
85
+ >>> gamma_recursive(172)
86
+ Traceback (most recent call last):
87
+ ...
88
+ OverflowError: math range error
89
+ >>> gamma_recursive(1.1)
90
+ Traceback (most recent call last):
91
+ ...
92
+ NotImplementedError: num must be an integer or a half-integer
93
+ """
94
+ if num <= 0 :
95
+ raise ValueError ("math domain error" )
96
+ if num > 171.5 :
97
+ raise OverflowError ("math range error" )
98
+ elif num - int (num ) not in (0 , 0.5 ):
99
+ raise NotImplementedError ("num must be an integer or a half-integer" )
100
+ elif num == 0.5 :
101
+ return math .sqrt (math .pi )
102
+ else :
103
+ return 1.0 if num == 1 else (num - 1 ) * gamma_recursive (num - 1 )
104
+
105
+
45
106
if __name__ == "__main__" :
46
107
from doctest import testmod
47
108
48
109
testmod ()
110
+ num = 1.0
111
+ while num :
112
+ num = float (input ("Gamma of: " ))
113
+ print (f"gamma_iterative({ num } ) = { gamma_iterative (num )} " )
114
+ print (f"gamma_recursive({ num } ) = { gamma_recursive (num )} " )
115
+ print ("\n Enter 0 to exit..." )
0 commit comments