From f9f3ca5504aa1afd494ca8f6e07b35c479e799db Mon Sep 17 00:00:00 2001 From: David Widmann Date: Wed, 2 Jun 2021 15:01:27 +0200 Subject: [PATCH 1/2] Fix `logabsbeta` for BigFloat --- src/beta_inc.jl | 3 ++- test/beta_inc.jl | 1 + test/gamma.jl | 18 ++++++++++++++++++ 3 files changed, 21 insertions(+), 1 deletion(-) diff --git a/src/beta_inc.jl b/src/beta_inc.jl index 6c6b87ff..89682b0e 100644 --- a/src/beta_inc.jl +++ b/src/beta_inc.jl @@ -9,8 +9,9 @@ const exparg_p = log(prevfloat(floatmax(Float64))) Computes ``log(\\Gamma(b)/\\Gamma(a+b))`` when b >= 8 """ -loggammadiv(a::Number, b::Number) = _loggammadiv(float(a), float(b)) +loggammadiv(a::Number, b::Number) = _loggammadiv(promote(float(a), float(b))...) +_loggammadiv(a::Number, b::Number) = loggamma(b) - loggamma(a + b) _loggammadiv(a::T, b::T) where {T<:Base.IEEEFloat} = T(_loggammadiv(Float64(a), Float64(b))) # handle Float16, Float32 function _loggammadiv(a::Float64, b::Float64) @assert b >= 8 diff --git a/test/beta_inc.jl b/test/beta_inc.jl index 591b11a5..30685099 100644 --- a/test/beta_inc.jl +++ b/test/beta_inc.jl @@ -227,6 +227,7 @@ @test beta_inc(a, b, x, 1 - x) == beta_inc(a, Float64(b), Float64(x), 1 - x) @test SpecialFunctions.loggammadiv(13.89, 21.0001) ≈ log(gamma(big(21.0001))/gamma(big(21.0001)+big(13.89))) + @test SpecialFunctions.loggammadiv(big(13.89), big(21.0001)) ≈ log(gamma(big(21.0001))/gamma(big(21.0001)+big(13.89))) @test SpecialFunctions.stirling_corr(11.99, 100.1) ≈ SpecialFunctions.stirling_error(11.99) + SpecialFunctions.stirling_error(100.1) - SpecialFunctions.stirling_error(11.99 + 100.1) end diff --git a/test/gamma.jl b/test/gamma.jl index 790b6718..26f287f6 100644 --- a/test/gamma.jl +++ b/test/gamma.jl @@ -274,6 +274,24 @@ end @test logabsbeta(1e8, 0.5)[1] ≈ log(0.00017724538531210809) rtol=1e-14 end + @testset "BigFloat" begin + @test beta(big(3), big(5)) ≈ inv(big(105)) + @test beta(big(3//2), big(7//2)) ≈ 5 * big(π) / 128 + @test beta(big(1e4), big(3//2)) ≈ 8.86193693673874630607029e-7 rtol=1e-14 + @test beta(big(1e8), big(1//2)) ≈ 0.00017724538531210809 rtol=1e-14 + + @test logbeta(big(5), big(4)) ≈ log(beta(big(5), big(4))) + @test logbeta(big(5.0), big(4)) ≈ log(beta(big(5), big(4))) + @test logbeta(big(1e4), big(3//2)) ≈ log(8.86193693673874630607029e-7) rtol=1e-14 + @test logbeta(big(1e8), big(1//2)) ≈ log(0.00017724538531210809) rtol=1e-14 + + @test logabsbeta(-big(2), big(2))[1] ≈ -log(big(2)) + @test logabsbeta(-big(2.0), big(2))[1] ≈ -log(big(2)) + @test logabsbeta(-big(2), big(2//1))[1] ≈ -log(big(2)) + @test logabsbeta(big(1e4), big(3//2))[1] ≈ log(8.86193693673874630607029e-7) rtol=1e-14 + @test logabsbeta(big(1e8), big(1//2))[1] ≈ log(0.00017724538531210809) rtol=1e-14 + end + @test beta(-1/2, 3) ≈ beta(-1/2 + 0im, 3 + 0im) ≈ -16/3 @test logabsbeta(-1/2, 3)[1] ≈ log(16/3) @test beta(Float32(5), Float32(4)) == beta(Float32(4), Float32(5)) From 9e7461fb9551ead20e05e03de4b5a11133ae78b3 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Fri, 4 Jun 2021 11:09:01 +0200 Subject: [PATCH 2/2] Add a comment --- src/beta_inc.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/beta_inc.jl b/src/beta_inc.jl index 89682b0e..62b6e6ba 100644 --- a/src/beta_inc.jl +++ b/src/beta_inc.jl @@ -11,7 +11,8 @@ Computes ``log(\\Gamma(b)/\\Gamma(a+b))`` when b >= 8 """ loggammadiv(a::Number, b::Number) = _loggammadiv(promote(float(a), float(b))...) -_loggammadiv(a::Number, b::Number) = loggamma(b) - loggamma(a + b) +# TODO: Add a proper implementation +_loggammadiv(a::Number, b::Number) = loggamma(b) - loggamma(a + b) # handle e.g. BigFloat (used by `logabsbeta`) _loggammadiv(a::T, b::T) where {T<:Base.IEEEFloat} = T(_loggammadiv(Float64(a), Float64(b))) # handle Float16, Float32 function _loggammadiv(a::Float64, b::Float64) @assert b >= 8