add lbinomial (#99)

Author: freeboson

src/SpecialFunctions.jl

@@ -68,7 +68,8 @@ else
 end
 
 export sinint,
-       cosint
+       cosint,
+       lbinomial
 
 include("bessel.jl")
 include("erf.jl")

src/gamma.jl

@@ -769,3 +769,25 @@ const lfactorial = lfact
 export lfactorial
 
 end # @static if
+
+"""
+    lbinomial(n, k) = log(abs(binomial(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`.
+"""
+function lbinomial(n::T, k::T) where {T<:Integer}
+    S = float(T)
+    (k < 0) && return typemin(S)
+    if n < 0
+        n = -n + k - 1
+    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) - lbeta(n - k + one(T), k + one(T))
+end
+lbinomial(n::Integer, k::Integer) = lbinomial(promote(n, k)...)

test/runtests.jl

@@ -642,6 +642,18 @@ 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)))
+end
+
 @static if isdefined(Base, :missing)
     @testset "missing data" begin
         for f in (digamma, erf, erfc, erfcinv, erfcx, erfi, erfinv, eta, gamma,