# Numerical evaluation of Hurwitz zeta function

Is there a way to evaluate numerically the Hurwitz zeta function

$$\zeta(s,a) = \sum_{n=0}^\infty \frac{1}{(n+a)^{s}}$$

that is more efficient (i.e., quick and precise) than simply explicitly adding the terms one by one?

Note that your variables names are uncommon, I will use the standard names from DLMF $$\zeta(s,a) = \sum\limits_{n=0}^{\infty}\frac{1}{(n+a)^s}$$ A good numerical algorithm uses the Euler-Maclaurin summation (DLMF or Wolfram) in the form $$\zeta(s,a) = \sum\limits_{k=0}^{n}\frac{1}{(a+k)^s} + \frac{(a+n)^{1-s}}{s-1} + \sum\limits_{k=0}^{\infty} \frac{B_{k+1}}{(k+1)!}\: \frac{(s+k-1)_k}{(a+n)^{s+k}}$$ where $(\cdot)_k$ is the Pochhammer symbol. Using the properties of the Bernoulli numbers $B_1=-1/2\;$ and $B_{2k+1}=0\;$ this can be rewritten as $$\zeta(s,a) = \sum\limits_{k=0}^{n}\frac{1}{(a+k)^s} + \frac{(a+n)^{1-s}}{s-1} - \frac{1}{2(a+n)^{s}} + \sum\limits_{k=1}^{\infty} \frac{B_{2k}}{(2k)!}\: \frac{(s+2k-2)_{2k-1}}{(a+n)^{s+2k-1}} \cdot$$ The upper bound of the first sum should be selected depending on your arithmetic (e.g $n=8$ or $n=9$ for double etc). The second (infinite) sum can computed by summing the terms until there is no change (i.e. the current term is smaller in magnitude than current sum).

This algorithm is efficient for the primary range $s>0, a>0:\;$ For 80-bit floating point arithmetic I get $\zeta(2.5,0.75)\approx 2.49154238551193522\;$ with $10$ terms of the second sum, and $\zeta(10,0.25)\approx 1048576.10768311475\;$ you only need two terms.

For negative $s$ you can use Bernoulli polynomials (for $s=-1,-2\dots$) or the reflection formula from DLMF.

Edit: The number of terms depends mainly on $s$ and $n$. Here some intermediate values from my implementation for 53-bit IEEE double. For the given $(s,a)\;$ pairs $n=9.\;$ s1 is the first sum from $k=0\dots n,\;$ s2 = s1 + the next two single terms, and s2 the total sum after finishing the last sum.

ZetaH(10, 0.25)
s1 = 1048576.10768311471
s2 = 1048576.10768311494
s3 = 1048576.10768311494       1 term

ZetaH(5.0, 0.25)
s1 = 1024.34894710160870
s2 = 1024.34897386673833
s3 = 1024.34897452658106       6 terms

ZetaH(2.5, 0.75)
s1 = 2.47125718157242114
s2 = 2.49147059743736321
s3 = 2.49154238551193474       8 terms

ZetaH(0.75, 0.25)
s1 =  6.13122575577889695
s2 = -0.938864666031481665
s3 = -0.937591964187151961     9 terms


In actual implementations, you can adjust the $n$ value, and of course there are other speed-up methods: Terminate the first sum, if there are no changes, pre-compute a certain numbers of the complicated coefficients of the second sum etc.

You find a basic C implementation in the file zeta.c of the well known Cephes library http://www.moshier.net/double.zip

Edit 2: Here the data for the two sums. The first column shows n (the upper limit of the first sum), j is the index of the second sum, for which the convergence criterion is achieved. The second sum is also terminated, if the terms are non-decreasing, this is the case for n=4,5 and the last three s,a pairs, in the table it is shown as ---.

        s,a=     s,a=     s,a=       s,a=
10,0.25  5,0.25  2.5,0.75   0.75,0.25
n=4      j=7      ---      ---       ---
n=5      j=3      ---      ---       ---
n=6      j=1      j=8      j=11      j=12
n=7      j=1      j=6      j=9       j=10
n=8      j=0      j=5      j=8       j=8
n=9      j=0      j=5      j=7       j=8
n=10              j=4      j=7       j=7
n=11              j=4      j=6       j=7
n=12              j=4      j=6       j=6
n=13              j=3      j=6       j=6
n=14              j=3      j=5       j=6
n=15              j=3      j=5       j=6
n=16              j=3      j=5       j=5
n=17              j=3      j=5       j=5
n=18              j=3      j=5       j=5
n=19              j=2      j=4       j=5
n=20              j=2      j=4       j=5
n=21              j=2      j=4       j=5
n=22              j=2      j=4       j=5
n=23              j=2      j=4       j=5
n=24              j=2      j=4       j=5
n=25              j=2      j=4       j=4
n=26              j=2      j=4       j=4

• Thanks @gammatester for the quick and detailed answer. I have implemented the algorithm in C++ using double precision. In my case, it seems that the precision of the answer I get depends mostly on the value of $n$ in the first sum (which corresponds basically to a direct evaluation of the Hurwitz zeta function). The contribution of the second sum falls rapidly to zero and seems to never be able to compensate for the truncation of the first sum. Any idea on what I may be doing wrong? Sep 3, 2014 at 14:20
• @allard: I just edited the answer to give some intermediate values (too long for a comment). What value of $n$ do you use? Sep 4, 2014 at 7:01
• Just a question: Does this formula give correct results for complex numbers too? If so, in what range? Aug 13, 2015 at 7:43
• @Anonymous Pi: DLMF lists the ranges $a>0$ and $\Re s >-1$ or $\Re s > -2n$. I would interpret this as $a$ should be a positive real. The deeper question may be the error estimate for complex $s$ (there are no error bounds in DLMF, but the situation may be similar to the Gamma function dlmf.nist.gov/5.11.ii) Aug 13, 2015 at 9:50
• @gammatester: Can you give a hint concerning the connection between the choice of $n$ and the involved arithmetic ? Dec 15, 2017 at 12:56

The first formula given by @gammatester agrees with that given on the Wolfram site, but the second, rewritten, form seems to be in error.

The first, simple, error is that the exponent for the $(a+n)$ term in the denominator of Bernoulli series is wrong by 2 - should be $s+2k-1$, but a more involved error is in interpreting the Pochhammer symbol when summing over the even indices, $2k$. As given, in the rewritten formula it amounts to $(s)2k+1$. This is incorrect: consider first the Pochhammer form in the original form as $(s + k-1)k$. This means the first factor in the Pochhammer symbol advances by one each time, and the last one by 2. In the even sum, it should therefore be $(s + 2k-2)2k-1$. The first term is $(s)1 =s$, the second is $(s+2)(s+3)(s+4)$, etc. The last factor is $s+4k-4$.

Numerically the quickest way to evaluate in the loop is to increment $k$ by one, two times. $(s+k-1)k$, is advanced one step by dividing by the first factor then adding two more.

I have implemented the correct algorithm in Fortran and Mathematica (to higher than double precision), and am in agreement with @allard that the error term after the first explicit summation - i.e. with the series including Bernoulli numbers - does not decrease as expected. The error saturates long before the series terms begin to diverge, and is not a precision issue.

• Your answer/comment is correct. Unfortunately I saw it only by accident, the @user notification works only from comments, and I guess your reputation was not sufficient to comment on my answer. The wrong formula is from the Cephes library documentation, AFAIK the actual code is correct. I updated the formula. Dec 12, 2017 at 10:58