From ef1e530b6c5b9247b07cd1ab6361295b80329a03 Mon Sep 17 00:00:00 2001
From: CarloLucibello <carlo.lucibello@gmail.com>
Date: Thu, 26 Nov 2015 18:30:51 -0500
Subject: [PATCH 01/14] generalize binomial to real and complex arguments

---
 base/special/gamma.jl | 4 ++++
 test/math.jl          | 5 +++++
 2 files changed, 9 insertions(+)

diff --git a/base/special/gamma.jl b/base/special/gamma.jl
index 0678287e1b064..5b23cfa879dda 100644
--- a/base/special/gamma.jl
+++ b/base/special/gamma.jl
@@ -482,3 +482,7 @@ eta(x::Integer) = eta(Float64(x))
 eta(x::Real)    = oftype(float(x),eta(Float64(x)))
 eta(z::Complex) = oftype(float(z),eta(Complex128(z)))
 @vectorize_1arg Number eta
+
+Base.binomial{S<:Number,T<:Number}(n::S, k::T) = ( (isinteger(n) && isinteger(k))
+                                ? convert(promote_type(T,S), binomial(round(Int64, n), round(Int64, k)))
+                                : gamma(1+n) / (gamma(1+n-k) * gamma(1+k)) )
diff --git a/test/math.jl b/test/math.jl
index b3039e8decf2c..713753ce98a03 100644
--- a/test/math.jl
+++ b/test/math.jl
@@ -712,3 +712,8 @@ end
 let A = [1 2; 3 4]; B = [5 6; 7 8]; C = [9 10; 11 12]
     @test muladd(A,B,C) == A*B + C
 end
+
+@test binomial(2.,1) == 2
+@test binomial(-2.,1) == -2
+@test typeof(binomial(-2.,1)) == Float64
+@test_approx_eq binomial(-2.1,1.1) -4.21013511374008

From d0055a1a5777d52c937da85f8216ff8b58294683 Mon Sep 17 00:00:00 2001
From: CarloLucibello <carlo.lucibello@gmail.com>
Date: Wed, 9 Dec 2015 13:24:22 +0100
Subject: [PATCH 02/14] port of CEPHES implementation of Beta function

---
 base/special/gamma.jl | 124 ++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 118 insertions(+), 6 deletions(-)

diff --git a/base/special/gamma.jl b/base/special/gamma.jl
index 99c77b7a2e6b3..70c073dc461a9 100644
--- a/base/special/gamma.jl
+++ b/base/special/gamma.jl
@@ -404,13 +404,125 @@ invdigamma(x::Float32) = Float32(invdigamma(Float64(x)))
 invdigamma(x::Real) = invdigamma(Float64(x))
 @vectorize_1arg Real invdigamma
 
-function beta(x::Number, w::Number)
-    yx, sx = lgamma_r(x)
-    yw, sw = lgamma_r(w)
-    yxw, sxw = lgamma_r(x+w)
-    return exp(yx + yw - yxw) * (sx*sw*sxw)
+
+const MAXGAM = 171.624376956302725
+const ASYMP_FACTOR = 1e6
+
+function beta(a::Number, b::Number)
+    a <= 0.0 && isinteger(a) && return beta_negint(a, b)
+    b <= 0.0 && isinteger(b) && return beta_negint(b, a)
+
+    if abs(a) < abs(b)
+        a, b = b, a
+    end
+
+    if abs(a) > ASYMP_FACTOR * abs(b) && a > ASYMP_FACTOR)
+        # Avoid loss of precision in lgamma(a + b) - lgamma(a)
+        y, s = lbeta_asymp(a, b)
+        return s * exp(y)
+    end
+
+    y = a + b
+    if abs(y) > MAXGAM || abs(a) > MAXGAM || abs(b) > MAXGAM)
+    	y, s = lgamma_r(y)
+    	yb, sb = lgamma_r(b)
+        y = yb - y
+    	ya, sa = lgamma_r(a)
+        y = ya + y
+    	# if (y > MAXLOG) {
+    	#     goto overflow;
+    	# }
+    	return s*sa*sb * exp(y)
+    end
+
+    y = gamma(y)
+    a = gamma(a)
+    b = gamma(b)
+    y == 0.0 && return Inf
+
+    if abs(abs(a) - abs(y)) > abs(abs(b) - abs(y)))
+        y = b / y
+        y *= a
+    else
+        y = a / y
+        y *= b
+    end
+
+    return y
 end
-lbeta(x::Number, w::Number) = lgamma(x)+lgamma(w)-lgamma(x+w)
+
+function lbeta(a::Number, b::Number)
+    a <= 0.0 && isinteger(a) && return lbeta_negint(a, b)
+    b <= 0.0 && isinteger(b) && return lbeta_negint(b, a)
+
+    if abs(a) < abs(b)
+        a, b = b, a
+    end
+
+    if abs(a) > ASYMP_FACTOR * abs(b) && a > ASYMP_FACTOR)
+        # Avoid loss of precision in lgamma(a + b) - lgamma(a)
+        y, s = lbeta_asymp(a, b)
+        return y
+    end
+
+    y = a + b
+    if abs(y) > MAXGAM || abs(a) > MAXGAM || abs(b) > MAXGAM)
+    	y, s = lgamma_r(y)
+    	yb, sb = lgamma_r(b)
+        y = yb - y
+    	ya, sa = lgamma_r(a)
+        y = ya + y
+    	return y
+    end
+
+    y = gamma(y)
+    a = gamma(a)
+    b = gamma(b)
+    y == 0.0 && return Inf
+
+    if abs(abs(a) - abs(y)) > abs(abs(b) - abs(y)))
+        y = b / y
+        y *= a
+    else
+        y = a / y
+        y *= b
+    end
+
+    return y < 0 ? log(-y) : log(y)
+end
+
+# assuming isinteger(x) and x < 0.
+function beta_negint(x::Number, w::Number)
+    if isinteger(w) && 1 - x - w > 0
+        s = ifelse(w % 2 , -1, 1)
+        return s * beta(1 - x - w, w)
+    else
+        return Inf
+    end
+end
+
+# assuming isinteger(x) and x < 0.
+function lbeta_negint(x::Number, w::Number)
+    if isinteger(w) && 1 - x - w > 0
+        s = ifelse(w % 2 , -1, 1)
+        return s * lbeta(1 - x - w, w)
+    else
+        return Inf
+    end
+end
+
+ # Asymptotic expansion for  ln(|B(a, b)|) for a > ASYMP_FACTOR*max(|b|, 1).
+function lbeta_asymp(a::Number, b::Number)
+    r, s = lgammma_r(b)
+    r -= b * log(a)
+
+    r += b*(1-b)/(2*a);
+    r += b*(1-b)*(1-2*b)/(12*a*a)
+    r += - b*b*(1-b)*(1-b)/(12*a*a*a)
+
+    return r, s
+end
+
 @vectorize_2arg Number beta
 @vectorize_2arg Number lbeta
 

