令\[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
1 #include <iostream> 2 #include <cstdio> 3 #include <cstring> 4 #include <algorithm> 5 using namespace std; 6 typedef long long LL; 7 const LL mod=1000000007; 8 inline LL sqr(LL num){return num*num;} 9 LL mexp(LL a,LL b,LL p){return b?sqr(mexp(a,b>>1,p))%p*((b&1)?a:1)%p:1;} 10 LL fac[1100],facinv[1100],powm1[1100]; 11 LL C(LL n,LL r){return fac[n]*facinv[r]%mod*facinv[n-r]%mod;} 12 LL n,m,s[1100]; 13 int main(int argc, char *argv[]) 14 { 15 cin>>n>>m; 16 if(m==1){cout<<n*(n+1)/2%mod<<endl;return 0;} 17 fac[0]=1;for(int i=1;i<=m;i++)fac[i]=fac[i-1]*i%mod; 18 facinv[m]=mexp(fac[m],mod-2,mod);facinv[0]=1; 19 for(int i=m-1;i>=1;i--)facinv[i]=facinv[i+1]*(i+1)%mod; 20 powm1[0]=1;for(int i=1;i<=m;i++)powm1[i]=powm1[i-1]*-1; 21 s[0]=((mexp(m,n+1,mod)-m)%mod+mod)%mod*mexp(m-1,mod-2,mod); 22 for(int i=1;i<=m;i++) 23 { 24 s[i]=mexp(n,i,mod)*mexp(m,n+1,mod)%mod; 25 for(int j=0;j<i;j++)s[i]=((s[i]+powm1[i-j]*C(i,j)*s[j])%mod+mod)%mod; 26 s[i]=s[i]*mexp(m-1,mod-2,mod)%mod; 27 } 28 cout<<s[m]<<endl; 29 return 0; 30 }