The gamma function extends the factorial

Confirms Gamma(n) = (n-1)! at integers, Gamma(1/2) = sqrt(pi), and the beta function’s relationship to gamma.

import math

import matplotlib.pyplot as plt
import numpy as np

from mathematicskit.special_functions import beta_function, gamma_function, log_gamma_function

Gamma at integers and at 1/2

for n in range(1, 7):
    print(f"Gamma({n}) = {gamma_function(float(n)):.1f}, (n-1)! = {math.factorial(n - 1)}")

x = np.linspace(-3.95, 5.0, 2000)
g = gamma_function(x)
g[np.abs(g) > 30] = np.nan
fig, ax = plt.subplots()
ax.plot(x, g, label=r"$\Gamma(x)$")
ns = np.arange(1, 6)
ax.plot(ns, [math.factorial(n - 1) for n in ns], "ko", label=r"$(n-1)!$")
ax.axhline(0, color="gray", lw=0.5)
ax.set_ylim(-10, 25)
ax.set_xlabel("x")
ax.set_title("The gamma function interpolates the factorial")
ax.legend()

print(f"\nGamma(0.5) = {gamma_function(0.5):.6f}, sqrt(pi) = {np.sqrt(np.pi):.6f}")

Gamma(1) = 1.0, (n-1)! = 1
Gamma(2) = 1.0, (n-1)! = 1
Gamma(3) = 2.0, (n-1)! = 2
Gamma(4) = 6.0, (n-1)! = 6
Gamma(5) = 24.0, (n-1)! = 24
Gamma(6) = 120.0, (n-1)! = 120

Gamma(0.5) = 1.772454, sqrt(pi) = 1.772454

Log-gamma for large arguments where gamma itself would overflow

print(f"\nlog(Gamma(500)) = {log_gamma_function(500.0):.4f}")

log(Gamma(500)) = 2605.1159

Beta function via the gamma-ratio identity

a, b = 3.0, 5.0
print(f"\nB({a},{b}) = {beta_function(a, b):.6f}")
print(f"Gamma({a})Gamma({b})/Gamma({a + b}) = {gamma_function(a) * gamma_function(b) / gamma_function(a + b):.6f}")

B(3.0,5.0) = 0.009524
Gamma(3.0)Gamma(5.0)/Gamma(8.0) = 0.009524