Fix beta/logabsbeta for negative args (#169)

src/gamma.jl

@@ -719,18 +719,38 @@ Euler integral of the first kind ``\\operatorname{B}(x,y) = \\Gamma(x)\\Gamma(y)
 """
 function beta end
 
-function beta(x::Real, w::Real)
-    yx, sx = logabsgamma(x)
-    yw, sw = logabsgamma(w)
-    yxw, sxw = logabsgamma(x+w)
-    return exp(yx + yw - yxw) * (sx*sw*sxw)
+function beta(a::Real, b::Real)
+    #Special case for negative integer argument
+    if a <= 0.0
+        if isinteger(a) && 1-a-b > 0
+            sgn = isinteger(b/2) ? 1 : -1
+            return sgn* beta(1-a-b,b)
+        end
+    end
+    if b <= 0.0
+        if isinteger(b) && 1-a-b > 0
+            sgn = isinteger(a/2) ? 1 : -1
+            return sgn* beta(1-a-b,a)
+        end
+    end
+    if a < b
+        a,b = b,a
+    end
+    #asymptotic expansion for log(B(a,b)) for |a| >> |b|
+    if abs(a) > 1e5*abs(b) && abs(a) > 1e5
+        return exp(loggammadiv(b,a) + loggamma(b))
+    end
+    ya, sa = logabsgamma(a)
+    yb, sb = logabsgamma(b)
+    yab, sab = logabsgamma(a+b)
+    return exp(ya + yb - yab) * (sa*sb*sab)
 end
 
-function beta(x::Number, w::Number)
-    yx = loggamma(x)
-    yw = loggamma(w)
-    yxw = loggamma(x+w)
-    return exp(yx + yw - yxw)
+function beta(a::Number, b::Number)
+    ya = loggamma(a)
+    yb = loggamma(b)
+    yab = loggamma(a+b)
+    return exp(ya + yb - yab)
 end
 
 """
@@ -740,7 +760,7 @@ Natural logarithm of the [`beta`](@ref) function ``\\log(|\\operatorname{B}(x,y)
 
 See also [`logabsbeta`](@ref).
 """
-logbeta(x::Number, w::Number) = loggamma(x)+loggamma(w)-loggamma(x+w)
+logbeta(a::Number, b::Number) = loggamma(a)+loggamma(b)-loggamma(a+b)
 
 """
     logabsbeta(x, y)
@@ -749,11 +769,30 @@ Compute the natural logarithm of the absolute value of the [`beta`](@ref) functi
 
 See also [`logbeta`](@ref).
 """
-function logabsbeta(x::Real, w::Real)
-    yx, sx = logabsgamma(x)
-    yw, sw = logabsgamma(w)
-    yxw, sxw = logabsgamma(x+w)
-    (yx + yw - yxw), (sx*sw*sxw)
+function logabsbeta(a::Real, b::Real)
+    if a <= 0.0
+        if isinteger(a) && 1-a-b > 0
+            sgn = isinteger(b/2) ? 1 : -1
+            return logabsbeta(1-a-b,b)
+        end
+    end
+    if b <= 0.0
+        if isinteger(b) && 1-a-b > 0
+            sgn = isinteger(a/2) ? 1 : -1
+            return logabsbeta(1-a-b,a)
+        end
+    end
+    if a < b
+        a,b = b,a
+    end
+    #asymptotic expansion for log(B(a,b)) for |a| >> |b|
+    if abs(a) > 1e5*abs(b) && abs(a) > 1e5
+        return (loggammadiv(b,a) + loggamma(b)), 1
+    end
+    ya, sa = logabsgamma(a)
+    yb, sb = logabsgamma(b)
+    yab, sab = logabsgamma(a+b)
+    (ya + yb - yab), (sa*sb*sab)
 end
 ## from base/mpfr.jl
 

test/runtests.jl

@@ -850,6 +850,10 @@ end
     @test beta(3,5) ≈ 1/105
     @test logbeta(5,4) ≈ log(beta(5,4))
     @test beta(5,4) ≈ beta(4,5)
+    @test beta(-2.0,2.0) ≈ 0.5
+    @test logabsbeta(-2.0,2.0)[1] ≈ -0.69314718055994529
+    @test beta(1e8,0.5) ≈ 0.00017724538531210809
+    @test logabsbeta(1e8,0.5)[1] ≈ -8.637975427801484
     @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))