Pure Julia implementation of OpenLibm's lgamma, lgamma_r

src/SpecialFunctions.jl

@@ -85,6 +85,8 @@ include("erf.jl")
 include("ellip.jl")
 include("expint.jl")
 include("sincosint.jl")
+include("logabsgamma/e_lgammaf_r.jl")
+include("logabsgamma/e_lgamma_r.jl")
 include("gamma.jl")
 include("gamma_inc.jl")
 include("betanc.jl")

src/gamma.jl

@@ -598,16 +598,6 @@ See also [`loggamma`](@ref).
 """
 logabsgamma(x::Real) = _logabsgamma(float(x))
 
-function _logabsgamma(x::Float64)
-    signp = Ref{Int32}()
-    y = ccall((:lgamma_r,libopenlibm),  Float64, (Float64, Ptr{Int32}), x, signp)
-    return y, Int(signp[])
-end
-function _logabsgamma(x::Float32)
-    signp = Ref{Int32}()
-    y = ccall((:lgammaf_r,libopenlibm),  Float32, (Float32, Ptr{Int32}), x, signp)
-    return y, Int(signp[])
-end
 function _logabsgamma(x::Float16)
     y, s = _logabsgamma(Float32(x))
     return Float16(y), s
@@ -649,6 +639,7 @@ function _loggamma(x::Real)
     return y
 end
 
+
 function _loggamma(x::BigFloat)
     isnan(x) && return x
     y = BigFloat()

src/logabsgamma/e_lgamma_r.jl

@@ -0,0 +1,176 @@
+#=
+/* @(#)e_lgamma_r.c 1.3 95/01/18 */
+/*
+* ====================================================
+* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+*
+* Developed at SunSoft, a Sun Microsystems, Inc. business.
+* Permission to use, copy, modify, and distribute this
+* software is freely granted, provided that this notice
+* is preserved.
+* ====================================================
+*
+*/
+=#
+
+#=
+
+/* __ieee754_lgamma_r(x, signgamp)
+ * Reentrant version of the logarithm of the Gamma function
+ * with user provide pointer for the sign of Gamma(x).
+ *
+ * Method:
+ *   1. Argument Reduction for 0 < x <= 8
+ * 	Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
+ * 	reduce x to a number in [1.5,2.5] by
+ * 		lgamma(1+s) = log(s) + lgamma(s)
+ *	for example,
+ *		lgamma(7.3) = log(6.3) + lgamma(6.3)
+ *			    = log(6.3*5.3) + lgamma(5.3)
+ *			    = log(6.3*5.3*4.3*3.3*2.3) + lgamma(2.3)
+ *   2. Polynomial approximation of lgamma around its
+ *	minimun ymin=1.461632144968362245 to maintain monotonicity.
+ *	On [ymin-0.23, ymin+0.27] (i.e., [1.23164,1.73163]), use
+ *		Let z = x-ymin;
+ *		lgamma(x) = -1.214862905358496078218 + z^2*poly(z)
+ *	where
+ *		poly(z) is a 14 degree polynomial.
+ *   2. Rational approximation in the primary interval [2,3]
+ *	We use the following approximation:
+ *		s = x-2.0;
+ *		lgamma(x) = 0.5*s + s*P(s)/Q(s)
+ *	with accuracy
+ *		|P/Q - (lgamma(x)-0.5s)| < 2**-61.71
+ *	Our algorithms are based on the following observation
+ *
+ *                             zeta(2)-1    2    zeta(3)-1    3
+ * lgamma(2+s) = s*(1-Euler) + --------- * s  -  --------- * s  + ...
+ *                                 2                 3
+ *
+ *	where Euler = 0.5771... is the Euler constant, which is very
+ *	close to 0.5.
+ *
+ *   3. For x>=8, we have
+ *	lgamma(x)~(x-0.5)log(x)-x+0.5*log(2pi)+1/(12x)-1/(360x**3)+....
+ *	(better formula:
+ *	   lgamma(x)~(x-0.5)*(log(x)-1)-.5*(log(2pi)-1) + ...)
+ *	Let z = 1/x, then we approximation
+ *		f(z) = lgamma(x) - (x-0.5)(log(x)-1)
+ *	by
+ *	  			    3       5             11
+ *		w = w0 + w1*z + w2*z  + w3*z  + ... + w6*z
+ *	where
+ *		|w - f(z)| < 2**-58.74
+ *
+ *   4. For negative x, since (G is gamma function)
+ *		-x*G(-x)*G(x) = pi/sin(pi*x),
+ * 	we have
+ * 		G(x) = pi/(sin(pi*x)*(-x)*G(-x))
+ *	since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0
+ *	Hence, for x<0, signgam = sign(sin(pi*x)) and
+ *		lgamma(x) = log(|Gamma(x)|)
+ *			  = log(pi/(|x*sin(pi*x)|)) - lgamma(-x);
+ *	Note: one should avoid compute pi*(-x) directly in the
+ *	      computation of sin(pi*(-x)).
+ *
+ *   5. Special Cases
+ *		lgamma(2+s) ~ s*(1-Euler) for tiny s
+ *		lgamma(1) = lgamma(2) = 0
+ *		lgamma(x) ~ -log(|x|) for tiny x
+ *		lgamma(0) = lgamma(neg.integer) = inf and raise divide-by-zero
+ *		lgamma(inf) = inf
+ *		lgamma(-inf) = inf (bug for bug compatible with C99!?)
+ *
+ */
+
+=#
+
+# Matches OpenLibm behavior (except commented out |x|≥2^52 early exit)
+function _logabsgamma(x::Float64)
+    ux = reinterpret(UInt64, x)
+    hx = ux >>> 32 % Int32
+    lx = ux % UInt32
+
+    #= purge off +-inf, NaN, +-0, tiny and negative arguments =#
+    signgam = 1
+    isneg = hx < Int32(0)
+    ix = hx & 0x7fffffff
+    ix ≥ 0x7ff00000 && return x * x, signgam
+    ix | lx == 0x00000000 && return Inf, signgam
+    if ix < 0x3b900000 #= |x|<2**-70, return -log(|x|) =#
+        if isneg
+            signgam = -1
+            return -log(-x), signgam
+        else
+            return -log(x), signgam
+        end
+    end
+    if isneg
+        # ix ≥ 0x43300000 && return Inf, signgam #= |x|>=2**52, must be -integer =#
+        t = sinpi(x)
+        iszero(t) && return Inf, signgam #= -integer =#
+        nadj = logπ - log(abs(t * x))
+        if t < 0.0; signgam = -1; end
+        x = -x
+    end
+    if ix < 0x40000000     #= x < 2.0 =#
+        i = round(x, RoundToZero)
+        f = x - i
+        if f == 0.0 #= purge off 1; 2 handled by x < 8.0 branch =#
+            return 0.0, signgam
+        elseif i == 1.0
+            r = 0.0
+            c = 1.0
+        else
+            r = -log(x)
+            c = 0.0
+        end
+        if f ≥ 0.7315998077392578
+            y = 1.0 + c - x
+            z = y * y
+            p1 = @evalpoly(z, 7.72156649015328655494e-02, 6.73523010531292681824e-02, 7.38555086081402883957e-03, 1.19270763183362067845e-03, 2.20862790713908385557e-04, 2.52144565451257326939e-05)
+            p2 = z * @evalpoly(z, 3.22467033424113591611e-01, 2.05808084325167332806e-02, 2.89051383673415629091e-03, 5.10069792153511336608e-04, 1.08011567247583939954e-04, 4.48640949618915160150e-05)
+            p = muladd(p1, y, p2)
+            r += muladd(y, -0.5, p)
+        elseif f ≥ 0.2316399812698364 # or, the lb? 0.2316322326660156
+            y = x - 0.46163214496836225 - c
+            z = y * y
+            w = z * y
+            p1 = @evalpoly(w, 4.83836122723810047042e-01, -3.27885410759859649565e-02, 6.10053870246291332635e-03, -1.40346469989232843813e-03, 3.15632070903625950361e-04)
+            p2 = @evalpoly(w, -1.47587722994593911752e-01, 1.79706750811820387126e-02, -3.68452016781138256760e-03, 8.81081882437654011382e-04, -3.12754168375120860518e-04)
+            p3 = @evalpoly(w, 6.46249402391333854778e-02, -1.03142241298341437450e-02, 2.25964780900612472250e-03, -5.38595305356740546715e-04, 3.35529192635519073543e-04)
+            p = muladd(z, p1, -muladd(w, -muladd(p3, y, p2), -3.63867699703950536541e-18))
+            r += p - 1.21486290535849611461e-1
+        else
+            y = x - c
+            p1 = y * @evalpoly(y, -7.72156649015328655494e-02, 6.32827064025093366517e-01, 1.45492250137234768737, 9.77717527963372745603e-01, 2.28963728064692451092e-01, 1.33810918536787660377e-02)
+            p2 = @evalpoly(y, 1.0, 2.45597793713041134822, 2.12848976379893395361, 7.69285150456672783825e-01, 1.04222645593369134254e-01, 3.21709242282423911810e-03)
+		    r += muladd(y, -0.5, p1 / p2)
+        end
+    elseif ix < 0x40200000              #= x < 8.0 =#
+        i = round(x, RoundToZero)
+        y = x - i
+	    z = 1.0
+        p = 0.0
+        u = x
+        while u ≥ 3.0
+            p -= 1.0
+            u = x + p
+            z *= u
+        end
+        p = y * @evalpoly(y, -7.72156649015328655494e-2, 2.14982415960608852501e-1, 3.25778796408930981787e-1, 1.46350472652464452805e-1, 2.66422703033638609560e-2, 1.84028451407337715652e-3, 3.19475326584100867617e-5)
+        q = @evalpoly(y, 1.0, 1.39200533467621045958, 7.21935547567138069525e-1, 1.71933865632803078993e-1, 1.86459191715652901344e-2, 7.77942496381893596434e-4, 7.32668430744625636189e-6)
+        r = log(z) + muladd(0.5, y, p / q)
+    elseif ix < 0x43900000              #= 8.0 ≤ x < 2^58 =#
+        z = 1.0 / x
+        y = z * z
+        w = muladd(z, @evalpoly(y, 8.33333333333329678849e-2, -2.77777777728775536470e-3, 7.93650558643019558500e-4, -5.95187557450339963135e-4, 8.36339918996282139126e-4, -1.63092934096575273989e-3), 4.18938533204672725052e-1)
+	    r = muladd(x - 0.5, log(x) - 1.0, w)
+    else #= 2^58 ≤ x ≤ Inf =#
+        r = muladd(x, log(x), -x)
+    end
+    if isneg
+        r = nadj - r
+    end
+    return r, signgam
+end

src/logabsgamma/e_lgammaf_r.jl

@@ -0,0 +1,106 @@
+#=
+/* e_lgammaf_r.c -- float version of e_lgamma_r.c.
+* Conversion to float by Ian Lance Taylor, Cygnus Support, [email protected].
+*/
+
+/*
+* ====================================================
+* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+*
+* Developed at SunPro, a Sun Microsystems, Inc. business.
+* Permission to use, copy, modify, and distribute this
+* software is freely granted, provided that this notice
+* is preserved.
+* ====================================================
+*/
+=#
+
+# Matches OpenLibm behavior exactly, including return of sign
+function _logabsgamma(x::Float32)
+    hx = reinterpret(Int32, x)
+
+    #= purge off +-inf, NaN, +-0, tiny and negative arguments =#
+    signgam = 1
+    isneg = hx < Int32(0)
+    ix = hx & 0x7fffffff
+    ix ≥ 0x7f800000 && return x * x, signgam
+    ix == 0x00000000 && return Inf32, signgam
+    if ix < 0x35000000 #= |x|<2**-21, return -log(|x|) =#
+        if isneg
+            signgam = -1
+            return -log(-x), signgam
+        else
+            return -log(x), signgam
+        end
+    end
+    if isneg
+        # ix ≥ 0x4b000000 && return Inf32, signgam #= |x|>=2**23, must be -integer =#
+        t = sinpi(x)
+        t == 0.0f0 && return Inf32, signgam #= -integer =#
+        nadj = logπ - log(abs(t * x))
+        if t < 0.0f0; signgam = -1; end
+        x = -x
+    end
+
+    if ix < 0x40000000 #= x < 2.0 =#
+        i = round(x, RoundToZero)
+        f = x - i
+        if f == 0.0f0 #= purge off 1; 2 handled by x < 8.0 branch =#
+            return 0.0f0, signgam
+        elseif i == 1.0f0
+            r = 0.0f0
+            c = 1.0f0
+        else
+            r = -log(x)
+            c = 0.0f0
+        end
+        if f ≥ 0.7315998f0
+            y = 1.0f0 + c - x
+            z = y * y
+            p1 = @evalpoly(z, 7.7215664089f-2, 6.7352302372f-2, 7.3855509982f-3, 1.1927076848f-3, 2.2086278477f-4, 2.5214456400f-5)
+            p2 = z * @evalpoly(z, 3.2246702909f-1, 2.0580807701f-2, 2.8905137442f-3, 5.1006977446f-4, 1.0801156895f-4, 4.4864096708f-5)
+            p = muladd(p1, y, p2)
+            r += muladd(y, -0.5f0, p)
+        elseif f ≥ 0.23163998f0 # or, the lb? 0.2316322f0
+            y = x - 0.46163213f0 - c
+            z = y * y
+            w = z * y
+            p1 = @evalpoly(w, 4.8383611441f-1, -3.2788541168f-2, 6.1005386524f-3, -1.4034647029f-3, 3.1563205994f-4)
+            p2 = @evalpoly(w, -1.4758771658f-1, 1.7970675603f-2, -3.6845202558f-3, 8.8108185446f-4, -3.1275415677f-4)
+            p3 = @evalpoly(w, 6.4624942839f-2, -1.0314224288f-2, 2.2596477065f-3, -5.3859531181f-4, 3.3552918467f-4)
+            p = muladd(z, p1, -muladd(w, -muladd(p3, y, p2), 6.6971006518f-9))
+            r += p - 1.2148628384f-1
+        else
+            y = x - c
+            p1 = y * @evalpoly(y, -7.7215664089f-2, 6.3282704353f-1, 1.4549225569f0, 9.7771751881f-1, 2.2896373272f-1, 1.3381091878f-2)
+            p2 = @evalpoly(y, 1.0f0, 2.4559779167f0, 2.1284897327f0, 7.6928514242f-1, 1.0422264785f-1, 3.2170924824f-3)
+            r += muladd(y, -0.5f0, p1 / p2)
+        end
+    elseif ix < 0x41000000 #= x < 8.0 =#
+        i = round(x, RoundToZero)
+        y = x - i
+        z = 1.0f0
+        p = 0.0f0
+        u = x
+        while u ≥ 3.0f0
+            p -= 1.0f0
+            u = x + p
+            z *= u
+        end
+        p = y * @evalpoly(y, -7.7215664089f-2, 2.1498242021f-1, 3.2577878237f-1, 1.4635047317f-1, 2.6642270386f-2, 1.8402845599f-3, 3.1947532989f-5)
+        q = @evalpoly(y, 1.0f0, 1.3920053244f0, 7.2193557024f-1, 1.7193385959f-1, 1.8645919859f-2, 7.7794247773f-4, 7.3266842264f-6)
+	    r = log(z) + muladd(0.5f0, y, p / q)
+    elseif ix < 0x5c800000 #= 8.0 ≤ x < 2^58 =#
+        z = 1.0f0 / x
+        y = z * z
+        w = muladd(z, @evalpoly(y, 8.3333335817f-2, -2.7777778450f-3, 7.9365057172f-4, -5.9518753551f-4, 8.3633989561f-4, -1.6309292987f-3), 4.1893854737f-1)
+        r = muladd(x - 0.5f0, log(x) - 1.0f0, w)
+    else
+        #= 2^58 ≤ x ≤ Inf =#
+        r = muladd(x, log(x), -x)
+    end
+    if isneg
+        r = nadj - r
+    end
+    return r, signgam
+end