The Generalization of Faulhaber’s Formula to Sums of Arbitrary Complex Powers

Raphael Schumacher Department of Mathematics ETH Zurich Rämistrasse 101 8092 Zurich Switzerland

                raphael.schumacher@math.ethz.ch

Abstract In this paper we present a generalization of Faulhaber’s formula to sums of arbitrary complex powers m ∈ ℂ. These summation formulas for sums of the form ∑_{k=1}^{⌊x⌋} k^{m} and ∑_{k=1}^{n} k^{m}, where x ∈ ℝ^{+} and n ∈ ℕ, are based on a series acceleration involving Stirling numbers of the first kind. While it is well-known that the corresponding expressions obtained from the Euler-Maclaurin summation formula diverge, our summation formulas are all very rapidly convergent.

1 Introduction

For two natural numbers m, n ∈ ℕ_{0}, the Faulhaber formula [1], given by

∑_{k=0}^{n} k^{m} = \frac{1}{m + 1} ∑_{k=0}^{m}(−1)^{k}\binom{m + 1}{k}B_{k}n^{m−k+1},  (1)

where the B_{k}’s are the Bernoulli numbers, provides a very efficient way to compute the sum of the m-th powers of the first n natural numbers. This formula was found by Jacob Bernoulli around 1700 and was first proved by Carl Gustav Jacobi in 1834.

We prove a rapidly convergent exact generalization of Faulhaber’s formula to finite sums of the form ∑_{k=1}^{⌊x⌋} k^{m} and ∑_{k=1}^{n} k^{m} for all exponents m ∈ ℂ. Our key tool is the so-called Weniger transformation [2, (4.1)] found by J. Weniger, transforming an inverse power series into an inverse factorial series [2, (1.1)]. This transformation of inverse power series was first found by Oskar Schlömilch around 1850 [3, 4, 5] based on earlier works of James Stirling in 1730 [6].

In an expanded form, one of our summation formulas for the sum ∑_{k=1}^{n} \sqrt{k}, where n ∈ ℕ, looks like

\begin{aligned}\sqrt[∑_{k=1}^{n}]{k} & = \frac{2}{3}n^{3/2} + \frac{1}{2}\sqrt{n} − \frac{1}{4π} ζ (\frac{3}{2}) + \sqrt{n} ∑_{k=1}^{∞}(−1)^{k+1} \frac{∑_{l=1}^{k} \frac{(2l−3)!!}{2^{l}(l+1)!}B_{l+1}S_{k}^{(1)}(l)}{(n + 1)(n + 2) ⋯ (n + k)} \\ & = \frac{2}{3}n^{3/2} + \frac{1}{2}\sqrt{n} − \frac{1}{4π} ζ (\frac{3}{2}) + \frac{\sqrt{n}}{24(n + 1)} + \frac{\sqrt{n}}{24(n + 1)(n + 2)} \\ & + \frac{53\sqrt{n}}{640(n + 1)(n + 2)(n + 3)} + \frac{79\sqrt{n}}{320(n + 1)(n + 2)(n + 3)(n + 4)} + ⋯, \end{aligned} (2)

where the B_{l}’s are the Bernoulli numbers and S_{k}^{(1)}(l) denotes the Stirling numbers of the first kind.

The above identity (2) is deduced by setting the variable x := n ∈ ℕ into the more general formula

\begin{matrix}∑_{k=1}^{⌊x⌋} & \sqrt{k} \\ = \frac{2}{3}x^{3/2} − \frac{1}{4π} ζ (\frac{3}{2}) − \sqrt{x}B_{1}(\{x\}) + \sqrt{x} ∑_{k=1}^{∞}(−1)^{k} \frac{∑_{l=1}^{k}(−1)^{l} \frac{(2l−3)!!}{2^{l}(l+1)!}S_{k}^{(1)}(l)B_{l+1}(\{x\})}{(x + 1)(x + 2) ⋯ (x + k)} \\ = \frac{2}{3}x^{3/2} − \frac{1}{4π} ζ (\frac{3}{2}) + (\frac{1}{2} − \{x\}) \sqrt{x} + \frac{(\frac{1}{4}\{x\}^{2} − \frac{1}{4}\{x\} + \frac{1}{24}) \sqrt{x}}{(x + 1)} \\ + \frac{(\frac{1}{24}\{x\}^{3} + \frac{3}{16}\{x\}^{2} − \frac{11}{48}\{x\} + \frac{1}{24}) \sqrt{x}}{(x + 1)(x + 2)} + \frac{(\frac{1}{64}\{x\}^{4} + \frac{3}{32}\{x\}^{3} + \frac{21}{64}\{x\}^{2} − \frac{7}{16}\{x\} + \frac{53}{640}) \sqrt{x}}{(x + 1)(x + 2)(x + 3)} \\ + \frac{(\frac{1}{128}\{x\}^{5} + \frac{19}{256}\{x\}^{4} + \frac{109}{384}\{x\}^{3} + \frac{29}{32}\{x\}^{2} − \frac{977}{768}\{x\} + \frac{79}{320}) \sqrt{x}}{(x + 1)(x + 2)(x + 3)(x + 4)} + ⋯, \end{matrix} ; (3)

where this time the B_{l}({x})’s are the fractional Bernoulli polynomials and x ∈ ℝ^{+} is a positive real number. All other formulas in this article have a similar shape, when we expand them.

We have searched our resulting formulas in the literature and on the internet. We could find only two of them, namely equation (46) and its analogues for the sums ∑_{k=1}^{n} \frac{1}{k^{m}} with m ∈ ℕ_{≥2}, which were already known to Stirling in 1730 [3, 7], and equation (45), which was obtained by Gregorio Fontana around 1780 [3, 8]. Both of these formulas were originally found in another form without the use of Bernoulli numbers.

We believe that all other generalized Faulhaber formulas presented in this article are new and that our method to obtain them has not been recognized before.

2 Definitions and basic facts

As usual, we denote the floor of x by ⌊x⌋ and the fractional part of x by {x}. The symbol ℕ := {1, 2, 3, 4, …} denotes the set of natural numbers and ℝ^{+} := {x ∈ ℝ : x > 0} represents the set of positive real numbers. We also set ℕ_{0} := ℕ ∪ {0} and ℝ_{0}^{+} := ℝ^{+} ∪ {0}. Moreover, we define ℍ^{+} := {z ∈ ℂ : Re(z) > 0} and let ζ(s) denote the Riemann zeta function at the point s ∈ ℂ \ {1}. For a complex number z = re^{iφ} ∈ ℂ, we let |z| = r ∈ ℝ_{0}^{+} denote its absolute value and by φ = arg(z) ∈ (−π, π] its argument or phase. We define for all θ ∈ (−\frac{π}{2}, \frac{π}{2}) the secant function by sec(θ) := \frac{1}{cos(θ)}. For z ∈ ℂ, the notation z → ∞ means that |z| → ∞.

The double factorial function for n ∈ ℕ_{0} is defined by n!! := ∏_{k=0}^{⌊(n−1)/2⌋}(n − 2k).

Definition 1. (Pochhammer symbol [2, p. 1429]) We define the Pochhammer symbol (or rising factorial function) (z)_{k} by

(z)_{k} := z(z + 1)(z + 2)(z + 3) ⋯ (z + k − 1) = \frac{Γ(z + k)}{Γ(z)},  (4)

where Γ(z) is the gamma function [9, (5.2.1), p. 136] defined as the meromorphic continuation of the integral

Γ(z) := ∫_{0}^{∞} e^{−t}t^{z−1}dt for all z ∈ ℂ with Re(z) > 0  (5)

to the whole complex plane ℂ.

Definition 2. (Stirling numbers of the first kind [2, (A.2), p. 1437] and [10, A008275]) Let k, l ∈ ℕ_{0} be two non-negative integers such that k ≥ l ≥ 0. We set the Stirling numbers of the first kind S_{k}^{(1)}(l) as the connecting coefficients in the identity

(z)_{k} = (−1)^{k} ∑_{l=0}^{k}(−1)^{l}S_{k}^{(1)}(l)z^{l},  (6)

where (z)_{k} is the rising factorial function. Furthermore, we set S_{k}^{(1)}(l) = 0 if k, l ∈ ℕ_{0} with l > k.

