import math

import numpy as np
from scipy.special import gamma, gammaln, rgamma


np.set_printoptions(precision=8, suppress=True)

n = np.arange(1, 6, dtype=float)
gamma_values = gamma(n)
expected = np.array([math.factorial(int(value) - 1) for value in n], dtype=float)

print("gamma(n):", gamma_values)
print("expected:", expected)
print("factorial identity:", np.allclose(gamma_values, expected))

half = gamma(0.5)
print(f"gamma(0.5): {half:.8f}")
print(f"sqrt(pi): {np.sqrt(np.pi):.8f}")

print(f"gammaln(50): {gammaln(50.0):.8f}")
print(f"exp(gammaln(5)): {np.exp(gammaln(5.0)):.8f}")

print("rgamma([0, -1, 4]):", rgamma([0.0, -1.0, 4.0]))
print("gamma([0, -1, -2]):", gamma([0.0, -1.0, -2.0]))