lgamma() fails for large negative imaginary arguments

[hsgg]

The lgamma() function in julia-0.4.6 has trouble with large negative imaginary arguments, e.g.

julia> lgamma(-300im)
-Inf + NaN*im

scipy.special.gammaln() gives a more reasonable answer:

>>> import scipy.special
>>> scipy.special.gammaln(-300.0j)
(-473.17185074259248-1410.3490664555823j)

[stevengj]

The problem is the log(sinpi(z)) call in gamma.jl, not the underlying clgamma_lanczos call. This code was written by @ViralBShah in #1466. by @JeffBezanson in 4675b10 ... Jeff, care to comment on the origins of this?

[stevengj]

For reference, here is how SciPy does it, albeit only for double precision.

SciPy's answer is correct except for the last digit — here is Mathematica's answer: -473.17185074259241355733179182866544204963885920016823743... - 1410.3490664555822107569308046418321236643870840962522425... i

[stevengj]

Note also that our lgamma function claims to be for Complex, but actually it is only double precision because of the const logpi = 1.14472988584940017. In general, we probably need to revisit some of these "ancient" special-function routines from 2011.

[stevengj]

The following implementation seems to avoid the numerous accuracy problems of our current lgamma, though it needs more testing (and checks for NaN/Inf arguments):

@inline function lgamma_asymptotic(z)
    zinv = inv(z)
    t = zinv*zinv
    # coefficients are bernoulli[2:n+1] .// (2*(1:n).*(2*(1:n) - 1))
    return (z-0.5)*log(z) - z + 9.1893853320467274178032927e-01 + zinv*@evalpoly(t, 8.3333333333333333333333368e-02,-2.7777777777777777777777776e-03,7.9365079365079365079365075e-04,-5.9523809523809523809523806e-04,8.4175084175084175084175104e-04,-1.9175269175269175269175262e-03,6.4102564102564102564102561e-03,-2.9550653594771241830065352e-02)
end

function my_lgamma(z::Complex{Float64})
    x = real(z)
    y = imag(z)
    yabs = abs(y)
    if x > 7 || yabs > 7 # use the Stirling asymptotic series for sufficiently large x or |y|
        return lgamma_asymptotic(z)
    end
    if x < 0.1 # use reflection formula to transform to x > 0
        if x == 0 && y == 0 # return Inf with the correct imaginary part for z == 0
            return Complex(Inf, signbit(x) ? copysign(oftype(x, pi), -y) : -y)
        end
        return Complex(1.1447298858494001741434262,
                       copysign(6.2831853071795864769252842, y)
                       * floor(0.5*x+0.25)) -
               log(sinpi(z)) - my_lgamma(1-z)
    end
    if abs(x - 1) + yabs < 0.1
        # taylor series around zero at z=1
        # ... coefficients are [-eulergamma; [(-1)^k * zeta(k)/k for k in 2:15]]
        w = Complex(x - 1, y)
        return w * @evalpoly(w, -5.7721566490153286060651188e-01,8.2246703342411321823620794e-01,-4.0068563438653142846657956e-01,2.705808084277845478790009e-01,-2.0738555102867398526627303e-01,1.6955717699740818995241986e-01,-1.4404989676884611811997107e-01,1.2550966952474304242233559e-01,-1.1133426586956469049087244e-01,1.000994575127818085337147e-01,-9.0954017145829042232609344e-02,8.3353840546109004024886499e-02,-7.6932516411352191472827157e-02,7.1432946295361336059232779e-02,-6.6668705882420468032903454e-02)
    end
    # use recurrence relation lgamma(z) = lgamma(z+1) - log(z) to shift to x > 7 for asymptotic series
    shiftprod = z
    x += 1
    while x <= 7
        shiftprod *= Complex(x,y)
        x += 1
    end
    return lgamma_asymptotic(Complex(x,y)) - log(shiftprod)
end
my_lgamma(z::Number) = my_lgamma(complex(float(z)))

It also seems to be about 3x faster than our current implementation (for zero-mean Gaussian random complex numbers with standard deviation 10). It's also about 40% faster (at 1/5 the lines of code) than the SciPy version (although I was hoping to beat SciPy by an even bigger margin).

[hsgg]

Thank you!