Definition 3. (Binomial coefficients [11]) We introduce the binomial coefficient \binom{z}{s} for all z ∈ ℂ and s ∈ ℂ by [11, (5) and (11), pp. 8–9]

\binom{z}{s} : = \frac{Γ(z + 1)}{Γ(s + 1)Γ(z − s + 1)}.  (7)

Moreover, for z ∈ ℂ \ {0, −1, −2, −3, …} we have the following asymptotic expansion as k → ∞ [11, (18), p. 2 and p. 35]

\binom{z}{k} = \frac{(−1)^{k}}{Γ(−z)k^{z+1}} + O (\frac{1}{k^{z+2}} ) for all z ∈ ℂ and k ∈ ℕ.  (8)

Definition 4. (Bernoulli polynomials and Bernoulli numbers [1, 12, 13, 14]) We define for n ∈ ℕ_{0} the n-th Bernoulli polynomial B_{n}(x) via the following exponential generating function [1] as

\frac{te^{xt}}{e^{t} − 1} = ∑_{n=0}^{∞} \frac{B_{n}(x)}{n!}t^{n} ∀t ∈ ℂ with |t| < 2π.  (9)

We also define the n-th Bernoulli number B_{n} as the value of the n-th Bernoulli polynomial B_{n}(x) at the point x = 0, that is

B_{n} := B_{n}(0).  (10)

For all n ∈ ℕ_{0} we have the explicit formula [12, Proposition 23.2, p. 86]

B_{n}(x) = ∑_{k=0}^{n} \binom{n}{k}B_{k}x^{n−k}.  (11)

For all 0 ≤ y ≤ 1 it is valid that [13, Corollary B.4, (B.21), p. 500]

|B_{1}(y)| ≤ \frac{1}{2} and that |B_{n}(y)| ≤ \frac{2ζ(n)n!}{(2π)^{n}} for all n ∈ ℕ_{≥2}.  (12)

We have [14, (1.10), p. 282]

(−1)^{k}B_{k}(1 − y) = B_{k}(y) for all k ∈ ℕ_{0} and 0 ≤ y ≤ 1.  (13)

Definition 5. (Digamma function [9, pp. 136–138]) We set the digamma function ψ(z) to

ψ(z) : = \frac{Γ′(z)}{Γ(z)} for all z ∈ ℂ \ \{0, −1, −2, −3, …\}.  (14)

Therefore, ψ(z) is an analytic function for all z ∈ ℂ\(−∞, 0]. For all z ∈ ℂ\{0, −1, −2, −3, …}, we have the identity [9, (5.5.2), p. 138]

ψ(z + 1) = ψ(z) + \frac{1}{z}  (15)

and for all n ∈ ℕ we have the formula [9, (5.4.14), p. 137]

∑_{k=1}^{n} \frac{1}{k} = ψ(n + 1) + γ.  (16)

Definition 6. (Hurwitz zeta function [9, p. 607]) We define the Hurwitz zeta function ζ(s, z) for all complex numbers s ∈ ℂ with Re(s) > 1 and all z ∈ ℂ \ {0, −1, −2, −3, …} by

ζ(s, z) : = ∑_{k=0}^{∞} \frac{1}{(k + z)^{s}} .  (17)

The function ζ(s, z) extends to an analytic function on ℂ \ {0, −1, −2, −3, …} in the zvariable and for every z ∉ {0, −1, −2, −3, …} to a meromorphic function in s ∈ ℂ\{1} with a simple pole at s = 1. For all s ∈ ℂ \ {1} and all z ∈ ℂ \ {0, −1, −2, −3, …} it satisfies the identity

ζ(s, z + 1) = ζ(s, z) − \frac{1}{z^{s}} .  (18)

For all m ∈ ℂ \ {1} and all n ∈ ℕ we have the formula [15, (1.2), p. 2]

∑_{k=1}^{n} k^{m} = ζ(−m) − ζ(−m, n + 1).  (19)

3 The structure of inverse factorial series expansions

In this section we study the structure of inverse factorial series expansions for analytic functions possessing an asymptotic series expansion by applying a theorem of G. N. Watson [16, Theorem 2, p. 45]. The main result of this section is Theorem 10, from which we later deduce convergent inverse factorial series expansions for the functions ζ(s, z + 1 − y) and ψ(z + 1 − y), where 0 ≤ y ≤ 1.

For this procedure, we need the following variant of a result found by J. Weniger [2, (4.1), p. 1433].

Lemma 7. (finite Weniger transformation [2]) For every finite inverse power series ∑_{k=1}^{n} \frac{a_{k}}{z^{k}}, where the a_{k}’s are any complex numbers and n ∈ ℕ, the following transformation formula holds:

∑_{k=1}^{n} \frac{a_{k}}{z^{k}} = ∑_{k=1}^{n} \frac{(−1)^{k} ∑_{l=1}^{k}(−1)^{l}S_{k}^{(1)}(l)a_{l}}{(z + 1)(z + 2) ⋯ (z + k)} + O (\frac{1}{z^{n+1}} ) as z → ∞.  (20)

Moreover, we have

∑_{k=1}^{n} \frac{a_{k}}{z^{k}} = ∑_{k=1}^{∞} \frac{(−1)^{k} ∑_{l=1}^{n}(−1)^{l}S_{k}^{(1)}(l)a_{l}}{(z + 1)(z + 2) ⋯ (z + k)} .  (21)

Proof. For l ∈ ℕ we have ([2, (A.14), p. 1438], [17, (6), p. 78])

\frac{1}{z^{l}} = ∑_{k=0}^{∞} \frac{(−1)^{k}S_{k+l}^{(1)}(l)}{(z + 1)(z + 2) ⋯ (z + k + l)}

= ∑_{k=l}^{∞} \frac{(−1)^{k−l}S_{k}^{(1)}(l)}{(z + 1)(z + 2) ⋯ (z + k)}

= ∑_{k=1}^{∞} \frac{(−1)^{k−l}S_{k}^{(1)}(l)}{(z + 1)(z + 2) ⋯ (z + k)}

= ∑_{k=1}^{n} \frac{(−1)^{k−l}S_{k}^{(1)}(l)}{(z + 1)(z + 2) ⋯ (z + k)} + O (\frac{1}{z^{n+1}} ) as z → ∞,

where we used in the third step that S_{k}^{(1)}(l) = 0 for k < l.

Therefore, we obtain that

∑_{k=1}^{n} \frac{a_{k}}{z^{k}} = ∑_{l=1}^{n} \frac{a_{l}}{z^{l}} = ∑_{l=1}^{n} a_{l} ∑_{k=1}^{n} \frac{(−1)^{k−l}S_{k}^{(1)}(l)}{(z + 1)(z + 2) ⋯ (z + k)} + O (\frac{1}{z^{n+1}} )

= ∑_{k=1}^{n} \frac{(−1)^{k} ∑_{l=1}^{n}(−1)^{l}S_{k}^{(1)}(l)a_{l}}{(z + 1)(z + 2) ⋯ (z + k)} + O (\frac{1}{z^{n+1}} )

= ∑_{k=1}^{n} \frac{(−1)^{k} ∑_{l=1}^{k}(−1)^{l}S_{k}^{(1)}(l)a_{l}}{(z + 1)(z + 2) ⋯ (z + k)} + O (\frac{1}{z^{n+1}} ) as z → ∞,

which is the first claimed formula (20).

The second formula (21) follows from the calculation

\begin{aligned}∑_{k=1}^{n} \frac{a_{k}}{z^{k}} = ∑_{l=1}^{n} \frac{a_{l}}{z^{l}} & = ∑_{l=1}^{n} a_{l} ∑_{k=1}^{∞} \frac{(−1)^{k−l}S_{k}^{(1)}(l)}{(z + 1)(z + 2) ⋯ (z + k)} \\ & = ∑_{k=1}^{∞} \frac{(−1)^{k} ∑_{l=1}^{n}(−1)^{l}S_{k}^{(1)}(l)a_{l}}{(z + 1)(z + 2) ⋯ (z + k)} , \end{aligned}

because we can always interchange a finite summation with an infinite summation.

