For some exercises with (divergent) summation of the Stieltjes constants I'm trying a formula, which involves derivatives of the $\zeta()$ -function at negative integers; perhaps better formulated as $$ f(x) = \zeta(1-x) + 1/x \qquad x \in \mathbb N $$ and then the scaled derivatives $$ g(x,d) = x^d \cdot {f^{(d)}(x)\over d!} $$.
In Pari/GP there is a very nice procedure for the computation of $\zeta$ implemented, but I don't see anything for the derivatives except to implement numerical differentiation along the line ${f(x+\epsilon/2) - f(x-\epsilon/2)\over \epsilon}$ and so on, but which becomes inaccurate for higher derivatives.
For small x I can use the representation of the function $f(x)$ as a power series involving the Stieltjes-constants themselves, but a suitable computation is then possible only for a small range of x-values if I use only, say 100 or 200 terms of the power series.
Q: What would be a more efficient representation of $f(x)$ and/or $g(x)$ for the actual computation, where -say- x and d both can go up to 100 or 200 ?
More context: what I'm finally after is to be able to explore the partial sums and its single terms in
$$ S \underset{\mathfrak E}{=} \lim_{c\to \infty} S(c) = \sum_{k=0}^{c} (-1)^k g( c+1-k ,k) $$ which should be a representation for the (divergent) sum of the Stieltjes constants, taken by some "Eulerian" summation procedure $\mathfrak E$.