extend binomial to real and complex numbers#14165
extend binomial to real and complex numbers#14165CarloLucibello wants to merge 2 commits intoJuliaLang:masterfrom
Conversation
There was a problem hiding this comment.
this would be better off using dispatch rather than runtime type checks
There was a problem hiding this comment.
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)
There was a problem hiding this comment.
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
There was a problem hiding this comment.
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
There was a problem hiding this comment.
Would be ideal to have more test cases for all these corner cases. What should be done for values too large for Int64 ?
There was a problem hiding this comment.
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
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
I agree, so better stay with Int64 conversion. In any case one could always resort to the BigInt version of binomial if needs to
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
This is how we compute the Beta function so I think you should just call beta here.
There was a problem hiding this comment.
@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).
There was a problem hiding this comment.
I think binomial(x,y) = inv((x+1)*beta(x - y + 1, y + 1)) should do the job.
There was a problem hiding this comment.
True. Probably the cost of the inv is negligible, and maybe someday we will have a more optimized beta function.
There was a problem hiding this comment.
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.)
NaNI could just do the residues calculation without casting to integers....
There was a problem hiding this comment.
@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.
There was a problem hiding this comment.
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.
|
@CarloLucibello It would be great to have this PR completed. |
|
@andreasnoack Right, I'll try to complete it during this week |
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). |
No description provided.