General formulas for a class of Euler sums

1. Problem Statement

The paper addresses the analytical evaluation of Euler sums (also known as harmonic sums), which are infinite series involving harmonic numbers $H_k = \sum_{j=1}^k \frac{1}{j}$. Specifically, the authors focus on sums of the general form:
$$S = \sum_{k=1}^{\infty} R(k)H_k$$
where $R(k)$ is a rational function (a quotient of two polynomials) such that the degree of the denominator is at least two greater than the degree of the numerator (ensuring convergence).

The authors aim to derive explicit closed-form evaluations for these sums, particularly for cases where the denominator consists of linear factors raised to powers, such as:
$$\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}}$$
and related variations involving powers of $k$ in the numerator (e.g., $\sum k^q H_k / (mk+n)^p$). While specific cases have been studied historically, a general, algorithmic approach for arbitrary rational functions $R(k)$ with mixed denominators was lacking.

2. Methodology

The core methodology relies on complex analysis, specifically the calculus of residues, combined with properties of the digamma ($\psi$) and polygamma functions.

A. Theoretical Foundation

The authors utilize a key proposition derived from the behavior of the digamma function near negative integers. Let $R(z)$ be a rational function with no poles at positive integers. The sum can be expressed as:
$$\sum_{k=1}^{\infty} R(k)H_k = \frac{1}{2} \sum_{\alpha \in \text{poles of } R} \text{Res}_{s=\alpha} \left[ R'(s)(\psi(-s) + \gamma) - R(s)(\psi(-s) + \gamma)^2 \right]$$
where $\gamma$ is the Euler-Mascheroni constant. This formula transforms the infinite summation problem into a finite calculation of residues at the poles of $R(s)$.

B. Algorithmic Implementation

The authors developed an algorithm (named HarmonicSum) to automate this process:

1. Factorization: The denominator of $R(s)$ is factored into irreducible polynomials.
2. Residue Calculation: For each pole $\alpha$ (which may be a root of a polynomial), the algorithm computes the Taylor and Laurent expansions of $R(s)$, $R'(s)$, and the term $(\psi(-s) + \gamma)$ around $\alpha$.
3. Symbolic Summation: Instead of finding numerical values for the poles, the algorithm performs symbolic summation over the conjugate roots of the irreducible factors. This avoids the need for explicit root finding and maintains exactness.
4. Output: The result is expressed as a linear combination of digamma and polygamma functions evaluated at rational arguments.

C. Numerical Verification

To ensure the correctness of the derived formulas, the authors employed a rigorous numerical verification scheme:

• They computed the Euler sums directly to extremely high precision (approx. 280–290 digits) using the Euler-Maclaurin summation formula to handle the tail of the series.
• They used the PSLQ integer relation algorithm to identify rational coefficients linking the high-precision numerical sum to hypothesized combinations of constants (like $\zeta$ values, $\pi$, and polygamma functions).
• This "numerical guessing" confirmed the analytic formulas derived via residues in hundreds of specific cases.

3. Key Contributions

A. General Closed-Form Theorems

The paper establishes three main theorems that provide general formulas for these sums:

• Theorem 1 (Simple Poles): Provides a formula for sums where the denominator has distinct linear factors (simple poles). The result involves the function $T(t, 1)$, which combines $\psi'(t)$, $\psi(t)$, and $\gamma$.
• Theorem 2 (Higher Order Poles): Generalizes the result to poles of order $p \geq 2$. It defines a function $T(t, p)$ involving higher-order polygamma functions $\psi^{(k)}(t)$ and their products.
• Theorem 3 (Mixed Denominators): Combines the previous results to handle arbitrary rational functions via partial fraction decomposition. The final sum is a weighted sum of $T(t_j, i)$ terms corresponding to the partial fraction coefficients.

B. Explicit Catalog of Formulas

The authors generated a comprehensive catalog of explicit formulas for:

• Denominators with one or two terms (e.g., $(mk+n)^p$ or $(m_1k+n_1)^{p_1}(m_2k+n_2)^{p_2}$).
• Powers $p_i$ ranging up to 3 (and up to 8 in the appendix for single-term denominators).
• Multipliers $m_i$ ranging up to 4.
• Variations including $k^q H_k$ in the numerator.

These formulas are expressed entirely in terms of $\gamma$, $\psi^{(n)}(rational)$, and products thereof.

C. Computational Tools

• Algorithm HarmonicSum: A Maple implementation (provided in the appendix) that takes a rational function as input and outputs the closed-form sum.
• High-Precision Evaluation: The paper details efficient algorithms for computing digamma and polygamma functions to arbitrary precision, noting that evaluating the closed-form formula is significantly faster (e.g., 30,000x) than direct summation of the series.

4. Results

• Explicit Formulas: The paper presents over 100 specific formulas in the appendices, covering cases like $\sum \frac{H_k}{(2k+1)^p}$, $\sum \frac{H_k}{(3k+1)^p}$, and mixed cases like $\sum \frac{H_k}{(2k+1)^2(3k+1)}$.
• Complex Roots: The method successfully handles cases where the denominator has complex roots (e.g., $k^3+1$), expressing the result in terms of $\psi$ evaluated at complex rational points, which can be further reduced using reflection and multiplication formulas.
• Efficiency: The authors demonstrate that while the closed-form expressions appear complex (involving sums of polygamma functions), they are computationally superior to direct summation for high-precision requirements.

5. Significance

• Unification: The work provides a unified framework for evaluating a broad class of Euler sums that previously required ad-hoc methods or were only known for specific small integer parameters.
• Algorithmic Automation: It demonstrates that these sums can be computed algorithmically by computer algebra systems, removing the need for manual residue calculus for each new case.
• Computational Efficiency: The paper highlights a practical application: using these analytic formulas to accelerate numerical computations. Evaluating a sum via the derived digamma/polygamma formula is orders of magnitude faster than summing the series term-by-term, even with high-precision arithmetic.
• Foundation for Future Work: The authors suggest that the residue-based approach can be extended to more complex sums, such as those involving $H_k^2$, alternating signs $(-1)^k H_k$, or even-indexed harmonic numbers $H_{2k}$, opening avenues for further research in mathematical physics and number theory.

In summary, this paper bridges the gap between theoretical residue calculus and practical computation, providing a robust, automated method to derive and verify closed-form solutions for a wide variety of Euler sums.