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?
 A: 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.
