伯努利数应用

上一次讲了一下什么是自然数幂和以及求法，但是个算法的时间复杂度是O(n2) 的，所以有的时候

不能解决一些题目，那么我就借助一下例题说一下伯努利数。

51NOD 1228 序列求和

T(n) = n^k，S(n) = T(1) + T(2) + …… T(n)。给出n和k，求S(n)。

例如k = 2，n = 5，S(n) = 1^2 + 2^2 + 3^2 + 4^2 + 5^2 = 55。

由于结果很大，输出S(n) Mod 1000000007的结果即可。

Input

第1行：一个数T，表示后面用作输入测试的数的数量。（1 <= T <= 5000)

第2 - T + 1行：每行2个数，N, K中间用空格分割。(1 <= N <= 10^18, 1 <= K <= 2000)

Output

共T行，对应S(n) Mod 1000000007的结果。

Input示例

3
5 3
4 2
4 1

Output示例

225
30
10

解题思路分析：

这个题目因为 T 非常大，所以不能直接用上一个递归的公式了，那么就引入了一个新的名词——伯努利数。

伯努利数是18世纪瑞士数学家雅各布·伯努利引入的一个数。设伯努利数为 Bn ，它的定义为： tet−1=∑∞n=0Bnn!∗tn这里|t|<2 。由计算知：B0=1,B1=−12 —— 摘自百度百科。

一般地，当 n≥2 时，我们可以通过 ∑ni=0C(n+1,i)∗Bk=0 通过这个公式我们就得到了 Bn 的公式：

Bn=−1C(n+1,n)∗(C(n+1,0)∗B0+C(n+1,1)∗B1+...+C(n+1,n−1)∗Bn−1)

=−1n+1∗(C(n+1,0)∗B0+C(n+1,1)∗B1+...+C(n+1,n−1)∗Bn−1)

那么 我现在给出 求 1k+2k+...+nk 的关于伯努利的公式：

∑i=1kik=1k+1∗C(k+1,i)∗Bk+1−i∗(n+1)i−−−(1)

这个公式的复杂度是 O(k) ，可以解决一些问题了。

其实这个公式中的大多数都是可以通过初始化来得到的，比如说 组合数 (n+1)i 和 1k+1 逆元，

都可以初始化一个数组来得到，然后就是比较复杂的 ”伯努利数“ 了，那么刚才介绍 伯努利数

的时候已经把公式给出了：

Bn=−1n+1∗(C(n+1,0)∗B0+C(n+1,1)∗B1+...+C(n+1,n−1)∗Bn−1)

这个公式很多东西也是可以初始化得到的， n+1 的逆元，组合数取模，等等。。

那么现在就可以写了，需要注意的是 在 (1) 式中 (n+1)i 需要先对 n 进行取模 ，否则 会爆

long long 的：

51NOD题目的AC代码：——可以作为模板

/**2016 - 08 - 07 下午Author: ITAKMotto:今日的我要超越昨日的我，明日的我要胜过今日的我，以创作出更好的代码为目标，不断地超越自己。**/#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;}