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 .)
 A: $$\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.
A: 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. 
A: 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.
