Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add lbinomial function based on lbeta #27419

Closed
wants to merge 1 commit into from
Closed

Add lbinomial function based on lbeta #27419

wants to merge 1 commit into from

Conversation

freeboson
Copy link
Contributor

@freeboson freeboson commented Jun 4, 2018 •
edited
Loading

This implementation of lbinomial uses the fact that binomial(n,k) = 1/( (n+1) * B(n - k + 1, k + 1) ). I have added tests for some simple cases, and then also comparing for very large n with binomial(::BigInt, ::BigInt).

The main motivation is evaluating log(binomial(n,k)) for large n with k near n/2 without approximating with Poisson:

julia> log(binomial(70, 25)) # OK
43.311504703669215

julia> log(binomial(70, 26)) # :(
ERROR: InexactError()
Stacktrace:
 [1] binomial(::Int64, ::Int64) at ./intfuncs.jl:810

julia> log(binomial(BigInt(70), BigInt(26))) # OK, if you can afford it
4.386007065541805521142055198415251198714576775237582336220765505130378175238536e+01

julia> ans - lbinomial(70, 26) # Close enough?
-8.056495743460440321531836165841374176637792344948696218247614638531207882439305e-15

In R, this function is available as lchoose(n, k) and is a builtin.

vtjnash
vtjnash reviewed Jun 5, 2018
View reviewed changes
base/special/gamma.jl Outdated
Natural logarithm of the [`binomial`](@ref) coefficient.
"""
function lbinomial(n::T, k::T) where {T<:Integer}
const S = float(T)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

┌ Warning: Deprecated syntax ``const declaration on local variable around special/gamma.jl:165.

base/special/gamma.jl Outdated
"""
function lbinomial(n::T, k::T) where {T<:Integer}
const S = float(T)
(k < 0) && return typemin(S)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

extra space between < 0

base/special/gamma.jl Outdated
(k == 0 || k == n) && return zero(S)
(k == 1) && return log(abs(n))
if k > (n>>1)
k = n -k
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

space after - (n - k)

base/special/gamma.jl Outdated

Natural logarithm of the [`binomial`](@ref) coefficient.
"""
function lbinomial(n::T, k::T) where {T<:Integer}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also define:

lbinomial(n::Integer, k::Integer) = lbinomial(promote(n, k)...)

?

vtjnash
vtjnash reviewed Jun 5, 2018
View reviewed changes
test/math.jl Outdated
@test lbinomial(10, 0) == 0.0
@test lbinomial(10, 10) == 0.0
@test lbinomial(10, 1) ≈ log(10)
@test lbfixed(200).(0:200) ≈ blbfixed(200).(0:200)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

maybe add a test for the behavior when n < 0?

vtjnash
vtjnash reviewed Jun 5, 2018
View reviewed changes
base/special/gamma.jl
end
k > n && return typemin(S)
(k == 0 || k == n) && return zero(S)
(k == 1) && return log(abs(n))
Copy link
Member

@vtjnash vtjnash Jun 5, 2018 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In the code for binomial, it looks like this is log(n < 0 && isodd(k) ? -n : n)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually, I guess this is equivalent, it just needs isodd(k) && throw(DomainError()) in the n < 0 branch?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well, it's typical for lbinomial to actually mean log(abs(binomial)). I didn't note that in the docstr however.

@freeboson
Copy link
Contributor Author

I've rebased with:

  • fixed ws
  • tests for n<0
  • docstring notes lbinomial(n, k) = log(abs(binomial(n, k)))
  • method that promotes different n and k of different Integer types
  • added lbinomial to doc/src/base/math.md

@freeboson
Copy link
Contributor Author

Any interest in lbinomial(::Real, ::Integer), as in R?

@dpsanders
Copy link
Contributor

Maybe the Combinatorics.jl package would be a good place for this?

@StefanKarpinski
Copy link
Member

StefanKarpinski commented Jun 6, 2018 •
edited
Loading

This is a great function to have somewhere but we are generally trying to slim down Base to a minimum. I also think that lgamma, gamma, lfact, beta and lbeta should all go elsewhere as well, see #27459. Perhaps lbinomial should also live there since it's defined in terms of the lbeta function rather than being a combinatorial definition?

@freeboson
Copy link
Contributor Author

I didn't know where lbinomial should go. I saw lgamma/lfact and lbeta all in one place, so it made sense to shove it there. If you are planning on moving those other functions, it makes sense to put lbinomial with them.

@freeboson
Copy link
Contributor Author

Since #27473 has been merged, I'll just close this and open in SpecialFunctions.jl.

@freeboson freeboson closed this Jun 19, 2018
@freeboson freeboson mentioned this pull request Jun 19, 2018
Merged
@freeboson freeboson deleted the lbinomial branch June 19, 2018 03:03
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Reviewers

@vtjnash vtjnash vtjnash left review comments

Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

Successfully merging this pull request may close these issues.
4 participants
@freeboson @dpsanders @StefanKarpinski @vtjnash