Handle arguments of 1 in beta and logabsbeta

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"

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)

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