2019.04.02【51nod1258】序列求和 V4（伯努利数）（生成函数）（多项式求逆）

传送门

---

解析：

伯努利数初探（还没写，几天后还看见这个说明我忘了，记得催更）

由于$k$变大了，$O(n^2)$，所以需要用伯努利数的生成函数通过多项式求逆来得到数列。

由于毒瘤模数需要用MTT。

---

代码：

#include<bits/stdc++.h>
#define ll long long
#define re register
#define gc get_char
#define cs const

namespace IO{
	inline char get_char(){
		static cs int Rlen=1<<20|1;
		static char buf[Rlen],*p1,*p2;
		return (p1==p2)&&(p2=(p1=buf)+fread(buf,1,Rlen,stdin),p1==p2)?EOF:*p1++;
	}
	
	inline ll getint(){
		re char c;
		while(!isdigit(c=gc()));re ll num=c^48;
		while(isdigit(c=gc()))num=(num+(num<<2)<<1)+(c^48);
		return num;
	}
}
using namespace IO;

using std::cout;
using std::cerr;

cs int mod=1e9+7;
inline int add(int a,int b){return a+b>=mod?a+b-mod:a+b;}
inline int dec(int a,int b){return a<b?a-b+mod:a-b;}
inline int mul(int a,int b){return (ll)a*b%mod;}
inline int quickpow(int a,int b,int res=1){
	while(b){
		if(b&1)res=mul(res,a);
		a=mul(a,a);
		b>>=1;
	}
	return res;
}
template<typename ...Args>
inline int mul(int a,Args ...args){return mul(a,mul(args...));}
template<typename ...Args>
inline int add(int a,Args ...args){return add(a,add(args...));}

cs double PI=acos(-1);
struct Complex{
	double x,y;
	Complex(){}
	Complex(cs double &_x,cs double &_y):x(_x),y(_y){}
	friend Complex operator+(cs Complex &a,cs Complex &b){return Complex(a.x+b.x,a.y+b.y);}
	friend Complex operator-(cs Complex &a,cs Complex &b){return Complex(a.x-b.x,a.y-b.y);}
	friend Complex operator*(cs Complex &a,cs Complex &b){return Complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
	friend Complex operator/(cs Complex &a,cs double &b){return Complex(a.x/b,a.y/b);}
}; 

cs int N=1<<17|7;
int r[N];
Complex w[N];
inline void FFT(Complex *A,int len,int typ){
	if(typ==-1)std::reverse(A+1,A+len);
	for(int re i=0;i<len;++i)if(i<r[i])std::swap(A[i],A[r[i]]);
	for(int re i=1;i<len;i<<=1)
	for(int re j=0;j<len;j+=i<<1)
	for(int re k=0;k<i;++k){
		static Complex x,y;
		x=A[j+k],y=A[j+k+i]*w[len/i/2*k];
		A[j+k]=x+y;
		A[j+k+i]=x-y;
	}
	if(typ==-1)for(int re i=0;i<len;++i)A[i]=A[i]/len;
}

int mxlen;
inline void init_rev(int len){
	mxlen=std::max(mxlen,len);
	for(int re i=0;i<len;++i){
		r[i]=r[i>>1]>>1|((i&1)*(len>>1));
		w[i]=Complex(cos(2*PI*i/len),sin(2*PI*i/len));
	}
}

typedef std::vector<int> Poly;

inline Poly operator*(cs Poly &a,cs Poly &b){
	int deg=a.size()+b.size()-1,len=1;
	while(len<deg)len<<=1;
	init_rev(len);
	static Complex A[N],B[N],C[N],D[N];
	for(int re i=0;i<len;++i){
		A[i]=i<a.size()?Complex(a[i]&0x7fff,a[i]>>15):Complex(0,0);
		B[i]=i<b.size()?Complex(b[i]&0x7fff,b[i]>>15):Complex(0,0);
	}
	FFT(A,len,1),FFT(B,len,1);
	for(int re i=0,u;i<len;++i){
		u=(len-i)%len;
		C[i]=Complex((A[i].x+A[u].x)/2,(A[i].y-A[u].y)/2)*B[i];
		D[i]=Complex((A[i].y+A[u].y)/2,(A[u].x-A[i].x)/2)*B[i];
	}
	FFT(C,len,-1),FFT(D,len,-1);
	Poly c;c.resize(deg);
	for(int re i=0;i<deg;++i){
		ll x=round(C[i].x),y=round(C[i].y)+round(D[i].x),z=round(D[i].y);
		x%=mod,y%=mod,z%=mod;
		x=(x+mod)%mod,y=(y+mod)%mod,z=(z+mod)%mod;
		c[i]=(x+(y<<15)+(z<<30))%mod;
	}
	return c;
}

inline Poly Inv(cs Poly &a){
	Poly c,b(1,quickpow(a[0],mod-2));int lim=a.size();
	for(int re l=4;(l>>2)<lim;l<<=1){
		init_rev(l);
		c=a,c.resize(l>>1);
		c=c*b,c.resize(l>>1);
		c=c*b,c.resize(l>>1);
		b.resize(l>>1);
		for(int re i=0;i<(l>>1);++i)
		b[i]=add(b[i],b[i],mod-c[i]);
	}
	b.resize(lim);
	return b;
}

cs int P=5e4+4;
int fac[P],ifac[P],inv[P];
int B[P];Poly a;
inline int C(int n,int m){return mul(fac[n],ifac[m],ifac[n-m]);}

signed main(){
	fac[0]=fac[1]=ifac[0]=ifac[1]=inv[0]=inv[1]=1;
	for(int re i=2;i<P;++i){
		fac[i]=mul(fac[i-1],i);
		inv[i]=mul(inv[mod%i],mod-mod/i);
		ifac[i]=mul(ifac[i-1],inv[i]);
	}
	for(int re i=0;i+1<P;++i)a.push_back(ifac[i+1]);
	a=Inv(a);
	for(int re i=0;i<P;++i)B[i]=mul(a[i],fac[i]);
	int q=getint();
	while(q--){
		int n=(getint()+1)%mod,k=getint(),ans=0,now=n;
		for(int re i=k;~i;--i,now=mul(now,n))ans=add(ans,mul(C(k+1,i),B[i],now));
		cout<<mul(ans,inv[k+1])<<"\n";
	} 
	return 0;
}