From 67b47cf21ea005aff542f5209fb227dcdd8cfe8f Mon Sep 17 00:00:00 2001
From: CarloLucibello <carlo.lucibello@gmail.com>
Date: Thu, 10 Dec 2015 14:46:58 +0100
Subject: [PATCH 03/14] add some tests

---
 test/math.jl | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/test/math.jl b/test/math.jl
index cf219e39cad6a..4ab597b18d9e0 100644
--- a/test/math.jl
+++ b/test/math.jl
@@ -459,6 +459,10 @@ end
       0.00634645247782269506319336871208405439180447035257028310080 -
       0.00169495384841964531409376316336552555952269360134349446910im)
 
+ for i=1:10
+     @test beta(-i, i) == (-1)^i / i
+ end
+
 # gamma, lgamma (complex argument)
 if Base.Math.libm == "libopenlibm"
     @test gamma(Float64[1:25;]) == gamma(1:25)

From 088508c593d6287c4ed537c29ddf8062a9bdc159 Mon Sep 17 00:00:00 2001
From: CarloLucibello <carlo.lucibello@gmail.com>
Date: Thu, 10 Dec 2015 15:23:11 +0100
Subject: [PATCH 04/14] some fix for complexes, still TODO though

---
 base/special/gamma.jl | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/base/special/gamma.jl b/base/special/gamma.jl
index 70c073dc461a9..5250effa42a88 100644
--- a/base/special/gamma.jl
+++ b/base/special/gamma.jl
@@ -409,8 +409,8 @@ const MAXGAM = 171.624376956302725
 const ASYMP_FACTOR = 1e6
 
 function beta(a::Number, b::Number)
-    a <= 0.0 && isinteger(a) && return beta_negint(a, b)
-    b <= 0.0 && isinteger(b) && return beta_negint(b, a)
+    real(a) <= 0.0 && isinteger(a) && return beta_negint(a, b)
+    real(b) <= 0.0 && isinteger(b) && return beta_negint(b, a)
 
     if abs(a) < abs(b)
         a, b = b, a
@@ -452,8 +452,8 @@ function beta(a::Number, b::Number)
 end
 
 function lbeta(a::Number, b::Number)
-    a <= 0.0 && isinteger(a) && return lbeta_negint(a, b)
-    b <= 0.0 && isinteger(b) && return lbeta_negint(b, a)
+    real(a) <= 0.0 && isinteger(a) && return lbeta_negint(a, b)
+    real(b) <= 0.0 && isinteger(b) && return lbeta_negint(b, a)
 
     if abs(a) < abs(b)
         a, b = b, a
@@ -488,7 +488,7 @@ function lbeta(a::Number, b::Number)
         y *= b
     end
 
-    return y < 0 ? log(-y) : log(y)
+    return real(y) < 0 ? log(-y) : log(y)
 end
 
 # assuming isinteger(x) and x < 0.

From 7fbfaff3c58e79d9bb52d48aa153a237568d4889 Mon Sep 17 00:00:00 2001
From: CarloLucibello <carlo.lucibello@gmail.com>
Date: Thu, 10 Dec 2015 15:38:09 +0100
Subject: [PATCH 05/14] fixes

---
 base/special/gamma.jl | 17 ++++++++---------
 1 file changed, 8 insertions(+), 9 deletions(-)

diff --git a/base/special/gamma.jl b/base/special/gamma.jl
index 5250effa42a88..5038df79ba3df 100644
--- a/base/special/gamma.jl
+++ b/base/special/gamma.jl
@@ -404,7 +404,6 @@ invdigamma(x::Float32) = Float32(invdigamma(Float64(x)))
 invdigamma(x::Real) = invdigamma(Float64(x))
 @vectorize_1arg Real invdigamma
 
-
 const MAXGAM = 171.624376956302725
 const ASYMP_FACTOR = 1e6
 
@@ -416,14 +415,14 @@ function beta(a::Number, b::Number)
         a, b = b, a
     end
 
-    if abs(a) > ASYMP_FACTOR * abs(b) && a > ASYMP_FACTOR)
+    if abs(a) > ASYMP_FACTOR * abs(b) && a > ASYMP_FACTOR
         # Avoid loss of precision in lgamma(a + b) - lgamma(a)
         y, s = lbeta_asymp(a, b)
         return s * exp(y)
     end
 
     y = a + b
-    if abs(y) > MAXGAM || abs(a) > MAXGAM || abs(b) > MAXGAM)
+    if abs(y) > MAXGAM || abs(a) > MAXGAM || abs(b) > MAXGAM
     	y, s = lgamma_r(y)
     	yb, sb = lgamma_r(b)
         y = yb - y
@@ -440,7 +439,7 @@ function beta(a::Number, b::Number)
     b = gamma(b)
     y == 0.0 && return Inf
 
-    if abs(abs(a) - abs(y)) > abs(abs(b) - abs(y)))
+    if abs(abs(a) - abs(y)) > abs(abs(b) - abs(y))
         y = b / y
         y *= a
     else
