Merge pull request #325 from devmotion/dw/logabsbeta

Fix `logabsbeta` for BigFloat

src/beta_inc.jl

@@ -9,8 +9,10 @@ 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))...)
 
+# 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

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
 

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