add logabsgamma (#72)

* add logabsgamma

* Update src/specfun.jl

Co-authored-by: RalphAS <[email protected]>

* handle Inf and NaN

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

* match sign also for negative integers

* remove deprecated lgamma (SpecialFunctions will take care of it)

---------

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

Author: 3 people

src/specfun.jl

@@ -1,6 +1,6 @@
 import .SpecialFunctions
 import .SpecialFunctions: erf, erfc, besselj0, besselj1, bessely0, bessely1,
-    besselj, bessely, gamma, lgamma
+    besselj, bessely, gamma, logabsgamma
 
 erf(x::Float128) =
     Float128(@ccall(libquadmath.erfq(x::Cfloat128)::Cfloat128))
@@ -24,5 +24,9 @@ bessely(n::Integer, x::Float128) =
 
 gamma(x::Float128) =
     Float128(@ccall(libquadmath.tgammaq(x::Cfloat128)::Cfloat128))
-lgamma(x::Float128) =
-    Float128(@ccall(libquadmath.lgammaq(x::Cfloat128)::Cfloat128))
+
+function logabsgamma(x::Float128)
+    result = Float128(@ccall(libquadmath.lgammaq(x::Cfloat128)::Cfloat128))
+    sign = !isfinite(result) || x >= 0 || !iszero(mod(ceil(x), 2)) ? 1 : -1
+    return result, sign
+end

test/specfun.jl

@@ -6,10 +6,13 @@ using SpecialFunctions
     piq = Float128(pi)
     halfq = Float128(0.5)
     @test gamma(halfq) ≈ sqrt(piq)
-    for func in (erf, erfc, besselj0, besselj1, bessely0, bessely1, lgamma)
+    for func in (erf, erfc, besselj0, besselj1, bessely0, bessely1, loggamma)
         @test func(halfq) ≈ func(0.5)
     end
     for func in (bessely, besselj)
         @test func(3,halfq) ≈ func(3,0.5)
     end
+    @test gamma(Float128(2),3) ≈ gamma(2,3)
+    @test all(logabsgamma(Float128(-0.5)) .≈ logabsgamma(-0.5))
+    @test all(logabsgamma(Float128(-1.5)) .≈ logabsgamma(-1.5))
 end