题目链接：

http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1228&judgeId=550328
就是求自然数幂的和

伯努利数来求

早就听说过伯努利数了，但是感觉太难了，一直没管，但是今天因为队友的一道题，突然推出一个递推式，貌似阔以用 $O\left(n^2\right)$O(n2)解决。。。结果别人的测试数据有几千组，每组都是 $O\left(n^2\right)$O(n2)的话就超时了，而伯努利数的这种方法是先 $O\left(n^2\right)$O(n2)求出伯努利数，然后再用 $O\left(n\right)$O(n)求出每组测试数据，这样就不得超时了

伯努利数的递推式： $\sum _{i=0}^k{C}_{k+1}^i{B}_i=0$i=0∑k​Ck+1i​Bi​=0每个 $B_k$Bk​都用这个算一下，就是 $O\left(n^2\right)了$O(n2)了

然后计算答案用这个公式：
$\frac{1}{k+1}\sum _{i=1}^{k+1}{C}_{k+1}^i{B}_{k+1-i}(n+1{)}^i$k+11​i=1∑k+1​Ck+1i​Bk+1−i​(n+1)i

伯努利这个公式怎么来的我也不知道。。。希望有一天能够懂怎么来的

#include"bits/stdc++.h"
using namespace std;
const int maxn=2e3+5;
const long long MOD=1e9+7;
long long fac[maxn],inv[maxn];
long long N,K,T;
long long B[maxn];//伯努利数
long long ksm(long long a,long long b)
{
    long long res=1,base=a;
    while(b)
    {
        if(b&1)res=(res*base)%MOD;
        base=(base*base)%MOD;
        b>>=1;
    }
    return res;
}
void Init()
{
    fac[0]=fac[1]=inv[0]=inv[1]=1;
    for(long long i=2;i<=2002;i++)
    {
        fac[i]=(fac[i-1]*i)%MOD;
        inv[i]=(inv[i-1]*ksm(i,MOD-2))%MOD;
    }
}
long long C(int n,int m)
{
    if(m==0||n==m)return 1;
    long long res=(fac[n]*inv[m])%MOD;
    res=(res*inv[n-m])%MOD;
    return res;
}
void Bernoulli()//O(n^2)求伯努利数
{
    B[0]=1;
    long long sum=1;
    for(int i=1;i<=2001;i++)
    {
        long long tp=0;
        for(int j=0;j<i;j++)tp=(tp+C(i+1,j)*B[j])%MOD;
        B[i]=((MOD-tp)*ksm(i+1,MOD-2))%MOD;
    }
}
long long f()//O(n)求答案
{
    long long p[maxn]={1};
    for(int i=1;i<=K+1;i++)p[i]=(p[i-1]*(N+1))%MOD;
    long long res=0;
    for(int i=1;i<=K+1;i++)
    {
        long long tp=(C(K+1,i)*B[K+1-i])%MOD;
        tp=(tp*p[i])%MOD;
        res=(res+tp)%MOD;
    }
    res=(res*ksm(K+1,MOD-2))%MOD;
    return res;
}
int main()
{
    Init();
    Bernoulli();
    cin>>T;
    while(T--)
    {
        cin>>N>>K;
        N%=MOD;//开始忘了加上这句话。。。。。
        cout<<f()<<endl;
    }
}

递推式来求

以平方和来举例：

