RFC: add complex polygamma and Hurwitz zeta functions

[stevengj]

This replaces the old polygamma(m,z) implementation for #7033. The new one seems to be around the same speed around 10-20x faster for real arguments z, but the same code also works for complex z, and it also handles negative z. I also added specialized digamma and trigamma implementations supporting real and complex arguments. The new trigamma is much faster, almost 10x more than 50x faster, than the old generic one based on polygamma. The new digamma seems to be about the same speed as the old digamma for real arguments, but is simpler and supports both real and complex z.

However, I still seem to be hitting #7060. My @chorner macro is supposed to fall back to @horner automatically for real arguments, but manually replacing @chorner with @horner for real arguments makes the digamma function 5x faster. cc: @carnaval

As a bonus, the way I implemented polygamma gives us the Hurwitz zeta function zeta(s,z) = ζ(s, z) for free, at least for real(s) ≥ 1. This is a common generalization of the Riemann zeta function zeta(s) = ζ(s) = ζ(s, 1), which we already export.

(In total, I deleted almost as much code as I added.)

[stevengj]

This version of digamma, which is exactly the same as the one in my patch except that it uses @horner instead of @chorner, is way faster for real arguments (and way slower for complex arguments) as mentioned above:

import Base.Math.@horner

function newdg(z::Union(Float64,Complex{Float64}))
    # Based on eq. (12), without looking at the accompanying source
    # code, of: K. S. Kölbig, "Programs for computing the logarithm of
    # the gamma function, and the digamma function, for complex
    # argument," Computer Phys. Commun.  vol. 4, pp. 221–226 (1972).
    x = real(z)
    if x <= 0 # reflection formula
        ψ = -π * cot(π*z)
        z = 1 - z
        x = real(z)
    else
        ψ = zero(z)
    end
    if x < 7
        # shift using recurrence formula
        n = 7 - ifloor(x)
        for ν = 1:n-1
            ψ -= inv(z + ν)
        end
        ψ -= inv(z)
        z += n
    end
    t = inv(z)
    ψ += log(z) - 0.5*t
    t *= t # 1/z^2
    # the coefficients here are float64(bernoulli[2:9] .// (2*(1:8)))
    ψ -= t * @horner(t,0.08333333333333333,-0.008333333333333333,0.003968253968253968,-0.004166666666666667,0.007575757575757576,-0.021092796092796094,0.08333333333333333,-0.4432598039215686)
end

Is there any way to fix @chorner (or Julia) so that it inlines to @horner for real z as intended?

[StefanKarpinski]

This is incredibly good work – I should expect no less, @stevengj, but still, you've outdone yourself. Are you comfortable merging this before addressing the constant propagation issue, or do you want to wait?

[stevengj]

It can be merged before that is fixed, since it only affects performance and not correctness. Will cause a performance regression in digamma, though.

[stevengj]

I've found a workaround for the type-inference problem (the problem was that Julia was performing type inference before dead-code elimination, so it thought the function wasn't type-stable), and will update the patch shortly.

[stevengj]

Compared to the implementations in SciPy (for large arrays so that I am just benchmarking their underlying C/Fortran code):

• digamma is almost exactly the same speed for real arguments, and about twice as fast for complex arguments.
• polygamma is about twice as fast (sometimes more, sometimes less) for real arguments. Their implementation does not support complex arguments.
• They don't have a specialized trigamma function. Compared to their generic polygamma(1, x) for real x, our trigamma(x) is 10x faster.

Also, it looks like their polygamma function is extremely inefficient for large negative arguments x, requiring time proportional to abs(x).

[stevengj]

Compared to Matlab (with -singleCompThread), our digamma function is about 3x faster for real arguments, and polygamma is about the same speed.

However, Matlab's psi function (digamma and polygamma) does not support complex or even negative arguments.

[stevengj]

cc @andreasnoackjensen

[andreasnoack]

Looks great.

[StefanKarpinski]

Phenomenal.

[StefanKarpinski]

Please merge whenever you like.

[StefanKarpinski]

Eh, let's leave it as is then. It can always be deprecated in the future.

[0joshuaolson1]

If I understand correctly, the current polygamma handles complex types now? It was mentioned in #7033 that there are several ways to extend it. Is there a relevant catch-all issue about linking function docs to the source code?

[stevengj]

@mrmormon, yes, the current polygamma(m, z) handles complex z. The docs should probably be clarified on this point.

Function documentation was overhauled in Julia 0.4; see the HISTORY.md file and #8791.