Lemma 8. (Uniqueness of inverse factorial series expansions) If for all z ∈ ℂ with Re(z) > 0 a function f has the absolutely convergent series expansion

f(z) = ∑_{k=1}^{∞} \frac{b_{k}}{(z + 1)(z + 2) ⋯ (z + k)}

and the asymptotic expansion

f(z) = ∑_{k=1}^{n} \frac{c_{k}}{(z + 1)(z + 2) ⋯ (z + k)} + O (\frac{1}{z^{n+1}} ) as z → ∞,

then we have c_{k} = b_{k} for all k ∈ ℕ and the absolutely convergent series expansion

f(z) = ∑_{k=1}^{∞} \frac{c_{k}}{(z + 1)(z + 2) ⋯ (z + k)}.

Proof. For all n ∈ ℕ we have ∑_{k=n+1}^{∞} \frac{|b_{k}|·|z|^{n}}{|(z+1)||(z+2)|···|(z+k)|} → 0 as z → ∞, because this infinite series is convergent for all z ∈ ℂ with Re(z) > 0 by assumption and consists of the monotone decreasing positive terms \frac{|b_{k}|·|z|^{n}}{|(z+1)||(z+2)|···|(z+k)|} → 0 as z → ∞.

From the above observation and the given absolutely convergent inverse factorial series expansion of f(z), we deduce for all n ∈ ℕ that

|∑_{k=n}^{∞} \frac{b_{k}}{(z + 1)(z + 2) ⋯ (z + k)}|

≤ ∑_{k=n}^{∞} \frac{|b_{k}|}{|(z + 1)||(z + 2)| ⋯ |(z + k)|}

= \frac{|b_{n}|}{|(z + 1)||(z + 2)| ⋯ |(z + n)|} + \frac{1}{|z|^{n}} ∑_{k=n+1}^{∞} \frac{|b_{k}| · |z|^{n}}{|(z + 1)||(z + 2)| ⋯ |(z + k)|}

= O (\frac{1}{z^{n}} ) as z → ∞,

which means that lim_{z→∞}(f(z)) = 0 and we have

z^{m} ∑_{k=n}^{∞} \frac{b_{k}}{(z + 1)(z + 2) ⋯ (z + k)} ⟶ 0 as z → ∞ for all m ∈ \{0, 1, 2, …, n − 1\}.

The result now follows by induction on n ∈ ℕ via a repeated application of the above limit. For n = 1, we get

f(z) = \frac{b_{1}}{z + 1} + ∑_{k=2}^{∞} \frac{b_{k}}{(z + 1)(z + 2) ⋯ (z + k)} = \frac{c_{1}}{z + 1} + O (\frac{1}{z^{2}} ) as z → ∞,

which implies by multiplying both sides with z + 1 and letting z → ∞ that c_{1} = b_{1}. Similarly, for n = 2, we get using c_{1} = b_{1} that

f(z) − \frac{c_{1}}{z + 1} = \frac{b_{2}}{(z + 1)(z + 2)} + ∑_{k=3}^{∞} \frac{b_{k}}{(z + 1)(z + 2) ⋯ (z + k)} = \frac{c_{2}}{(z + 1)(z + 2)} + O (\frac{1}{z^{3}} ) ,

which implies by multiplying both sides with (z + 1)(z + 2) and letting z → ∞ that c_{2} = b_{2}.

In general, we can induct from n−1 to n using that c_{k} = b_{k} for all k ∈ {1, 2, 3, …, n−1} by the identity

f(z) − ∑_{k=1}^{n−1} \frac{c_{k}}{(z + 1)(z + 2) ⋯ (z + k)}

= \frac{b_{n}}{(z + 1)(z + 2) ⋯ (z + n)} + ∑_{k=n+1}^{∞} \frac{b_{k}}{(z + 1)(z + 2) ⋯ (z + k)}

= \frac{c_{n}}{(z + 1)(z + 2) ⋯ (z + n)} + O (\frac{1}{z^{n+1}} ) as z → ∞,

again by multiplying both sides with (z + 1)(z + 2) ⋯ (z + n) and letting z → ∞ to conclude that c_{k} = b_{k} holds also for k = n. This proves that c_{k} = b_{k} for all k ∈ ℕ.

The key to our generalized Faulhaber formulas is the following theorem.

Theorem 9. (Watson’s Transformation Theorem [16, Theorem 2, p. 45]) Let f(z) be a function of z ∈ ℂ which is analytic when Re(z) > 0; and let f(z) be also analytic in the region D of the complex plane defined by

D : = \{z ∈ ℂ : |z| > γ and | arg(z)| ≤ \frac{π}{2} + α + 3δ\} ,

where γ ≥ 0 is a finite number, α > 0, δ > 0 and α + 3δ < \frac{π}{2}.

In the region D let f(z) possess the asymptotic expansion

f(z) = ∑_{k=0}^{n} \frac{a_{k}}{z^{k}} + R_{n}(z) = a_{0} + \frac{a_{1}}{z} + \frac{a_{2}}{z^{2}} + \frac{a_{3}}{z^{3}} + \frac{a_{4}}{z^{4}} + ⋯ + \frac{a_{n}}{z^{n}} + R_{n}(z),

where

|a_{n}| < Aρ^{n}n! and |R_{n}(z)z^{n+1}| < Bσ^{n}n!

with some constants A, B, ρ and σ, which are independent of n.

Let M ≤ M_{0} be any positive real number, where M_{0} is the largest positive root of the equation

e^{− \frac{2 cos(α)}{ρM_{0}}} − 2 cos (\frac{sin(α)}{ρM_{0}}) · e^{− \frac{cos(α)}{ρM_{0}}} + 1 − p^{2} = 0,

where

1 < p < 1 + e^{−π cot(α)}.

Then the function f(z) can be expanded into the absolutely convergent series

f(z) = b_{0} + ∑_{k=1}^{∞} \frac{b_{k}}{(Mz + w + 1)(Mz + w + 2) ⋯ (Mz + w + k)},

when Re(z) > 0 and w ∈ ℂ with Re(w) ≥ 0.

Proof. The proof of this theorem is given in Watson’s paper [16].

Theorem 10. (Structure of inverse factorial series expansions) Let f(z) be a function of z ∈ ℂ which is analytic when Re(z) > 0; and let f(z) be also analytic in the region D of the complex plane defined by

D : = {z ∈ ℂ : |z| > 0 and | arg(z)| ≤ π − ε, where ε > 0 is arbitrarily small} .

In the region D let f(z) possess the asymptotic expansion

f(z) = ∑_{k=1}^{n} \frac{a_{k}}{z^{k}} + R_{n}(z) = \frac{a_{1}}{z} + \frac{a_{2}}{z^{2}} + \frac{a_{3}}{z^{3}} + ⋯ + \frac{a_{n}}{z^{n}} + R_{n}(z),

where

|a_{n}| < Aρ^{n}n! and |R_{n}(z)z^{n+1}| < Bσ^{n}n!

with some constants A, B, ρ < \frac{3}{π} and σ, which are independent of n.

Then the function f(z) is equal to the absolutely convergent inverse factorial series

f(z) = ∑_{k=1}^{∞} \frac{(−1)^{k} ∑_{l=1}^{k}(−1)^{l}S_{k}^{(1)}(l)a_{l}}{(z + 1)(z + 2) ⋯ (z + k)}  (22)

for Re(z) > 0.

Proof. Let f(z) and the region D be as described in the above Theorem 10. Because of the conditions on the function f(z) and the region D, we can choose in Theorem 9 the variables γ := 0, α := \frac{π}{2} −4ε, δ := ε and p := 1+ε for ε > 0 arbitrarily small by [16, beginning of p. 85]. We have then that M_{0} = \frac{3}{πρ} − ε for some arbitrarily small number ε > 0 and because ρ < \frac{3}{π}, we obtain that M_{0} > 1. According to Watson’s Theorem 9 with M := 1 < M_{0}, w := 0 and a_{0} = b_{0} = 0, we know that we can expand the function f(z) into an absolutely convergent series of the form

f(z) = ∑_{k=1}^{∞} \frac{b_{k}}{(z + 1)(z + 2) ⋯ (z + k)}

for some constants b_{k} ∈ ℂ and all z ∈ ℂ with Re(z) > 0.

