rename lfactorial/lbinomial (#161)

Ensure consistency with #156.

Author: simonbyrne

Project.toml

@@ -7,12 +7,12 @@ BinDeps = "9e28174c-4ba2-5203-b857-d8d62c4213ee"
 BinaryProvider = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
 Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
 
+[compat]
+BinaryProvider = "≥ 0.5.0"
+julia = "≥ 0.7.0"
+
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
 test = ["Test"]
-
-[compat]
-BinaryProvider = "≥ 0.5.0"
-julia = "≥ 0.7.0"

docs/src/index.md

@@ -49,10 +49,11 @@ libraries.
 | [`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                                                                   |
+| [`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`     |
+| [`logabsbinomial(x,y)`](@ref SpecialFunctions.logabsbinomial)                                          | 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

@@ -51,8 +51,9 @@ SpecialFunctions.zeta
 SpecialFunctions.gamma
 SpecialFunctions.loggamma
 SpecialFunctions.logabsgamma
-SpecialFunctions.lfactorial
+SpecialFunctions.logfactorial
 SpecialFunctions.beta
 SpecialFunctions.logbeta
 SpecialFunctions.logabsbeta
+SpecialFunctions.logabsbinomial
 ```

src/SpecialFunctions.jl

@@ -65,7 +65,7 @@ include("gamma.jl")
 include("deprecated.jl")
 
 for f in (:digamma, :erf, :erfc, :erfcinv, :erfcx, :erfi, :erfinv,
-          :eta, :gamma, :invdigamma, :lfactorial, :lgamma, :trigamma, :ellipk, :ellipe)
+          :eta, :gamma, :invdigamma, :logfactorial, :lgamma, :trigamma, :ellipk, :ellipe)
     @eval $(f)(::Missing) = missing
 end
 for f in (:beta, :lbeta)

src/deprecated.jl

@@ -9,6 +9,11 @@
 @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)
+@deprecate lfact(x) logfactorial(x)
+@deprecate lfactorial(x) logfactorial(x)
+@deprecate lbinomial(x) logabsbinomial(x)[1]
+
+
 
 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

@@ -551,7 +551,7 @@ function polygamma(m::Integer, z::Number)
     polygamma(m, x)
 end
 
-export gamma, loggamma, logabsgamma, beta, logbeta, logabsbeta, lfactorial
+export gamma, loggamma, logabsgamma, beta, logbeta, logabsbeta, logfactorial, logabsbinomial
 
 ## from base/special/gamma.jl
 
@@ -584,27 +584,31 @@ logabsgamma(x::AbstractFloat) = throw(MethodError(logabsgamma, x))
 
 
 """
-    lfactorial(x)
+    logfactorial(x)
 
 Compute the logarithmic factorial of a nonnegative integer `x`.
 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.")) : loggamma(x + oneunit(x))
-Base.@deprecate lfact lfactorial
+logfactorial(x::Integer) = x < 0 ? throw(DomainError(x, "`x` must be non-negative.")) : loggamma(x + oneunit(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))).
+[`Real`](@ref) `x`and returns a tuple `(log(abs(gamma(x))), sign(gamma(x)))`.
+
+See also [`loggamma`](@ref).
 """
 function logabsgamma end
 
 """
     loggamma(x)
 
-Computes the logarithm of [`gamma`](@ref) for given `x`
+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).
 """
 function loggamma end
 
@@ -732,15 +736,18 @@ end
 """
     logbeta(x, y)
 
-Natural logarithm of the [`beta`](@ref)
-function ``\\log(|\\operatorname{B}(x,y)|)``.
+Natural logarithm of the [`beta`](@ref) function ``\\log(|\\operatorname{B}(x,y)|)``.
+
+See also [`logabsbeta`](@ref).
 """
 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 .
+Compute the natural logarithm of the absolute value of the [`beta`](@ref) function, returning a tuple `(log(abs(beta(x,y))), sign(beta(x,y)))`
+
+See also [`logbeta`](@ref).
 """
 function logabsbeta(x::Real, w::Real)
     yx, sx = logabsgamma(x)
@@ -802,23 +809,37 @@ factorial(x) = Base.factorial(x) # to make SpecialFunctions.factorial work uncon
 factorial(x::Number) = gamma(x + 1) # fallback for x not Integer
 
 """
-    lbinomial(n, k) = log(abs(binomial(n, k)))
+    logabsbinomial(n, k)
 
 Accurate natural logarithm of the absolute value of the [`binomial`](@ref)
 coefficient `binomial(n, k)` for large `n` and `k` near `n/2`.
