Gamma Functions Example

Demonstrates Gamma, log-Gamma, and digamma functions using specialfunctions.gamma(), specialfunctions.log_gamma(), and specialfunctions.digamma().

Example Code

"""Example: IMSL Special Functions — Gamma and related functions."""
from __future__ import annotations

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from specialfunctions import gamma, log_gamma, digamma

# Print some known values
print(f"Gamma(0.5) = {gamma(0.5):.10f}  (expected: {np.sqrt(np.pi):.10f})")
print(f"Gamma(1)   = {gamma(1.0):.10f}  (expected: 1.0000000000)")
print(f"Gamma(5)   = {gamma(5.0):.10f}  (expected: 24.0000000000)")
print(f"Gamma(1.5) = {gamma(1.5):.10f}  (expected: {0.5 * np.sqrt(np.pi):.10f})")
print(f"log_gamma(5) = {log_gamma(5.0):.10f}  (expected: {np.log(24.0):.10f})")
print(f"digamma(1)   = {digamma(1.0):.10f}  (expected: -0.5772156649 [Euler-Mascheroni])")

fig, axes = plt.subplots(1, 3, figsize=(15, 5))
fig.suptitle("Gamma and Related Functions", fontsize=14)

# Plot gamma on multiple segments to avoid poles
segments = [(-3.9, -2.1), (-1.9, -0.1), (0.1, 5.0)]
for xlo, xhi in segments:
    x = np.linspace(xlo, xhi, 400)
    y = gamma(x)
    y = np.where(np.abs(y) > 20, np.nan, y)
    axes[0].plot(x, y, "b-", lw=1.5)

axes[0].axhline(0, color="k", lw=0.5, ls="--")
axes[0].set_xlim(-4, 5)
axes[0].set_ylim(-15, 15)
axes[0].set_title("Γ(x)")
axes[0].set_xlabel("x")
axes[0].grid(True, alpha=0.3)

x2 = np.linspace(0.1, 10, 400)
axes[1].plot(x2, log_gamma(x2), "g-", lw=1.5, label="log Γ(x)")
axes[1].set_title("log Γ(x)")
axes[1].set_xlabel("x")
axes[1].grid(True, alpha=0.3)

x3 = np.linspace(0.5, 5, 400)
axes[2].plot(x3, digamma(x3), "r-", lw=1.5)
axes[2].set_title("ψ(x) = digamma")
axes[2].set_xlabel("x")
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("test_output/example_imsl_gamma_functions.svg", bbox_inches="tight")
print("Saved: test_output/example_imsl_gamma_functions.svg")

Plot Output

Gamma(x), log-Gamma(x), and digamma(x) plotted over their domains.

Console Output

Gamma(0.5) = 1.7724538509  (expected: 1.7724538509)
Gamma(1)   = 1.0000000000  (expected: 1.0000000000)
Gamma(5)   = 24.0000000000  (expected: 24.0000000000)
Gamma(1.5) = 0.8862269255  (expected: 0.8862269255)
log_gamma(5) = 3.1780538303  (expected: 3.1780538303)
digamma(1)   = -0.5772156649  (expected: -0.5772156649 [Euler-Mascheroni])
Saved: test_output/example_imsl_gamma_functions.svg