Generalization of Faulhaber's formula

Is there a way to calculate a sum of non-integer positive powers of integers?

$\sum_{k=1}^nk^p: n \in \mathbb{N}, p \in \mathbb{R^+}$

There's a Faulhaber's formula, but as far as I can see, it is applicable only to integer powers. Is it possible to generalize it w/o getting too complex computation?

The thing is the formula needs to be calculated on a computer and if a solution involves calculating integrals up to infinity it might be simpler to calculate the sum directly.

$n$ might be up to $10^{12}$, approximations are also an option.

This is related to the (Hurwitz) Zeta function or generalized harmonic numbers $$\sum_{k=1}^nk^p = \zeta(-p) -\zeta(-p, 1 + n) = H_n^{(-p)}$$

Example: for $n=10^9, p=0.1$ I get with double precision functions the value $7.2211657737752\cdot 10^9$ or for $n=10^{12},\; p=1/4$ you have $2.8109134757068\cdot 10^{13}$

These numbers are computed with my Pascal routines, you can find C functions in the GSL or Cephes libraries, with Python there is mpmath.

• I feel like this just comes from the definition of those special functions though. Commented Nov 2, 2017 at 16:34
• @Ethan: Yes of course, but you can use them to actually compute the sums without adding $10^9$ terms, see my example. Commented Nov 2, 2017 at 16:38

For integer $p$, the first Faulhaber terms are

$$\frac{n^{p+1}}{p+1},\frac{n^p}2,\frac{p{n^{p-1}}}{12},0,\cdots$$

For large $n$, these terms are quickly decreasing and I wouldn't be surprised that you can simply plug fractional values of $p$ to get precise estimates.

(I have no serious justification.)

• If you add the three terms you get an error of $0.417$ for my $10^p, 0.1$ example, i.e. relative error of $0.58\cdot10^{-10}.$ A justificatiion might come from the representations by the Euler–Maclaurin formula dlmf.nist.gov/25.2#iii and dlmf.nist.gov/25.11#iii Commented Nov 2, 2017 at 19:05