I want to evaluate the sum:
$$E(a,n) = \sum_{i=a}^n S(i,a)$$
with $S(n,k)$ denoting the Stirling numbers of the second kind. My problem is that as $n$ gets bigger this calculation gets quite computationally-expensive with either using the recursive definition
$$S(n+1,k) = kS(n,k)+S(n, k-1)$$
or the explicit formula
$$S(n,k)=\frac{1}{k!} \sum_{j=0}^k (-1)^{k-j} {k \choose j} j^n$$
for $S$.
So my question arises: is there a closed-form expression for $E$? If not, are there some sorts of other improvements to reduce the computational complexity of this sum?
For fixed $a$ the function $S(i, a)$ is a sum of various powers, so the partial sum over $i$ can be computed with a termwise sum of geometric series. That is, we have
$$E(n, a) = \sum_{i=a}^n \frac{1}{a!} \sum_{j=0}^a (-1)^{a-j} {a \choose j} j^i$$
(I am switching $a$ and $n$ in the input for $E$ to maintain consistency with $a$ being on the right for the Stirling numbers) and exchanging the order of the sum, as well as removing the $\frac{1}{a!}$, gives
$$E(n, a) = \frac{1}{a!} \sum_{j=0}^a (-1)^{a-j} {a \choose j} \sum_{i=a}^n j^i.$$
The sum $\sum_{i=a}^n j^i$ is a geometric series $\frac{j^{n+1} - j^a}{j - 1}$, so this gives
$$E(n, a) = \frac{1}{a!} \sum_{j=0}^a (-1)^{a-j} {a \choose j} \frac{j^{n+1} - j^a}{j - 1}.$$
For fixed $a$ most of this expression doesn't depend on $n$ and so can be precomputed (I don't know if that's the case in your application though), and then we only need to compute the powers $j^{n+1}$. If you're working mod a prime this can be done very efficiently using binary exponentiation and Fermat's little theorem.
Furthermore it sounds like in your application the prime is bigger than $a$, which allows us to make use of the observation that we can cancel the $a!$ factor from the binomial coefficients:
$$E(n, a) = \sum_{j=0}^a (-1)^{a-j} \frac{1}{j! (a - j)!} \frac{j^{n+1} - j^a}{j - 1}.$$
(The point being that we can compute $j! (a - j)! \bmod p$ and then invert it, without needing to keep track of a potentially large denominator that eventually cancels.) This should be a big save overall; we no longer have to perform any operation $n$ times.
Formally, to compute $j^{n+1} \bmod p$ we first reduce the exponent $\bmod (p-1)$; I'm not familiar with fast algorithms for dividing with remainder so I don't know how long this takes but probably not very long. This gives us a new exponent $m < p-1$ and then binary exponentiation will let us compute $j^m \bmod p$ in about $O(\log m) \le O(\log p)$ steps, and since $p = 10^9 + 7$ here this is not only constant but a pretty small constant!
Similar improvements are probably not available for the binomial coefficients / factorials but we still only have to do about $O(a)$ multiplications and we can work recursively and so forth, and computing inverses $\bmod p$ can be done efficiently using the Euclidean algorithm.