[51nod1228][伯努利数][自然数k次幂和]序列求和

题意

---

给定n,k，求$\sum_{i=1}^{n} i^k$

---

因为绍一模拟考的20分部分分，所以来学一发伯努利数

首先暴力求这个式子是要nlogk的，然而n是10^18……

然后了解到一个叫伯努利数的东西，这个k^2次预处理后O(k)求出解

伯努利数是一个乱七八糟的有理数数列，

定义$B_0=1$

根据关系式$\sum_{i=0}^{n} C_{n+1}^{i}B_i=0$可以得到递推式

$B_n=-\frac{1}{n+1}\sum_{i=0}^{n-1}C_{n+1}^{i}B_i$

这样就可以O(k^2)预处理出伯努利数

因为又知道一个神奇的公式

$\sum_{i=1}^n i^k=\frac{1}{k+1} \sum_{i=1}^{k+1}C_{k+1}^{i}B_{k+1-i}(n+1)^i$

那么就可以O(k)求出自然数k次幂和啦

学习自http://blog.csdn.net/acdreamers/article/details/38929067

#include <cstdio>
#include <iostream>
#include <algorithm>
#define N 2010
#define p 1000000007

using namespace std;

typedef long long ll;

int t;
ll n,k;
ll inv[N],fac[N],invf[N],B[N];

inline ll C(ll x,ll y){
  return fac[x]*invf[y]%p*invf[x-y]%p;
}

inline void reaD(int &x){
  char c=getchar();x=0;
  for(;c>57||c<48;c=getchar());for(;c>=48&&c<=57;x=x*10+c-48,c=getchar());
}

inline void reaD(ll &x){
  char c=getchar();x=0;
  for(;c>57||c<48;c=getchar());for(;c>=48&&c<=57;x=x*10+c-48,c=getchar());
}

ll pow_(ll x,ll y){
  x%=p;
  ll r=1;
  while(y){
    if(y&1) r=r*x%p;
    x=x*x%p;
    y>>=1;
  }
  return r;
}

int main(){
  freopen("1.in","r",stdin);
  freopen("1.out","w",stdout);
  invf[0]=inv[0]=inv[1]=fac[0]=fac[1]=1;
  for(ll i=1;i<=2005;i++) fac[i]=fac[i-1]*i%p;
  for(int i=2;i<=2005;i++) inv[i]=(p-p/i)*inv[p%i]%p;
  for(int i=1;i<=2005;i++) invf[i]=invf[i-1]*inv[i]%p;
  B[0]=1;
  for(int i=1;i<=2001;i++){
    B[i]=0;
    for(int j=0;j<i;j++) B[i]=(B[i]+C(i+1,j)*B[j])%p;
    B[i]=((B[i]*-inv[i+1])%p+p)%p;
  }
  reaD(t);
  while(t--){
    reaD(n); reaD(k);
    ll Ans=0;
    for(int i=1;i<=k+1;i++) Ans=(Ans+C(k+1,i)*B[k+1-i]%p*pow_(n+1,i))%p;
    Ans=Ans*inv[k+1]%p;
    printf("%lld\n",Ans);
  }
  return 0;
}

查了下维基百科

伯努利数有两类，两类伯努利数的区别在于$B_1$的正负

第一类伯努利数，由美国国家标准技术研究所 (NIST)制定，在这标准下$B_1=-\frac{1}{2}$

第二类伯努利数，又被称为是“原始的白努利数“，在这标准下$B_1=\frac{1}{2}$

上面的求自然数k次幂和是根据第一类伯努利数得出，可以看出这个做法得到的多项式是关于(n+1)的，那么在某些题目中，要求多项式关于n，这时候就需要用第二类伯努利数。

第二类伯努利数求自然数k次幂和：

$\sum_{i=1}^{n} i^k=\frac{1}{k+1}\sum_{i=0}^{k}C_{k+1}^{i} B_i n^{k+1-i}$

很神奇啊

#include <cstdio>
#include <iostream>
#include <algorithm>
#define N 2010
#define p 1000000007

using namespace std;

typedef long long ll;

int t;
ll n,k;
ll inv[N],fac[N],invf[N],B[N];

inline ll C(ll x,ll y){
  return fac[x]*invf[y]%p*invf[x-y]%p;
}

inline void reaD(int &x){
  char c=getchar();x=0;
  for(;c>57||c<48;c=getchar());for(;c>=48&&c<=57;x=x*10+c-48,c=getchar());
}

inline void reaD(ll &x){
  char c=getchar();x=0;
  for(;c>57||c<48;c=getchar());for(;c>=48&&c<=57;x=x*10+c-48,c=getchar());
}

ll pow_(ll x,ll y){
  x%=p;
  ll r=1;
  while(y){
    if(y&1) r=r*x%p;
    x=x*x%p;
    y>>=1;
  }
  return r;
}

int main(){
  invf[0]=inv[0]=inv[1]=fac[0]=fac[1]=1;
  for(ll i=1;i<=2005;i++) fac[i]=fac[i-1]*i%p;
  for(int i=2;i<=2005;i++) inv[i]=(p-p/i)*inv[p%i]%p;
  for(int i=1;i<=2005;i++) invf[i]=invf[i-1]*inv[i]%p;
  B[0]=1;
  for(int i=1;i<=2001;i++){
    B[i]=0;
    for(int j=0;j<i;j++) B[i]=(B[i]+C(i+1,j)*B[j])%p;
    B[i]=((B[i]*-inv[i+1])%p+p)%p;
  }
  B[1]=(p-B[1])%p;
  reaD(t);
  while(t--){
    reaD(n); reaD(k);
    ll Ans=0;
    for(int i=0;i<=k;i++) Ans=(Ans+C(k+1,i)*B[i]%p*pow_(n,k+1-i))%p;
    Ans=Ans*inv[k+1]%p;
    printf("%lld\n",Ans);
  }
  return 0;
}