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.

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

Handle arguments of 1 in beta and logabsbeta #349

Open
wants to merge 5 commits into
base: master
Choose a base branch
from

Conversation

devmotion
Copy link
Member

@devmotion devmotion commented Sep 28, 2021

This PR adds special handling for arguments of 1 in beta and logabsbeta. This fixes one of the issues in JuliaStats/StatsFuns.jl#126 (I checked locally that the problem is fixed) but probably it is useful more generally (see JuliaStats/StatsFuns.jl#126 (comment) for a problem in the linked PR).

Comparison (copied from JuliaStats/StatsFuns.jl#126 (comment) and JuliaStats/StatsFuns.jl#126 (comment)):

On the master branch we have

julia> beta(1.0, 200.0)
0.005000000000000002

julia> logabsbeta(1.0, 4.0)
(-1.3862943611198908, 1)

julia> -log(4.0)
-1.3862943611198906

julia> @btime beta(1.0, 200.0);
  47.907 ns (0 allocations: 0 bytes)

julia> @btime logabsbeta(1.0, 4.0);
  60.584 ns (0 allocations: 0 bytes)

julia> @btime beta(3.4, 200.0);
  66.080 ns (0 allocations: 0 bytes)

julia> @btime logabsbeta(3.4, 4.0);
  78.036 ns (0 allocations: 0 bytes)

whereas with this PR one gets

julia> beta(1.0, 200.0)
0.005

julia> logabsbeta(1.0, 4.0)
(-1.3862943611198906, 1)

julia> -log(4.0)
-1.3862943611198906

julia> @btime beta(1.0, 200.0); # seems it is compiled away completely...
  0.019 ns (0 allocations: 0 bytes)

julia> @btime logabsbeta(1.0, 4.0);
  1.451 ns (0 allocations: 0 bytes)

julia> @btime beta(3.4, 200.0);
  67.098 ns (0 allocations: 0 bytes)

julia> @btime logabsbeta(3.4, 4.0);
  79.259 ns (0 allocations: 0 bytes)

The additional dispatches of beta that are introduced to ensure type stability also fix a bug where e.g. beta(::BigInt, ::BigFloat) did not use and differed from the MPFR implementation.

src/gamma.jl Outdated
Comment on lines 745 to 747
# here `T` is a floating point type but we don't want to restrict the implementation
# to `AbstractFloat` or `Float64`
function _beta(a::T, b::T) where {T<:Number}
Copy link
Contributor

Choose a reason for hiding this comment

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

How do you know that it's a floating point type then? :-) It could also be a complex.

Copy link
Member Author

Choose a reason for hiding this comment

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

I guess I was not completely precise here. What I meant was that the types are the result of calling float (this is what I called slightly incorrectly "floating point type"), so also in the case of complex numbers it has to be Complex{T} where T is a floating point number.

@devmotion
Copy link
Member Author

It seems there is some problem with the MPFR library on Windows: beta(big(1.0e8), big(1 // 2)) returns NaN only on Windows. The problem did only show up now since previously for these arguments the Julia implementation was used.

@devmotion
Copy link
Member Author

The documentation of mpfr_beta says:

Note: the current code does not try to avoid internal overflow or underflow, and might use a huge internal precision in some cases.

Maybe on Windows the default precision can't handle large arguments such as big(1e8) yet? Maybe we should not use MPFR on Windows?

@ViralBShah
Copy link
Member

This is an old one, but good to merge if we can - @oscardssmith?

Copy link

codecov bot commented May 7, 2024

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 92.96%. Comparing base (124915f) to head (7d4a576).
Report is 1 commits behind head on master.

❗ Current head 7d4a576 differs from pull request most recent head c655873. Consider uploading reports for the commit c655873 to get more accurate results

Additional details and impacted files
@@            Coverage Diff             @@
##           master     #349      +/-   ##
==========================================
- Coverage   94.10%   92.96%   -1.15%     
==========================================
  Files          14       12       -2     
  Lines        2935     2715     -220     
==========================================
- Hits         2762     2524     -238     
- Misses        173      191      +18     
Flag Coverage Δ
unittests 92.96% <100.00%> (-1.15%) ⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants