WIP: add loggamma(a,z) (#283)

* add loggamma(a,z)

* note branch cut mismatch

* test large args

* more tests, fix a domainerror

Author: stevengj

docs/src/functions_list.md

@@ -57,13 +57,14 @@ SpecialFunctions.ellipe
 SpecialFunctions.eta
 SpecialFunctions.zeta
 SpecialFunctions.gamma(::Number)
+SpecialFunctions.loggamma(::Number)
+SpecialFunctions.logabsgamma
+SpecialFunctions.logfactorial
 SpecialFunctions.gamma(::Number,::Number)
+SpecialFunctions.loggamma(::Number,::Number)
 SpecialFunctions.gamma_inc
 SpecialFunctions.gamma_inc_inv
 SpecialFunctions.beta_inc
-SpecialFunctions.loggamma
-SpecialFunctions.logabsgamma
-SpecialFunctions.logfactorial
 SpecialFunctions.beta
 SpecialFunctions.logbeta
 SpecialFunctions.logabsbeta

docs/src/functions_overview.md

@@ -6,18 +6,18 @@ Here the *Special Functions* are listed according to the structure of [NIST Digi
 | Function | Description |
 |:-------- |:----------- |
 | [`gamma(z)`](@ref SpecialFunctions.gamma(::Number)) | [gamma function](https://en.wikipedia.org/wiki/Gamma_function) ``\Gamma(z)`` |
-| [`digamma(x)`](@ref SpecialFunctions.digamma)                 | [digamma function](https://en.wikipedia.org/wiki/Digamma_function) (i.e. the derivative of `lgamma` at `x`)                                                     |
+| [`loggamma(x)`](@ref SpecialFunctions.loggamma(::Number))                                           | accurate `log(gamma(x))` for large `x`                 |
+| [`logabsgamma(x)`](@ref SpecialFunctions.logabsgamma)                                           | accurate `log(abs(gamma(x)))` for large `x`                                                                                  |
+| [`logfactorial(x)`](@ref SpecialFunctions.logfactorial)                                            | accurate `log(factorial(x))` for large `x`; same as `loggamma(x+1)` for `x > 1`, zero otherwise                                                                   |
+| [`digamma(x)`](@ref SpecialFunctions.digamma)                 | [digamma function](https://en.wikipedia.org/wiki/Digamma_function) (i.e. the derivative of `loggamma` at `x`)                                                     |
 | [`invdigamma(x)`](@ref SpecialFunctions.invdigamma)   | [invdigamma function](http://bariskurt.com/calculating-the-inverse-of-digamma-function/) (i.e. inverse of `digamma` function at `x` using fixed-point iteration algorithm) |
 | [`trigamma(x)`](@ref SpecialFunctions.trigamma)     | [trigamma function](https://en.wikipedia.org/wiki/Trigamma_function) (i.e the logarithmic second derivative of `gamma` at `x`) |
-| [`polygamma(m,x)`](@ref SpecialFunctions.polygamma)  | [polygamma function](https://en.wikipedia.org/wiki/Polygamma_function) (i.e the (m+1)-th derivative of the `lgamma` function at `x`) |
+| [`polygamma(m,x)`](@ref SpecialFunctions.polygamma)  | [polygamma function](https://en.wikipedia.org/wiki/Polygamma_function) (i.e the (m+1)-th derivative of the `loggamma` function at `x`) |
 | [`gamma(a,z)`](@ref SpecialFunctions.gamma(::Number,::Number))  | [upper incomplete gamma function ``\Gamma(a,z)``](https://en.wikipedia.org/wiki/Incomplete_gamma_function) |
+| [`loggamma(a,z)`](@ref SpecialFunctions.loggamma(::Number,::Number))                                           | accurate `log(gamma(a,x))` for large arguments                 |
 | [`gamma_inc(a,x,IND)`](@ref SpecialFunctions.gamma_inc)  | [incomplete gamma function ratio P(a,x) and Q(a,x)](https://en.wikipedia.org/wiki/Incomplete_gamma_function) (i.e evaluates P(a,x) and Q(a,x)for accuracy specified by IND and returns tuple (p,q)) |
 | [`beta_inc(a,b,x,y)`](@ref SpecialFunctions.beta_inc)  | [incomplete beta function ratio Ix(a,b) and Iy(a,b)](https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function) (i.e evaluates Ix(a,b) and Iy(a,b) and returns tuple (p,q)) |
 | [`gamma_inc_inv(a,p,q)`](@ref SpecialFunctions.gamma_inc_inv)  | [inverse of incomplete gamma function ratio P(a,x) and Q(a,x)](https://en.wikipedia.org/wiki/Incomplete_gamma_function) (i.e evaluates x given P(a,x)=p and Q(a,x)=q  |
-| [`loggamma(x)`](@ref SpecialFunctions.loggamma)                                           | accurate `log(gamma(x))` for large `x`                 |
-| [`logabsgamma(x)`](@ref SpecialFunctions.logabsgamma)                                           | accurate `log(abs(gamma(x)))` for large `x`                                                                                  |
-| [`lgamma(x)`](@ref SpecialFunctions.lgamma)                                           | accurate `log(gamma(x))` for large `x`                                                                                                                          |
-| [`logfactorial(x)`](@ref SpecialFunctions.logfactorial)                                            | accurate `log(factorial(x))` for large `x`; same as `lgamma(x+1)` for `x > 1`, zero otherwise                                                                   |
 | [`beta(x,y)`](@ref SpecialFunctions.beta)                                           | [beta function](https://en.wikipedia.org/wiki/Beta_function) at `x,y`                                                                                           |
 | [`logbeta(x,y)`](@ref SpecialFunctions.logbeta)                                          | accurate `log(beta(x,y))` for large `x` or `y`      |
 | [`logabsbeta(x,y)`](@ref SpecialFunctions.logabsbeta)                                          | accurate `log(abs(beta(x,y)))` for large `x` or `y`     |

src/gamma.jl

@@ -658,16 +658,19 @@ function logabsgamma end
 Computes the logarithm of [`gamma`](@ref) for given `x`. If `x` is a `Real`, then it
 throws a `DomainError` if `gamma(x)` is negative.
 
-See also [`logabsgamma`](@ref).
+If `x` is complex, then `exp(loggamma(x))` matches `gamma(x)` (up to floating-point error),
+but `loggamma(x)` may differ from `log(gamma(x))` by an integer multiple of ``2\\pi i``
+(i.e. it may employ a different branch cut).
+
+See also [`logabsgamma`](@ref) for real `x`.
 """
-function loggamma end
+loggamma(x::Number) = loggamma(float(x))
 
 function loggamma(x::Real)
     (y, s) = logabsgamma(x)
     s < 0.0 && throw(DomainError(x, "`gamma(x)` must be non-negative"))
     return y
 end
-loggamma(x::Number) = loggamma(float(x))
 
 # asymptotic series for log(gamma(z)), valid for sufficiently large real(z) or |imag(z)|
 function loggamma_asymptotic(z::Complex{Float64})

src/gamma_inc.jl

@@ -1079,7 +1079,7 @@ function _gamma(a::Number,x::Number)
             end
         end
     end
-    return iszero(x) ? one(x)*gamma(a) : x^a * expint(1-a, x)
+    return iszero(x) ? gamma(one(x)*a) : x^a * expint(1-a, x)
 end
 
 _gamma(a::Integer, x::BigFloat) = _gamma_big(a, x)
@@ -1099,3 +1099,45 @@ function _gamma_big(a::Real,x::Real)
     ccall((:mpfr_gamma_inc, :libmpfr), Int32 , (Ref{BigFloat} , Ref{BigFloat} , Ref{BigFloat} , Int32) , z , a , x , ROUNDING_MODE[])
     return z
 end
+
+"""
+    loggamma(a,x)
+
+Returns the log of the upper incomplete gamma function [`gamma(a,x)`](@ref):
+```math
+\\log \\Gamma(a,x) = \\log \\int_x^\\infty t^{a-1} e^{-t} dt \\, ,
+```
+supporting arbitrary real or complex `a` and `x`.
+
+If `a` and/or `x` is complex, then `exp(loggamma(a,x))` matches `gamma(a,x)` (up to floating-point error),
+but `loggamma(a,x)` may differ from `log(gamma(a,x))` by an integer multiple of ``2\\pi i``
+(i.e. it may employ a different branch cut).
+
+See also [`loggamma(x)`](@ref).
+"""
+loggamma(a::Number,x::Number) = _loggamma(promotereal(float(a),float(x))...)
+loggamma(a::Integer,x::Number) = _loggamma(a, float(x))
+function _loggamma(a::Number,x::Number)
+    if a isa Real && x isa Real && !isfinite(a*x)
+        if isinf(x) && isfinite(a)
+            if x > 0 # == +Inf
+                return -one(a)*x
+            elseif x < 0
+                throw(DomainError((a,x), "loggamma will only return a complex result if called with a complex argument"))
+            end
+        elseif isinf(a) && isfinite(x)
+            if a > 0 && x ≥ 0
+                return a*one(x) # +Inf
+            elseif a < 0
+                return a*one(x) # -Inf
+            end
+        end
+    end
+    # from gamma(a,x) = x^a * expintx(1-a, x) * exp(-x):
+    iszero(x) && return loggamma(one(x)*a)
+    if x isa Real && x < 0 && a isa Integer && isodd(a)
+        # minus signs in expintx and x^a may cancel
+        return a*log(-x) + log(-expintx(1-a, x)) - x
+    end
+    return a*log(x) + log(expintx(1-a, x)) - x
+end

test/gamma_inc.jl

@@ -100,3 +100,23 @@ double(x::Complex) = ComplexF64(x)
     @test gamma(2, -Inf) == -Inf
     @test_throws DomainError gamma(2.2, -Inf)
 end
+
+@testset "upper incomplete gamma function logarithm" begin
+    for (a,z) in ((3,5), (3,-5), (3,5+4im), (3,-5+4im), (3,5-4im), (-3,5+4im), (-3,-5+4im))
+        @test exp(loggamma(a,z)) ≈ gamma(a,z) rtol=1e-13
+    end
+    @test loggamma(50, 1e7) ≅ -9.9992102133082030308153473164868727954977460876571275797855e6
+    @test real(loggamma(50, 1e7 + 1e8im)) ≅ -9.999097142860392e6
+    @test cis(imag(loggamma(50, 1e7 + 1e8im))) ≈ cis(1.0275220422549918) rtol=1e-8
+    @test real(loggamma(10+20im, -1e5 + 1e8im)) ≈ 100134.3502048662475864409896160245625409130538724704329203542339
+    @test cis(imag(loggamma(10+20im, -1e5 + 1e8im))) ≈ cis(-2.6572071454623347) rtol=1e-8
+    @test loggamma(-1e8, 1e9) ≅ -3.0723266045132171331933746054573197040165846554476396719312e9
+    @test loggamma(3, Inf) == -Inf
+    @test_throws DomainError loggamma(3, -Inf)
+    @test loggamma(Inf, 3.2) == Inf
+    @test loggamma(-Inf, 3.2) == -Inf
+    @test_throws DomainError loggamma(Inf, -3.2)
+    @test loggamma(117.3, 0) == loggamma(117.3)
+    @test loggamma(7, -300.2) ≅ log(gamma(7, -300.2))
+    @test_throws DomainError loggamma(6, -3.2)
+end