@@ -459,14 +458,14 @@ function lbeta(a::Number, b::Number)
         a, b = b, a
     end
 
-    if abs(a) > ASYMP_FACTOR * abs(b) && a > ASYMP_FACTOR)
+    if abs(a) > ASYMP_FACTOR * abs(b) && a > ASYMP_FACTOR
         # Avoid loss of precision in lgamma(a + b) - lgamma(a)
         y, s = lbeta_asymp(a, b)
         return y
     end
 
     y = a + b
-    if abs(y) > MAXGAM || abs(a) > MAXGAM || abs(b) > MAXGAM)
+    if abs(y) > MAXGAM || abs(a) > MAXGAM || abs(b) > MAXGAM
     	y, s = lgamma_r(y)
     	yb, sb = lgamma_r(b)
         y = yb - y
@@ -480,7 +479,7 @@ function lbeta(a::Number, b::Number)
     b = gamma(b)
     y == 0.0 && return Inf
 
-    if abs(abs(a) - abs(y)) > abs(abs(b) - abs(y)))
+    if abs(abs(a) - abs(y)) > abs(abs(b) - abs(y))
         y = b / y
         y *= a
     else
@@ -494,7 +493,7 @@ end
 # assuming isinteger(x) and x < 0.
 function beta_negint(x::Number, w::Number)
     if isinteger(w) && 1 - x - w > 0
-        s = ifelse(w % 2 , -1, 1)
+        s = ifelse(w % 2 == 0, 1., -1.)
         return s * beta(1 - x - w, w)
     else
         return Inf
@@ -504,7 +503,7 @@ end
 # assuming isinteger(x) and x < 0.
 function lbeta_negint(x::Number, w::Number)
     if isinteger(w) && 1 - x - w > 0
-        s = ifelse(w % 2 , -1, 1)
+        s = ifelse(w % 2 == 0 , 1., -1.)
         return s * lbeta(1 - x - w, w)
     else
         return Inf

From 86017108f88736c075193985dcc61687869221f5 Mon Sep 17 00:00:00 2001
From: CarloLucibello <carlo.lucibello@gmail.com>
Date: Thu, 10 Dec 2015 15:47:05 +0100
Subject: [PATCH 06/14] fix lgammma

---
 base/special/gamma.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/base/special/gamma.jl b/base/special/gamma.jl
index 5038df79ba3df..087627b104abb 100644
--- a/base/special/gamma.jl
+++ b/base/special/gamma.jl
@@ -512,7 +512,7 @@ end
 
  # Asymptotic expansion for  ln(|B(a, b)|) for a > ASYMP_FACTOR*max(|b|, 1).
 function lbeta_asymp(a::Number, b::Number)
-    r, s = lgammma_r(b)
+    r, s = lgamma_r(b)
     r -= b * log(a)
 
     r += b*(1-b)/(2*a);

From efce6cc1e7481d2c4080e2289729799d28b67591 Mon Sep 17 00:00:00 2001
From: CarloLucibello <carlo.lucibello@gmail.com>
Date: Wed, 9 Dec 2015 13:24:22 +0100
Subject: [PATCH 07/14] port of CEPHES implementation of Beta function

---
 base/special/gamma.jl | 124 ++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 118 insertions(+), 6 deletions(-)

diff --git a/base/special/gamma.jl b/base/special/gamma.jl
index 5884ba418585a..7745020241662 100644
--- a/base/special/gamma.jl
+++ b/base/special/gamma.jl
@@ -404,13 +404,125 @@ invdigamma(x::Float32) = Float32(invdigamma(Float64(x)))
 invdigamma(x::Real) = invdigamma(Float64(x))
 @vectorize_1arg Real invdigamma
 
-function beta(x::Number, w::Number)
-    yx, sx = lgamma_r(x)
-    yw, sw = lgamma_r(w)
-    yxw, sxw = lgamma_r(x+w)
-    return exp(yx + yw - yxw) * (sx*sw*sxw)
+
+const MAXGAM = 171.624376956302725
+const ASYMP_FACTOR = 1e6
+
+function beta(a::Number, b::Number)
+    a <= 0.0 && isinteger(a) && return beta_negint(a, b)
+    b <= 0.0 && isinteger(b) && return beta_negint(b, a)
+
+    if abs(a) < abs(b)
+        a, b = b, a
+    end
+
+    if abs(a) > ASYMP_FACTOR * abs(b) && a > ASYMP_FACTOR)
+        # Avoid loss of precision in lgamma(a + b) - lgamma(a)
+        y, s = lbeta_asymp(a, b)
+        return s * exp(y)
+    end
+
+    y = a + b
+    if abs(y) > MAXGAM || abs(a) > MAXGAM || abs(b) > MAXGAM)
+    	y, s = lgamma_r(y)
+    	yb, sb = lgamma_r(b)
+        y = yb - y
+    	ya, sa = lgamma_r(a)
+        y = ya + y
+    	# if (y > MAXLOG) {
+    	#     goto overflow;
+    	# }
+    	return s*sa*sb * exp(y)
+    end
+
+    y = gamma(y)
+    a = gamma(a)
+    b = gamma(b)
+    y == 0.0 && return Inf
+
+    if abs(abs(a) - abs(y)) > abs(abs(b) - abs(y)))
+        y = b / y
+        y *= a
+    else
+        y = a / y
+        y *= b
+    end
+
+    return y
 end
