How to compute gamma functions with SciPy

The gamma function appears in probability distributions, factorial extensions, normalization constants, and scientific models. SciPy provides gamma-family functions in scipy.special so Python code can evaluate those formulas with arrays, real values, and complex values without hand-coding numerical approximations.

The direct gamma() ufunc returns scalar or array values and follows the identity gamma(n) = (n - 1)! for positive integers. For real calculations where products become too large, gammaln() keeps the logarithm of abs(gamma(x)) in floating-point range, and rgamma() evaluates reciprocal gamma terms without dividing by a pole.

Zero and negative integers are poles of the gamma function. SciPy 1.15 and newer returns inf at positive zero and nan at negative integer poles, so formulas with gamma terms in a denominator should use rgamma() rather than dividing by gamma() after the fact.

1. Create a Python script named gamma_function_demo.py.

   gamma_function_demo.py

   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]))

   gamma(), gammaln(), and rgamma() are vectorized ufunc objects, so the same call works for scalars and NumPy arrays.

2. Run the script and compare the printed identities.

   $ python3 gamma_function_demo.py
   gamma(n): [ 1.  1.  2.  6. 24.]
   expected: [ 1.  1.  2.  6. 24.]
   factorial identity: True
   gamma(0.5): 1.77245385
   sqrt(pi): 1.77245385
   gammaln(50): 144.56574395
   exp(gammaln(5)): 24.00000000
   rgamma([0, -1, 4]): [0.         0.         0.16666667]
   gamma([0, -1, -2]): [inf nan nan]

3. Check positive integer inputs against the factorial identity.

   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)

   For positive integer n, gamma(n) equals (n - 1)!. Matching the array result catches accidental off-by-one use when replacing factorial formulas.

4. Use gamma(0.5) to verify the square-root-of-pi case.

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

   gamma(0.5) equals sqrt(pi), which gives a compact scalar check before using gamma values in a larger formula.

5. Use gammaln() when direct gamma values would overflow or when products should stay in log space.

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

   Only exponentiate a gammaln() result when the final value fits the numeric range needed by the calculation.

6. Use rgamma() for reciprocal gamma factors and denominator terms.

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

   rgamma(x) evaluates 1 / gamma(x) directly and returns zero at nonpositive integer poles, which avoids nan propagation from denominator expressions.

7. Check pole handling before reusing values from user-supplied inputs.

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

   In SciPy 1.15 and newer, gamma(0.0) returns inf and gamma() at negative integer poles returns nan. Filter or handle those inputs before passing results into later arithmetic.

8. Remove the demo script after adapting the calculation.

   $ rm gamma_function_demo.py