LGP5282题解

Min_25牛逼！Min_25牛逼！Min_25牛逼！

设有 \(f_d(x)=\prod_{i=1}^{d}(x+i)=f_{d-1}(x)(x+d)\)，设 \(B=\lfloor\sqrt{n}\rfloor\)，那么计算阶乘就相当于需要计算 \(f_B(iB)\)（\(i\in[0,B]\)）

考虑倍增。如果我们知道了 \(f_d(iB)\)（\(i\in[0,d]\)），考虑应该如何计算 \(f_{d+1}(iB))\)（\(i\in[0,d+1]\)） 和 \(f_{2d}(iB)\)（\(i\in[0,2d]\)）

第一个是较为显然的，只需要对 \(x+d+1\) 求出 \(x\in[0,dB],B\mid x\) 时的点值即可，对于 \(f_{d+1}(dB)\) 直接暴力计算即可。

第二个比较麻烦。

根据定义显然有 \(f_{2d}(x)=f_d(x)f_d(x+d)\)。也就是说我们只要知道了 \(f_d(iB),f_d(iB+d)\)（\(i\in[0,2d]\)） 的点值即可。

令 \(g_d(x)=f_d(xB)\)，上述问题变为给定 \(g_d(x)\)（\(x\in[0,2d]\)），求 \(g(x+\frac{d}{B})\)（\(x\in[0,2d]\)）。

很显然在模意义下有 \(\frac{d}{B}>B\)。

证明：因为 \(n<p\)，所以对于 \(k\in[0,B]\) 都有 \(kB<n<p\)，所以 \(kB\ne d\)，那么 \([0,B]\) 中一定不存在 \(\frac{d}{B}\)。

你会发现这个和 拉格朗日插值2 很像，但是我们并不知道 \(g_d(i)\)（\(i\in[d,2d]\)）。

可以考虑通过 \(g_d(i)\)（\(i\in[0,d]\)）求出 \(g_d(i)\)（\(i\in[d,2d]\)），再求上面的那个东西。

这里可以直接套用 LOJ166 的做法了。复杂度 \(T(n)=T(\frac{n}{2})+O(\log n)=T(n\log n)\)，这里 \(n=\sqrt N\) 所以复杂度是 \(O(\sqrt{n}\log n)\)。

第二部分的流程：\(f_d(iB)\)（\(i\in[0,d]\)）->\(f_d(iB)\)（\(i\in[0,2d]\)）->\(f_d(iB+d)\)（\(i\in[2d]\)）->\(f_{2d}(iB)\)（\(i\in[0,2d]\)）