-lbeta(x::Number, w::Number) = lgamma(x)+lgamma(w)-lgamma(x+w)
+
+function lbeta(a::Number, b::Number)
+    a <= 0.0 && isinteger(a) && return lbeta_negint(a, b)
+    b <= 0.0 && isinteger(b) && return lbeta_negint(b, a)
+
+    if abs(a) < abs(b)
+        a, b = b, a
+    end
+
+    if abs(a) > ASYMP_FACTOR * abs(b) && a > ASYMP_FACTOR)
+        # Avoid loss of precision in lgamma(a + b) - lgamma(a)
+        y, s = lbeta_asymp(a, b)
+        return y
+    end
+
+    y = a + b
+    if abs(y) > MAXGAM || abs(a) > MAXGAM || abs(b) > MAXGAM)
+    	y, s = lgamma_r(y)
+    	yb, sb = lgamma_r(b)
+        y = yb - y
+    	ya, sa = lgamma_r(a)
+        y = ya + y
+    	return y
+    end
+
+    y = gamma(y)
+    a = gamma(a)
+    b = gamma(b)
+    y == 0.0 && return Inf
+
+    if abs(abs(a) - abs(y)) > abs(abs(b) - abs(y)))
+        y = b / y
+        y *= a
+    else
+        y = a / y
+        y *= b
+    end
+
+    return y < 0 ? log(-y) : log(y)
+end
+
+# assuming isinteger(x) and x < 0.
+function beta_negint(x::Number, w::Number)
+    if isinteger(w) && 1 - x - w > 0
+        s = ifelse(w % 2 , -1, 1)
+        return s * beta(1 - x - w, w)
+    else
+        return Inf
+    end
+end
+
+# assuming isinteger(x) and x < 0.
+function lbeta_negint(x::Number, w::Number)
+    if isinteger(w) && 1 - x - w > 0
+        s = ifelse(w % 2 , -1, 1)
+        return s * lbeta(1 - x - w, w)
+    else
+        return Inf
+    end
+end
+
+ # Asymptotic expansion for  ln(|B(a, b)|) for a > ASYMP_FACTOR*max(|b|, 1).
+function lbeta_asymp(a::Number, b::Number)
+    r, s = lgammma_r(b)
+    r -= b * log(a)
+
+    r += b*(1-b)/(2*a);
+    r += b*(1-b)*(1-2*b)/(12*a*a)
+    r += - b*b*(1-b)*(1-b)/(12*a*a*a)
+
+    return r, s
+end
+
 @vectorize_2arg Number beta
 @vectorize_2arg Number lbeta
 

From d2de1b3354cb6011fb432e06e7f77f836160681a Mon Sep 17 00:00:00 2001
From: CarloLucibello <carlo.lucibello@gmail.com>
Date: Thu, 10 Dec 2015 14:46:58 +0100
Subject: [PATCH 08/14] add some tests

---
 test/math.jl | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/test/math.jl b/test/math.jl
index 09a6beb330c35..bcfabb927b16d 100644
--- a/test/math.jl
+++ b/test/math.jl
@@ -459,6 +459,10 @@ end
       0.00634645247782269506319336871208405439180447035257028310080 -
       0.00169495384841964531409376316336552555952269360134349446910im)
 
+ for i=1:10
+     @test beta(-i, i) == (-1)^i / i
+ end
+
 # gamma, lgamma (complex argument)
 if Base.Math.libm == "libopenlibm"
     @test gamma(Float64[1:25;]) == gamma(1:25)

From dab853de8035eda1e753b04748ea1c36b2edadf4 Mon Sep 17 00:00:00 2001
From: CarloLucibello <carlo.lucibello@gmail.com>
Date: Thu, 10 Dec 2015 15:23:11 +0100
Subject: [PATCH 09/14] some fix for complexes, still TODO though

---
 base/special/gamma.jl | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/base/special/gamma.jl b/base/special/gamma.jl
index 7745020241662..b9c6d3b8f830e 100644
--- a/base/special/gamma.jl
+++ b/base/special/gamma.jl
@@ -409,8 +409,8 @@ const MAXGAM = 171.624376956302725
 const ASYMP_FACTOR = 1e6
 
 function beta(a::Number, b::Number)
-    a <= 0.0 && isinteger(a) && return beta_negint(a, b)
-    b <= 0.0 && isinteger(b) && return beta_negint(b, a)
+    real(a) <= 0.0 && isinteger(a) && return beta_negint(a, b)
+    real(b) <= 0.0 && isinteger(b) && return beta_negint(b, a)
 
     if abs(a) < abs(b)
         a, b = b, a
@@ -452,8 +452,8 @@ function beta(a::Number, b::Number)
 end
 
 function lbeta(a::Number, b::Number)
-    a <= 0.0 && isinteger(a) && return lbeta_negint(a, b)
-    b <= 0.0 && isinteger(b) && return lbeta_negint(b, a)
+    real(a) <= 0.0 && isinteger(a) && return lbeta_negint(a, b)
+    real(b) <= 0.0 && isinteger(b) && return lbeta_negint(b, a)
 
     if abs(a) < abs(b)
         a, b = b, a
@@ -488,7 +488,7 @@ function lbeta(a::Number, b::Number)
         y *= b
     end
 
-    return y < 0 ? log(-y) : log(y)
+    return real(y) < 0 ? log(-y) : log(y)
 end
 
 # assuming isinteger(x) and x < 0.

From 0b93ab94c0006cd0e7375148f36d10aa275c8ef2 Mon Sep 17 00:00:00 2001
From: CarloLucibello <carlo.lucibello@gmail.com>
Date: Thu, 10 Dec 2015 15:38:09 +0100
Subject: [PATCH 10/14] fixes

---
 base/special/gamma.jl | 17 ++++++++---------
 1 file changed, 8 insertions(+), 9 deletions(-)

diff --git a/base/special/gamma.jl b/base/special/gamma.jl
index b9c6d3b8f830e..712e2e3559540 100644
--- a/base/special/gamma.jl
+++ b/base/special/gamma.jl
@@ -404,7 +404,6 @@ invdigamma(x::Float32) = Float32(invdigamma(Float64(x)))
 invdigamma(x::Real) = invdigamma(Float64(x))
 @vectorize_1arg Real invdigamma
 
-
 const MAXGAM = 171.624376956302725
 const ASYMP_FACTOR = 1e6
 
@@ -416,14 +415,14 @@ function beta(a::Number, b::Number)
         a, b = b, a
     end
 
