How to compute gamma functions with SciPy

The gamma function extends the factorial relationship to non-integer and complex inputs, which makes it useful in probability distributions, normalization constants, and scientific models. SciPy exposes the function as a vectorized operation that can evaluate one value or a whole NumPy array.

For positive integers, gamma(n) equals (n - 1)!, while gamma(0.5) equals sqrt(pi). Those identities provide known reference values for checking a calculation before its results feed a larger model.

Zero and negative integers are poles rather than finite results. SciPy 1.15 and newer returns positive infinity at 0.0 and nan at negative integer poles, so calculations should identify those values instead of treating every array element as finite. Use rgamma() for reciprocal gamma factors and gammaln() when large magnitudes need to stay in logarithmic form.

1. Create gamma_function_demo.py with the SciPy import and representative finite and pole inputs.

   gamma_function_demo.py

   import numpy as np
   from scipy.special import gamma
    
    
   x = np.array([0.5, 1.0, 5.0, 0.0, -1.0])

2. Append the gamma calculation and known finite reference values after the input array.

   values = gamma(x)
   expected = np.array([np.sqrt(np.pi), 1.0, 24.0])

3. Append fail-capable checks for the three defined values and both pole inputs.

   finite_values_match = np.allclose(values[:3], expected)
   poles_detected = np.isinf(values[3]) and np.isnan(values[4])
    
   if not finite_values_match or not poles_detected:
       raise RuntimeError("gamma result checks failed")

   The first comparison covers gamma(0.5), gamma(1.0), and gamma(5.0); the second check keeps the pole results separate from ordinary floating-point values.

4. Append concise output formatting after the checks.

   np.set_printoptions(precision=8)
   print("x:", x)
   print("gamma(x):", values)
   print("finite values match:", finite_values_match)
   print("poles detected:", poles_detected)

5. Run gamma_function_demo.py to execute the value and pole checks.

   $ python3 gamma_function_demo.py
   x: [ 0.5  1.   5.   0.  -1. ]
   gamma(x): [ 1.77245385  1.         24.                 inf         nan]
   finite values match: True
   poles detected: True

   The script exits with an error if a finite reference value does not match or if either pole is misclassified.