# Gamma function at a negative whole number: a cautionary tale

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.

1 Like

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

7 Likes

OK. I see so

``````julia> gamma(0.0)
Inf

julia> gamma(-0.0)
-Inf
``````