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.