ζ function works incorrectly on some complex inputs

[rominf]

For example,

julia> ζ(0 + 99.69im)
2.685576248606084 - 21.972880545101557im

But this answer is not correct.

python> sympy.zeta(0 + 99.69j).n()
4.67192766128947 + 3.89448062985268*I

and

wolframalpha> zeta(0 + 99.69i)
4.67193... + 3.89448... i

[staticfloat]

This seems to be because the eta function is giving us wonky answers.

[stevengj]

Matlab gives:

>> zeta(99.69i)

ans =
   4.6719 + 3.8945i

so it may be possible to compute this accurately without arbitrary precision. (This was not entirely clear to me since the zeta function is highly oscillatory along the imaginary axis.)

(Scipy's zeta function does not support complex arguments.)

[stevengj]

@ViralBShah, where did this eta implementation come from? It looks like the eta implementation might be using Borwein's method as described in Wikipedia, which has an error bound that gets exponentially worse with abs(imag(z)) and hence is not a good idea far from the real axis.

[JeffBezanson]

Yes I think you're right. That's an implementation I threw together several
years ago, originally in C then I ported it. We should use something better.
On Jun 10, 2014 5:22 PM, "Steven G. Johnson" notifications@github.com
wrote:

@ViralBShah https://github.com/ViralBShah, where did this eta
implementation come from? It looks like the eta implementation might be
using
https://en.wikipedia.org/wiki/Dirichlet_eta_function#Borwein.27s_method
http://Borwein's%20method, which has an error bound that gets
exponentially worse with abs(imag(z)) and hence is not a good idea far
from the real axis.

—
Reply to this email directly or view it on GitHub
#7169 (comment).

[stevengj]

The following seems to do a better job. It is based on my Hurwitz ζ function, which itself is based on differentiating the Stirling asymptotic series for the digamma function, specialized for the case of z==1. It gives 4.671927662426274 + 3.894480629951745im, which is at least reasonably accurate (10 significant digits).

import Base.Math: @pg_horner, inv_oftype, zeta_returntype

function myzeta(s::Union(Float64,Complex{Float64}))
    if real(s) < 0.5
        return myzeta(1 - s) * gamma(1 - s) * sinpi(s*0.5) * (2π)^s / π
    end

    m = s - 1

    # shift using recurrence formula                                        
    n = iceil(6 + real(m) + abs(imag(m))) # FIXME: work out more accurate imag(m) dependence
    ζ = one(s)
    for ν = 2:n
        ζₒ= ζ
        ζ += inv_oftype(ζ, ν)^s
        ζ == ζₒ && break # prevent long loop for large m 
    end
    z = 1 + n
    t = inv(z)
    w = isa(t, Real) ? conj(oftype(ζ, t))^m : oftype(ζ, t)^m
    ζ += w * (inv(m) + 0.5*t)

    t *= t # 1/z^2                                                              
    ζ += w*t * @pg_horner(t,m,0.08333333333333333,-0.008333333333333333,0.003968253968253968,-0.004166666666666667,0.007575757575757576,-0.021092796092796094,0.08333333333333333,-0.4432598039215686,3.0539543302701198)

    return ζ
end

The computation of the shift n needs work (basically it should be an empirical fit, but I haven't had time to wok out the dependence on imag(s) carefully). And I've only eyeballed the accuracy at a few points; I haven't done any careful checks yet.

As a bonus, however, it seems to be nearly 10x faster than our current zeta function for both real and complex arguments. (And of course we can compute eta from this zeta.)

[ViralBShah]

We will have to look at the git history. I did not write the eta implementation.

[ViralBShah]

Just saw Jeff's comment that he wrote it.

[stevengj]

It might also be nice to have a zetac(z) = zeta(z) - 1 function, similar to Scipy. This is potentially useful, much like erfc, because the zeta function approaches 1 for large real(z) > 0.

It could use almost the same implementation, but the cutoff for the recurrence would need to change.

(The Scipy zetac function does not accept complex arguments. I'm guessing this allows them to use much simpler rational approximants, and is probably why they are about 5x faster than our zeta right now for real arguments.)

[StefanKarpinski]

Go for it. Seems like a good thing to have.

[stevengj]

Although now that I look at it further, I'm less convinced that zetac makes sense. Yes, zeta(z)-1 loses all of its significant digits for z ≥ 55 or so. But zetac(z) loses all of its significant digits for z > 1080 or so, because of underflow; worse, it becomes trivially equal to 0.5^z for z > 90 or so. You could, of course, define a zetacx(z) function that multiplies by 2^z in order to eliminate the underflow, but it still effectively loses all useful digits for z > 90 or so: it just becomes 1.

The basic problem is that zetac(z) is simply a sum of exponentially decaying terms, each term decaying exponentially faster than the previous one (this is obvious by inspecting the definition). It makes no difference how many terms you subtract off; it becomes "boring" exponentially quickly.

Because of this, I'm inclined to wait until there is evidence of a real-world application for zetac before including it.