51nod 1228、1258 序列求和
这里一次讲两题...貌似都是板子？
所以两题其实可以一起做 【雾
noteskey
总之就是伯努利数的两道入门题啦，就是第二道有点鬼畜了，居然要任意模数的！（好吧是 1e9+7 但也没什么区别了）
伯努利数学习可以看这里
第二题的式子其实学习笔记里已经有写了...这里就再推一遍吧（不然没什么好写的呢~）
注意下面的伯努利数 \(B\) 用的是 \(B^+\) 伯努利
\[\begin{aligned} ANS=&S(n,k) \\=& {1\over k+1} \sum_{i=0}^k \begin{pmatrix} k+1\\i \end{pmatrix} B_i ~· n^{k+1-i} \\=& k! \sum_{i=0}^k {B_i\over i!} ~· {n^{k+1-i}\over (k+1-i)! }\end{aligned}
\]
然后后面的式子一眼卷积，一遍...咳咳...一遍 MTT 出解...
code
然后第一题其实可以暴力做的...但是反正第二题也是要打的，还打个 P 的题 1 暴力，当然是两题一起肝
所以代码只给第二道题的（虽说两道都能 A ）... MTT 就抄个板子好了，我是不来再打一遍的，打错还得调半天 XD
//by Judge
#include<bits/stdc++.h>
#define Rg register
#define fp(i,a,b) for(Rg int i=(a),I=(b)+1;i<I;++i)
#define fd(i,a,b) for(Rg int i=(a),I=(b)-1;i>I;--i)
#define ll long long
#define db double
using namespace std;
const db PI=acos(-1);
const int mod=1e9+7;
const int bl=32768;
const int M=5e5+3;
typedef int arr[M];
#ifndef Judge
#define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++)
#endif
char buf[1<<21],*p1=buf,*p2=buf;
inline bool cmax(int& a,int b){return a<b?a=b,1:0;}
inline bool cmin(int& a,int b){return a>b?a=b,1:0;}
inline int mul(int x,int y){return 1ll*x*y%mod;}
inline int dec(int x,int y){return x<y?x-y+mod:x-y;}
inline int inc(int x,int y){return x+y>=mod?x+y-mod:x+y;}
inline ll read(){ ll x=0,f=1; char c=getchar();
for(;!isdigit(c);c=getchar()) if(c=='-') f=-1;
for(;isdigit(c);c=getchar()) x=x*10+c-'0'; return x*f;
} char sr[1<<21],z[20];int CCF=-1,Z;
inline void Ot(){fwrite(sr,1,CCF+1,stdout),CCF=-1;}
inline void print(int x,char chr='\n'){
if(CCF>1<<20)Ot();if(x<0)sr[++CCF]=45,x=-x;
while(z[++Z]=x%10+48,x/=10);
while(sr[++CCF]=z[Z],--Z);sr[++CCF]=chr;
} int n,k,now,ans,d,limit; arr fac,inv,ifac,B,r,lg; db iv[M];
inline int qpow(int x,int p=mod-2,int s=1){
for(;p;p>>=1,x=mul(x,x)) if(p&1) s=mul(s,x); return s;
}
struct cp{ db x,y; cp(db _x,db _y=0){x=_x,y=_y;} cp(){}
cp operator +(const cp& b)const{return cp(x+b.x,y+b.y);}
cp operator -(const cp& b)const{return cp(x-b.x,y-b.y);}
cp operator *(const cp& b)const{return cp(x*b.x-y*b.y,x*b.y+y*b.x);}
cp operator *(const db& b)const{return cp(x*b,y*b);}
cp operator ~()const{return cp(x,-y);}
}w[2][M],x,y;
inline void prep(){ for(;1<<d<limit;++d);
fp(i,0,d) lg[1<<i]=i,iv[i]=(db)1.0/(1<<i);
fp(i,1,limit-1) r[i]=(r[i>>1]>>1)|((i&1)<<(d-1));
for(int i=1,d=0;i<limit;i<<=1,++d) fp(k,0,i-1){
w[1][i+k]=cp(cos(PI*k*iv[d]),sin(PI*k*iv[d]));
w[0][i+k]=cp(cos(PI*k*iv[d]),-sin(PI*k*iv[d]));
}
}
inline void FFT(cp* a,int tp){
fp(i,0,limit-1) if(i<r[i]) swap(a[i],a[r[i]]);
for(int mid=1;mid<limit;mid<<=1)
for(int j=0;j<limit;j+=mid<<1) for(int k=0;k<mid;++k)
x=a[j+k],y=w[tp][mid+k]*a[j+k+mid],a[j+k]=x+y,a[j+k+mid]=x-y;
if(!tp) fp(i,0,limit-1) a[i]=a[i]*iv[d];
}
inline void MTT(int* a,int* b,int* c,int n){ static cp f[M],g[M],p[M],q[M];
fp(i,0,n) f[i]=cp(a[i]>>15,a[i]&32767),g[i]=cp(b[i]>>15,b[i]&32767);
fp(i,n+1,limit-1) f[i]=g[i]=cp(0,0); prep(),FFT(f,1),FFT(g,1);
fp(i,0,limit-1){ cp t,f0,f1,g0,g1;
t=~f[i?limit-i:0],f0=(f[i]-t)*cp(0,-0.5),f1=(f[i]+t)*0.5,
t=~g[i?limit-i:0],g0=(g[i]-t)*cp(0,-0.5),g1=(g[i]+t)*0.5;
p[i]=f1*g1,q[i]=f1*g0+f0*g1+f0*g0*cp(0,1);
} FFT(p,0),FFT(q,0);
fp(i,0,n) c[i]=(((ll)(p[i].x+0.5)%mod<<30)+((ll)(q[i].x+0.5)<<15)+((ll)(q[i].y+0.5)))%mod;
}
inline void get_inv(int* a,int* b,int n){ static arr c,d;
if(n==1) return b[0]=qpow(a[0]),void();
get_inv(a,b,n>>1); int len=-1;
for(limit=1;limit<=n;limit<<=1) ++len;
fp(i,0,limit-1) r[i]=(r[i>>1]>>1)|((i&1)<<len);
MTT(a,b,c,n),MTT(c,b,d,n);
fp(i,0,n-1) b[i]=dec(inc(b[i],b[i]),d[i]);
}
inline void prep(int len){
inv[0]=inv[1]=ifac[0]=ifac[1]=fac[0]=fac[1]=1;
fp(i,2,len) fac[i]=mul(fac[i-1l],i),
inv[i]=mul(mod-mod/i,inv[mod%i]),
ifac[i]=mul(ifac[i-1],inv[i]);
get_inv(ifac+1,B,len<<1); B[1]=mod-B[1];
fp(i,0,len) B[i]=mul(B[i],fac[i]);
}
inline int C(int n,int m){return mul(fac[n],mul(ifac[m],ifac[n-m]));}
int main(){ prep(65536);
for(int T=read();T;--T,print(mul(ans,inv[k+1]))){ n=read()%mod,k=read(),ans=0;
now=1; fd(i,k,0) now=mul(now,n),ans=inc(ans,mul(C(k+1,i),mul(B[i],now)));
} return Ot(),0;
}
```cpp
