In RFC: special: Behavior of `gamma` function and related functions at poles. · Issue #21810 · scipy/scipy · GitHub I propose we fix a very old bug in the gamma function (and related functions) where the behavior at the poles deviates from the C99 standard for tgamma (and also reasonable expectations). Basically gamma(x) for x a negative integer should return NaN, because the limit approaches plus or minus infinity depending on the direction in which x approaches the pole. Currently gamma(x) always returns +inf at all poles, (except that it returns -inf at x = -inf).
I’d been dragging my feet on fixing this because I was concerned it might break user code which relies on gamma(a)/gamma(b) evaluating to zero, when b is a pole. I was unsure if we might need to go through a painful deprecation process. Recently I saw in ENH: extend factorial{,2,k} to allow complex inputs by h-vetinari · Pull Request #21801 · scipy/scipy · GitHub that @h-vetinari had to go through an annoying workaround for the incorrect pole behavior and it convinced me that it’s time to fix this. The current behavior is clearly a bug in my opinion.
Feel free to join the discussion in #21810 if you have any thoughts.