On the other hand, we have by applying a finite Weniger transformation (20) to the asymptotic expansion of f(z) that

f(z) = ∑_{k=1}^{n} \frac{(−1)^{k} ∑_{l=1}^{k}(−1)^{l}S_{k}^{(1)}(l)a_{l}}{(z + 1)(z + 2) ⋯ (z + k)} + O (\frac{1}{z^{n+1}} ) as z → ∞

also holds.

Comparing the two expressions above for f(z) by using Lemma 8 with

c_{k} := (−1)^{k} ∑_{l=1}^{k}(−1)^{l}S_{k}^{(1)}(l)a_{l},

we conclude that there is an absolutely convergent series

f(z) = ∑_{k=1}^{∞} \frac{(−1)^{k} ∑_{l=1}^{k}(−1)^{l}S_{k}^{(1)}(l)a_{l}}{(z + 1)(z + 2) ⋯ (z + k)}

for all z ∈ ℂ with Re(z) > 0.

4 The convergent inverse factorial series expansions for ζ(s, z + 1 − y) and ψ(z + 1 − y)

In this section, in Theorem 13 we deduce the convergent inverse factorial series expansions for the functions ζ(s, z + 1 − y) and ψ(z + 1 − y), where 0 ≤ y ≤ 1.

For this, we need the following lemma.

Lemma 11. (Euler-Maclaurin summation formula [13, Theorem B.5, pp. 500–501]) Suppose that n ∈ ℕ is a positive integer and that the function f(t) has continuous derivatives through the n-th order on the interval [a, b] where a and b are real numbers with a < b. Then we have

\begin{aligned}∑_{a<k≤b} f(k) & = ∫_{a}^{b} f(t)dt + ∑_{k=1}^{n}(−1)^{k} \frac{B_{k}(\{b\})}{k!}f^{(k−1)}(b) − ∑_{k=1}^{n}(−1)^{k} \frac{B_{k}(\{a\})}{k!}f^{(k−1)}(a) \\ & + \frac{(−1)^{n+1}}{n!} ∫_{a}^{b} f^{(n)}(t)B_{n}(\{t\})dt. \end{aligned} (23)

Proof. The proof of this Lemma 11 is given in [13, p. 501].

From the above Lemma 11, it follows the next lemma.

Lemma 12. (Asymptotic series expansions for ζ(s, z + h) and ψ(z + h) with 0 ≤ h ≤ 1) Let n ∈ ℕ_{0} and let s ∈ ℂ \ {1} such that Re(s) > −n.We for z ∈ ℂ with | arg(z)| < π and 0 ≤ h ≤ 1 the asymptotic series expansions

ζ(s, z + h) = \frac{z^{1−s}}{s − 1} + \frac{z^{1−s}}{s − 1} ∑_{k=1}^{n} \binom{1 − s}{k}\frac{B_{k}(h)}{z^{k}} + O_{n}(z)  (24)

with

|O_{n}(z)| = |\binom{1 − s}{n + 2}\frac{n + 2}{s − 1} ∫_{0}^{∞} \frac{B_{n+1} (\{x − h\}) − (−1)^{n+1}B_{n+1}(h)}{(x + z)^{n+s+1}}dx|

≤ \frac{2(n + 2)}{|s − 1|} |\binom{1 − s}{n + 2}| \frac{|B_{n+1}| sec^{n+Re(s)+1} (\frac{1}{2} arg(z))}{(n + Re(s))|z|^{n+Re(s)}} max \{1, e^{Im(s) arg(z)}\}

(25)

and

ψ(z + h) = log(z) − ∑_{k=1}^{n} \frac{(−1)^{k}B_{k}(h)}{kz^{k}} + U_{n}(z)  (26)

with

|U_{n}(z)| = |∫_{0}^{∞} \frac{(−1)^{n+1}B_{n+1}(h) − B_{n+1}(\{x − h\})}{(x + z)^{n+2}}dx| ≤ \frac{2 |B_{n+1}| sec^{n+2} (\frac{1}{2} arg(z))}{(n + 1)|z|^{n+1}}. (27)

Proof. Let 0 < h ≤ 1 and let z ∈ ℂ with |arg(z)| < π. Setting a := −h, b := N and f(x) := \frac{1}{(z+h+x)^{s}} with \frac{d^{n}f(x)}{dx^{n}} = \frac{d^{n}}{dx^{n}} (\frac{1}{(z+h+x)^{s}} ) = −\frac{(n+1)!}{s−1} \binom{1−s}{n+1} \frac{1}{(z+h+x)^{n+s}} into Lemma 11, we obtain for s ∈ ℂ with Re(s) > 1 that

ζ(s, z + h) = ∑_{k=0}^{∞} \frac{1}{(z + h + k)^{s}} = lim_{N→∞} ( ∑_{−h<k≤N} \frac{1}{(z + h + k)^{s}} )

= ∫_{−h}^{∞} \frac{1}{(z + h + x)^{s}} dx + lim_{N→∞} (∑_{k=1}^{n}(−1)^{k} \frac{B_{k}(\{N\})}{k!} \frac{d^{k−1} (\frac{1}{(z+h+x)^{s}} )}{dx^{k−1}} |_{x=N} )

− ∑_{k=1}^{n}(−1)^{k} \frac{B_{k}(\{−h\})}{k!} \frac{d^{k−1} (\frac{1}{(z+h+x)^{s}} )}{dx^{k−1}} |_{x=−h} + \frac{(−1)^{n+1}}{n!} ∫_{−h}^{∞} \frac{d^{n} (\frac{1}{(z+h+x)^{s}} )}{dx^{n}} |_{x=t}B_{n}(\{t\})dt

= \frac{z^{1−s}}{s − 1} + \frac{z^{1−s}}{s − 1} ∑_{k=1}^{n}(−1)^{k}\binom{1 − s}{k}\frac{B_{k}(\{1 − h\})}{z^{k}} + (−1)^{n}\binom{1 − s}{n + 1}\frac{n + 1}{s − 1} ∫_{0}^{∞} \frac{B_{n}(\{x − h\})}{(x + z)^{n+s}}dx.

Now replacing n by n + 1, we get

\begin{aligned}ζ(s, z + h) & = \frac{z^{1−s}}{s − 1} + \frac{z^{1−s}}{s − 1} ∑_{k=1}^{n+1}(−1)^{k}\binom{1 − s}{k}\frac{B_{k}(\{1 − h\})}{z^{k}} \\ & + (−1)^{n+1}\binom{1 − s}{n + 2}\frac{n + 2}{s − 1} ∫_{0}^{∞} \frac{B_{n+1}(\{x − h\})}{(x + z)^{n+s+1}}dx. \end{aligned}

This is equivalent to

\begin{aligned}ζ(s, z + h) & = \frac{z^{1−s}}{s − 1} + \frac{z^{1−s}}{s − 1} ∑_{k=1}^{n}(−1)^{k}\binom{1 − s}{k}\frac{B_{k}(1 − h)}{z^{k}} \\ & + (−1)^{n+1}\binom{1 − s}{n + 2}\frac{n + 2}{s − 1} ∫_{0}^{∞} \frac{B_{n+1}(\{x − h\}) − B_{n+1}(1 − h)}{(x + z)^{n+s+1}}dx. \end{aligned}

We use the relation (13) and deduce equation (24), which extends ζ(s, z + h) analytically to the whole punctured complex s-plane ℂ \ {1}. Therefore, the equation (24) is also true for all s ∈ ℂ \ {1}. By using the identity (18), we see that the formula (24) is also true for h = 0. The bound (25) for the error term O_{n}(z) follows from [14, p. 294] and [15, p. 6]. This proves the first part about ζ(s, z + h) of the above Lemma 12.

Now, we prove the second part for ψ(z + h). For z ∈ ℂ with |arg(z)| < π, n ≥ 2 and 0 ≤ h ≤ 1 we have the following series expansion [14, Ex. 4.4, p. 295]:

ln (Γ(z + h)) = (z + h − \frac{1}{2}) ln(z) − z + \frac{1}{2} ln(2π) + ∑_{k=2}^{n} \frac{(−1)^{k}B_{k}(h)}{k(k − 1)z^{k−1}} − \frac{1}{n} ∫_{0}^{∞} \frac{B_{n} (\{x − h\})}{(x + z)^{n}}dx.