+
+Returns a tuple `(log(abs(binomial(n,k))), sign(binomial(n,k)))`.
 """
-function lbinomial(n::T, k::T) where {T<:Integer}
+function logabsbinomial(n::T, k::T) where {T<:Integer}
     S = float(T)
-    (k < 0) && return typemin(S)
+    s = one(S)
     if n < 0
+        # reflection formula
+        #  binomial(n,k) = (-1)^k * binomial(-n+k-1,k)
         n = -n + k - 1
+        s = isodd(k) ? -s : s
+    end
+    if k < 0 || k > n
+        # binomial(n,k) == 0
+        return (typemin(S), zero(S))
+    elseif k == 0 || k == n
+        # binomial(n,k) == ±1
+        return (zero(S), s)
     end
-    k > n && return typemin(S)
-    (k == 0 || k == n) && return zero(S)
-    (k == 1) && return log(abs(n))
     if k > (n>>1)
         k = n - k
     end
-    -log1p(n) - logabsbeta(n - k + one(T), k + one(T))[1]
+    if k == 1
+        # binomial(n,k) == ±n
+        return (log(abs(n)), s)
+    else
+        return (-log1p(n) - logabsbeta(n - k + one(T), k + one(T))[1], s)
+    end
 end
-lbinomial(n::Integer, k::Integer) = lbinomial(promote(n, k)...)
+logabsbinomial(n::Integer, k::Integer) = logabsbinomial(promote(n, k)...)

test/runtests.jl

@@ -634,12 +634,12 @@ end
         for x in (3.2, 2+1im, 3//2, 3.2+0.1im)
             @test SpecialFunctions.factorial(x) == gamma(1+x)
         end
-        @test lfactorial(0) == lfactorial(1) == 0
-        @test lfactorial(2) == loggamma(3)
-        # Ensure that the domain of lfactorial matches that of factorial (issue #21318)
-        @test_throws DomainError lfactorial(-3)
+        @test logfactorial(0) == logfactorial(1) == 0
+        @test logfactorial(2) == loggamma(3)
+        # Ensure that the domain of logfactorial matches that of factorial (issue #21318)
+        @test_throws DomainError logfactorial(-3)
         @test_throws DomainError loggamma(-4.2)
-        @test_throws MethodError lfactorial(1.0)
+        @test_throws MethodError logfactorial(1.0)
     end
 
     # loggamma & logabsgamma test cases (from Wolfram Alpha)
@@ -708,21 +708,24 @@ end
     @test beta(big(1.0),big(1.2)) ≈ beta(1.0,1.2) rtol=4*eps()
 end
 
-@testset "lbinomial" begin
-    @test lbinomial(10, -1) == -Inf
-    @test lbinomial(10, 11) == -Inf
-    @test lbinomial(10,  0) == 0.0
-    @test lbinomial(10, 10) == 0.0
-
-    @test lbinomial(10,  1)      ≈ log(10)
-    @test lbinomial(-6, 10)      ≈ log(binomial(-6, 10))
-    @test lbinomial(-6, 11)      ≈ log(abs(binomial(-6, 11)))
-    @test lbinomial.(200, 0:200) ≈ log.(binomial.(BigInt(200), (0:200)))
+@testset "logabsbinomial" begin
+    @test logabsbinomial(10, -1) == (-Inf, 0.0)
+    @test logabsbinomial(10, 11) == (-Inf, 0.0)
+    @test logabsbinomial(10,  0) == ( 0.0, 1.0)
+    @test logabsbinomial(10, 10) == ( 0.0, 1.0)
+
+    @test logabsbinomial(10,  1)[1]   ≈ log(10)
+    @test logabsbinomial(10,  1)[2]   == 1.0
+    @test logabsbinomial(-6, 10)[1]   ≈ log(binomial(-6, 10))
+    @test logabsbinomial(-6, 10)[2]   == 1.0
+    @test logabsbinomial(-6, 11)[1]   ≈ log(abs(binomial(-6, 11)))
+    @test logabsbinomial(-6, 11)[2]   == -1.0
+    @test first.(logabsbinomial.(200, 0:200)) ≈ log.(binomial.(BigInt(200), (0:200)))
 end
 
 @testset "missing data" begin
     for f in (digamma, erf, erfc, erfcinv, erfcx, erfi, erfinv, eta, gamma,
-              invdigamma, lfactorial, trigamma)
+              invdigamma, logfactorial, trigamma)
         @test f(missing) === missing
     end
     @test beta(1.0, missing) === missing