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

**URL:** <https://discourse.julialang.org/t/gamma-function-at-a-negative-whole-number-a-cautionary-tale/41307>\
**Category:** General Usage\
**Created:** [June 13, 2020, 12:58am UTC](https://discourse.julialang.org/t/gamma-function-at-a-negative-whole-number-a-cautionary-tale/41307 "2020-06-13T00:58:42Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![billmclean](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/billmclean/32/3790_2.png) [@billmclean](https://discourse.julialang.org/u/billmclean)\
**Post date:** [June 13, 2020, 12:58am UTC](https://discourse.julialang.org/t/gamma-function-at-a-negative-whole-number-a-cautionary-tale/41307/1 "2020-06-13T00:58:42Z")

</div>

Consider the following results:

```julia
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
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
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
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
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
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,

```julia
>> 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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [June 13, 2020, 2:26am UTC](https://discourse.julialang.org/t/gamma-function-at-a-negative-whole-number-a-cautionary-tale/41307/2 "2020-06-13T02:26:28Z")

</div>

> [@billmclean](#):
>
> Why doesn’t `gamma(0.0)` throw a domain error?

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

---

<div class="post-metadata">

**Author:** ![billmclean](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/billmclean/32/3790_2.png) [@billmclean](https://discourse.julialang.org/u/billmclean)\
**Post date:** [June 13, 2020, 3:09am UTC](https://discourse.julialang.org/t/gamma-function-at-a-negative-whole-number-a-cautionary-tale/41307/3 "2020-06-13T03:09:00Z")

</div>

OK. I see so

```julia
julia> gamma(0.0)
Inf

julia> gamma(-0.0)
-Inf

```