-    if abs(a) > ASYMP_FACTOR * abs(b) && a > ASYMP_FACTOR)
+    if abs(a) > ASYMP_FACTOR * abs(b) && a > ASYMP_FACTOR
         # Avoid loss of precision in lgamma(a + b) - lgamma(a)
         y, s = lbeta_asymp(a, b)
         return s * exp(y)
     end
 
     y = a + b
-    if abs(y) > MAXGAM || abs(a) > MAXGAM || abs(b) > MAXGAM)
+    if abs(y) > MAXGAM || abs(a) > MAXGAM || abs(b) > MAXGAM
     	y, s = lgamma_r(y)
     	yb, sb = lgamma_r(b)
         y = yb - y
@@ -440,7 +439,7 @@ function beta(a::Number, b::Number)
     b = gamma(b)
     y == 0.0 && return Inf
 
-    if abs(abs(a) - abs(y)) > abs(abs(b) - abs(y)))
+    if abs(abs(a) - abs(y)) > abs(abs(b) - abs(y))
         y = b / y
         y *= a
     else
@@ -459,14 +458,14 @@ function lbeta(a::Number, b::Number)
         a, b = b, a
     end
 
-    if abs(a) > ASYMP_FACTOR * abs(b) && a > ASYMP_FACTOR)
+    if abs(a) > ASYMP_FACTOR * abs(b) && a > ASYMP_FACTOR
         # Avoid loss of precision in lgamma(a + b) - lgamma(a)
         y, s = lbeta_asymp(a, b)
         return y
     end
 
     y = a + b
-    if abs(y) > MAXGAM || abs(a) > MAXGAM || abs(b) > MAXGAM)
+    if abs(y) > MAXGAM || abs(a) > MAXGAM || abs(b) > MAXGAM
     	y, s = lgamma_r(y)
     	yb, sb = lgamma_r(b)
         y = yb - y
@@ -480,7 +479,7 @@ function lbeta(a::Number, b::Number)
     b = gamma(b)
     y == 0.0 && return Inf
 
-    if abs(abs(a) - abs(y)) > abs(abs(b) - abs(y)))
+    if abs(abs(a) - abs(y)) > abs(abs(b) - abs(y))
         y = b / y
         y *= a
     else
@@ -494,7 +493,7 @@ end
 # assuming isinteger(x) and x < 0.
 function beta_negint(x::Number, w::Number)
     if isinteger(w) && 1 - x - w > 0
-        s = ifelse(w % 2 , -1, 1)
+        s = ifelse(w % 2 == 0, 1., -1.)
         return s * beta(1 - x - w, w)
     else
         return Inf
@@ -504,7 +503,7 @@ end
 # assuming isinteger(x) and x < 0.
 function lbeta_negint(x::Number, w::Number)
     if isinteger(w) && 1 - x - w > 0
-        s = ifelse(w % 2 , -1, 1)
+        s = ifelse(w % 2 == 0 , 1., -1.)
         return s * lbeta(1 - x - w, w)
     else
         return Inf

From 8065f51e8f5ad0d11895df499b6c25236afbb135 Mon Sep 17 00:00:00 2001
From: CarloLucibello <carlo.lucibello@gmail.com>
Date: Thu, 10 Dec 2015 15:47:05 +0100
Subject: [PATCH 11/14] fix lgammma

---
 base/special/gamma.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/base/special/gamma.jl b/base/special/gamma.jl
index 712e2e3559540..5e64a06cc6c64 100644
--- a/base/special/gamma.jl
+++ b/base/special/gamma.jl
@@ -512,7 +512,7 @@ end
 
  # Asymptotic expansion for  ln(|B(a, b)|) for a > ASYMP_FACTOR*max(|b|, 1).
 function lbeta_asymp(a::Number, b::Number)
-    r, s = lgammma_r(b)
+    r, s = lgamma_r(b)
     r -= b * log(a)
 
     r += b*(1-b)/(2*a);

From 5e2e09f72dc9885547a480c250af8341221ffef5 Mon Sep 17 00:00:00 2001
From: CarloLucibello <carlo.lucibello@gmail.com>
Date: Tue, 19 Jan 2016 13:43:27 +0100
Subject: [PATCH 12/14] - Add credits to Cephes - use fastabs for complex
 numbers - add a test for asymptotic behaviour from issue #4301

---
 base/special/beta.jl  | 134 ++++++++++++++++++++++++++++++++++++++++++
 base/special/gamma.jl | 121 --------------------------------------
 test/math.jl          |   2 +
 3 files changed, 136 insertions(+), 121 deletions(-)
 create mode 100644 base/special/beta.jl