#include<cassert>
#include<cstdio>
#include<cmath>
#define IMP(lim,act) for(int qwq=(lim),i=0;i^qwq;++i)act
typedef double db;
const int M=3e5+5;
const db Pi=acos(-1);
int n,P,F[M],ifac[M];
struct Barrett{
	typedef unsigned long long ull;
	typedef __uint128_t LL;
	ull m,B;
	Barrett(const ull&m=2):m(m),B((LL(1)<<64)/m){}
	friend inline ull operator%(const ull&a,const Barrett&mod){
		ull r=a-mod.m*(LL(mod.B)*a>>64);return r>=mod.m?r-mod.m:r;
	}
}mod;
struct complex{
	db x,y;
	complex(const db&x=0,const db&y=0):x(x),y(y){}
	inline complex operator+(const complex&it)const{
		return complex(x+it.x,y+it.y);
	}
	inline complex operator-(const complex&it)const{
		return complex(x-it.x,y-it.y);
	}
	inline complex operator*(const complex&it)const{
		return complex(x*it.x-y*it.y,x*it.y+y*it.x);
	}
}buf[M],*w[20];
inline int Getlen(const int&n){
	int len(0);while((1<<len)<n)++len;return len;
}
inline void swap(complex&a,complex&b){
	complex c=a;a=b;b=c;
}
inline int Get(const db&x){
	return((long long)(x+.5))%mod;
}
inline int qpow(int a,int b=P-2){
	int ans(1);for(;b;b>>=1,a=1ll*a*a%mod)if(b&1)ans=1ll*ans*a%mod;return ans;
}
inline void NTT_init(const int&n){
	const int&m=Getlen(n)-1;complex*now=buf;w[m]=now;now+=1<<m;
	IMP(1<<m,w[m][i]=complex(std::cos(i*Pi/(1<<m)),std::sin(i*Pi/(1<<m))));
	for(int k=m-1;k>=0&&(w[k]=now,now+=1<<k);--k)IMP(1<<k,w[k][i]=w[k+1][i<<1]);
}
inline void DFT(complex*f,const int&M){
	const int&n=1<<M;
	for(int len=n>>1,d=M-1;d>=0;--d,len>>=1)for(int k=0;k^n;k+=len<<1){
		complex*W=w[d],*L=f+(k),*R=f+(k|len),x,y;IMP(len,(x=*L,y=*R)),*L++=(x+y),*R++=*W++*(x-y);
	}
}
inline void IDFT(complex*f,const int&M){
	const int&n=1<<M;
	for(int len=1,d=0;d^M;++d,len<<=1)for(int k=0;k^n;k+=len<<1){
		complex*W=w[d],*L=f+(k),*R=f+(k|len),x,y;IMP(len,(x=*L,y=*W++**R)),*L++=(x+y),*R++=(x-y);
	}
	IMP(n,(f[i].x/=n,f[i].y/=n));for(int i=1;(i<<1)<n;++i)swap(f[i],f[n-i]);
}
inline void MTT(int*f,int*g,int*h,const int&n,const int&m,const int&t,const int&Len=-1){
	static complex Q[M],P[M],T[M];const int&len=Getlen(!~Len?n+m-1:Len);
	IMP(n,(Q[i].x=f[i]&32767,P[i].x=f[i]>>15));IMP(m,T[i]=complex(g[i]&32767,g[i]>>15));
	DFT(Q,len);DFT(P,len);DFT(T,len);IMP(1<<len,(Q[i]=Q[i]*T[i],P[i]=P[i]*T[i]));IDFT(Q,len);IDFT(P,len);
	IMP(t,h[i]=(Get(Q[i].x)+(1ll*Get(Q[i].y+P[i].x)<<15)+(1ll*Get(P[i].y)<<30))%mod);IMP(1<<len,Q[i]=P[i]=T[i]=0);
}
inline void Getinv(int*f,const int&n){
	static int pre[M];pre[0]=f[0];for(int i=1;i<n;++i)pre[i]=1ll*pre[i-1]*f[i]%mod;
	int t,c=qpow(pre[n-1]);for(int i=n-1;i>=1;--i)t=f[i],f[i]=1ll*pre[i-1]*c%mod,c=1ll*c*t%mod;f[0]=c;
}
inline void PT(int*f,int*g,const int&n,const int&m){
	static int F[M],G[M],H[M];H[0]=1;IMP(n,H[0]=1ll*H[0]*(m-i)%mod);IMP(n+n,G[i]=m+i-n);G[0]=1;Getinv(G,n+n);
	IMP(n,F[i]=1ll*ifac[i]*(n-i-1&1?P-ifac[n-i-1]:ifac[n-i-1])%mod*f[i]%mod);
	const int&len=Getlen(n+n);for(int i=1;i<n;++i)H[i]=1ll*(m+i)*G[i]%mod*H[i-1]%mod;MTT(F,G,F,n,n+n,n+n,n+n);
	IMP(n,g[i]=1ll*F[n+i]*H[i]%mod);IMP(1<<len,F[i]=G[i]=H[i]=0);
}
inline int Getfac(const int&n,const int&p){
	static int F[M],G[M];const int&B=sqrt(n),m=Getlen(B-1)-1;
	mod=Barrett(P=p);ifac[0]=ifac[1]=1;for(int i=2;i<=B;++i)ifac[i]=1ll*(P-P/i)*ifac[P%i]%mod;
	for(int i=1;i<=B;++i)ifac[i]=1ll*ifac[i-1]*ifac[i]%mod;F[0]=1;F[1]=B+1;
	for(int i=m;i>=0;--i){
		const int&q=B>>i+1,p=B>>i;if(!q)continue;
		PT(F,G,q+1,q+1);IMP(q,F[i+q+1]=G[i]),G[i]=0;G[q]=0;
		PT(F,G,q*2+1,1ll*q*qpow(B)%mod);IMP(q*2+1,F[i]=1ll*F[i]*G[i]%mod),G[i]=0;
		if(q*2+1==p){
			IMP(q*2+1,F[i]=1ll*F[i]*(i*B+p)%mod);F[p]=1;IMP(p,F[p]=1ll*F[p]*(p*B+i+1)%mod);
		}
	}
	int ans(1);IMP(B,ans=1ll*ans*F[i]%mod),F[i]=0;F[B]=0;for(int i=B*B+1;i<=n;++i)ans=1ll*ans*i%mod;return ans;
}
signed main(){
	int T,N,P;NTT_init(1<<17);scanf("%d",&T);while(T--)scanf("%d%d",&N,&P),printf("%d\n",Getfac(N,P));
}