51nod1228：序列求和（伯努利数）

题面

题意：求幂和$\sum_{i=0}^ni^k$ 模1e9+7
n≤1e18，k≤2e3，5000组数据

曾记得L指导koi后教我们杜教筛，普及了下幂和怎么求
一本正经地推试子
原来那就是伯努利数
还说实在不会就写高斯消元

由于我还很菜，还在学，就随便吹一些。

当某个数列${B_i}​$指数生成函数为$\frac{x}{e^x-1}​$，
我们定义它叫伯努利数
即$\sum_{n=0}^\infty \frac{B_n}{n!}x^n=\frac{x}{e^x-1}​$

根据常识（泰勒展开），数列{${t^i}$}的指数生成函数为$e^{tx}$

定义伯努利多项式$\beta_n(t)=\sum_{k=0}^n C_n^k B_{n-k}t^k$
为${t^i}$与${B_i}​$的指数生成多项式的积

设$S_k(x)=\sum_{i=1}^x i^k$
根据DYH的ptt上的推导，得到了一条很重要的试子

$$S_k(x)=\frac{1}{k+1}\sum_{g=0}^k C_{k+1}^g B_g x^{k+1-g}$$

把x取0，得到B的递推试

$$B_n=-\frac{1}{n+1}\sum_{k=0}^{n-1}C_{n+1}^{k}B_k$$

根据b的定义，有

$$\sum_{n=0}^\infty \frac{B_n}{n!}x^n=\frac{x}{e^x-1}​=\frac{1}{\sum_{i=0}^\infty \frac{x^i}{(i+1)!}}$$

多项式求逆就可以了。

#include <iostream>
#include <fstream>
#include <algorithm>
#include <cmath>
#include <ctime>
#include <cstdio>
#include <cstdlib>
#include <cstring>

using namespace std;
#define mmst(a, b) memset(a, b, sizeof(a))
#define mmcp(a, b) memcpy(a, b, sizeof(b))

typedef long long LL;

const int N=5050;
const LL mo=1e9+7;

int T;
LL B[N],jc[N],I[N],Ijc[N];
LL n,k,ans;
LL f[N];

int main()
{
    jc[0]=Ijc[0]=I[1]=1;
    for(int i=2;i<=3000;i++)
    I[i]=(mo-mo/i)*I[mo%i]%mo;

    for(int i=1;i<=3000;i++)
    jc[i]=jc[i-1]*i%mo,Ijc[i]=Ijc[i-1]*I[i]%mo;

    B[0]=1;
    for(int i=1;i<=2000;i++)
    {
        for(int j=0;j<=i-1;j++)
        B[i]=(B[i]+jc[i+1]*Ijc[j]%mo*Ijc[i+1-j]%mo*B[j]%mo)%mo;
        B[i]=mo-B[i];
        B[i]=B[i]*I[i+1]%mo;
    }

    cin>>T;

    while(T--)
    {

    scanf("%lld%lld",&n,&k);
    ans=0;
    f[0]=1;

    for(int i=1;i<=k+1;i++)
    f[i]=f[i-1]*((n+1)%mo)%mo;

    for(int i=1;i<=k+1;i++)
    ans=(ans+jc[k+1]*Ijc[i]%mo*Ijc[k+1-i]%mo*B[k+1-i]%mo*f[i]%mo)%mo;
    ans=ans*I[k+1]%mo;

    printf("%lld\n",ans);
    }

    return 0;
}