diff --git a/base/special/beta.jl b/base/special/beta.jl
new file mode 100644
index 0000000000000..c1153575c1d5e
--- /dev/null
+++ b/base/special/beta.jl
@@ -0,0 +1,134 @@
+import Base.Math.lgamma_r
+import Base.Math.lbeta
+import Base.Math.beta
+using Base.Test
+#=
+Port of Chephes' Beta function implementation as found in
+scipy libray https://github.com/scipy/scipy/blob/0cff7a5fe6226668729fc2551105692ce114c2b3/scipy/special/cephes/beta.c
+
+* Cephes Math Library Release 2.0:  April, 1987
+* Copyright 1984, 1987 by Stephen L. Moshier
+* Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+=#
+
+const MAXGAM = 171.624376956302725
+const ASYMP_FACTOR = 1e6
+
+fastabs(x::Real) = abs(x)
+fastabs(x::Complex) = abs(real(x)) + abs(imag(x))
+
+function beta(a::Number, b::Number)
+    real(a) <= 0.0 && isinteger(a) && return beta_negint(a, b)
+    real(b) <= 0.0 && isinteger(b) && return beta_negint(b, a)
+
+    if fastabs(a) < fastabs(b)
+        a, b = b, a
+    end
+
+    if fastabs(a) > ASYMP_FACTOR * fastabs(b) && a > ASYMP_FACTOR
+        # Avoid loss of precision in lgamma(a + b) - lgamma(a)
+        y, s = lbeta_asymp(a, b)
+        return s * exp(y)
+    end
+
+    y = a + b
+    if fastabs(y) > MAXGAM || fastabs(a) > MAXGAM || fastabs(b) > MAXGAM
+    	y, s = lgamma_r(y)
+    	yb, sb = lgamma_r(b)
+        y = yb - y
+    	ya, sa = lgamma_r(a)
+        y = ya + y
+    	# if (y > MAXLOG) {
+    	#     goto overflow;
+    	# }
+    	return s*sa*sb * exp(y)
+    end
+
+    y = gamma(y)
+    a = gamma(a)
+    b = gamma(b)
+    y == 0.0 && return Inf
+
+    if fastabs(fastabs(a) - fastabs(y)) > fastabs(fastabs(b) - fastabs(y))
+        y = b / y
+        y *= a
+    else
+        y = a / y
+        y *= b
+    end
+
+    return y
+end
+
+function lbeta(a::Number, b::Number)
+    real(a) <= 0.0 && isinteger(a) && return lbeta_negint(a, b)
+    real(b) <= 0.0 && isinteger(b) && return lbeta_negint(b, a)
+
+    if fastabs(a) < fastabs(b)
+        a, b = b, a
+    end
+
+    if fastabs(a) > ASYMP_FACTOR * fastabs(b) && a > ASYMP_FACTOR
+        # Avoid loss of precision in lgamma(a + b) - lgamma(a)
+        y, s = lbeta_asymp(a, b)
+        return y
+    end
+
+    y = a + b
+    if fastabs(y) > MAXGAM || fastabs(a) > MAXGAM || fastabs(b) > MAXGAM
+    	y, s = lgamma_r(y)
+    	yb, sb = lgamma_r(b)
+        y = yb - y
+    	ya, sa = lgamma_r(a)
+        y = ya + y
+    	return y
+    end
+
+    y = gamma(y)
+    a = gamma(a)
+    b = gamma(b)
+    y == 0.0 && return Inf
+
+    if fastabs(fastabs(a) - fastabs(y)) > fastabs(fastabs(b) - fastabs(y))
+        y = b / y
+        y *= a
+    else
+        y = a / y
+        y *= b
+    end
+
+    return real(y) < 0 ? log(-y) : log(y)
+end
+
+# assuming isinteger(x) and x < 0.
+function beta_negint(x::Number, w::Number)
+    if isinteger(w) && 1 - x - w > 0
+        s = ifelse(w % 2 == 0, 1., -1.)
+        return s * beta(1 - x - w, w)
+    else
+        return Inf
+    end
+end
+
+# assuming isinteger(x) and x < 0.
+function lbeta_negint(x::Number, w::Number)
+    if isinteger(w) && 1 - x - w > 0
+        s = ifelse(w % 2 == 0 , 1., -1.)
+        return s * lbeta(1 - x - w, w)
+    else
+        return Inf
+    end
+end
+
+ # Asymptotic expansion for  ln(|B(a, b)|) for a > ASYMP_FACTOR*max(|b|, 1).
+function lbeta_asymp(a::Number, b::Number)
+    r, s = lgamma_r(b)
+    r -= b * log(a)
+    r += b*(1-b)/(2*a);
+    r += b*(1-b)*(1-2*b)/(12*a*a)
+    r += - b*b*(1-b)*(1-b)/(12*a*a*a)
+    return r, s
+end
+
+@vectorize_2arg Number beta
+@vectorize_2arg Number lbeta
diff --git a/base/special/gamma.jl b/base/special/gamma.jl
index 5e64a06cc6c64..cbc709e7ec1d1 100644
--- a/base/special/gamma.jl
+++ b/base/special/gamma.jl
@@ -404,127 +404,6 @@ invdigamma(x::Float32) = Float32(invdigamma(Float64(x)))
 invdigamma(x::Real) = invdigamma(Float64(x))
 @vectorize_1arg Real invdigamma
 
-const MAXGAM = 171.624376956302725
-const ASYMP_FACTOR = 1e6
-
-function beta(a::Number, b::Number)
-    real(a) <= 0.0 && isinteger(a) && return beta_negint(a, b)
-    real(b) <= 0.0 && isinteger(b) && return beta_negint(b, a)
-
-    if abs(a) < abs(b)
-        a, b = b, a
-    end
-
-    if abs(a) > ASYMP_FACTOR * abs(b) && a > ASYMP_FACTOR
-        # Avoid loss of precision in lgamma(a + b) - lgamma(a)
-        y, s = lbeta_asymp(a, b)
-        return s * exp(y)
-    end
-
-    y = a + b
-    if abs(y) > MAXGAM || abs(a) > MAXGAM || abs(b) > MAXGAM
-    	y, s = lgamma_r(y)
-    	yb, sb = lgamma_r(b)
-        y = yb - y
-    	ya, sa = lgamma_r(a)
-        y = ya + y
-    	# if (y > MAXLOG) {
-    	#     goto overflow;
-    	# }
-    	return s*sa*sb * exp(y)
-    end
-
-    y = gamma(y)
-    a = gamma(a)
-    b = gamma(b)
-    y == 0.0 && return Inf
-
-    if abs(abs(a) - abs(y)) > abs(abs(b) - abs(y))
-        y = b / y
-        y *= a
-    else
-        y = a / y
-        y *= b
-    end
-
-    return y
-end
-
-function lbeta(a::Number, b::Number)
-    real(a) <= 0.0 && isinteger(a) && return lbeta_negint(a, b)
-    real(b) <= 0.0 && isinteger(b) && return lbeta_negint(b, a)
-
-    if abs(a) < abs(b)
-        a, b = b, a
-    end
-
-    if abs(a) > ASYMP_FACTOR * abs(b) && a > ASYMP_FACTOR
-        # Avoid loss of precision in lgamma(a + b) - lgamma(a)
-        y, s = lbeta_asymp(a, b)
-        return y
-    end
-
-    y = a + b
-    if abs(y) > MAXGAM || abs(a) > MAXGAM || abs(b) > MAXGAM
-    	y, s = lgamma_r(y)
-    	yb, sb = lgamma_r(b)
-        y = yb - y
-    	ya, sa = lgamma_r(a)
-        y = ya + y
-    	return y
-    end
-
-    y = gamma(y)
-    a = gamma(a)
-    b = gamma(b)
-    y == 0.0 && return Inf
-
-    if abs(abs(a) - abs(y)) > abs(abs(b) - abs(y))
-        y = b / y
-        y *= a
-    else
-        y = a / y
-        y *= b
-    end
-
-    return real(y) < 0 ? log(-y) : log(y)
-end
-
-# assuming isinteger(x) and x < 0.
-function beta_negint(x::Number, w::Number)
-    if isinteger(w) && 1 - x - w > 0
-        s = ifelse(w % 2 == 0, 1., -1.)
-        return s * beta(1 - x - w, w)
-    else
-        return Inf
-    end
-end
-
-# assuming isinteger(x) and x < 0.
-function lbeta_negint(x::Number, w::Number)
-    if isinteger(w) && 1 - x - w > 0
-        s = ifelse(w % 2 == 0 , 1., -1.)
-        return s * lbeta(1 - x - w, w)
-    else
-        return Inf
-    end
-end
-
- # Asymptotic expansion for  ln(|B(a, b)|) for a > ASYMP_FACTOR*max(|b|, 1).
-function lbeta_asymp(a::Number, b::Number)
-    r, s = lgamma_r(b)
-    r -= b * log(a)
-
-    r += b*(1-b)/(2*a);
-    r += b*(1-b)*(1-2*b)/(12*a*a)
-    r += - b*b*(1-b)*(1-b)/(12*a*a*a)
-
-    return r, s
-end
-
-@vectorize_2arg Number beta
-@vectorize_2arg Number lbeta
-
 # Riemann zeta function; algorithm is based on specializing the Hurwitz
 # zeta function above for z==1.
 function zeta(s::Union{Float64,Complex{Float64}})
