Python/maths/gamma.py
pre-commit-ci[bot] bc8df6de31
[pre-commit.ci] pre-commit autoupdate (#11322)
* [pre-commit.ci] pre-commit autoupdate

updates:
- [github.com/astral-sh/ruff-pre-commit: v0.2.2 → v0.3.2](https://github.com/astral-sh/ruff-pre-commit/compare/v0.2.2...v0.3.2)
- [github.com/pre-commit/mirrors-mypy: v1.8.0 → v1.9.0](https://github.com/pre-commit/mirrors-mypy/compare/v1.8.0...v1.9.0)

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
2024-03-13 07:52:41 +01:00

117 lines
3.4 KiB
Python

"""
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...")