iframe-proxy | Sunbelt Computer Software

History

115 lines (103 loc) · 3.44 KB

Raw

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

"""

Gamma function is a very useful tool in math and physics.

It helps calculating complex integral in a convenient way.

for more info: https://en.wikipedia.org/wiki/Gamma_function

In mathematics, the gamma function is one commonly

used extension of the factorial function to complex numbers.

The gamma function is defined for all complex numbers except

the non-positive integers

Python's Standard Library math.gamma() function overflows around gamma(171.624).

"""

import math

from numpy import inf

from scipy.integrate import quad

def gamma_iterative(num: float) -> float:

"""

Calculates the value of Gamma function of num

where num is either an integer (1, 2, 3..) or a half-integer (0.5, 1.5, 2.5 ...).

>>> gamma_iterative(-1)

Traceback (most recent call last):

...

ValueError: math domain error

>>> gamma_iterative(0)

Traceback (most recent call last):

...

ValueError: math domain error

>>> gamma_iterative(9)

40320.0

>>> from math import gamma as math_gamma

>>> all(.99999999 < gamma_iterative(i) / math_gamma(i) <= 1.000000001

... for i in range(1, 50))

True

>>> gamma_iterative(-1)/math_gamma(-1) <= 1.000000001

Traceback (most recent call last):

...

ValueError: math domain error

>>> gamma_iterative(3.3) - math_gamma(3.3) <= 0.00000001

True

"""

if num <= 0:

raise ValueError("math domain error")

return quad(integrand, 0, inf, args=(num))[0]

def integrand(x: float, z: float) -> float:

return math.pow(x, z - 1) * math.exp(-x)

def gamma_recursive(num: float) -> float:

"""

Calculates the value of Gamma function of num

where num is either an integer (1, 2, 3..) or a half-integer (0.5, 1.5, 2.5 ...).

Implemented using recursion

Examples:

>>> from math import isclose, gamma as math_gamma

>>> gamma_recursive(0.5)

1.7724538509055159

>>> gamma_recursive(1)

1.0

>>> gamma_recursive(2)

1.0

>>> gamma_recursive(3.5)

3.3233509704478426

>>> gamma_recursive(171.5)

9.483367566824795e+307

>>> all(isclose(gamma_recursive(num), math_gamma(num))

... for num in (0.5, 2, 3.5, 171.5))

True

>>> gamma_recursive(0)

Traceback (most recent call last):

...

ValueError: math domain error

>>> gamma_recursive(-1.1)

Traceback (most recent call last):

...

ValueError: math domain error

>>> gamma_recursive(-4)

Traceback (most recent call last):

...

ValueError: math domain error

>>> gamma_recursive(172)

Traceback (most recent call last):

...

OverflowError: math range error

>>> gamma_recursive(1.1)

Traceback (most recent call last):

...

NotImplementedError: num must be an integer or a half-integer

"""

if num <= 0:

raise ValueError("math domain error")

if num > 171.5:

raise OverflowError("math range error")

elif num - int(num) not in (0, 0.5):

raise NotImplementedError("num must be an integer or a half-integer")

elif num == 0.5:

return math.sqrt(math.pi)

else:

return 1.0 if num == 1 else (num - 1) * gamma_recursive(num - 1)

if __name__ == "__main__":

from doctest import testmod

testmod()

num = 1.0

while num:

num = float(input("Gamma of: "))

print(f"gamma_iterative({num}) = {gamma_iterative(num)}")

print(f"gamma_recursive({num}) = {gamma_recursive(num)}")

print("\nEnter 0 to exit...")

Sunbelt Computer Software

PL/B Language Development and Support

Navigation Menu

Sunbelt Computer Software

PL/B Language Development and Support

FilesExpand file tree

gamma.py

Latest commit

History

gamma.py

File metadata and controls