diff --git a/test/math.jl b/test/math.jl
index bcfabb927b16d..95a5c1d40275e 100644
--- a/test/math.jl
+++ b/test/math.jl
@@ -459,6 +459,8 @@ end
       0.00634645247782269506319336871208405439180447035257028310080 -
       0.00169495384841964531409376316336552555952269360134349446910im)
 
+@test_approx_eq lbeta(1e10, 1e-10) 23.02585092758015 # comment in issue #4301
+
  for i=1:10
      @test beta(-i, i) == (-1)^i / i
  end

From c5e346d4c4054d030a72df7507c1e706ec227e6b Mon Sep 17 00:00:00 2001
From: CarloLucibello <carlo.lucibello@gmail.com>
Date: Sun, 24 Jan 2016 12:50:43 +0100
Subject: [PATCH 13/14] fix typewarns for real types. Still to address
 complexes.

---
 base/special/beta.jl | 33 +++++++++++++++++++--------------
 1 file changed, 19 insertions(+), 14 deletions(-)

diff --git a/base/special/beta.jl b/base/special/beta.jl
index c1153575c1d5e..54b9c44915d32 100644
--- a/base/special/beta.jl
+++ b/base/special/beta.jl
@@ -17,14 +17,16 @@ const ASYMP_FACTOR = 1e6
 fastabs(x::Real) = abs(x)
 fastabs(x::Complex) = abs(real(x)) + abs(imag(x))
 
-function beta(a::Number, b::Number)
+beta(a::Number, b::Number) = beta(promote(float(a), float(b))...)
+
+function beta{T<:Number}(a::T, b::T)
     real(a) <= 0.0 && isinteger(a) && return beta_negint(a, b)
     real(b) <= 0.0 && isinteger(b) && return beta_negint(b, a)
-
+    #
     if fastabs(a) < fastabs(b)
         a, b = b, a
     end
-
+    #
     if fastabs(a) > ASYMP_FACTOR * fastabs(b) && a > ASYMP_FACTOR
         # Avoid loss of precision in lgamma(a + b) - lgamma(a)
         y, s = lbeta_asymp(a, b)
@@ -47,7 +49,7 @@ function beta(a::Number, b::Number)
     y = gamma(y)
     a = gamma(a)
     b = gamma(b)
-    y == 0.0 && return Inf
+    y == 0.0 && return convert(T, Inf)
 
     if fastabs(fastabs(a) - fastabs(y)) > fastabs(fastabs(b) - fastabs(y))
         y = b / y
@@ -60,7 +62,9 @@ function beta(a::Number, b::Number)
     return y
 end
 
-function lbeta(a::Number, b::Number)
+lbeta(a::Number, b::Number) = lbeta(promote(float(a), float(b))...)
+
+function lbeta{T<:Number}(a::T, b::T)
     real(a) <= 0.0 && isinteger(a) && return lbeta_negint(a, b)
     real(b) <= 0.0 && isinteger(b) && return lbeta_negint(b, a)
 
@@ -87,7 +91,7 @@ function lbeta(a::Number, b::Number)
     y = gamma(y)
     a = gamma(a)
     b = gamma(b)
-    y == 0.0 && return Inf
+    y == 0.0 && return inf(T)
 
     if fastabs(fastabs(a) - fastabs(y)) > fastabs(fastabs(b) - fastabs(y))
         y = b / y
@@ -101,33 +105,34 @@ function lbeta(a::Number, b::Number)
 end
 
 # assuming isinteger(x) and x < 0.
-function beta_negint(x::Number, w::Number)
+function beta_negint{T}(x::T, w::T)
     if isinteger(w) && 1 - x - w > 0
-        s = ifelse(w % 2 == 0, 1., -1.)
+        s = ifelse(Int(w) % 2 == 0, T(1), T(-1))
         return s * beta(1 - x - w, w)
     else
-        return Inf
+        return convert(T, Inf)
     end
 end
 
 # assuming isinteger(x) and x < 0.
-function lbeta_negint(x::Number, w::Number)
+function lbeta_negint{T}(x::T, w::T)
     if isinteger(w) && 1 - x - w > 0
