Yo veía esto en OEIS A008292 y
De Wikipedia.
Al parecer
$$k^m = \sum_{q=0}^{m-1} {k+q\elegir m}
\langle {m\cima de q}\rangle$$
con $\langle {m\atop q}\rangle$ la Euleriano números que han
la generación de la función
$$\frac{t-1}{t-\exp((t-1)z)}.$$
Ahora podemos probar esta suma fórmula.
Tenemos los siguientes exponencial de generación de función
$$G_k(z) = \sum_{m\ge 0} k^m \frac{z^m}{m!} = \exp(kz).$$
Por otro lado, la fórmula de la suma da el EGF
$$H_k(z) = \sum_{m\ge 0}
\frac{z^m}{m.} \sum_{q=0}^{m-1} {k+q\elegir m}
\langle {m\cima de q}\rangle.$$
Podemos extender este a $q=m$ debido a que el número de Euler
$\langle {m\atop m}\rangle$ es cero para obtener
$$\sum_{m\ge 0}
\frac{z^m}{m.} \sum_{q=0}^{m} {k+q\elegir m}
\langle {m\cima de q}\rangle
\\ = \sum_{q\ge 0} \sum_{m\ge p}
{k+q\elegir m} \frac{z^m}{m.}
\langle {m\cima de q}\rangle.$$
Ahora introducir
$${k+q\elegir m}
= \frac{1}{2\pi i}
\int_{|w|=\epsilon} \frac{(1+w)^{k+p}}{w^{m+1}} \; ps.$$
Sustituir esto en la suma para obtener
$$\frac{1}{2\pi i}
\int_{|w|=\epsilon}
\sum_{q\ge 0} \frac{(1+w)^{k+p}}{w}
\sum_{m\ge q} \frac{1}{w^m} \frac{z^m}{m.}
\langle {m\cima de q}\rangle \; dw
\\ = \frac{1}{2\pi i}
\int_{|w|=\epsilon} \frac{(1+w)^{k}}{w}
\sum_{q\ge 0} (1+w)^q
\sum_{m\ge q} \frac{1}{w^m} \frac{z^m}{m.}
\langle {m\cima de q}\rangle \; ps.$$
Ahora lo que tenemos aquí es un doble
aniquilado el coeficiente de extractor
la que ahora nos colapso:
$$\frac{1}{2\pi i}
\int_{|w|=\epsilon} \frac{(1+w)^{k}}{w}
\sum_{q\ge 0} (1+w)^q
[t^q] \frac{t-1}{t\exp((t-1)z/w)} \; dw
\\ = \frac{1}{2\pi i}
\int_{|w|=\epsilon} \frac{(1+w)^{k}}{w}
\frac{w}{1+w-\exp(z)} \; dw
\\ = \frac{1}{2\pi i}
\int_{|w|=\epsilon} \frac{(1+w)^{k}}{1+w-\exp(z)} \; ps.$$
La contribución de la pole en $w=\exp(z)-1$ que es simple es
precisamente
$$(1+(\exp(z)-1))^k = \exp(kz) = G_k(z),$$
QED.
Hay otro aniquilado coeficiente de extractor en este
MSE enlace y en este
MSE vínculo II y también aquí, en este MSE enlace III.
El Arce código para la inicial de búsqueda en la OEIS fue como sigue:
Q :=
proc(m)
local s, sys, sol;
s := expand(añadir(a[q]*binomial(k+q,m), q=0..m-1));
sys := [coef(s, k, m)=1];
sys :=
[op(sys), seq(coef(s, k, q)=0, q=0..m-1)];
sol := solve(sys, [seq ([p], p=0..m-1)]);
subs(sol[1], [seq ([p], p=0..m-1)]);
end;