From c075b6389bb4d904781004a94158dcd9bf7a1646 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Tue, 28 Sep 2021 19:08:14 +0200 Subject: [PATCH 1/3] Handle arguments of 1 in `beta` and `logabsbeta` --- src/gamma.jl | 23 +++++++++++++++++++++-- test/gamma.jl | 37 +++++++++++++++++++++++++++++++++++++ 2 files changed, 58 insertions(+), 2 deletions(-) diff --git a/src/gamma.jl b/src/gamma.jl index ab8c4eed..b0d960b8 100644 --- a/src/gamma.jl +++ b/src/gamma.jl @@ -740,13 +740,21 @@ end Euler integral of the first kind ``\\operatorname{B}(x,y) = \\Gamma(x)\\Gamma(y)/\\Gamma(x+y)``. """ -function beta(a::Number, b::Number) +beta(a::Number, b::Number) = _beta(map(float, promote(a, b))...) + +# 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} + # two special cases + a == 1 && return inv(b) + b == 1 && return inv(a) + lab, sign = logabsbeta(a, b) return sign*exp(lab) end if Base.MPFR.version() >= v"4.0.0" - function beta(y::BigFloat, x::BigFloat) + function _beta(y::BigFloat, x::BigFloat) z = BigFloat() ccall((:mpfr_beta, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Int32), z, y, x, ROUNDING_MODE[]) return z @@ -781,6 +789,17 @@ function logabsbeta(a::T, b::T) where T<:Real return logabsbeta(b, a) end + # minimum of arguments is 1 + if a == 1 + # b >= a >= 1, so `abs` is not needed and the sign is always 1 + return -log(b), 1 + end + # maximum of arguments is 1 + if b == 1 + sgn = a >= 0 ? 1 : -1 + return -log(abs(a)), sgn + end + if a <= 0 && isinteger(a) if a + b <= 0 && isinteger(b) r = logbeta(1 - a - b, b) diff --git a/test/gamma.jl b/test/gamma.jl index 8f49290f..5ee9544b 100644 --- a/test/gamma.jl +++ b/test/gamma.jl @@ -268,6 +268,32 @@ end @test logbeta(-2.0, 2.1) == Inf end + @testset "one of arguments is 1" begin + x = rand() + y = 100 + rand() + for (a, b) in ( + (1.0, 1.0), + (3.0, 1.0), + (1.0, 4.0), + (1.0, x), + (x, 1.0), + (1.0, y), + (y, 1.0), + (1.0, -x), + (-x, 1.0), + (1.0, -y), + (-y, 1.0), + ) + z = a == 1 ? b : a + @test beta(a, b) == inv(z) + @test logabsbeta(a, b) == (-log(abs(z)), z > 0 ? 1 : -1) + + if z > 0 + @test logbeta(a, b) == -log(z) + end + end + end + @testset "large difference in magnitude" begin @test beta(1e4, 1.5) ≈ 8.86193693673874630607029e-7 rtol=1e-14 @test logabsbeta(1e4, 1.5)[1] ≈ log(8.86193693673874630607029e-7) rtol=1e-14 @@ -281,6 +307,17 @@ end @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 + # check that results match the ones we get with MPFR + if Base.MPFR.version() >= v"4.0.0" + function beta_mpfr(x::BigFloat, y::BigFloat) + return invoke(SpecialFunctions._beta, Tuple{BigFloat,BigFloat}, x, y) + end + @test beta(big(3), big(5)) == beta_mpfr(big(3.0), big(5.0)) + @test beta(big(3//2), big(7//2)) == beta_mpfr(big(1.5), big(3.5)) + @test beta(big(1e4), big(3//2)) == beta_mpfr(big(1e4), big(1.5)) + @test beta(big(1e8), big(1//2)) == beta_mpfr(big(1e8), big(0.5)) + end + @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 From 30174c5a3130593c0bcfc2bec67e395fe4fea5d1 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Tue, 28 Sep 2021 19:08:52 +0200 Subject: [PATCH 2/3] Bump version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index e0d9275d..322b245f 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "SpecialFunctions" uuid = "276daf66-3868-5448-9aa4-cd146d93841b" -version = "1.7.0" +version = "1.7.1" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" From 7d4a57610b5824b6e8aa078be6dd2636fc48339b Mon Sep 17 00:00:00 2001 From: David Widmann Date: Mon, 11 Oct 2021 23:10:32 +0200 Subject: [PATCH 3/3] Add comment --- src/gamma.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/gamma.jl b/src/gamma.jl index 515a5ce5..c93f0027 100644 --- a/src/gamma.jl +++ b/src/gamma.jl @@ -750,8 +750,9 @@ Euler integral of the first kind ``\\operatorname{B}(x,y) = \\Gamma(x)\\Gamma(y) """ beta(a::Number, b::Number) = _beta(map(float, promote(a, b))...) -# here `T` is a floating point type but we don't want to restrict the implementation -# to `AbstractFloat` or `Float64` +# here `T` is a floating point type (e.g., `Float64` or `ComplexF64`) since +# it is called by `beta` above with arguments converted with `float` +# we don't want to restrict the implementation to `AbstractFloat` or `Float64` though function _beta(a::T, b::T) where {T<:Number} # two special cases a == 1 && return inv(b)