extend binomial to real and complex numbers

[CarloLucibello]

No description provided.

[tkelman]

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

[CarloLucibello]

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)

[tkelman]

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

[CarloLucibello]

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

[tkelman]

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

[CarloLucibello]

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

[nalimilan]

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.

[CarloLucibello]

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

[tkelman]

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.

[stevengj]

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.

[andreasnoack]

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

[stevengj]

@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).

[andreasnoack]

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

[stevengj]

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

[CarloLucibello]

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....

[stevengj]

@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.

[stevengj]

See #14256.

[stevengj]

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.

[andreasnoack]

@CarloLucibello It would be great to have this PR completed.

[CarloLucibello]

@andreasnoack Right, I'll try to complete it during this week

[musm]

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.

Since the beta special function has been moved to SpecialFunctions, addressing the issue here requires the method to be added in that package (moved to JuliaMath/SpecialFunctions.jl#282).