-        s = ifelse(w % 2 == 0 , 1., -1.)
+        s = ifelse(Int(w) % 2 == 0 , T(1), T(-1))
+        return s
         return s * lbeta(1 - x - w, w)
     else
-        return Inf
+        return convert(T, Inf)
     end
 end
 
  # Asymptotic expansion for  ln(|B(a, b)|) for a > ASYMP_FACTOR*max(|b|, 1).
-function lbeta_asymp(a::Number, b::Number)
+function lbeta_asymp{T}(a::T, b::T)
     r, s = lgamma_r(b)
     r -= b * log(a)
     r += b*(1-b)/(2*a);
     r += b*(1-b)*(1-2*b)/(12*a*a)
     r += - b*b*(1-b)*(1-b)/(12*a*a*a)
-    return r, s
+    return r, T(s)
 end
 
 @vectorize_2arg Number beta

From 8bf6196ce460ae6fe2d21b979c0cc4c3a68770ef Mon Sep 17 00:00:00 2001
From: CarloLucibello <carlo.lucibello@gmail.com>
Date: Mon, 25 Jan 2016 01:33:06 +0100
Subject: [PATCH 14/14] type stability for complex numbers

---
 base/special/beta.jl | 39 +++++++++++++++++++++------------------
 1 file changed, 21 insertions(+), 18 deletions(-)

diff --git a/base/special/beta.jl b/base/special/beta.jl
index 54b9c44915d32..0312ae2d24a3a 100644
--- a/base/special/beta.jl
+++ b/base/special/beta.jl
@@ -15,19 +15,19 @@ const MAXGAM = 171.624376956302725
 const ASYMP_FACTOR = 1e6
 
 fastabs(x::Real) = abs(x)
-fastabs(x::Complex) = abs(real(x)) + abs(imag(x))
+fastabs(x::Complex) = abs(real(x))
 
 beta(a::Number, b::Number) = beta(promote(float(a), float(b))...)
 
 function beta{T<:Number}(a::T, b::T)
     real(a) <= 0.0 && isinteger(a) && return beta_negint(a, b)
     real(b) <= 0.0 && isinteger(b) && return beta_negint(b, a)
-    #
+
     if fastabs(a) < fastabs(b)
         a, b = b, a
     end
-    #
-    if fastabs(a) > ASYMP_FACTOR * fastabs(b) && a > ASYMP_FACTOR
+
+    if fastabs(a) > ASYMP_FACTOR * fastabs(b) && real(a) > ASYMP_FACTOR
         # Avoid loss of precision in lgamma(a + b) - lgamma(a)
         y, s = lbeta_asymp(a, b)
         return s * exp(y)
@@ -35,10 +35,10 @@ function beta{T<:Number}(a::T, b::T)
 
     y = a + b
     if fastabs(y) > MAXGAM || fastabs(a) > MAXGAM || fastabs(b) > MAXGAM
-    	y, s = lgamma_r(y)
-    	yb, sb = lgamma_r(b)
+    	y, s::T = lgamma_r(y)
+    	yb, sb::T = lgamma_r(b)
         y = yb - y
-    	ya, sa = lgamma_r(a)
+    	ya, sa::T = lgamma_r(a)
         y = ya + y
     	# if (y > MAXLOG) {
     	#     goto overflow;
@@ -65,33 +65,34 @@ end
 lbeta(a::Number, b::Number) = lbeta(promote(float(a), float(b))...)
 
 function lbeta{T<:Number}(a::T, b::T)
-    real(a) <= 0.0 && isinteger(a) && return lbeta_negint(a, b)
-    real(b) <= 0.0 && isinteger(b) && return lbeta_negint(b, a)
+    # real(a) <= 0.0 && isinteger(a) && return lbeta_negint(a, b)
+    # real(b) <= 0.0 && isinteger(b) && return lbeta_negint(b, a)
 
     if fastabs(a) < fastabs(b)
         a, b = b, a
     end
 
-    if fastabs(a) > ASYMP_FACTOR * fastabs(b) && a > ASYMP_FACTOR
+    if fastabs(a) > ASYMP_FACTOR * fastabs(b) && real(a) > ASYMP_FACTOR
         # Avoid loss of precision in lgamma(a + b) - lgamma(a)
-        y, s = lbeta_asymp(a, b)
+        y,  = lbeta_asymp(a, b)
         return y
     end
 
     y = a + b
     if fastabs(y) > MAXGAM || fastabs(a) > MAXGAM || fastabs(b) > MAXGAM
-    	y, s = lgamma_r(y)
-    	yb, sb = lgamma_r(b)
+    	y, = lgamma_r(y)
+    	yb, = lgamma_r(b)
         y = yb - y
-    	ya, sa = lgamma_r(a)
+    	ya, = lgamma_r(a)
         y = ya + y
-    	return y
+    	return y #don't need the s' since the are 1 for complex (add test for this)
+                 #and for real we are computing log|B(a,b)|
     end
 
     y = gamma(y)
     a = gamma(a)
     b = gamma(b)
-    y == 0.0 && return inf(T)
+    y == 0.0 && return convert(T, Inf)
 
     if fastabs(fastabs(a) - fastabs(y)) > fastabs(fastabs(b) - fastabs(y))
         y = b / y
@@ -101,9 +102,12 @@ function lbeta{T<:Number}(a::T, b::T)
         y *= b
     end
 
-    return real(y) < 0 ? log(-y) : log(y)
+    return _log(y)
 end
 
+_log(y::Complex) = log(y)
+_log(y::Real) = log(abs(y))
+
 # assuming isinteger(x) and x < 0.
 function beta_negint{T}(x::T, w::T)
     if isinteger(w) && 1 - x - w > 0
@@ -118,7 +122,6 @@ end
 function lbeta_negint{T}(x::T, w::T)
     if isinteger(w) && 1 - x - w > 0
         s = ifelse(Int(w) % 2 == 0 , T(1), T(-1))
-        return s
         return s * lbeta(1 - x - w, w)
     else
         return convert(T, Inf)