Use mpfr_lngamma for loggamma(::BigFloat)

Project.toml

@@ -1,6 +1,6 @@
 name = "SpecialFunctions"
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "1.7.0"
+version = "1.8.0"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"

src/gamma.jl

@@ -1,6 +1,6 @@
 # This file contains code that was formerly a part of Julia. License is MIT: http://julialang.org/license
 
-using Base.MPFR: ROUNDING_MODE
+using Base.MPFR: MPFRRoundingMode, ROUNDING_MODE
 
 export gamma, loggamma, logabsgamma, beta, logbeta, logabsbeta, logfactorial, logabsbinomial
 
@@ -643,12 +643,18 @@ See also [`logabsgamma`](@ref) for real `x`.
 """
 loggamma(x::Number) = _loggamma(float(x))
 
-function _loggamma(x::Union{Real,BigFloat})
+function _loggamma(x::Real)
     (y, s) = logabsgamma(x)
     s < 0 && throw(DomainError(x, "`gamma(x)` must be non-negative"))
     return y
 end
 
+function _loggamma(x::BigFloat)
+    y = BigFloat()
+    ccall((:mpfr_lngamma, :libmpfr), Cint, (Ref{BigFloat}, Ref{BigFloat}, MPFRRoundingMode), y, x, ROUNDING_MODE[])
+    return y
+end
+
 # Compute the logΓ(z) function using a combination of the asymptotic series,
 # the Taylor series around z=1 and z=2, the reflection formula, and the shift formula.
 # Many details of these techniques are discussed in D. E. G. Hare,

test/gamma.jl

@@ -80,6 +80,17 @@
         @test_throws MethodError logfactorial(1.0)
     end
 
+    @testset "BigFloat" begin
+        loggamma_mpfr(x::BigFloat) = invoke(SpecialFunctions._loggamma, Tuple{BigFloat}, x)
+        logabsgamma_mpfr(x::BigFloat) = invoke(SpecialFunctions._logabsgamma, Tuple{BigFloat}, x)
+        @test loggamma(big(3124)) == loggamma_mpfr(big(3124.0))
+        @test loggamma(big(3124)) ≈ loggamma(3124)
+        @test loggamma(big(1e6)) == loggamma_mpfr(big(1e6))
+        @test logabsgamma(big(1e7)) == logabsgamma_mpfr(big(1e7))
+        @test logabsgamma(big(1e8)) == (loggamma(big(1e8)), 1)
+        @test logabsgamma(big(1e8))[1] ≈ logabsgamma(1e8)[1]
+    end
+
     # loggamma & logabsgamma test cases (from Wolfram Alpha)
     @test loggamma(-300im) ≅ -473.17185074259241355733179182866544204963885920016823743 - 1410.3490664555822107569308046418321236643870840962522425im
     @test loggamma(3.099) ≅ loggamma(3.099+0im) ≅ 0.786413746900558058720665860178923603134125854451168869796