Replacing n by n + 1 again, we get

ln (Γ(z + h)) = (z + h − \frac{1}{2}) ln(z) − z + \frac{1}{2} ln(2π) + ∑_{k=2}^{n+1} \frac{(−1)^{k}B_{k}(h)}{k(k − 1)z^{k−1}}

− \frac{1}{n + 1} ∫_{0}^{∞} \frac{B_{n+1} (\{x − h\})}{(x + z)^{n+1}}dx

= (z + h − \frac{1}{2}) ln(z) − z + \frac{1}{2} ln(2π) + ∑_{k=2}^{n} \frac{(−1)^{k}B_{k}(h)}{k(k − 1)z^{k−1}}

+ \frac{1}{n + 1} ∫_{0}^{∞} \frac{(−1)^{n+1}B_{n+1}(h) − B_{n+1} (\{x − h\})}{(x + z)^{n+1}}dx.

Differentiating this identity with respect to the variable z, we get equation (26). The estimate (27) for the error term U_{n}(z) follows from [14, p. 294 and Ex. 4.2, p. 295].

We get the following theorem.

Theorem 13. (Inverse factorial series expansions for ζ(s, z + 1 − y) and ψ(z + 1 − y)) Let 0 ≤ y ≤ 1 and let s ∈ ℂ \ {1}. For z ∈ ℍ^{+} and all a ∈ ℕ_{0} we have the absolutely convergent inverse factorial series expansions

\begin{aligned}ζ(s, z + 1 − y) & = \frac{z^{1−s}}{s − 1} + \frac{z^{1−s}}{s − 1} ∑_{k=1}^{a}(−1)^{k}\binom{1 − s}{k}\frac{B_{k}(y)}{z^{k}} \\ & + \frac{z^{1−s−a}}{s − 1} ∑_{k=1}^{∞}(−1)^{k+a} \frac{∑_{l=1}^{k} \binom{1−s}{l+a}S_{k}^{(1)}(l)B_{l+a}(y)}{(z + 1)(z + 2) ⋯ (z + k)} \end{aligned}  (28)

and

ψ(z + 1 − y) = log(z) − ∑_{k=1}^{a} \frac{B_{k}(y)}{kz^{k}} + \frac{1}{z^{a}} ∑_{k=1}^{∞}(−1)^{k+1} \frac{∑_{l=1}^{k} \frac{(−1)^{l}}{l+a} S_{k}^{(1)}(l)B_{l+a}(y)}{(z + 1)(z + 2) ⋯ (z + k)}.  (29)

Proof. Let s ∈ ℂ \ {1, 0, −1, −2, −3, …} be a fixed complex number and let z ∈ ℂ \ (−∞, 0] with | arg(z)| ≤ π − ε < π for some arbitrarily small, but fixed ε > 0. Setting h := 1 − y for 0 ≤ y ≤ 1 into the identities (24) and (26), by using the relation (13) and by exchanging n with n + a we deduce that

ζ(s, z + 1 − y) = \frac{z^{1−s}}{s − 1} + \frac{z^{1−s}}{s − 1} ∑_{k=1}^{n+a}(−1)^{k}\binom{1 − s}{k}\frac{B_{k}(y)}{z^{k}} + O_{n+a}(z)  (30)

and that

ψ(z + 1 − y) = log(z) − ∑_{k=1}^{n+a} \frac{B_{k}(y)}{kz^{k}} + U_{n+a}(z),  (31)

where O_{n+a}(z) and U_{n+a}(z) are as in the previous Lemma 12 with h := 1 − y. We can write the equations (30) and (31) in the form

\begin{aligned}ζ(s, z + 1 − y) & = \frac{z^{1−s}}{s − 1} + \frac{z^{1−s}}{s − 1} ∑_{k=1}^{a}(−1)^{k}\binom{1 − s}{k}\frac{B_{k}(y)}{z^{k}} \\ & + \frac{z^{1−s−a}}{s − 1} ∑_{k=1}^{n}(−1)^{k+a}\binom{1 − s}{k + a}\frac{B_{k+a}(y)}{z^{k}} + O_{n+a}(z) \end{aligned}  (32)

and

ψ(z + 1 − y) = log(z) − ∑_{k=1}^{a} \frac{B_{k}(y)}{kz^{k}} − \frac{1}{z^{a}} ∑_{k=1}^{n} \frac{B_{k+a}(y)}{(k + a)z^{k}} + U_{n+a}(z).  (33)

In the following calculations, we use the fact that the function g(k) := 2^{k} grows faster than any polynomial p(k) as k → ∞.

From equation (32) we get (28) for s ∈ ℂ \ {1, 0, −1, −2, −3, …} by applying Theorem 10 with R_{n}(z) := (s − 1)z^{s+a−1}O_{n+a}(z) to the analytic function f_{1}(z) defined by

\begin{aligned}f_{1}(z) : & = (s − 1)z^{s+a−1} [ζ(s, z + 1 − y) − \frac{z^{1−s}}{s − 1} − \frac{z^{1−s}}{s − 1} ∑_{k=1}^{a}(−1)^{k}\binom{1 − s}{k}\frac{B_{k}(y)}{z^{k}}] \\ & = ∑_{k=1}^{n}(−1)^{k+a}\binom{1 − s}{k + a}\frac{B_{k+a}(y)}{z^{k}} + (s − 1)z^{s+a−1}O_{n+a}(z) \end{aligned}

on z ∈ ℂ \ (−∞, 0] with | arg(z)| ≤ π − ε < π. Because of the relation |x^{s}| = x^{Re(s)} for x ∈ ℝ_{0}^{+} and the use of the identity (8), we have

\begin{aligned}|\binom{1 − s}{k + a}B_{k+a}(y)| & ≤ \frac{2ζ(k + a)(k + a)!(k + a)^{Re(s)−2}}{(2π)^{k+a} |Γ(s − 1)|} + O (\frac{2ζ(k + a)(k + a)!(k + a)^{Re(s)−3}}{(2π)^{k+a}} ) \\ & < \frac{C_{1}(a)k!}{π^{k}} \end{aligned}

and with A_{1} := max {1, e^{Im(s) arg(z)}}, A_{2} := max {1, e^{− Im(s) arg(z)}}, as well as Re(s) > −n, we get that

|(s − 1)z^{s+a−1}O_{n+a}(z)|

≤ 2(n + a + 2)A_{1} |\binom{1 − s}{n + a + 2}| \frac{|B_{n+a+1}| e^{− Im(s) arg(z)} sec^{n+a+Re(s)+1} (\frac{1}{2} arg(z))}{|n + a + Re(s)| · |z|^{n+1}}

≤ \frac{4(n + a + 2)A_{2}}{|n + a + Re(s)|} · \frac{ζ(n + a + 1)(n + a + 1)!(n + a + 2)^{Re(s)−2} sec^{n+a+Re(s)+1} (\frac{1}{2} arg(z))}{(2π)^{n+a+1} |Γ(s − 1)| |z|^{n+1}}

+ O (\frac{4(n + a + 2)A_{2}}{|n + a + Re(s)|} · \frac{ζ(n + a + 1)(n + a + 1)!(n + a + 2)^{Re(s)−3} sec^{n+a+Re(s)+1} (\frac{1}{2} arg(z))}{(2π)^{n+a+1}|z|^{n+1}})

< \frac{C_{2}(a) sec^{n} (\frac{1}{2} arg(z)) n!}{π^{n}|z|^{n+1}}

for some positive constants C_{1}(a), C_{2}(a) depending on a and independent of n. In the last computation above, we have used the relation |z^{s}| = |z|^{Re(s)}e^{− Im(s) arg(z)}.

The above bound for |(s − 1)z^{s+a−1}O_{n+a}(z)| also holds if Re(s) ≤ −n by taking C_{2}(a) large enough, because Re(s) ≤ −n is only possible for finitely many n’s and in each case we have |(s − 1)z^{s+a−1}O_{n+a}(z)| ≤ \frac{C(n)}{|z|^{n+1}} for all n ∈ ℕ and some positive constants C(n).

