首先我来介绍一下什么是自然数幂和：
1+2+3+...+i+...+n=？
12+22+32+...+i2+...+n2=？
1k+2k+3k+...+ik+...+nk=？
类似上述式子的就是自然数幂和了，那么具体怎么求呢，这就是今天的重点了：
1+2+3+...+i+...+n=？ 这个式子大家肯定都会求，这就是一个等差数列求和，然后套公式：
S=(n?n+1)2
那么 1k+2k+3k+...+ik+...+nk这样的式子怎么求呢。
首先求 (n+1)k+1?nk+1????????????????????????(1)
根据牛顿二项式展开定理：
(n+1)k+1=C(k+1,0)?nk+1+C(k+1,1)?nk+...+C(k+1,k+1)?n0
所以：
(1)式=C(k+1,1)?nk+...+C(k+1,k+1)?n0
(n+1)k+1?nk+1=C(k+1,1)?nk+...+1??????????????(2)
那么我们现在消去了最高次幂的项 xk+1 那么我们现在对 (2)式从1?n求和 得到：
(n+1)k+1?nk+1=C(k+1,1)?nk+...+1
nk+1?(n?1)k+1=C(k+1,1)?(n?1)k+...+1
…
2k+1?1k+1=C(k+1,1)?1k+...+1
左边加左边，右边加右边，得到：
(n+1)k+1?1=C(k+1,1)?∑ni=1ik+...+n
根据上式我们得到幂指数是 k 的：
∑ni=1ik=1k+1?[(n+1)k+1?...?C(k+1,i)?∑ni=0ik+1?i?...?n?1]
也可以写为：
∑ni=1ik=1k+1?((n+1)k+1?(C(k+1,i)?∑ni=0ik+1?i+...+n+1))
现在 我们将 ∑ni=1ik 记作 S(n,k)
那么我们现在得到一个递推式：
S(n,k)=1k+1?((n+1)k+1?(C(k+1,2)?S(n,k?1)+C(k+1,i)?S(n,k+1?i)+C(k+1,k)?S(n,1)+n+1))
当 k==1 的时候是递归出口：
S(n,1)=n?(n+1)2
我现在来举个例子：
12+22+32+...+n2
我们先对 (n+1)3?n3 求和得到：
(n+1)3?n3=3?n2+3?n+1
n3?(n?1)3=3?(n?1)2+3?(n?1)+1
...
(1+1)3?b3=3?12+3?1+1
求和得：
(n+1)3?13=3?(12+22+...+n2)+3?(1+2+...+n)+n
那么(12+22+...+n2)=(n+1)3?(3?(1+2+...n)+n+1)3
又因为：
1+2+...+n=n?(n+1)2
所以：
(12+22+...+n2)=(n+1)3?(3?n?(n+1)2+n+1)3 =n?(2?n+1)?(n+1)6
模板：
/**
2016 - 08 - 06 上午
Author: ITAK
Motto:
今日的我要超越昨日的我，明日的我要胜过今日的我，
以创作出更好的代码为目标，不断地超越自己。
**/
#include <iostream>
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <vector>
#include <queue>
#include <algorithm>
#include <set>
using namespace std;
typedef long long LL;
typedef unsigned long long ULL;
const int INF = 1e9+5;
const int MAXN = 2e3+5;
const LL MOD = 1e9+7;
const double eps = 1e-7;
const double PI = acos(-1);
using namespace std;
LL c[MAXN][MAXN];
void Get_Fac()
{
for(int i=0; i<MAXN; i++)
c[i][0] = 1;
for(int i=1; i<MAXN; i++)
for(int j=1; j<=i; j++)
c[i][j] = (c[i-1][j]+c[i-1][j-1])%MOD;
}
LL ans[MAXN];
LL Solve(LL n, LL k)///1^k+2^k+...+n^k
{
Get_Fac();
if(ans[k] != -1)
return ans[k];
if(k == 1)
return (n+1)*n/2;
LL tmp = 1;
for(int i=0; i<=k; i++)
tmp = tmp*(n+1);
tmp = tmp - (n+1);
LL sum = 0;
for(int i=1; i<k; i++)
sum += c[k+1][i+1]*Solve(n,k-i);
ans[k] = (tmp-sum)/(k+1);
return ans[k];
}
int main()
{
LL n, k;
while(cin>>n>>k)
{
memset(ans, -1, sizeof(ans));
cout<<Solve(n,k)<<endl;
}
return 0;
}