Consider the following results:
julia> using SpecialFunctions
julia> gamma(-1.999)
500.4623295054461
julia> gamma(-2.0)
ERROR: DomainError with -2.0:
NaN result for non-NaN input.
julia> gamma(-2.001)
-499.5395437293685
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)
999.4237724845955
julia> gamma(0.0)
Inf
julia> gamma(-0.001)
-1000.5782056293585
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)
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.