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

extend binomial to real and complex numbers #14165

Closed
wants to merge 2 commits into from
Closed

extend binomial to real and complex numbers #14165

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ef1e530
generalize binomial to real and complex arguments
CarloLucibello Nov 26, 2015
360e4d1
Merge branch 'master' of https://github.com/CarloLucibello/julia
CarloLucibello Nov 27, 2015
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions base/special/gamma.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -482,3 +482,7 @@ eta(x::Integer) = eta(Float64(x))
eta(x::Real) = oftype(float(x),eta(Float64(x)))
eta(z::Complex) = oftype(float(z),eta(Complex128(z)))
@vectorize_1arg Number eta

Base.binomial{S<:Number,T<:Number}(n::S, k::T) = ( (isinteger(n) && isinteger(k))
Copy link
Contributor

Choose a reason for hiding this comment

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

this would be better off using dispatch rather than runtime type checks

Copy link
Contributor Author

Choose a reason for hiding this comment

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

here with isinteger I don't check on type but only if I can cast to some integer type without loss. This is to be able to compute binomial(2.0, 1)

Copy link
Contributor

Choose a reason for hiding this comment

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

ah, right. the generic isinteger for floating point numbers is comparing whether the value is equal to trunc of itself and is finite, which doesn't mean it'll necessarily fit in an Int64 - for type stability I wonder whether it would be any simpler to just always use the definition in terms of gamma

Copy link
Contributor Author

Choose a reason for hiding this comment

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

it would be easier but we also have to address the case where we have infinities both in the numerator and in the denominator, that result in a finite number for the fraction. This can happen only for negative integer arguments of the gamma functions, and it is handed nicely by the integer version of binomial

Copy link
Contributor

Choose a reason for hiding this comment

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

Would be ideal to have more test cases for all these corner cases. What should be done for values too large for Int64 ?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

What about converting always to BigInt instead of Int64?

? convert(promote_type(T,S), binomial(round(BigInt, n), round(BigInt, k)))

it would be a little slower but will address the corner case of those crazy enough to evaluate binomials with such big numbers

Copy link
Member

Choose a reason for hiding this comment

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

The general practice in Julia is to use native Int by default, to avoid imposing an overhead on all use cases just to support very specific cases. Always using BigInt goes against this idea.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I agree, so better stay with Int64 conversion. In any case one could always resort to the BigInt version of binomial if needs to

Copy link
Contributor

Choose a reason for hiding this comment

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

If you know when it's appropriate to overflow I think it's better to do so rather than throw an InexactError in round(Int64, n). This also needs quite a few more tests.

? convert(promote_type(T,S), binomial(round(Int64, n), round(Int64, k)))
: gamma(1+n) / (gamma(1+n-k) * gamma(1+k)) )
Copy link
Member

Choose a reason for hiding this comment

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

This seems like it is susceptible to spurious overflow. It might be more reliable to do exp(lgamma(1+n) - lgamma(1+n-k) - lgamma(1+k))?

I wouldn't bother with the integer check and binomial calls, which are even more susceptible to overflow. If the user has integer values and wants specialized code, they should convert to integers themselves.

For example, try n=200.0 and k=100.0. The correct answer is representable, 9.054851465609245e58, but both your integer calls and your gamma calls will overflow.

Copy link
Member

Choose a reason for hiding this comment

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

This is how we compute the Beta function so I think you should just call beta here.

Copy link
Member

Choose a reason for hiding this comment

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

@andreasnoack, the beta function is similar but does not seem identical. If x = 1 + n - k and w = 1 + k, then x + w = 2 + n and hence the function here does not seem to match 1/beta(x,w).

Copy link
Member

Choose a reason for hiding this comment

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

I think binomial(x,y) = inv((x+1)*beta(x - y + 1, y + 1)) should do the job.

Copy link
Member

Choose a reason for hiding this comment

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

True. Probably the cost of the inv is negligible, and maybe someday we will have a more optimized beta function.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I think binomial(x,y) = inv((x+1)*beta(x - y + 1, y + 1)) should do the job

This looks like a nice solution.

I wouldn't bother with the integer check and binomial calls, which are even more susceptible to overflow. If the user has integer values and wants specialized code, they should convert to integers themselves.

I would like to cover cases such as

julia> f(x,y) = inv((x+1)*beta(x - y + 1, y + 1))
f (generic function with 1 method)

julia> binomial(-2.,4.)
5.0

julia> f(-2.,4.)
NaN

I could just do the residues calculation without casting to integers....

Copy link
Member

Choose a reason for hiding this comment

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

@CarloLucibello, I think that is a bug in our beta function. beta(-n,n) should return -1/n, not NaN. We should fix the problem there.

Copy link
Member

Choose a reason for hiding this comment

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

See #14256.

Copy link
Member

Choose a reason for hiding this comment

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

In fact, the kinds of subtleties handled by the Cephes implementation I linked from 14256 make it even clearer that the binomial function should call beta for non-Integer arguments — we don't want to implement those subtleties twice.

5 changes: 5 additions & 0 deletions test/math.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -712,3 +712,8 @@ end
let A = [1 2; 3 4]; B = [5 6; 7 8]; C = [9 10; 11 12]
@test muladd(A,B,C) == A*B + C
end

@test binomial(2.,1) == 2
@test binomial(-2.,1) == -2
@test typeof(binomial(-2.,1)) == Float64
@test_approx_eq binomial(-2.1,1.1) -4.21013511374008