Fix beta for negative args

src/gamma.jl

@@ -719,18 +719,61 @@ 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 a == floor(a)
+            if isinteger(a)
+                if isinteger(b) && 1-a-b > 0
+                    sgn = (isinteger(b/2)) ? 1 : -1
+                    return sgn* beta(1-a-b,b)
+                else
+                    return NaN
+                end
+            else
+                return NaN
+            end
+        end
+    end
+    if b <= 0.0
+        if b == floor(b)
+            if isinteger(b)
+                if isinteger(b) && 1-a-b > 0
+                    sgn = (isinteger(b/2)) ? 1 : -1
+                    return sgn* beta(1-a-b,b)
+                else
+                    return NaN
+                end
+            else
+                return NaN
+            end
+        end
+    end
+    if abs(a) < abs(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
+        r, sgn = logabsgamma(b)
+
+        r-= b*log(a)
+        r+= b*(1-b)/(2*a)
+        r+= b*(1-b)*(1-2*b)/(12*a^2)
+        r+= -b^2*(1-b)^2/(12*a^3)
+
+        return sgn*exp(r)
+    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
 
 """

test/runtests.jl

@@ -802,6 +802,8 @@ 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 beta(1e8,0.5) ≈ 0.00017724538531210809
     @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))