From 15a599ce27afe72ad052fd7738e89232c6b506c9 Mon Sep 17 00:00:00 2001 From: Sumegh-git Date: Sun, 23 Jun 2019 22:33:21 +0530 Subject: [PATCH 01/10] fix beta for negative args --- src/gamma.jl | 72 ++++++++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 1 + 2 files changed, 73 insertions(+) diff --git a/src/gamma.jl b/src/gamma.jl index f61dc439..63ea2718 100644 --- a/src/gamma.jl +++ b/src/gamma.jl @@ -712,6 +712,28 @@ loggamma(z::Complex{T}) where {T<:Union{Float32,Float16}} = Complex{T}(loggamma( gamma(z::Complex) = exp(loggamma(z)) +#Special case for negative integer argument +function beta_neg(a::Integer, b::Float64) + if b == trunc(Int,b) && 1-a-b > 0 + sgn = ((trunc(Int,b))%2 == 0) ? 1 : -1 + return sgn* beta(1-a-b,b) + else + return error("overflow") + end +end + +#asymptotic expansion for log(B(a,b)) for |a| >> |b| +function beta_asymp(a::Float64, b::Float64) + r, sgn = loggamma(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 + """ beta(x, y) @@ -720,6 +742,31 @@ Euler integral of the first kind ``\\operatorname{B}(x,y) = \\Gamma(x)\\Gamma(y) function beta end function beta(x::Real, w::Real) + if x <= 0.0 + if x == floor(x) + if x == trunc(Int,x) + return beta_neg(trunc(Int,x), w) + else + return error("overflow") + end + end + end + if w <= 0.0 + if w == floor(w) + if w == trunc(Int,w) + return beta_neg(trunc(Int,w), x) + else + return error("overflow") + end + end + end + if abs(x) < abs(w) + x,w = w,x + end + + if abs(x) > 1e5*abs(w) && abs(x) > 1e5 + return beta_asymp(x,w) + end yx, sx = logabsgamma(x) yw, sw = logabsgamma(w) yxw, sxw = logabsgamma(x+w) @@ -727,6 +774,31 @@ function beta(x::Real, w::Real) end function beta(x::Number, w::Number) + if x <= 0.0 + if x == floor(x) + if x == trunc(Int,x) + return beta_neg(trunc(Int,x), w) + else + return error("overflow") + end + end + end + if w <= 0.0 + if w == floor(w) + if w == trunc(Int,w) + return beta_neg(trunc(Int,w), x) + else + return error("overflow") + end + end + end + if abs(x) < abs(w) + x,w = w,x + end + + if abs(x) > 1e4*abs(w) && abs(x) > 1e4 + return beta_asymp(x,w) + end yx = loggamma(x) yw = loggamma(w) yxw = loggamma(x+w) diff --git a/test/runtests.jl b/test/runtests.jl index b2ad1ac1..fe480cbd 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -802,6 +802,7 @@ 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(-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)) From 28df26d8d8e6461f6f247651448839f95b423112 Mon Sep 17 00:00:00 2001 From: Sumegh-git Date: Sun, 23 Jun 2019 22:38:52 +0530 Subject: [PATCH 02/10] remove for complex --- src/gamma.jl | 7 ------- 1 file changed, 7 deletions(-) diff --git a/src/gamma.jl b/src/gamma.jl index 63ea2718..18be30f0 100644 --- a/src/gamma.jl +++ b/src/gamma.jl @@ -792,13 +792,6 @@ function beta(x::Number, w::Number) end end end - if abs(x) < abs(w) - x,w = w,x - end - - if abs(x) > 1e4*abs(w) && abs(x) > 1e4 - return beta_asymp(x,w) - end yx = loggamma(x) yw = loggamma(w) yxw = loggamma(x+w) From 20229026f3a039aea8b9d622c1ca545e8d81e5b0 Mon Sep 17 00:00:00 2001 From: Sumegh-git Date: Sun, 23 Jun 2019 22:40:07 +0530 Subject: [PATCH 03/10] fix --- src/gamma.jl | 18 ------------------ 1 file changed, 18 deletions(-) diff --git a/src/gamma.jl b/src/gamma.jl index 18be30f0..b171dde8 100644 --- a/src/gamma.jl +++ b/src/gamma.jl @@ -774,24 +774,6 @@ function beta(x::Real, w::Real) end function beta(x::Number, w::Number) - if x <= 0.0 - if x == floor(x) - if x == trunc(Int,x) - return beta_neg(trunc(Int,x), w) - else - return error("overflow") - end - end - end - if w <= 0.0 - if w == floor(w) - if w == trunc(Int,w) - return beta_neg(trunc(Int,w), x) - else - return error("overflow") - end - end - end yx = loggamma(x) yw = loggamma(w) yxw = loggamma(x+w) From 0a5af4b49f43192c8f3d73c8aae8460284180008 Mon Sep 17 00:00:00 2001 From: Sumegh-git Date: Tue, 25 Jun 2019 02:36:45 +0530 Subject: [PATCH 04/10] changes as reviewed --- src/gamma.jl | 88 +++++++++++++++++++++++------------------------- test/runtests.jl | 1 + 2 files changed, 43 insertions(+), 46 deletions(-) diff --git a/src/gamma.jl b/src/gamma.jl index b171dde8..4e83c2aa 100644 --- a/src/gamma.jl +++ b/src/gamma.jl @@ -712,28 +712,6 @@ loggamma(z::Complex{T}) where {T<:Union{Float32,Float16}} = Complex{T}(loggamma( gamma(z::Complex) = exp(loggamma(z)) -#Special case for negative integer argument -function beta_neg(a::Integer, b::Float64) - if b == trunc(Int,b) && 1-a-b > 0 - sgn = ((trunc(Int,b))%2 == 0) ? 1 : -1 - return sgn* beta(1-a-b,b) - else - return error("overflow") - end -end - -#asymptotic expansion for log(B(a,b)) for |a| >> |b| -function beta_asymp(a::Float64, b::Float64) - r, sgn = loggamma(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 - """ beta(x, y) @@ -741,43 +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) - if x <= 0.0 - if x == floor(x) - if x == trunc(Int,x) - return beta_neg(trunc(Int,x), w) +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 error("overflow") + return NaN end end end - if w <= 0.0 - if w == floor(w) - if w == trunc(Int,w) - return beta_neg(trunc(Int,w), x) + 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 error("overflow") + return NaN end end end - if abs(x) < abs(w) - x,w = w,x + 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) - if abs(x) > 1e5*abs(w) && abs(x) > 1e5 - return beta_asymp(x,w) + return sgn*exp(r) end - yx, sx = logabsgamma(x) - yw, sw = logabsgamma(w) - yxw, sxw = logabsgamma(x+w) - return exp(yx + yw - yxw) * (sx*sw*sxw) + 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 """ diff --git a/test/runtests.jl b/test/runtests.jl index fe480cbd..27c81162 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -803,6 +803,7 @@ end @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)) From e858da3cbf1859c42cb05d9bd9131796a2634987 Mon Sep 17 00:00:00 2001 From: Sumegh-git Date: Sat, 29 Jun 2019 00:49:19 +0530 Subject: [PATCH 05/10] refactor code --- src/gamma.jl | 37 ++++++++++++++----------------------- 1 file changed, 14 insertions(+), 23 deletions(-) diff --git a/src/gamma.jl b/src/gamma.jl index 4e83c2aa..ec113bf6 100644 --- a/src/gamma.jl +++ b/src/gamma.jl @@ -722,38 +722,29 @@ function beta end 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 + if isinteger(a) && isinteger(b) && 1-a-b > 0 + sgn = isinteger(b/2) ? 1 : -1 + return sgn* beta(1-a-b,b) + else + return NaN 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 + if isinteger(a) && isinteger(b) && 1-a-b > 0 + sgn = isinteger(a/2) ? 1 : -1 + return sgn* beta(1-a-b,a) + else + return NaN end end - if abs(a) < abs(b) + 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 + if a > 8 + return exp(loggammadiv(b,a) + loggamma(b)) + end r, sgn = logabsgamma(b) r-= b*log(a) From 2cd508c14fef0277d7659f49950ff43d2c93dbc4 Mon Sep 17 00:00:00 2001 From: Sumegh-git Date: Sat, 29 Jun 2019 00:58:04 +0530 Subject: [PATCH 06/10] fix bug --- src/gamma.jl | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/src/gamma.jl b/src/gamma.jl index ec113bf6..ee84235c 100644 --- a/src/gamma.jl +++ b/src/gamma.jl @@ -722,19 +722,15 @@ function beta end function beta(a::Real, b::Real) #Special case for negative integer argument if a <= 0.0 - if isinteger(a) && isinteger(b) && 1-a-b > 0 + if isinteger(a) && 1-a-b > 0 sgn = isinteger(b/2) ? 1 : -1 return sgn* beta(1-a-b,b) - else - return NaN end end if b <= 0.0 - if isinteger(a) && isinteger(b) && 1-a-b > 0 + if isinteger(b) && 1-a-b > 0 sgn = isinteger(a/2) ? 1 : -1 return sgn* beta(1-a-b,a) - else - return NaN end end if a < b From 4f4da3d526321f8f6c61564f2bc18e3d5b659ac5 Mon Sep 17 00:00:00 2001 From: Sumegh-git Date: Mon, 1 Jul 2019 15:53:07 +0530 Subject: [PATCH 07/10] remove redundant condition --- src/gamma.jl | 12 +----------- 1 file changed, 1 insertion(+), 11 deletions(-) diff --git a/src/gamma.jl b/src/gamma.jl index ee84235c..c91a4bbe 100644 --- a/src/gamma.jl +++ b/src/gamma.jl @@ -738,17 +738,7 @@ function beta(a::Real, b::Real) end #asymptotic expansion for log(B(a,b)) for |a| >> |b| if abs(a) > 1e5*abs(b) && abs(a) > 1e5 - if a > 8 - return exp(loggammadiv(b,a) + loggamma(b)) - end - 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) + return exp(loggammadiv(b,a) + loggamma(b)) end ya, sa = logabsgamma(a) yb, sb = logabsgamma(b) From 2a828a4cee4b70aca3c06cddf73e313b580bbea4 Mon Sep 17 00:00:00 2001 From: Sumegh-git Date: Fri, 5 Jul 2019 13:10:03 +0530 Subject: [PATCH 08/10] changes for logabsbeta --- src/gamma.jl | 31 +++++++++++++++++++++++++------ test/runtests.jl | 2 ++ 2 files changed, 27 insertions(+), 6 deletions(-) diff --git a/src/gamma.jl b/src/gamma.jl index c91a4bbe..cd4af9ca 100644 --- a/src/gamma.jl +++ b/src/gamma.jl @@ -760,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) @@ -769,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 sgn* 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 sgn* 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) + 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 diff --git a/test/runtests.jl b/test/runtests.jl index 27c81162..bca5aed4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -803,7 +803,9 @@ end @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) ≈ -0.69314718055994529 @test beta(1e8,0.5) ≈ 0.00017724538531210809 + @test logabsbeta(1e8,0.5) ≈ -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)) From be0cd39316856384fa70a507baad05c9a45a8a19 Mon Sep 17 00:00:00 2001 From: Sumegh-git Date: Fri, 5 Jul 2019 13:16:52 +0530 Subject: [PATCH 09/10] return tuple for consistency --- src/gamma.jl | 6 +++--- test/runtests.jl | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/gamma.jl b/src/gamma.jl index cd4af9ca..3ab73c8c 100644 --- a/src/gamma.jl +++ b/src/gamma.jl @@ -773,13 +773,13 @@ function logabsbeta(a::Real, b::Real) if a <= 0.0 if isinteger(a) && 1-a-b > 0 sgn = isinteger(b/2) ? 1 : -1 - return sgn* logabsbeta(1-a-b,b) + return logabsbeta(1-a-b,b), sgn end end if b <= 0.0 if isinteger(b) && 1-a-b > 0 sgn = isinteger(a/2) ? 1 : -1 - return sgn* logabsbeta(1-a-b,a) + return logabsbeta(1-a-b,a), sgn end end if a < b @@ -787,7 +787,7 @@ function logabsbeta(a::Real, b::Real) 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) + return (loggammadiv(b,a) + loggamma(b)), 1 end ya, sa = logabsgamma(a) yb, sb = logabsgamma(b) diff --git a/test/runtests.jl b/test/runtests.jl index bca5aed4..fc95afc2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -803,9 +803,9 @@ end @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) ≈ -0.69314718055994529 + @test logabsbeta(-2.0,2.0)[1] ≈ -0.69314718055994529 @test beta(1e8,0.5) ≈ 0.00017724538531210809 - @test logabsbeta(1e8,0.5) ≈ -8.637975427801484 + @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)) From 2262c86758676b3eb56665ae5b46f2f346d7a28b Mon Sep 17 00:00:00 2001 From: Sumegh-git Date: Fri, 5 Jul 2019 13:24:34 +0530 Subject: [PATCH 10/10] fix --- src/gamma.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/gamma.jl b/src/gamma.jl index 3ab73c8c..dcb92103 100644 --- a/src/gamma.jl +++ b/src/gamma.jl @@ -773,13 +773,13 @@ 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), sgn + 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), sgn + return logabsbeta(1-a-b,a) end end if a < b