Newton-Cotes integration for approximating Stieltjes (generalized Euler) constants

In the Laurent expansion

of the Riemann-Hurwitz zeta function, the coefficients $\gamma _{k}(a)$ are known as Stieltjes, or generalized Euler, constants. [When $a=1$, $\zeta (s,1)=\zeta (s)$ (the Riemann zeta function), and $\gamma _{k}(1)=\gamma _{k}$.] We present a new approach to high-precision approximation of $\gamma _{k}(a)$. Plots of our results reveal much structure in the growth of the generalized Euler constants. Our results when $1\leq k\leq 3200$ for $\gamma _{k}$, and when $1\leq k\leq 600$ for $\gamma _{k}(a)$ (for $a$ such as 53/100, 1/2, etc.) suggest that published bounds on the growth of the Stieltjes constants can be much improved, and lead to several conjectures. Defining $g(k)=\sup _{0<a\leq 1}\vert\gamma _{k}(a)-\frac{\log ^{k} a}{a}\vert$, we conjecture that $g$ is attained: for any given $k$, $g(k)= \vert\gamma _{k}(a)-\frac{\log ^{k} a}{a}\vert$ for some $a$ (and similarly that, given $\epsilon$ and $a$, $g(k)$ is within $\epsilon$ of $\vert\gamma _{k}(a)-\frac{\log ^{k} a}{a}\vert$ for infinitely many $k$). In addition we conjecture that $g$ satisfies $[\log \big (g(k)\big )]/k <\log(\log(k))$ for $k>1$. We also conjecture that $\lim _{k\rightarrow \infty }\big ( \gamma _{k}(1/2)+\gamma _{k}\big )/\gamma _{k} = 0$, a special case of a more general conjecture relating the values of $\gamma _{k}(a)$ and $\gamma _{k}(a+\frac{1}{2})$ for $0<a\leq \frac{1}{2}$. Finally, it is known that $\gamma _{k} = \lim _{n\rightarrow \infty }\{\sum _{j=2}^{n} \frac{\log ^{k} j}{j}- \frac{\log ^{k+1} n}{k+1}\}$ for $k=1,2,\dots$. Using this to define $\gamma _{r}$ for all real $r>0$, we conjecture that for nonintegral $r$, $\gamma _{r}$ is precisely $(-1)^{r}$ times the $r$-th (Weyl) fractional derivative at $s=1$ of the entire function $\zeta (s)-1/(s-1)-1$. We also conjecture that $g$, now defined for all real arguments $r>0$, is smooth. Our numerical method uses Newton-Cotes integration formulae for very high-degree interpolating polynomials; it differs in implementation from, but compares in error bounding to, Euler-Maclaurin summation based methods.