果然我数学不行啊，题解君： http://www.cnblogs.com/zhuohan123/p/3726933.html
1 const 2 h=1000000007; 3 var 4 fac,facinv,powm,s:array[0..1100]of int64; 5 n,m:int64; 6 7 function mexp(a,b:int64):int64; 8 begin 9 if b=0 then exit(1); 10 mexp:=sqr(mexp(a,b>>1))mod h; 11 if b and 1=1 then mexp:=mexp*a mod h; 12 end; 13 14 function C(n,r:int64):int64; 15 begin 16 exit((fac[n]*facinv[r] mod h)*facinv[n-r] mod h); 17 end; 18 19 procedure main; 20 var 21 i,j:longint; 22 begin 23 read(n,m); 24 if m=1 then 25 begin 26 write((n*(n+1)div 2)mod h); 27 exit; 28 end; 29 fac[0]:=1; 30 for i:=1 to m do 31 fac[i]:=fac[i-1]*i mod h; 32 facinv[m]:=mexp(fac[m],h-2); 33 facinv[0]:=1; 34 for i:=m-1 downto 1 do 35 facinv[i]:=facinv[i+1]*(i+1)mod h; 36 powm[0]:=1; 37 for i:=1 to m do 38 powm[i]:=-powm[i-1]; 39 s[0]:=((mexp(m,n+1)-m+h)mod h)*mexp(m-1,h-2)mod h; 40 for i:=1 to m do 41 begin 42 s[i]:=mexp(n,i)*mexp(m,n+1)mod h; 43 for j:=0 to i-1 do 44 s[i]:=(s[i]+powm[i-j]*(C(i,j)*s[j]mod h)+h)mod h; 45 s[i]:=s[i]*mexp(m-1,h-2)mod h; 46 end; 47 write(s[m]); 48 end; 49 50 begin 51 main; 52 end.