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

题意

---

给定n,k，求∑ni=1ik

---

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

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

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

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

定义B0=1

根据关系式∑ni=0Cin+1Bi=0可以得到递推式

Bn=−1n+1∑n−1i=0Cin+1Bi

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

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

∑ni=1ik=1k+1∑k+1i=1Cik+1Bk+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 1000000007using 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;}

查了下维基百科

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

第一类伯努利数，由美国国家标准技术研究所 (NIST)制定，在这标准下B1=−12

第二类伯努利数，又被称为是“原始的白努利数“，在这标准下B1=12

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

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

∑ni=1ik=1k+1∑ki=0Cik+1Bink+1−i

很神奇啊

#include <cstdio>#include <iostream>#include <algorithm>#define N 2010#define p 1000000007using 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;}