令\[S_i=\sum_{k=1}^n k^i m^k\]我们有\[\begin{eqnarray*}(m-1)S_i & = & mS_i - S_i \\& = & \sum_{k=1}^n k^i m^{k+1} - \sum_{k=1}^n k^i m^k \\& = & \sum_{k=2}^{n+1} (k-1)^i m^k - \sum_{k=1}^n k^i m^k \\& = & n^i m^{n+1} + \sum_{k=1}^n m^k\big( (k-1)^i - k^i\big) \\& = & n^i m^{n+1} + \sum_{k=1}^n \bigg( \sum_{j=1}^{i-1} (-1)^{i-j}{i \choose j}k^j m^k \bigg) \\& = & n^i m^{n+1} + \sum_{j=0}^{i-1} (-1)^{i-j}{i \choose j} S_j\end{eqnarray*}\]直接按照这条式子\(O(m^2)\)递推即可。
P.S. 来自http://taorunz.logdown.com/posts/193397-bzoj3157
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long LL;
const LL mod=;
inline LL sqr(LL num){return num*num;}
LL mexp(LL a,LL b,LL p){return b?sqr(mexp(a,b>>,p))%p*((b&)?a:)%p:;}
LL fac[],facinv[],powm1[];
LL C(LL n,LL r){return fac[n]*facinv[r]%mod*facinv[n-r]%mod;}
LL n,m,s[];
int main(int argc, char *argv[])
{
cin>>n>>m;
if(m==){cout<<n*(n+)/%mod<<endl;return ;}
fac[]=;for(int i=;i<=m;i++)fac[i]=fac[i-]*i%mod;
facinv[m]=mexp(fac[m],mod-,mod);facinv[]=;
for(int i=m-;i>=;i--)facinv[i]=facinv[i+]*(i+)%mod;
powm1[]=;for(int i=;i<=m;i++)powm1[i]=powm1[i-]*-;
s[]=((mexp(m,n+,mod)-m)%mod+mod)%mod*mexp(m-,mod-,mod);
for(int i=;i<=m;i++)
{
s[i]=mexp(n,i,mod)*mexp(m,n+,mod)%mod;
for(int j=;j<i;j++)s[i]=((s[i]+powm1[i-j]*C(i,j)*s[j])%mod+mod)%mod;
s[i]=s[i]*mexp(m-,mod-,mod)%mod;
}
cout<<s[m]<<endl;
return ;
}