The accepted answer takes $O(kn)$ multiplications and additions. There is also a method using $O(k\log n)$ additions, multiplications and divisions, using the below identity.
$$
S(n,k)=\frac1{k!}\sum_{j=0}^k(-1)^j\binom{k}j(k-j)^n
$$
We can use DLMF Equation 26.3.E6 as a recurrence relation to greatly reduce the amount of computation necessary to compute the binomial coefficient.
$$
\binom{m}{n} = \frac{m-n+1}{n}\binom{m}{n-1}
$$
Pseudocode:
accum ← 0
k_choose_j ← 1
for j = 0,1,...,k-1:
accum ← accum + (-1)^j * k_choose_j * (k - j)^n
k_choose_j ← k_choose_j * (k-j) / (j+1)
return accum / factorial(k)
Note that in the above pseudocode, we update k_choose_j
at the end of each iteration, so we need to substitute j + 1
for n
in the above relation.