[伯努利数] 51Nod 1258 序列求和 V4

自然数幂和的求法有很多
大概有

• 由二项式定理得递推式
  
• 拉格朗日插值法&牛顿插值法
• 第一类斯特林数
• 伯努利数

wiki传送门
ACdreamer的blog

伯努利数有两种$B^+$和$B^-$
唯一的争议就在$B_1$上 因为 对于$n>1$有$B_{2n+1}=0$
其实这不是什么大问题 一般个人习惯用$B^-$ 以下不注明$B$都表示$B^-$

然后我就只会伯努利数最裸的应用了

伯努利数的预处理
可以直接通过得到递推式
这样是$O(m^2)$

范围大时使用$B$的指数生成函数
膜拜pls

顺便 对于另一类伯努利数$B^+$

母函数据百度说应该是 $\sum_{i=0}^\infty {B_i^+\over i!}x^i={xe^x\over e^x-1}$

然后回到这个题 %1e9+7？ 我只会三模数NTT
然后一次倍增要15次NTT
数据范围不大 跑的竟然比我想象的要快

#include<cstdio>
#include<cstdlib>
#include<algorithm>
using namespace std;
typedef long long ll;

const int P=1e9+7;
const int M[]={998244353,1004535809,469762049};
const int G[]={3,3,3};
const ll _M=(ll)M[0]*M[1];

inline ll Pow(ll a,int b,int p){
  ll ret=1;
  for (;b;b>>=1,a=a*a%p) if (b&1) ret=ret*a%p;
  return ret;
}
inline ll mul(ll a,ll b,ll p){
  a%=p; b%=p;
  return ((a*b-(ll)((ll)((long double)a/p*b+1e-3)*p))%p+p)%p;
}

const int m1=M[0],m2=M[1],m3=M[2];
const int inv1=Pow(m1%m2,m2-2,m2),inv2=Pow(m2%m1,m1-2,m1),inv12=Pow(_M%m3,m3-2,m3);
inline int CRT(int a1,int a2,int a3){
  ll A=(mul((ll)a1*m2%_M,inv2,_M)+mul((ll)a2*m1%_M,inv1,_M))%_M;
  ll k=((ll)a3+m3-A%m3)*inv12%m3;
  return (k*(_M%P)+A)%P;
}

const int N=264000;
int R[N];
struct NTT{
  int P,G;
  int num,w[2][N];
  void Pre(int _P,int _G,int m){
    num=m; P=_P; G=_G;
    int g=Pow(G,(P-1)/num,P);
    w[1][0]=1; for (int i=1;i<num;i++) w[1][i]=(ll)w[1][i-1]*g%P;
    w[0][0]=1; for (int i=1;i<num;i++) w[0][i]=w[1][num-i];
  }
  void FFT(int *a,int n,int r){
    for (int i=0;i<n;i++) if (i<R[i]) swap(a[i],a[R[i]]);
    for (int i=1;i<n;i<<=1)
      for (int j=0;j<n;j+=(i<<1))
    for (int k=0;k<i;k++){
      int x=a[j+k],y=(ll)a[j+i+k]*w[r][num/(i<<1)*k]%P;
      a[j+k]=(x+y)%P; a[j+i+k]=(x+P-y)%P;
    }
    if (!r) for (int i=0,inv=Pow(n,P-2,P);i<n;i++) a[i]=(ll)a[i]*inv%P;
  }
}ntt[3];

int tmp[N],b2[N],b3[N],_b[3][N],c[N];

inline void GetInv(int *a,int *b,int n){
  if (n==1) return void(b[0]=CRT(Pow(a[0],m1,m1-2),Pow(a[0],m2,m2-2),Pow(a[0],m3,m3-2)));
  GetInv(a,b,n>>1);
  int L=0; while (!(n>>L&1)) L++;
  for (int i=1;i<(n<<1);i++) R[i]=(R[i>>1]>>1)|((i&1)<<L);

  for (int i=0;i<n;i++) b2[i]=b3[i]=b[i],tmp[i]=a[i],b2[i+n]=b3[i+n]=0;
  ntt[0].FFT(b,n<<1,1);
  ntt[1].FFT(b2,n<<1,1);
  ntt[2].FFT(b3,n<<1,1);

  for (int I=0;I<3;I++){
    int *d=I==0?b:(I==1?b2:b3);
    for (int i=0;i<n;i++) tmp[n+i]=0,tmp[i]=a[i];
    ntt[I].FFT(tmp,n<<1,1);
    for (int i=0;i<(n<<1);i++)
      tmp[i]=(ll)tmp[i]*d[i]%M[I];
    ntt[I].FFT(tmp,n<<1,0);
    for (int i=0;i<(n<<1);i++) _b[I][i]=tmp[i];
  }

  for (int i=0;i<(n<<1);i++) c[i]=CRT(_b[0][i],_b[1][i],_b[2][i]),c[i]=c[i]==0?0:P-c[i];
  (c[0]+=2)%=P;

  for (int I=0;I<3;I++){
    int *d=I==0?b:(I==1?b2:b3);
    for (int i=0;i<n;i++) tmp[n+i]=0,tmp[i]=c[i];
    ntt[I].FFT(tmp,n<<1,1);
    for (int i=0;i<(n<<1);i++)
      tmp[i]=(ll)tmp[i]*d[i]%M[I];
    ntt[I].FFT(tmp,n<<1,0);
    for (int i=0;i<(n<<1);i++) _b[I][i]=tmp[i];
  }

  for (int i=0;i<n;i++) b[i]=CRT(_b[0][i],_b[1][i],_b[2][i]),b[n+i]=0;
}

ll fac[N],inv[N];

inline void Init(int n){
  fac[0]=1; for (int i=1;i<=n;i++) fac[i]=fac[i-1]*i%P;
  inv[1]=1; for (int i=2;i<=n;i++) inv[i]=(ll)(P-P/i)*inv[P%i]%P;
  inv[0]=1; for (int i=1;i<=n;i++) inv[i]=inv[i]*inv[i-1]%P;
}
inline ll C(int n,int m){
  return fac[n]*inv[m]%P*inv[n-m]%P;
}

int n,m;
int a[N],b[N];
int _a[3][N],B[N];

inline int Solve(ll n,int m){
  static ll pw[N];
  pw[0]=1; for (int i=1;i<=m+1;i++) pw[i]=pw[i-1]*(n%P)%P;
  ll ret=((ll)m+1)*((P+1)>>1)%P*pw[m]%P;
  for (int k=0;k<=m;k+=2)
    ret+=C(m+1,k)*B[k]%P*pw[m+1-k]%P;
  return ret%P*Pow(m+1,P-2,P)%P;
}

int main(){
  freopen("t.in","r",stdin);
  freopen("t.out","w",stdout);
  n=50000;
  for (m=1;m<=n;m<<=1);
  for (int i=0;i<3;i++) ntt[i].Pre(M[i],G[i],m<<1);
  Init(m+1);
  for (int i=0;i<m;i++) a[i]=inv[i+1];
  GetInv(a,b,m);
  for (int i=0;i<m;i++) B[i]=(ll)b[i]*fac[i]%P;
  int T,_K; ll _n;  scanf("%d",&T);
  while (T--){
    scanf("%lld%d",&_n,&_K);
    printf("%d\n",Solve(_n,_K));
  }
  return 0;
}