51nod 1258 FFT 多项式求逆 伯努利数

#include<iostream>
#include<algorithm>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#define rd round
using namespace std;
typedef long long LL;
inline LL read()
{
	LL x=0;bool f=0;char c=getchar();
	for (;c<'0'||c>'9';c=getchar()) f=c=='-'?1:0;
	for (;c>='0'&&c<='9';c=getchar()) x=x*10+c-'0';
	return f?-x:x;
}
const unsigned int K=(1<<17)+5,P=1000000007,M=31627;
const double Pi=acos(-1.0);
struct na
{
	long double a,b;
	na(){}
	na(long double x,long double y) {a=x;b=y;}
	inline na operator + (const na &x) {return na(a+x.a,b+x.b);}
	inline na operator - (const na &x) {return na(a-x.a,b-x.b);}
	inline na operator * (const na &x) {return na(a*x.a-b*x.b,a*x.b+b*x.a);}
	inline na operator ~ () {return na(a,-b);}
}A[K],B[K],C[K],D[K],w[2][K];
int rev[K],n=50005,N,s;
LL fac[K],fr[K],f[K],g[K],h[K],_B[K];
inline LL _C(int a,int b){
	if (b<0||b>a) return 0;
	return fac[a]*fr[b]%P*fr[a-b]%P;
}
inline LL _pow(LL a,int k)
{
	LL rec=1;
	for (;k;a=a*a%P,k>>=1)
		(k&1)?(rec=rec*a%P):0;
	return rec;
}

void fft(na A[],int r)
{
	for (int i=0;i<N;i++)
		if (rev[i]<i) swap(A[i],A[rev[i]]);
	na x,y;
	for (int i=1,d=1<<16;i<N;i<<=1,d>>=1)
		for (int j=0;j<N;j+=i<<1)
			for (int k=0,p=0;k<i;k++,p+=d)
				x=A[j+k],y=A[j+i+k]*w[r][p],A[j+k]=x+y,A[j+i+k]=x-y;
	if (r) for (int i=0;i<N;i++) A[i].a/=N;
}

void calc(int n)
{
	if (n==1) {g[0]=1;return;}
	calc(n+1>>1);
	for (N=2,s=0;N<n<<1;N<<=1,s++);
	for (int i=1;i<N;i++) rev[i]=rev[i>>1]>>1|(i&1)<<s;

	for (int i=0;i<n+1>>1;i++) A[i]=na(g[i]%M,0),B[i]=na(g[i]/M,0);
	for (int i=n+1>>1;i<N;i++) A[i]=B[i]=na(0,0);
	fft(A,0);fft(B,0);
	for (int i=0;i<N;i++)
		C[i]=B[i]*B[i],B[i]=A[i]*B[i],A[i]=A[i]*A[i];
	fft(A,1);fft(B,1);fft(C,1);
	for (int i=0;i<N;i++)
		h[i]=((LL)rd(C[i].a)%P*M*M%P+(LL)rd(B[i].a)%P*2*M%P+(LL)rd(A[i].a)%P)%P;

	for (int i=0;i<n;i++)
		A[i]=na(f[i]%M,0),B[i]=na(f[i]/M,0),C[i]=na(h[i]%M,0),D[i]=na(h[i]/M,0);
	for (int i=n;i<N;i++) A[i]=B[i]=C[i]=D[i]=na(0,0);
	fft(A,0);fft(B,0);fft(C,0);fft(D,0);
	for (int i=0;i<N;i++)
	{
		na t=A[i];
		A[i]=A[i]*C[i],C[i]=C[i]*B[i],B[i]=B[i]*D[i],D[i]=D[i]*t;
	}
	fft(A,1);fft(B,1);fft(C,1);fft(D,1);
	for (int i=0;i<N;i++)
		h[i]=((LL)rd(B[i].a)%P*M*M%P+(LL)rd(C[i].a+D[i].a)%P*M%P+(LL)rd(A[i].a)%P)%P;
	for (int i=0;i<n;i++)
		g[i]=(g[i]+g[i]-h[i]+P)%P;
}

void solve()
{
	LL n=read()%P,p=read(),ans=0;
	for (int i=0,k=1;i<=p;i++,k=-k)
		ans=(ans+_B[i]*_C(p+1,i)%P*_pow(n,p+1-i)%P*k+P)%P;
	ans=ans*_pow(p+1,P-2)%P;
	printf("%lld\n",ans);
}

int main()
{
	for (int i=0;i<1<<16;i++) w[0][i]=na(cos(Pi*i/(1<<16)),sin(Pi*i/(1<<16)));
	for (int i=0;i<1<<16;i++) w[1][i]=~w[0][i];
	fac[0]=fr[0]=1;
	for (int i=1;i<n;i++) fac[i]=fac[i-1]*i%P;
	for (int i=1;i<n;i++) fr[i]=_pow(fac[i],P-2);
	for (int i=0;i<n;i++) f[i]=fr[i+1];
	calc(n);
	for (int i=0;i<n;i++) _B[i]=g[i]*fac[i]%P;
	for (int cas=read();cas--;) solve();
	return 0;
}