Update src/gamma.jl

Co-authored-by: David Widmann <[email protected]>

Author: andreasnoack

src/beta_inc.jl

@@ -630,7 +630,7 @@ function beta_inc_diff(a::Float64, b::Float64, x::Float64, y::Float64, n::Intege
         d = exp(-t)
     end
 
-    ans = beta_integrand(a,b,x,y,mu)/a
+    ans = beta_integrand(a, b, x, y, mu)/a
     if n == 1 || ans == 0.0
         return ans
     end
@@ -820,53 +820,53 @@ function beta_inc(a::Float64, b::Float64, x::Float64, y::Float64)
     end
 
     if b0 < min(epps, epps*a0)
-        p = beta_inc_power_series2(a0,b0,x0,epps)
+        p = beta_inc_power_series2(a0, b0, x0, epps)
         q = 1.0 - p
     elseif a0 < min(epps, epps*b0) && b0*x0 <= 1.0
-        q = beta_inc_power_series1(a0,b0,x0,epps)
+        q = beta_inc_power_series1(a0, b0, x0, epps)
         p = 1.0 - q
     elseif max(a0,b0) > 1.0
         if b0 <= 1.0
-            p = beta_inc_power_series(a0,b0,x0,epps)
+            p = beta_inc_power_series(a0, b0, x0, epps)
             q = 1.0 - p
         elseif x0 >= 0.3
-            q = beta_inc_power_series(b0,a0,y0,epps)
+            q = beta_inc_power_series(b0, a0, y0, epps)
             p = 1.0 - q
         elseif x0 >= 0.1
             if b0 > 15.0
-                q = beta_inc_asymptotic_asymmetric(b0,a0,y0,x0,q,15.0*eps())
+                q = beta_inc_asymptotic_asymmetric(b0, a0, y0, x0, q, 15.0*eps())
                 p = 1.0 - q
             else
                 n = 20
-                q = beta_inc_diff(b0,a0,y0,x0,n,epps)
+                q = beta_inc_diff(b0, a0, y0, x0, n, epps)
                 b0 += n
-                q = beta_inc_asymptotic_asymmetric(b0,a0,y0,x0,q,15.0*eps())
+                q = beta_inc_asymptotic_asymmetric(b0, a0, y0, x0, q, 15.0*eps())
                 p = 1.0 - q
             end
         elseif (x0*b0)^(a0) <= 0.7
-            p = beta_inc_power_series(a0,b0,x0,epps)
+            p = beta_inc_power_series(a0, b0, x0, epps)
             q = 1.0 - p
         else
             n = 20
-            q = beta_inc_diff(b0,a0,y0,x0,n,epps)
+            q = beta_inc_diff(b0, a0, y0, x0, n, epps)
             b0 += n
-            q = beta_inc_asymptotic_asymmetric(b0,a0,y0,x0,q,15.0*eps())
+            q = beta_inc_asymptotic_asymmetric(b0, a0, y0, x0, q, 15.0*eps())
             p = 1.0 - q
         end
     elseif a0 >= min(0.2, b0)
-        p = beta_inc_power_series(a0,b0,x0,epps)
+        p = beta_inc_power_series(a0, b0, x0, epps)
         q = 1.0 - p
     elseif x0^a0 <= 0.9
-        p = beta_inc_power_series(a0,b0,x0,epps)
+        p = beta_inc_power_series(a0, b0, x0, epps)
         q = 1.0 - p
     elseif x0 >= 0.3
-        q = beta_inc_power_series(b0,a0,y0,epps)
+        q = beta_inc_power_series(b0, a0, y0, epps)
         p = 1.0 - q
     else
         n = 20
-        q = beta_inc_diff(b0,a0,y0,x0,n,epps)
+        q = beta_inc_diff(b0, a0, y0, x0, n, epps)
         b0 += n
-        q = beta_inc_asymptotic_asymmetric(b0,a0,y0,x0,q,15.0*eps())
+        q = beta_inc_asymptotic_asymmetric(b0, a0, y0, x0, q, 15.0*eps())
         p = 1.0 - q
     end
 

src/gamma.jl

@@ -806,12 +806,12 @@ function logabsbeta(a::T, b::T) where T<:Real
     end
 
     if a <= 0 && isinteger(a)
-        if a + b <= 0
+        if a + b <= 0 && isinteger(b)
             r = logbeta(1 - a - b, b)
-            sgn = isinteger(b/2) ? 1 : -1
+            sgn = iseven(Int(b)) ? 1 : -1
             return r, sgn
         else
-            return -log(zero(first(promote(a, b)))), 1
+            return -log(zero(a)), 1
         end
     end
 

test/gamma.jl

@@ -255,11 +255,14 @@ end
     @test beta(5, 4)     ≈ beta(4,5)
 
     @testset "negative integer argument" begin
-        @test beta(-2.0, 2.0)          ≈ 1/2     rtol=1e-14
-        @test beta(-5.0, 2.0)          ≈ 1/20    rtol=1e-14
+        @test beta(-2.0, 1.0)          ≈ -1/2    rtol=1e-14
+        @test beta(-2.0, 2.0)          ≈  1/2    rtol=1e-14
+        @test beta(-5.0, 2.0)          ≈  1/20   rtol=1e-14
         @test logabsbeta(-2.0, 2.0)[1] ≈ -log(2) rtol=1e-14
         @test beta(-2.0, -2.0)         == Inf
         @test logbeta(-2.0, -2.0)      == Inf
+        @test beta(-2.0, 1.9)          == Inf
+        @test logbeta(-2.0, 1.9)       == Inf
         @test beta(-2.0, 2.1)          == Inf
         @test logbeta(-2.0, 2.1)       == Inf
     end