General formulas for a class of Euler sums

Let $H_k = 1 + 1/2 + 1/3 + \cdots + 1/k$ denote the $k$th harmonic number. We present an easy-to-implement algorithm for the computation of explicit closed-form evaluations, in terms of the digamma and polygamma functions, for Euler sums of the form \begin{align} \sum_{k=1}^\infty R(k) H_k, \end{align} where $R(k)$ is a rational function (quotient of two polynomials) whose denominator degree is at least two larger than the numerator degree. We apply the same method to show how the computation of a general formula for Euler sums of the form \begin{align*} \sum_{k=1}^\infty \frac{H_k}{(m_1 k + n_1)^{p_1} (m_2 k + n_2)^{p_2} \cdots (m_r k + n_r)^{p_r}} \end{align*} reduces to partial fraction decomposition. We present explicit formulae for sums with one or two terms in the denominator, with powers $p_i$ ranging up to 3, and with multipliers $m_i$ ranging up to 4. We also include results for related Euler sums such as \begin{align*} \sum_{k=1}^\infty \frac{k^q H_k}{(m k + n)^p}. \end{align*} Computation of Euler sums directly to very high precision enables us to rigorously check the above-mentioned formulas in many specific cases.