Rename lgamma/lbeta to log(abs)gamma/log(abs)beta (#156)

See behaviour of `logdet`/`logabsdet` in the LinearAlgebra stdlib

docs/src/index.md

@@ -46,11 +46,13 @@ libraries.
 | [`besselk(nu,z)`](@ref SpecialFunctions.besselk)              | modified [Bessel function](https://en.wikipedia.org/wiki/Bessel_function) of the second kind of order `nu` at `z`                                               |
 | [`besselkx(nu,z)`](@ref SpecialFunctions.besselkx)            | scaled modified Bessel function of the second kind of order `nu` at `z`                                                                                         |
 | [`gamma(x)`](@ref SpecialFunctions.gamma)                                            | [gamma function](https://en.wikipedia.org/wiki/Gamma_function) at `x`                                                                                           |
+| [`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`                                                                                                                          |
 | [`lfactorial(x)`](@ref SpecialFunctions.lfactorial)                                            | 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`                                                                                           |
-| [`lbeta(x,y)`](@ref SpecialFunctions.lbeta)                                          | accurate `log(beta(x,y))` for large `x` or `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`     |
 ## Installation
 
 The package is available for Julia versions 0.5 and up. To install it, run

docs/src/special.md

@@ -49,8 +49,10 @@ SpecialFunctions.ellipe
 SpecialFunctions.eta
 SpecialFunctions.zeta
 SpecialFunctions.gamma
-SpecialFunctions.lgamma
+SpecialFunctions.loggamma
+SpecialFunctions.logabsgamma
 SpecialFunctions.lfactorial
 SpecialFunctions.beta
-SpecialFunctions.lbeta
+SpecialFunctions.logbeta
+SpecialFunctions.logabsbeta
 ```

src/deprecated.jl

@@ -3,6 +3,12 @@
 @deprecate airy(z::Number) airyai(z)
 @deprecate airyx(z::Number) airyaix(z)
 @deprecate airyprime(z::Number) airyaiprime(z)
+@deprecate lgamma(x::Real) logabsgamma(x)[1]
+@deprecate lgamma(x::Number) loggamma(x)
+@deprecate lgamma_r(x::Real) logabsgamma(x)
+@deprecate lgamma_r(x::Number) (loggamma(x), 1)
+@deprecate lbeta(x::Real, w::Real) logabsbeta(x, w)[1]
+@deprecate lbeta(x::Number, w::Number) logbeta(x, w)
 
 function _airy(k::Integer, z::Complex{Float64})
     Base.depwarn("`airy(k,x)` is deprecated, use `airyai(x)`, `airyaiprime(x)`, `airybi(x)` or `airybiprime(x)` instead.",:airy)

src/gamma.jl

@@ -412,7 +412,7 @@ function zeta(s::ComplexOrReal{Float64})
         end
         if absim > 12 # amplitude of sinpi(s/2) ≈ exp(imag(s)*π/2)
             # avoid overflow/underflow (issue #128)
-            lg = lgamma(1 - s)
+            lg = loggamma(1 - s)
             ln2pi = 1.83787706640934548356 # log(2pi) to double precision
             rehalf = real(s)*0.5
             return zeta(1 - s) * exp(lg + absim*(pi/2) + s*ln2pi) * (0.5/π) *
@@ -551,7 +551,7 @@ function polygamma(m::Integer, z::Number)
     polygamma(m, x)
 end
 
-export gamma, lgamma, beta, lbeta, lfactorial
+export gamma, loggamma, logabsgamma, beta, logbeta, logabsbeta, lfactorial
 
 ## from base/special/gamma.jl
 
@@ -567,47 +567,56 @@ Compute the gamma function of `x`.
 """
 gamma(x::Real) = gamma(float(x))
 
-function lgamma_r(x::Float64)
+function logabsgamma(x::Float64)
     signp = Ref{Int32}()
     y = ccall((:lgamma_r,libm),  Float64, (Float64, Ptr{Int32}), x, signp)
     return y, signp[]
 end
-function lgamma_r(x::Float32)
+function logabsgamma(x::Float32)
     signp = Ref{Int32}()
     y = ccall((:lgammaf_r,libm),  Float32, (Float32, Ptr{Int32}), x, signp)
     return y, signp[]
 end
-lgamma_r(x::Real) = lgamma_r(float(x))
-lgamma_r(x::Number) = lgamma(x), 1 # lgamma does not take abs for non-real x
-"""
-    lgamma_r(x)
+logabsgamma(x::Real) = logabsgamma(float(x))
+logabsgamma(x::Float16) = Float16.(logabsgamma(Float32(x)))
+logabsgamma(x::Integer) = logabsgamma(float(x))
+logabsgamma(x::AbstractFloat) = throw(MethodError(logabsgamma, x))
 
-Return L,s such that `gamma(x) = s * exp(L)`
-"""
-lgamma_r
 
 """
     lfactorial(x)
 
 Compute the logarithmic factorial of a nonnegative integer `x`.
-Equivalent to [`lgamma`](@ref) of `x + 1`, but `lgamma` extends this function
+Equivalent to [`loggamma`](@ref) of `x + 1`, but `loggamma` extends this function
 to non-integer `x`.
 """
-lfactorial(x::Integer) = x < 0 ? throw(DomainError(x, "`x` must be non-negative.")) : lgamma(x + oneunit(x))
+lfactorial(x::Integer) = x < 0 ? throw(DomainError(x, "`x` must be non-negative.")) : loggamma(x + oneunit(x))
 Base.@deprecate lfact lfactorial
 
 """
-    lgamma(x)
+    logabsgamma(x)
+
+Compute the logarithm of absolute value of [`gamma`](@ref) for
+[`Real`](@ref) `x`and returns a tuple (logabsgamma(x), sign(gamma(x))).
+"""
+function logabsgamma end
+
+"""
+    loggamma(x)
 
-Compute the logarithm of the absolute value of [`gamma`](@ref) for
-[`Real`](@ref) `x`, while for [`Complex`](@ref) `x` compute the
-principal branch cut of the logarithm of `gamma(x)` (defined for negative `real(x)`
-by analytic continuation from positive `real(x)`).
+Computes the logarithm of [`gamma`](@ref) for given `x`
 """
-function lgamma end
+function loggamma end
+
+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)|
-@inline function lgamma_asymptotic(z::Complex{Float64})
+function loggamma_asymptotic(z::Complex{Float64})
     zinv = inv(z)
     t = zinv*zinv
     # coefficients are bernoulli[2:n+1] .// (2*(1:n).*(2*(1:n) - 1))
@@ -625,7 +634,7 @@ end
 # and similar techniques are used (in a somewhat different way) by the
 # SciPy loggamma function.  The key identities are also described
 # at http://functions.wolfram.com/GammaBetaErf/LogGamma/
-function lgamma(z::Complex{Float64})
+function loggamma(z::Complex{Float64})
     x = real(z)
     y = imag(z)
     yabs = abs(y)
@@ -638,7 +647,7 @@ function lgamma(z::Complex{Float64})
             return Complex(NaN, NaN)
         end
     elseif x > 7 || yabs > 7 # use the Stirling asymptotic series for sufficiently large x or |y|
-        return lgamma_asymptotic(z)
+        return loggamma_asymptotic(z)
     elseif x < 0.1 # use reflection formula to transform to x > 0
         if x == 0 && y == 0 # return Inf with the correct imaginary part for z == 0
             return Complex(Inf, signbit(x) ? copysign(oftype(x, pi), -y) : -y)
@@ -647,7 +656,7 @@ function lgamma(z::Complex{Float64})
         return Complex(1.1447298858494001741434262, # log(pi)
                        copysign(6.2831853071795864769252842, y) # 2pi
                        * floor(0.5*x+0.25)) -
-               log(sinpi(z)) - lgamma(1-z)
+               log(sinpi(z)) - loggamma(1-z)
     elseif abs(x - 1) + yabs < 0.1
         # taylor series around zero at z=1
         # ... coefficients are [-eulergamma; [(-1)^k * zeta(k)/k for k in 2:15]]
@@ -671,7 +680,7 @@ function lgamma(z::Complex{Float64})
                                -2.2315475845357937976132853e-04,9.9457512781808533714662972e-05,
                                -4.4926236738133141700224489e-05,2.0507212775670691553131246e-05)
     end
-    # use recurrence relation lgamma(z) = lgamma(z+1) - log(z) to shift to x > 7 for asymptotic series
+    # use recurrence relation loggamma(z) = loggamma(z+1) - log(z) to shift to x > 7 for asymptotic series
     shiftprod = Complex(x,yabs)
     x += 1
     sb = false # == signbit(imag(shiftprod)) == signbit(yabs)
@@ -692,33 +701,53 @@ function lgamma(z::Complex{Float64})
     else
         shift = Complex(real(shift), imag(shift) + signflips*6.2831853071795864769252842)
     end
-    return lgamma_asymptotic(Complex(x,y)) - shift
+    return loggamma_asymptotic(Complex(x,y)) - shift
 end
-lgamma(z::Complex{T}) where {T<:Union{Integer,Rational}} = lgamma(float(z))
-lgamma(z::Complex{T}) where {T<:Union{Float32,Float16}} = Complex{T}(lgamma(Complex{Float64}(z)))
+loggamma(z::Complex{T}) where {T<:Union{Integer,Rational}} = loggamma(float(z))
+loggamma(z::Complex{T}) where {T<:Union{Float32,Float16}} = Complex{T}(loggamma(Complex{Float64}(z)))
 
-gamma(z::Complex) = exp(lgamma(z))
+gamma(z::Complex) = exp(loggamma(z))
 
 """
     beta(x, y)
 
 Euler integral of the first kind ``\\operatorname{B}(x,y) = \\Gamma(x)\\Gamma(y)/\\Gamma(x+y)``.
 """
-function beta(x::Number, w::Number)
-    yx, sx = lgamma_r(x)
-    yw, sw = lgamma_r(w)
-    yxw, sxw = lgamma_r(x+w)
+function beta end
+
+function beta(x::Real, w::Real)
+    yx, sx = logabsgamma(x)
+    yw, sw = logabsgamma(w)
+    yxw, sxw = logabsgamma(x+w)
     return exp(yx + yw - yxw) * (sx*sw*sxw)
 end
 
+function beta(x::Number, w::Number)
+    yx = loggamma(x)
+    yw = loggamma(w)
+    yxw = loggamma(x+w)
+    return exp(yx + yw - yxw)
+end
+
 """
-    lbeta(x, y)
+    logbeta(x, y)
 
-Natural logarithm of the absolute value of the [`beta`](@ref)
+Natural logarithm of the [`beta`](@ref)
 function ``\\log(|\\operatorname{B}(x,y)|)``.
 """
-lbeta(x::Number, w::Number) = lgamma(x)+lgamma(w)-lgamma(x+w)
+logbeta(x::Number, w::Number) = loggamma(x)+loggamma(w)-loggamma(x+w)
+
+"""
+    logabsbeta(x, y)
 
+Returns a tuple of the natural logarithm of the absolute value of the [`beta`](@ref) function ``\\log(|\\operatorname{B}(x,y)|)`` and sign of the [`beta`](@ref) function .
+"""
+function logabsbeta(x::Real, w::Real)
+    yx, sx = logabsgamma(x)
+    yw, sw = logabsgamma(w)
+    yxw, sxw = logabsgamma(x+w)
+    (yx + yw - yxw), (sx*sw*sxw)
+end
 ## from base/mpfr.jl
 
 # Functions for which NaN results are converted to DomainError, following Base
@@ -731,15 +760,18 @@ function gamma(x::BigFloat)
 end
 
 # log of absolute value of gamma function
-function lgamma_r(x::BigFloat)
+function logabsgamma(x::BigFloat)
     z = BigFloat()
     lgamma_signp = Ref{Cint}()
     ccall((:mpfr_lgamma,:libmpfr), Cint, (Ref{BigFloat}, Ref{Cint}, Ref{BigFloat}, Int32), z, lgamma_signp, x, ROUNDING_MODE[])
     return z, lgamma_signp[]
 end
 
-lgamma(x::BigFloat) = lgamma_r(x)[1]
-
+function loggamma(x::BigFloat)
+    (y, s) = logabsgamma(x)
+    s < 0.0 && throw(DomainError(x, "`gamma(x)` must be non-negative"))
+    return y
+end
 if Base.MPFR.version() >= v"4.0.0"
     function beta(y::BigFloat, x::BigFloat)
         z = BigFloat()
@@ -760,12 +792,6 @@ end
 
 ## from base/math.jl
 
-@inline lgamma(x::Float64) = nan_dom_err(ccall((:lgamma, libm), Float64, (Float64,), x), x)
-@inline lgamma(x::Float32) = nan_dom_err(ccall((:lgammaf, libm), Float32, (Float32,), x), x)
-@inline lgamma(x::Float16) = Float16(lgamma(Float32(x)))
-@inline lgamma(x::Real) = lgamma(float(x))
-lgamma(x::AbstractFloat) = throw(MethodError(lgamma, x))
-
 ## from base/numbers.jl
 
 # this trickery is needed while the deprecated method in Base exists
@@ -793,6 +819,6 @@ function lbinomial(n::T, k::T) where {T<:Integer}
     if k > (n>>1)
         k = n - k
     end
-    -log1p(n) - lbeta(n - k + one(T), k + one(T))
+    -log1p(n) - logabsbeta(n - k + one(T), k + one(T))[1]
 end
 lbinomial(n::Integer, k::Integer) = lbinomial(promote(n, k)...)