[WIP] port of Cephes implementation of Beta function (fixes #14256)

base/special/beta.jl

@@ -0,0 +1,142 @@
+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))
+
+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) && real(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::T = lgamma_r(y)
+    	yb, sb::T = lgamma_r(b)
+        y = yb - y
+    	ya, sa::T = 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 convert(T, 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
+
+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)
+
+    if fastabs(a) < fastabs(b)
+        a, b = b, a
+    end
+
+    if fastabs(a) > ASYMP_FACTOR * fastabs(b) && real(a) > ASYMP_FACTOR
+        # Avoid loss of precision in lgamma(a + b) - lgamma(a)
+        y,  = lbeta_asymp(a, b)
+        return y
+    end
+
+    y = a + b
+    if fastabs(y) > MAXGAM || fastabs(a) > MAXGAM || fastabs(b) > MAXGAM
+    	y, = lgamma_r(y)
+    	yb, = lgamma_r(b)
+        y = yb - y
+    	ya, = lgamma_r(a)
+        y = ya + 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 convert(T, 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 _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
+        s = ifelse(Int(w) % 2 == 0, T(1), T(-1))
+        return s * beta(1 - x - w, w)
+    else
+        return convert(T, Inf)
+    end
+end
+
+# assuming isinteger(x) and x < 0.
+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 * lbeta(1 - x - w, w)
+    else
+        return convert(T, Inf)
+    end
+end
+
+ # Asymptotic expansion for  ln(|B(a, b)|) for a > ASYMP_FACTOR*max(|b|, 1).
+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, T(s)
+end
+
+@vectorize_2arg Number beta
+@vectorize_2arg Number lbeta

base/special/gamma.jl

@@ -404,16 +404,6 @@ 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)
-end
-lbeta(x::Number, w::Number) = lgamma(x)+lgamma(w)-lgamma(x+w)
-@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}})
@@ -484,3 +474,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)) )

test/math.jl

@@ -459,6 +459,12 @@ 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
+
 # gamma, lgamma (complex argument)
 if Base.Math.libm == "libopenlibm"
     @test gamma(Float64[1:25;]) == gamma(1:25)
@@ -728,3 +734,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