To get formula (28) also for all s ∈ {0, −1, −2, −3, …}, we apply the Weniger transformation formula (21) directly to the function f_{1}(z) with n := 1−s−a and O_{n+a}(z) = O_{1−s}(z) = 0.

Similarly, from equation (33) we obtain the formula (29) by applying Theorem 10 with R_{n}(z) := z^{a}U_{n+a}(z) to the analytic function f_{2}(z) defined by

f_{2}(z) : = z^{a} [log(z) − ψ(z + 1 − y) − ∑_{k=1}^{a} \frac{B_{k}(y)}{kz^{k}}] = ∑_{k=1}^{n} \frac{B_{k+a}(y)}{(k + a)z^{k}} + z^{a}U_{n+a}(z)

on z ∈ ℂ with | arg(z)| ≤ π − ε < π, because we have

|\frac{B_{k+a}(y)}{k + a}| ≤ \frac{2ζ(k + a)(k + a)!}{(2π)^{k+a}(k + a)} < \frac{C_{3}(a)k!}{π^{k}}

and

|z^{a}U_{n+a}(z)| ≤ \frac{2 |B_{n+a+1}| sec^{n+a+2} (\frac{1}{2} arg(z))}{(n + a + 1)|z|^{n+1}} ≤ \frac{4ζ(n + a + 1) sec^{n+a+2} (\frac{1}{2} arg(z)) (n + a + 1)!}{(2π)^{n+a+1}(n + a + 1)|z|^{n+1}}

< \frac{C_{4}(a) sec^{n} (\frac{1}{2} arg(z)) n!}{π^{n}|z|^{n+1}}

for some positive constants C_{3}(a), C_{4}(a) depending on a and independent of n.

5 The generalized Faulhaber formulas

In this section we prove our generalized versions of Faulhaber’s formula, which all converge very rapidly. For their proofs, we use the above Theorem 13.

Theorem 14. (extended generalized Faulhaber formulas) For every complex number m ∈ ℂ \ {−1} and every positive real number x ∈ ℝ^{+}, we have

∑_{k=1}^{⌊x⌋} k^{m} = \frac{1}{m + 1}x^{m+1} + ζ (−m) + \frac{x^{m+1}}{m + 1} ∑_{k=1}^{∞}(−1)^{k} \frac{∑_{l=1}^{k} \binom{m+1}{l}S_{k}^{(1)}(l)B_{l}(\{x\})}{(x + 1)(x + 2) ⋯ (x + k)}.  (34)

More generally, for every x ∈ ℝ^{+} and every a ∈ ℕ_{0}, we have

\begin{aligned}∑_{k=1}^{⌊x⌋} k^{m} & = \frac{1}{m + 1}x^{m+1} + ζ (−m) + \frac{1}{m + 1} ∑_{k=1}^{a}(−1)^{k}\binom{m + 1}{k}B_{k}(\{x\})x^{m−k+1} \\ & + \frac{x^{m−a+1}}{m + 1} ∑_{k=1}^{∞}(−1)^{k+a} \frac{∑_{l=1}^{k} \binom{m+1}{l+a}S_{k}^{(1)}(l)B_{l+a} (\{x\})}{(x + 1)(x + 2) ⋯ (x + k)} \end{aligned} (35)

and for m = m_{1} + im_{2} ∈ ℂ \ {−1} with m_{1} = Re(m) ≥ −1 the special case

\begin{aligned}∑_{k=1}^{⌊x⌋} k^{m} & = \frac{1}{m + 1}x^{m+1} + ζ (−m) + \frac{1}{m + 1} ∑_{k=1}^{⌊m_{1}+1⌋}(−1)^{k}\binom{m + 1}{k}B_{k}(\{x\})x^{m−k+1} \\ & + (−1)^{⌊m_{1}+1⌋} \frac{x^{m−⌊m_{1}+1⌋+1}}{m + 1} ∑_{k=1}^{∞}(−1)^{k} \frac{∑_{l=1}^{k} \binom{m+1}{l+⌊m_{1}+1⌋}S_{k}^{(1)}(l)B_{l+⌊m_{1}+1⌋}(\{x\})}{(x + 1)(x + 2) ⋯ (x + k)}. \end{aligned} (36)

Moreover, if m = −1, for every positive real number x ∈ ℝ^{+} and every a ∈ ℕ_{0} we have

∑_{k=1}^{⌊x⌋} \frac{1}{k} = log(x) + γ − ∑_{k=1}^{a} \frac{B_{k}(\{x\})}{kx^{k}} + \frac{1}{x^{a}} ∑_{k=1}^{∞}(−1)^{k+1} \frac{∑_{l=1}^{k} \frac{(−1)^{l}}{l+a} S_{k}^{(1)}(l)B_{l+a}(\{x\})}{(x + 1)(x + 2) ⋯ (x + k)}. (37)

In particular, for x ∈ ℝ^{+} we have

∑_{k=1}^{⌊x⌋} \frac{1}{k} = log(x) + γ + ∑_{k=1}^{∞}(−1)^{k+1} \frac{∑_{l=1}^{k} \frac{(−1)^{l}}{l} S_{k}^{(1)}(l)B_{l}(\{x\})}{(x + 1)(x + 2) ⋯ (x + k)}  (38)

and that

∑_{k=1}^{⌊x⌋} \frac{1}{k} = log(x) + γ − \frac{B_{1}(\{x\})}{x} + ∑_{k=1}^{∞}(−1)^{k+1} \frac{∑_{l=1}^{k} \frac{(−1)^{l}}{l+1} S_{k}^{(1)}(l)B_{l+1}(\{x\})}{x(x + 1)(x + 2) ⋯ (x + k)}.  (39)

Proof. From the formula (28) with the parameters s := −m, z := x and y := {x}, we get

∑_{k=1}^{⌊x⌋} k^{m} − ζ(−m) = −ζ(−m, x + 1 − {x})

= \frac{1}{m + 1}x^{m+1} + \frac{1}{m + 1} ∑_{k=1}^{a}(−1)^{k}\binom{m + 1}{k}B_{k}(\{x\})x^{m−k+1}

+ \frac{x^{m−a+1}}{m + 1} ∑_{k=1}^{∞}(−1)^{k+a} \frac{∑_{l=1}^{k} \binom{m+1}{l+a}S_{k}^{(1)}(l)B_{l+a} (\{x\})}{(x + 1)(x + 2) ⋯ (x + k)}

by using the formula (19) with n := ⌊x⌋ = x − {x} in the first step. This gives the above identity (35) with its special cases (34) and (36).

Similarly, we now use the formula (29) again with the variables z := x and y := {x}, and then we get

\begin{aligned}∑_{k=1}^{⌊x⌋} \frac{1}{k} − γ & = ψ(x + 1 − \{x\}) \\ & = log(x) − ∑_{k=1}^{a} \frac{B_{k}(\{x\})}{kx^{k}} + \frac{1}{x^{a}} ∑_{k=1}^{∞}(−1)^{k+1} \frac{∑_{l=1}^{k} \frac{(−1)^{l}}{l+a} S_{k}^{(1)}(l)B_{l+a}(\{x\})}{(x + 1)(x + 2) ⋯ (x + k)} \end{aligned}

by employing the formula (16) with n := ⌊x⌋ = x − {x} in the first line of the above calculation. This gives the above identity (37) with its special cases (38) and (39).

By setting x := n ∈ ℕ into Theorem 14, we obtain the following corollary.

Corollary 15. (generalized Faulhaber formulas) For every complex number m ∈ ℂ \ {−1} and every natural number n ∈ ℕ, we have

∑_{k=1}^{n} k^{m} = \frac{1}{m + 1}n^{m+1} + ζ (−m) + \frac{n^{m+1}}{m + 1} ∑_{k=1}^{∞} \frac{(−1)^{k} ∑_{l=1}^{k} \binom{m+1}{l}B_{l}S_{k}^{(1)}(l)}{(n + 1)(n + 2) ⋯ (n + k)}  (40)

and more generally when a ∈ ℕ_{0} that

