# Numerical evaluation of derivatives of Hurwitz zeta function

I am looking a an efficient algorithm to compute in C derivatives of the Hurwitz zeta function $$\zeta^{(n)}(s,a) = (-1)^{n} \sum_{k=0}^\infty \frac{\log^{n}(k+a)}{(k+a)^{s}} .$$ for $1<s$ and $0<q$.

(This question is related to the previous question Numerical evaluation of Hurwitz zeta function .)

$$\sum_{k\geq 0}\frac{1}{(k+a)^s} \stackrel{\mathcal{L}^{-1}}{=}\frac{1}{\Gamma(s)}\int_{0}^{+\infty}\frac{e^{(1-a)u}u^{s-1}}{e^u-1}\,du$$ so $$(-1)^n \sum_{k\geq 0}\frac{\log(k+a)^n}{(k+a)^s} = \frac{1}{\Gamma(s)}\int_{0}^{+\infty}\frac{e^{(1-a)u}u^{s-1}\log(u)^n}{e^u-1}\,du$$ and any accurate algorithm for numerical integration can do the job.

• did I get your answer wrong or is it really missing something? See my answer. Feb 6, 2019 at 15:38
• It seems to me that this answer is wrong, because we have to take multiple derivatives of a product $u^s/ \Gamma(s)$ and not just $u^s$. Could you please check? Jul 15, 2019 at 21:21

I tried commenting on @JackD'Aurizio's answer, but I'm three reputation points shy of being able to. While the idea of using the integral representation is great, I believe the answer given by @JackD'Aurizio is in error, at least for the first derivative.

In what follows I will provide what I believe to be the correct expression for the first derivative. I believe the correct expression is $$\zeta^{(1)}(s, a) = \frac{1}{\Gamma(s)}\int_0^\infty \frac{e^{(1-a)x}x^{s-1}\left(\log(x) - \psi_0(s) \right)}{e^x -1}dx.$$ where $$\psi_0$$ is the polygamma function. I obtained this by differentiating $$\zeta(s) = \frac{1}{\Gamma(s)}\int_{0}^\infty \frac{e^{(1-a)x}x^{s-1}}{e^x -1}dx$$ with respect to $$s$$ and then applying the product and reciprocal rules to get $$\zeta^{(1)}(s, a) = -\frac{\psi_0(s)}{\Gamma(s)}\int_0^\infty \frac{e^{(1-a)x}x^{s-1}}{e^x -1}dx + \frac{1}{\Gamma(s)}\int_0^\infty \frac{e^{(1-a)x}x^{s-1}\log(x)}{e^x -1}dx,$$ from which the result follows. Accompanying R code (for $$a = 1$$):

zeta_int <- function(s){
func <- function(x){
x^(s-1)/(exp(x) - 1)
}
func <- Vectorize(func)
return(integrate(func, 0, Inf)$value/gamma(s)) } zeta_deriv <- function(s){ k <- digamma(s) func <- function(x){ (x^(s-1) *(log(x) - k))/(exp(x) -1) } func <- Vectorize(func) ans <- integrate(func, 0, Inf)$value/gamma(s)
return(ans)
}
####
# Testing
ss <- 3/2
zeta_int(ss)
VGAM::zeta(ss)
zeta_deriv_alt(ss)
VGAM::zeta(ss, deriv = 1)


I know this does not answer the question as posed, but I think it's a valid contribution nonetheless.

If you can handle working in C++ rather than C, use the Chebyshev transform in Boost.Math. You will need a way of evaluate the Hurwitz zeta function to compute its value at the Chebyshev nodes, but once you do, you'll get a stable method of numerical differentiation via differentiating the associated Chebyshev series.