Gamma function at a negative whole number: a cautionary tale

Consider the following results:

julia> using SpecialFunctions

julia> gamma(-1.999)

julia> gamma(-2.0)
ERROR: DomainError with -2.0:
NaN result for non-NaN input.

julia> gamma(-2.001)

I see the logic here: gamma(x) tends to plus infinity if x tends to -2.0 from above, but to minus infinity if x tends to -2.0 from below, so there is no sensible result when x = -2.0.
However, for x tending to zero we have the following behaviour.

julia> gamma(0.001)

julia> gamma(0.0)

julia> gamma(-0.001)

Why doesn’t gamma(0.0) throw a domain error?

What if we make the argument complex?

julia> gamma(Complex(-2.0))
-Inf - Inf*im
julia> gamma(Complex(0.0))
Inf - 0.0im

Things become even trickier when we consider the reciprocal of the Gamma function (which is often more important in applications).

julia> 1/gamma(0.0)

julia> 1/gamma(-2.0)
ERROR: DomainError with -2.0:
NaN result for non-NaN input.

julia> 1/gamma(Complex(0.0))
0.0 + 0.0im

julia> 1/gamma(Complex(-2.0))
-0.0 + 0.0im

julia> 1im/gamma(Complex(-2.0))
NaN + NaN*im

Actually, the last result just reflects the way IEEE complex arithmetic works.

julia> 1.0/Complex(Inf,Inf)
0.0 - 0.0im

julia> 1.0im/Complex(Inf,Inf)
NaN + NaN*im

I guess the compiler figures that (a+ib)/(c+id) has real part
(ac+bd)/(c^2+d^2) and imaginary part (bc-ad)/(c^2+d^2), which both evaluate to Inf/Inf = NaN if c = d = Inf. However, I’m not sure what is going on with a real numerator since if b is zero we should still have Inf/Inf for both real and imaginary parts, but I notice

julia> (1.0+0.0im)/Complex(Inf,Inf)
NaN + NaN*im

The impetus for this post came from translating a complicated Matlab function into Julia. In Matlab,

>> 1/gamma(-2.0)
ans = 0.0

and it took quite some time to track down why the Julia code returned NaN when the Matlab code did not.

1 Like

Because zero has a sign in floating-point arithmetic, so gamma can return infinity with the correct sign.


OK. I see so

julia> gamma(0.0)

julia> gamma(-0.0)