\begin{aligned}∑_{k=1}^{n} k^{m} & = \frac{1}{m + 1}n^{m+1} + ζ (−m) + \frac{1}{m + 1} ∑_{k=1}^{a}(−1)^{k}\binom{m + 1}{k}B_{k}n^{m−k+1} \\ & + \frac{n^{m−a+1}}{m + 1} ∑_{k=1}^{∞} \frac{(−1)^{k+a} ∑_{l=1}^{k} \binom{m+1}{l+a}B_{l+a}S_{k}^{(1)}(l)}{(n + 1)(n + 2) ⋯ (n + k)}. \end{aligned}  (41)

We have again when m = m_{1} + im_{2} ∈ ℂ \ {−1} with m_{1} = Re(m) ≥ −1 the special case

\begin{aligned}∑_{k=1}^{n} k^{m} & = \frac{1}{m + 1}n^{m+1} + ζ (−m) + \frac{1}{m + 1} ∑_{k=1}^{⌊m_{1}+1⌋}(−1)^{k}\binom{m + 1}{k}B_{k}n^{m−k+1} \\ & + (−1)^{⌊m_{1}+1⌋} \frac{n^{m−⌊m_{1}+1⌋+1}}{m + 1} ∑_{k=1}^{∞}(−1)^{k} \frac{∑_{l=1}^{k} \binom{m+1}{l+⌊m_{1}+1⌋}B_{l+⌊m_{1}+1⌋}S_{k}^{(1)}(l)}{(n + 1)(n + 2) ⋯ (n + k)}. \end{aligned}  (42)

For m = −1, for every natural number n ∈ ℕ and every a ∈ ℕ_{0} we have

∑_{k=1}^{n} \frac{1}{k} = log(n) + γ − ∑_{k=1}^{a} \frac{B_{k}}{kn^{k}} + \frac{1}{n^{a}} ∑_{k=1}^{∞} \frac{(−1)^{k+1} ∑_{l=1}^{k} \frac{(−1)^{l}}{l+a} B_{l+a}S_{k}^{(1)}(l)}{(n + 1)(n + 2) ⋯ (n + k)}.  (43)

In particular, for every n ∈ ℕ we have

\begin{aligned}∑_{k=1}^{n} \frac{1}{k} & = log(n) + γ + ∑_{k=1}^{∞} \frac{(−1)^{k+1} ∑_{l=1}^{k} \frac{(−1)^{l}}{l} B_{l}S_{k}^{(1)}(l)}{(n + 1)(n + 2) ⋯ (n + k)} \\ & = log(n) + γ + \frac{1}{2(n + 1)} + \frac{5}{12(n + 1)(n + 2)} + \frac{3}{4(n + 1)(n + 2)(n + 3)} \\ & + \frac{251}{120(n + 1)(n + 2)(n + 3)(n + 4)} + ⋯ \end{aligned}  (44)

and

\begin{aligned}∑_{k=1}^{n} \frac{1}{k} & = log(n) + γ + \frac{1}{2n} + ∑_{k=1}^{∞} \frac{(−1)^{k} ∑_{l=1}^{k} \frac{B_{l+1}}{l+1} S_{k}^{(1)}(l)}{n(n + 1)(n + 2) ⋯ (n + k)} \\ & = log(n) + γ + \frac{1}{2n} − \frac{1}{12n(n + 1)} − \frac{1}{12n(n + 1)(n + 2)} − \frac{19}{120n(n + 1)(n + 2)(n + 3)} \\ & − \frac{9}{20n(n + 1)(n + 2)(n + 3)(n + 4)} − · · · . \end{aligned} ; (45)

For every positive real number x ∈ ℝ^{+} and for every natural number n ∈ ℕ, we list the following 8 most used generalized Faulhaber summation formulas:

1.) Generalized Faulhaber formula for the partial sums of ζ(2):

For every natural number n ∈ ℕ, we have

\begin{aligned}∑_{k=1}^{n} \frac{1}{k^{2}} & = ζ(2) − \frac{1}{n} + ∑_{k=1}^{∞} \frac{(−1)^{k+1} ∑_{l=1}^{k}(−1)^{l}B_{l}S_{k}^{(1)}(l)}{n(n + 1)(n + 2) ⋯ (n + k)} \\ & = ζ(2) − \frac{1}{n} + ∑_{k=1}^{∞} \frac{1}{k + 1} · \frac{(k − 1)!}{n(n + 1)(n + 2) ⋯ (n + k)} \\ & = ζ(2) − \frac{1}{n} + \frac{1}{2n(n + 1)} + \frac{1}{3n(n + 1)(n + 2)} + \frac{1}{2n(n + 1)(n + 2)(n + 3)} \\ & + \frac{6}{5n(n + 1)(n + 2)(n + 3)(n + 4)} + · · · . \end{aligned} (46)

2.) Extended generalized Faulhaber formula for the partial sums of ζ(3): For every real number x ∈ ℝ^{+}, we obtain

∑_{k=1}^{⌊x⌋} \frac{1}{k^{3}} = ζ(3) − \frac{1}{2x^{2}} + \frac{1}{2x} ∑_{k=1}^{∞}(−1)^{k+1} \frac{∑_{l=1}^{k}(−1)^{l}(l + 1)S_{k}^{(1)}(l)B_{l}(\{x\})}{x(x + 1)(x + 2) ⋯ (x + k)}.  (47)

3.) Extended generalized Faulhaber formula for the sum of the square roots: For every real number x ∈ ℝ^{+}, we get

\sqrt[∑_{k=1}^{⌊x⌋}]{k} = \frac{2}{3}x^{3/2} − \frac{1}{4π} ζ (\frac{3}{2}) + x\sqrt{x} ∑_{k=1}^{∞}(−1)^{k} \frac{∑_{l=1}^{k} \frac{(−1)^{l}(2l−5)!!}{2^{l−1}l!}S_{k}^{(1)}(l)B_{l}(\{x\})}{(x + 1)(x + 2) ⋯ (x + k)}. (48)

4.) Generalized Faulhaber formula for the partial sums of ζ(−3/2): For every natural number n ∈ ℕ, we have

\begin{matrix}∑_{k=1}^{n} k\sqrt{k} \\ = \frac{2}{5}n^{5/2} + \frac{1}{2}n^{3/2} + \frac{1}{8}\sqrt{n} − \frac{3}{16π^{2}} ζ (\frac{5}{2}) + 3\sqrt{n} ∑_{k=1}^{∞}(−1)^{k+1} \frac{∑_{l=1}^{k} \frac{(2l−3)!!}{2^{l+1}(l+2)!}B_{l+2}S_{k}^{(1)}(l)}{(n + 1)(n + 2) ⋯ (n + k)} \\ = \frac{2}{5}n^{5/2} + \frac{1}{2}n^{3/2} + \frac{1}{8}\sqrt{n} − \frac{3}{16π^{2}} ζ (\frac{5}{2}) + \frac{\sqrt{n}}{1920(n + 1)(n + 2)} \\ + \frac{\sqrt{n}}{640(n + 1)(n + 2)(n + 3)} + \frac{611\sqrt{n}}{107520(n + 1)(n + 2)(n + 3)(n + 4)} \\ + \frac{275\sqrt{n}}{10752(n + 1)(n + 2)(n + 3)(n + 4)(n + 5)} \\ + \frac{159157\sqrt{n}}{1146880(n + 1)(n + 2)(n + 3)(n + 4)(n + 5)(n + 6)} + · · · . \end{matrix} ; (49)

5.) Generalized Faulhaber formula for the partial sums of ζ(−5/2):

For every natural number n ∈ ℕ, we obtain that

\begin{matrix}∑_{k=1}^{n} k^{2}\sqrt{k} \\ = \frac{2}{7}n^{7/2} + \frac{1}{2}n^{5/2} + \frac{5}{24}n^{3/2} + \frac{15}{64π^{3}} ζ (\frac{7}{2}) + 15\sqrt{n} ∑_{k=1}^{∞}(−1)^{k+1} \frac{∑_{l=1}^{k} \frac{(2l−3)!!}{2^{l+2}(l+3)!}B_{l+3}S_{k}^{(1)}(l)}{(n + 1)(n + 2) ⋯ (n + k)} \\ = \frac{2}{7}n^{7/2} + \frac{1}{2}n^{5/2} + \frac{5}{24}n^{3/2} + \frac{15}{64π^{3}} ζ (\frac{7}{2}) − \frac{\sqrt{n}}{384(n + 1)} − \frac{\sqrt{n}}{384(n + 1)(n + 2)} \\ − \frac{37\sqrt{n}}{7168(n + 1)(n + 2)(n + 3)} − \frac{55\sqrt{n}}{3584(n + 1)(n + 2)(n + 3)(n + 4)} \\ − \frac{1995\sqrt{n}}{32768(n + 1)(n + 2)(n + 3)(n + 4)(n + 5)} − · · · . \end{matrix} ; (50)

