Evaluates the natural logarithm of the absolute value of gamma function.

The following was formatted as Nim-Favor Markdown from FreeBSD `lgamma` source with some minor amendment.

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. A polynomial approximation of lgamma

Create a 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.

3. 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.

4. For x>=8, ...

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

5. For negative x, ...

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);

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  # see below

• some implementations, like CPython's math, R

and C/C++ returns +inf; which is not suitable, as gamma(x) where x < about -200 is always truncated to 0.0 at ieee754 float domain. This behavior was said for bug compatible with C99, and is readlly documented by POSIX man and cppreference.com

• While others like scipy.special.gammaln, Go's math.Lgamma, returns

In my option, ln(|gamma(-oo)|) -[ieee754 float trunc]-> ln(0+) -> -inf

But in this function it returns +inf to keep consistent with Python,

Example:

from std/math import isNaN
assert lgamma(1.0) == 0.0
assert lgamma(Inf) == Inf
assert lgamma(NaN).isNaN