From 8136181a2f7b44d90f8f486e27c013cb747ceb32 Mon Sep 17 00:00:00 2001 From: Simon Byrne Date: Tue, 7 May 2019 10:20:46 -0700 Subject: [PATCH] rename lfactorial/lbinomial (#161) Ensure consistency with #156. --- Project.toml | 8 +++--- docs/src/index.md | 3 ++- docs/src/special.md | 3 ++- src/SpecialFunctions.jl | 2 +- src/deprecated.jl | 5 ++++ src/gamma.jl | 55 ++++++++++++++++++++++++++++------------- test/runtests.jl | 35 ++++++++++++++------------ 7 files changed, 71 insertions(+), 40 deletions(-) diff --git a/Project.toml b/Project.toml index ccda5a0a..a72b522a 100644 --- a/Project.toml +++ b/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" diff --git a/docs/src/index.md b/docs/src/index.md index 476f9e8f..e0f2320f 100644 --- a/docs/src/index.md +++ b/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 diff --git a/docs/src/special.md b/docs/src/special.md index b54b620d..342a4297 100644 --- a/docs/src/special.md +++ b/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 ``` diff --git a/src/SpecialFunctions.jl b/src/SpecialFunctions.jl index 9dd9200d..170c0791 100644 --- a/src/SpecialFunctions.jl +++ b/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) diff --git a/src/deprecated.jl b/src/deprecated.jl index 075074ff..a1a79c12 100644 --- a/src/deprecated.jl +++ b/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) diff --git a/src/gamma.jl b/src/gamma.jl index 3d9c7a2e..f61dc439 100644 --- a/src/gamma.jl +++ b/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)...) diff --git a/test/runtests.jl b/test/runtests.jl index 8a92a810..4562c997 100644 --- a/test/runtests.jl +++ b/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