Gamma function should return Inf at -1, -2, -3......

[Jinwei1987]

at -1, -2, -3... gamma should return Inf like gamma(0)

julia> gamma(0)
Inf

julia> gamma(-1.0)
ERROR: DomainError with -1.0:
NaN result for non-NaN input.
Stacktrace:
[1] nan_dom_err at ./math.jl:339 [inlined]
[2] gamma(::Float64) at /Users/zhengjinwei/.julia/packages/SpecialFunctions/fvheQ/src/gamma.jl:558
[3] top-level scope at none:0

[simonbyrne]

There isn't a well-defined limit:

julia> gamma(-0.99999999)
-9.999999992030846e7

julia> gamma(-1.00000001)
1.0000000018496278e8

(I am sympathetic to returning NaN instead)

[Jinwei1987]

NaN makes sense

If I wan to plot gamma function from [-5, 5]. I can't do it, because it throws an error. If NaN is returned then it will be fine.

right now I have to use python scipy.special.gamma instead, which return Inf at -1, -2, -3. I think it would be good not throwing any error on undefined points.

[simonbyrne]

I'd be in favour of this change: it's particularly annoying for the gamma function where the set of "valid" values doesn't make a single contiguous set. We could add a keyword option (say domainerror) to control this behaviour.

[dlfivefifty]

It should return Inf for consistency with other special functions:

julia> cot(0.0)
Inf

[simonbyrne]

That only works for because zeros are signed:

julia> cot(-0.0)
-Inf

[putianyi889]

Although NaN makes sense, it is not convenient in applications. For example, you would expect 1/gamma(-1) to be zero since -1 is a simple pole of gamma function, but it's NaN when gamma(-1) is NaN.
The disadvantage of assigning Inf to gamma(-1) is that gamma(0)+gamma(-1) should be NaN instead of Inf.
If there were distinct Inf, +Inf and -Inf, where Inf represents the infinity on the complex plane, +Inf and -Inf represent the ends of the real line,gamma(-1) should definitely be Inf.

[bobbielf2]

I agree with @putianyi889
The reciprocal function 1/gamma frequently appears in applications. Since 1/gamma is an entire function with simple zeros, it makes sense that gamma returns Inf at its simple poles, which simplifies programming when 1/gamma is needed.

This reciprocal argument also applies to @dlfivefifty 's example of cot.

And this of course simplifies plotting as @Jinwei1987 mentioned.

[bobbielf2]

An update:

I notice that exp(loggamma(n)) gives Inf for any negative integer n, is there a reason why this behavior is not used to compute gamma(n) == exp(loggamma(n))?

Also, I found the following inconsistency

julia> exp(loggamma(-1-0.01))
99.59128511327795

julia> exp(loggamma(-1+0.01))
ERROR: DomainError with -0.99:
`gamma(x)` must be non-negative

[stevengj]

@bobbielf2, the latter example is not inconsistent: gamma(-1-0.01) ≈ 99.59 but gamma(-1+0.01) ≈ -100.44, so the latter does not have a real-valued log. You need to use logabsgamma if you want to handle both cases.

This may be why loggamma(-1.0) gives Inf — only the limit from the positive side is defined in this function.

I would be in favor of returning NaN from gamma in these cases. The case for returning +Inf is more questionable.

[bobbielf2]

Hi @stevengj, thanks for the explanation. I thought that all the negative inputs to loggamma would be promoted to complex, but apparently they are not, i.e. I have to manually make it complex

julia> exp(loggamma(-0.99+0im))
-100.43695466580859 - 1.229997950475903e-14im

I guess each time I need to enforce the continuity of 1/gamma, I have to manually promote the input to complex:

julia> real(gamma(-1+0im))
-Inf

This behavior falls into "questionable" by your definition, but I honestly feel like it would be nice to have gamma(-n) = real(gamma(-n+0im)), which simplifies programming in a lot of applications. Otherwise, it would really take every newcomer a lot of hassle to get to this workaround, like what I am doing now.

[stevengj]

I have to manually make it complex

Yes, for the same reason that sqrt(-1) throws an error unless you make the input explicitly complex: this is type stability. Returning different types for different real inputs would hurt performance of every function calling gamma.

[bobbielf2]

Hi @stevengj, I see that type stability is important.
But back to the original question of giving a value to gamma(-n), I think defining gamma(-n)=Inf is good because:

1. It maintains the continuity of 1/gamma (which is needed in almost all relevant applications)
2. It is not breaking (because gamma(-n) was undefined anyway)
3. It is consistent with real(gamma(-n+0im)) or exp(prod(logabsgamma(-n)))
4. It maintains type stability
5. It makes plotting easier

If it ended up that gamma(-n)=NaN is preferred, then it would still be nice to document the alternative approach of using real(gamma(-n+0im)) in place of gamma(-n) whenever Inf is desired, so that there are less hiccup and more user-friendly for newcomers.

Thanks for the discussion.