6.) Generalized Faulhaber formula for the sum of the inverses of the square roots: For every natural number n ∈ ℕ, we get that

\begin{aligned}∑_{k=1}^{n} \frac{1}{\sqrt{k}} & = 2\sqrt{n} + ζ (\frac{1}{2}) + \frac{1}{2\sqrt{n}} + \frac{1}{\sqrt{n}} ∑_{k=1}^{∞}(−1)^{k} \frac{∑_{l=1}^{k} \frac{(2l−1)!!}{2^{l}(l+1)!}B_{l+1}S_{k}^{(1)}(l)}{(n + 1)(n + 2) ⋯ (n + k)} \\ & = 2\sqrt{n} + ζ (\frac{1}{2}) + \frac{1}{2\sqrt{n}} − \frac{1}{24\sqrt{n}(n + 1)} − \frac{1}{24\sqrt{n}(n + 1)(n + 2)} \\ & − \frac{31}{384\sqrt{n}(n + 1)(n + 2)(n + 3)} − \frac{15}{64\sqrt{n}(n + 1)(n + 2)(n + 3)(n + 4)} − · · · . \end{aligned} ; (51)

7.) Extended generalized Faulhaber formula for the partial sums of ζ(3/2): For every real number x ∈ ℝ^{+}, we have

\begin{aligned}∑_{k=1}^{⌊x⌋} \frac{1}{k\sqrt{k}} & = ζ (\frac{3}{2}) − \frac{2}{\sqrt{x}} − \frac{B_{1}(\{x\})}{x\sqrt{x}} + \frac{2}{\sqrt{x}} ∑_{k=1}^{∞}(−1)^{k+1} \frac{∑_{l=1}^{k} \frac{(−1)^{l}(2l+1)!!}{2^{l+1}(l+1)!}S_{k}^{(1)}(l)B_{l+1}(\{x\})}{x(x + 1)(x + 2) ⋯ (x + k)}. \\ & (52) \end{aligned}

8.) Extended generalized Faulhaber formula for the partial sums of ζ(5/2): For every real number x ∈ ℝ^{+}, we obtain

∑_{k=1}^{⌊x⌋} \frac{1}{k^{2}\sqrt{k}} = ζ (\frac{5}{2}) − \frac{2}{3x^{3/2}} + \frac{4}{3\sqrt{x}} ∑_{k=1}^{∞}(−1)^{k+1} \frac{∑_{l=1}^{k} \frac{(−1)^{l}(2l+1)!!}{2^{l+1}l!}S_{k}^{(1)}(l)B_{l}(\{x\})}{x(x + 1)(x + 2) ⋯ (x + k)}. (53)

In equations 3.), 4.) and 5.), we have used that ζ (−\frac{1}{2}) = −\frac{1}{4π}ζ (\frac{3}{2}), ζ (−\frac{3}{2}) = −\frac{3}{16π^{2}} ζ (\frac{5}{2}) and that ζ (−\frac{5}{2}) = \frac{15}{64π^{3}} ζ (\frac{7}{2}), which follow from the functional equation of the Riemann zeta function [14, (11.05), p. 63] and are also given in [23, (2), (5) and (9)].

6 Conclusion

We have proved a rapidly convergent generalization of Faulhaber’s formula to sums of arbitrary complex powers m ∈ ℂ. In our eyes, these formulas are useful because of their rapid convergence. We believe that they may also have applications in physics [18], such as the extended version of Faulhaber’s formula [19, 20]. With the universal technique, explained in this paper, one can obtain other summation formulas of this type [21, 22]. A generalization of Faulhaber’s formula for alternating sums can be found in [22].

7 Acknowledgment

The author thanks the referee for the valuable suggestions regarding the improvement of this paper.

This work was supported by SNSF (Swiss National Science Foundation) under grant 169247.

References

[1] Kevin J. McGown and Harold R. Parks, The generalization of Faulhaber’s formula to sums of non-integral powers, J. Math. Anal. Appl. 330 (2007), 571–575.

[2] Ernst Joachim Weniger, Summation of divergent power series by means of factorial series, Appl. Num. Math. 60 (2010), 1429–1441.

[3] Oskar Schlömilch, Ueber Facultätenreihen, Zeitschrift für Mathematik und Physik 4 (1859), 390–415.

[4] Oskar Schlömilch, Relationen zwischen den Facultätencoefficienten, Arch. für Math. und Phys. 9 (1847), 333–335.

[5] Oskar Schlömilch, Über die independente Bestimmung der Coëfficienten unendlicher Reihen und der Facultätencoefficienten insbesondere, Arch. für Math. und Phys. 18 (1852), 306–327.

[6] I. Tweedle, James Stirling’s Methodus Differentialis: An Annotated Translation of Stirling’s Text, Springer, 2003.

[7] A. C. Aitken, A note on inverse central factorial series, Proc. Edinburgh Math. Soc. 7 (1946), 168–170.

[8] Donatella Merlini, Renzo Sprugnoli, and M. Cecilia Verri, The Cauchy numbers, Discrete Math. 306 (2006), 1906–1920.

[9] Frank W. J. Olver, Daniel W. Lozier, Ronald F. Boisvert, and Charles W. Clark, NIST Handbook of Mathematical Functions, Cambridge University Press, 2010.

[10] N. J. A. Sloane et al., The Online Encyclopedia of Integer Sequences. Published electronically at https://oeis.org, 2022.

[11] Yudell L. Luke, The Special Functions and their Approximations, Academic Press, 1969.

[12] Victor Kac and Pokman Cheung, Quantum Calculus, Springer, 2002.

[13] Hugh L. Montgomery and Robert C. Vaughan, Multiplicative Number Theory. I. Classical Theory, Cambridge University Press, 2007.

[14] Frank W. J. Olver, Asymptotics and Special Functions, A. K. Peters Ltd., 1997.

[15] G. Nemes, Error bounds for the asymptotic expansion of the Hurwitz zeta function, Proc. A. 473 (2017), 1–16.

[16] G. N. Watson, The transformation of an asymptotic series into a convergent series of inverse factorials, Rend. Circ. Mat. Palermo 34 (1912), 41–88.

[17] Niels Nielsen, Handbuch der Theorie der Gammafunktion, Teubner, 1906.

[18] Hanno Sahlmann and Thomas Zilker, Quantum surface holonomies for loop quantum gravity and their application to black hole horizons, Phys. Rev. D 102 (2020), 026009: 1–33.

[19] Antônio F. Neto, A note on a theorem of Schumacher, J. Integer Sequences 19 (2016), Article 16.8.5.

[20] Raphael Schumacher, An extended version of Faulhaber’s formula, J. Integer Sequences 19 (2016), Article 16.4.2.

[21] Raphael Schumacher, Rapidly convergent summation formulas involving Stirling series,

   preprint, 2016. Available at https://arxiv.org/abs/1602.00336.

[22] Raphael Schumacher, Extension of summation formulas involving Stirling series,

   preprint, 2016. Available at https://arxiv.org/abs/1605.09204.

[23] Srinivasa Ramanujan, On the sum of the square roots of the first n natural numbers, J. Indian Math. Soc. 7 (1915), 173–175.

2020 Mathematics Subject Classification: Primary 65B15; Secondary 11B68. Keywords: generalization of Faulhaber’s formula, extended Faulhaber formula, finite Weniger transformation, Stirling number of the first kind, Bernoulli polynomial, Bernoulli number.

(Concerned with sequences \underline{A000142}, \underline{A008275}, \underline{A027641} and \underline{A027642}.)

Received March 11 2021; revised versions received March 14 2021; December 14 2021; March 6 2022; March 8 2022. Published in Journal of Integer Sequences, March 25 2022.

Return to Journal of Integer Sequences home page.