# Summation of an expression

Hi,

Does anyone have an idea about an exact or approximate formulae for the following summation? $$\sum_{j=1}^n \frac{j^k}{(j-1)!}$$ where k is a positive integer (the denominator of the j^th term is of course $\Gamma(j)$).

In retrospect, deriving the approximation is quite easy. Indeed, $$\sum\limits_{j = 1}^n {\frac{{j^k }}{{(j - 1)!}}} = e\sum\limits_{j = 1}^n {e^{ - 1} \frac{{j^{k + 1} }}{{j!}}} \approx e\sum\limits_{j = 0}^\infty {e^{ - 1} \frac{{j^{k + 1} }}{{j!}}} = eB_{k + 1}.$$ (For large $k$, you may consider Asymptotic limit and bounds.)

Assume that $n$ is sufficiently large. Since $$\sum\limits_{j = 1}^n {\frac{{j^k }}{{(j - 1)!}}} = \sum\limits_{j = 0}^{n - 1} {\frac{{(j + 1)^k }}{{j!}}} ,$$ the problem reduces to approximating $\sum\nolimits_{j = 0}^{n - 1} {\frac{{j^m }}{{j!}}}$, for $0 \leq m \leq k$. Now, $e^{ - 1} \sum\nolimits_{j = 0}^\infty {\frac{{j^m }}{{j!}}}$ is the $m$-th moment of the Poisson distribution with mean $1$. For the latter, see, for example, this.

EDIT: Specifically, the approximation (with the right-hand side being an upper bound) is $$\sum\limits_{j = 1}^n {\frac{{j^k }}{{(j - 1)!}}} \approx e\sum\limits_{m = 0}^k {{k \choose m}B_m },$$ where $B_m$ is the $m$-th Bell number (a list is given here). Numerical results indicate that this approximation is very accurate, even for moderate values of $n$. For example, the absolute error for $n=15$, $k=3$ is $\approx 3.4 \times 10^{-9}$; for $n=20$, $k=5$ is $\approx 1.8 \times 10^{-12}$; for $n=20$, $k=7$ is $\approx 7.9 \times 10^{-10}$; for $n=23$, $k=9$ is $\approx 5.8 \times 10^{-11}$.

EDIT: If only the relative error is concerned, then the approximation is quite accurate even for relatively small $n$ values (of course, depending on $k$). For example, for $n=8$, $k=4$: $\approx 141.156$ compared to $\approx 141.351$; for $n=9$, $k=5$: $\approx 551.484$ compared to $\approx 551.811$; for $n=10$, $k=6$: $\approx 2383.359$ compared to $\approx 2383.933$.

EDIT: As Mike Spivey observed, using the identity $$\sum\limits_{m = 0}^k {{k \choose m}B_m } = B_{k + 1},$$ the above approximation can be simplified greatly to $$\sum\limits_{j = 1}^n {\frac{{j^k }}{{(j - 1)!}}} \approx e B_{k+1}.$$

• Shai, $\sum_{m=0}^k \binom{k}{m} B_m = B_{k+1}$. See the first identity under en.wikipedia.org/wiki/Bell_number#Properties_of_Bell_numbers. So this will simplify your answer greatly. – Mike Spivey Feb 15 '11 at 6:25
• @Mike: Great observation, thanks! I'll edit the answer soon. – Shai Covo Feb 15 '11 at 6:37
• I think you can improve on your approximation by subtracting a shifted infinite sum from the unshifted one, then using the binomial theorem and the original approximation. – Steve Huntsman Feb 15 '11 at 11:28
• @Steve and all: You are welcome to edit my answer. – Shai Covo Feb 15 '11 at 15:47

There might be something in Hall and Knight, Higher Algebra. On page 333, they have the formula, $\sum_1^{n+2}{j^2+j-1\over(j+2)!}={1\over2}-{1\over(n+2)(n!)}$.

The infinite sum is a multiple of $e\$ and a Bell number. So find the discrepancy (which will be very small), maybe on order of $e^{-2n}$).

• I discovered this by experiment with mathematica (the sum from 0 to $\infty$) -- what's the reference? – Igor Rivin Feb 15 '11 at 0:47
• But for the discrepancy to be small, $n$ has to be very large (crudely, on the order of $k.$) – Igor Rivin Feb 15 '11 at 0:50
• A reference is en.wikipedia.org/wiki/Bell_number see Dobinski's formula. – Gerry Myerson Feb 15 '11 at 0:52

It should be possible to derive the asymptotics of the expression automatically:

1) let the computer guess and prove a recurrence for the expression, which is P-recursive in $n$.

2) use Doron Zeilberger's package asyrec for Maple to obtain the asymptotics (but not the constant term, unfortunately)

3) try to guess the constant term...

I did only part one, for your convenience I post the complete code. (in FriCAS, Bruno Salvy's gfun for Maple and Manuel Kauers' guess for Mathematica should work just as well. In fact, FriCAS is not particularly well suited because it does not yet support multivariate guessing...)

Setup:

)expose RECOP
f := operator 'f; ops := [N^l for l in 0..2]; vars := [f(n+l) for l in 0..2];


Guess the recurrences for each $k$ separately (and transform the results into operator notation for easier postprocessing:

g(n,k) == reduce(+, [j^k/factorial(j-1) for j in 1..n], 0)

r := [eval(getEq first guessPRec([g(n,k) for n in 0..100], maxShift==2), vars, ops)::UP(N, POLY INT)::UP(N, FR POLY INT) for k in 1..15]


The result is (showing only $k=1..4$)

           2 2     2
[(n + 1) N  - (n  + 3n + 3)N + n + 2,
3 2     3     2                     2
(n + 1) N  - (n  + 4n  + 7n + 5)N + (n + 2) ,
4 2     4     3      2                      3
(n + 1) N  - (n  + 5n  + 12n  + 16n + 9)N + (n + 2) ,
5 2     5     4      3      2                       4
(n + 1) N  - (n  + 6n  + 18n  + 34n  + 37n + 17)N + (n + 2) ,


So the coefficients of $N^2$ (i.e., $f(n+2)$) and $N^0$ (i.e., $f(n)$) are obvious. The coefficient of $N$ turns out to have a rational generating function:

guessPade([coefficient(t, 1) for t in r], indexName=='k)


resulting in

                   2      2
k (n + 1)(n + 2) x - (n  + 3n + 3)
[[x ]--------------------------------]
2
(n + 1)(n + 2)x  - (2n + 3)x + 1


You can use Wolfram Alpha to get exact solutions for different k, although the results are expressed in not-so-good functions

This allows to build smooth plots:

(both plots are for k=5)