$二项式展开:(n+1{)}^3=C_3^0{1}^0{n}^3+C_3^1{1}^1{n}^2+C_3^2{1}^2{n}^1+C_3^3{1}^3{n}^0$二项式展开:(n+1)3=C30​10n3+C31​11n2+C32​12n1+C33​13n0
把第一项移到左边去，并且系数是1,就省略了
$(n+1{)}^3-n^3=C_3^1{1}^1{n}^2+C_3^2{1}^2{n}^1+C_3^3{1}^3{n}^0$(n+1)3−n3=C31​11n2+C32​12n1+C33​13n0
然后就开始推
$(n+1{)}^3-n^3=C_3^1{1}^1{n}^2+C_3^2{1}^2{n}^1+C_3^3{1}^3{n}^0$(n+1)3−n3=C31​11n2+C32​12n1+C33​13n0
$(n{)}^3-(n-1{)}^3=C_3^1{1}^1(n-1{)}^2+C_3^2{1}^2(n-1{)}^1+C_3^3{1}^3(n-1{)}^0$(n)3−(n−1)3=C31​11(n−1)2+C32​12(n−1)1+C33​13(n−1)0
$(n-1{)}^3-(n-2{)}^3=C_3^1{1}^1(n-2{)}^2+C_3^2{1}^2(n-2{)}^1+C_3^3{1}^3(n-2{)}^0$(n−1)3−(n−2)3=C31​11(n−2)2+C32​12(n−2)1+C33​13(n−2)0
$...$...
$(2-1{)}^3-1^3=C_3^1{1}^1{1}^2+C_3^2{1}^2{1}^1+C_3^3{1}^3{1}^0$(2−1)3−13=C31​1112+C32​1211+C33​1310

然后把这 $n$n个式子求和，左边一加一减一加一减就只剩下 $(n+1{)}^3-1^3$(n+1)3−13，所以
$(n+1{)}^3-1^3=C_3^1{1}^1(1^2+2^2+...+n^2)+C_3^2{1}^2(1^1+2^1+...+n^1)+C_3^3{1}^3(1^0+2^0+...+n^0)$(n+1)3−13=C31​11(12+22+...+n2)+C32​12(11+21+...+n1)+C33​13(10+20+...+n0)
简写一哈就是这样：
$A=C_3^1{1}^1{I}_2+C_3^2{1}^2{I}_1+C_3^3{1}^3{I}_0$A=C31​11I2​+C32​12I1​+C33​13I0​
其中:
$A=(n+1{)}^3-1,只需要对数级别的快速幂计算一哈就行$A=(n+1)3−1,只需要对数级别的快速幂计算一哈就行
$I_2就是我们要求的答案平方和,I_1是一次方和,I_0是0次方和$I2​就是我们要求的答案平方和,I1​是一次方和,I0​是0次方和
所以递推一次差不多要 $k$k次运算，然后要差不多 $k$k次递推，所以是 $k^2$k2的复杂度

/*递推版，每一计算是O(k^2),过不了*/
#include"bits/stdc++.h"
using namespace std;
const int maxn=2e3+5;
const int MOD=1e9+7;
long long dp[maxn];
long long fac[maxn],inv[maxn];
long long N,K,T;
long long C[maxn][maxn];
long long ksm(long long a,long long b)
{
    long long res=1,base=a;
    while(b)
    {
        if(b&1)res=(res*base)%MOD;
        base=(base*base)%MOD;
        b>>=1;
    }
    return res;
}
void Init()
{
    fac[0]=fac[1]=inv[0]=inv[1]=1;
    for(long long i=2;i<=2000;i++)
    {
        fac[i]=(fac[i-1]*i)%MOD;
        inv[i]=(MOD-MOD/i)*inv[MOD%i]%MOD;
    }

    for(int i=0;i<=2000;i++)
    {
        C[i][0]=C[i][i]=1;
        for(int j=1;j<i;j++)
        {
            C[i][j]=(C[i-1][j-1]+C[i-1][j])%MOD;
        }
    }

}
long long f(int k)
{
    long long res=(ksm(N+1,k+1)-1+MOD)%MOD;
    long long tp=0;
    for(int i=2;i<=k+1;i++)
    {
        tp+=(C[k+1][i]*dp[k+1-i])%MOD;
        tp%=MOD;
    }
    res=(res-tp+MOD)%MOD;
    res=(res*ksm(k+1,MOD-2))%MOD;
    return dp[k]=res;
}
int main()
{
    Init();
    cin>>T;
    while(T--)
    {
        cin>>N>>K;
        dp[0]=N%MOD;
        for(int i=1;i<=K;i++)f(i);
        cout<<dp[K]<<endl;
    }
}