伯努利数

应用

计算$\Sigma(i^p)(1<=i<=n)$

代码

#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 LL 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], Inv[MAXN], B[MAXN];
void Exgcd(LL a, LL b, LL &x, LL &y)
{
    if(b == 0)
    {
        x = 1;
        y = 0;
        return ;
    }
    LL x1, y1;
    Exgcd(b, a%b, x1, y1);
    x = y1;
    y = x1 - (a/b)*y1;
}
void Get_Fac()
{
    for(int i=0; i<MAXN; i++)
    {
        c[i][0] = 1;
        c[i][i] = 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;
}

void Get_Inv()
{
    for(int i=1; i<MAXN; i++)
    {
        LL x, y;
        Exgcd(i, MOD, x, y);
        x = (x%MOD+MOD)%MOD;
        Inv[i] = x;
    }
}
LL quick_MOD(LL a, LL b)
{
    LL ans = 1;
    while(b)
    {
        if(b & 1)
            ans = (ans*a)%MOD;
        b>>=1;
        a = (a*a)%MOD;
    }
    return ans;
}
void Get_Bonuli()
{
    B[0] = 1;
    for(int i=1; i<MAXN-1; i++)
    {
        LL tmp = 0;
        for(int j=0; j<i; j++)
            tmp = (tmp+c[i+1][j]*B[j])%MOD;
        B[i] = tmp;
        B[i] = B[i]*(-Inv[i+1]);
        B[i] = (B[i]%MOD+MOD)%MOD;
    }
}
void Init()
{
    Get_Fac();
    Get_Inv();
    Get_Bonuli();
}
int main()
{
    Init();
    int T, k;
    LL n;
    scanf("%d",&T);
    while(T--)
    {
        scanf("%lld%d",&n,&k);
        n++;
        n %= MOD;
        LL ans = 0;
        for(int i=1; i<=k+1; i++)
        {
            ans = (ans+((c[k+1][i]*B[k+1-i])%MOD)*quick_MOD(n,(LL)i))%MOD;
            ans = (ans%MOD+MOD)%MOD;
        }
        ans = ans*Inv[k+1];
        ans = (ans%MOD+MOD)%MOD;
        printf("%lld\n",ans);
    }
    return 0;
}