【模板】快速阶乘算法/阶乘模大质数（拉格朗日插值）（MTT）

洛谷板子传送门

---

题解：

首先，很显然地，我们需要一个低于$O(n)$的做法。

考虑设多项式$g_d(x)=\prod _{i=1}^d(x+i)$

设$s=\lfloor \sqrt{n}\rfloor$，则我们求出$g_s(0),g_s(s),\dots ,g_s((s-1)s)$，再在最后暴力算不超过$O(\sqrt{n})$个数的乘积即可。

如果按照根号分块构建多项式然后计算的话可以用多点求值做到$O(\sqrt{n}{\log }^{2}n)$。

可能有人觉得可以通过改变块的大小来优化复杂度，实际上设块的大小为$S$，你前面求这个多项式还需要分治乘法，后面多点求值还有一个分治操作，也就是你实际上的复杂度是$O(S{\log }^{2}S+\frac{n}{S}{\log }^2\frac{n}{S})$，显然，复杂度最优当且仅当$S=O(\sqrt{n})$，且复杂度为$O(\sqrt{n}{\log }^{2}n)$

实际上还有更加优秀的$O(\sqrt{n}\log n)$的做法。

我们注意到让上面复杂度多一个$\log$ 的地方在于分治。

根据经验，能倍增的话会比分治少一个$\log$，考虑倍增。

注意我们最后要求的其实就是点值，我们考虑直接由点值得到点值，这个过程也就是插值。

$g_d(x)$是$d$次多项式需要$d+1$个点值。

假设现在我们有点值序列$P_d:g_d(0),g_d(s),g_d(2s),\dots ,g_d(ds)$，考虑支持两种操作：

1. $P_d$得到$P_{d+1}$
2. $P_d$得到$P_{2d}$

如果能够在较低的时间复杂度内支持这两个操作就可以比较方便地进行倍增。

首先考虑第一个操作：
容易注意到$g_{d+1}(x)=g_d(x)\cdot (x+d+1)$
然后$g_{d+1}((d+1)s)$直接暴力计算即可，反正不超过$O(\sqrt{n})$

考虑第二个操作：

容易注意到$g_{2d}(x)=g_d(x)g_d(x+d)$，所以我们实际上想要求的是两个部分：

$g_d((d+1)s),g_d((d+2)s),\dots ,g_d(2ds)$
$g_d(d),g_d(s+d),\dots ,g_d(2ds+d)$

给一点微调，我们发现实际上我们要解决的两个形式相同的问题：

1. 已知$g_d(0),\dots g_d(ds)$，求$g_d(0+(d+1)s),\dots g_d(ds+(d+1)s)$（注意这里比我们要用的多求了一项）
2. 已知$g_d(0),\dots g_d(2ds)$，求$g_d(0+d),\dots g_d(2ds+d)。$

设$h(i)=g_d(is)$，不难发现上述两个问题本质相同。

给定点值$h(0),h(1),\cdots h(n)$，求点值$h(\Delta +0),h(\Delta +1),\cdots h(\Delta +n)$

点值转点值，实际上这个操作就叫插值，你仔细想一想感觉牛顿插值不太行，考虑拉格朗日插值，注意我们现在要求的点值和原来的点值是肯定没有交集的，所以$x-i$肯定非0：

$$\begin{array}{rl}h(x)& =\sum _{i=0}^{d}h(i)\prod _{j\ne i}\frac{x-j}{i-j}\\ & =\sum _{i=0}^{d}h(i)\frac{\prod _{j=0}^d(x-j)}{(x-i)\prod _{k\ne i}(i-k)}\\ & =(\prod _{j=0}^d(x-j))\left(\sum _{i=0}^d\frac{h(i)}{(x-i)i!(d-i)!(-1{)}^{d-i}}\right)\end{array}$$

前面的那个$\prod$显然可以直接在算出后面的东西之后上一个双指针。

所以考虑后面的那个东西，容易发现把$\frac{1}{x-i}$分开之后就是一个卷积，由于每一项都有贡献所以需要平移一下数组，由于平移后卷出来只会用到中间项的系数，所以直接利用循环卷积的特性把后面的丢到前面去来优化常数。由于丢人模数不能用NTT所以写一个MTT。

目前是luogu上加强数据之后的rk1（看日期，我前面的代码全部都是改数据前的。

---

代码：

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

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

int mod;
inline int add(int a,int b){a+=b-mod;return a+(a>>31&mod);}
inline int dec(int a,int b){return a-b<0?a-b+mod:a-b;}
inline int mul(int a,int b){ll r=(ll)a*b;return r>=mod?r%mod:r;}
inline void Inc(int &a,int b){a+=b-mod;a+=a>>31&mod;}
inline void Dec(int &a,int b){a-=b;a+=a>>31&mod;}
inline void Mul(int &a,int b){a=mul(a,b);}
inline void ex_gcd(int a,int b,int &x,int &y){
	if(!b){x=1,y=0;return ;}ex_gcd(b,a%b,y,x);y-=a/b*x;
}
int po(int a,int b){
	int r=1;for(;b;b>>=1,a=mul(a,a))
	if(b&1)r=mul(r,a);return r;
}
inline int inv(int a){
	return po(a,mod-2);
}

namespace MTT{
	struct cp{
		double x,y;cp(){}cp(double _x,double _y=0):x(_x),y(_y){}
		friend cp operator+(cs cp &a,cs cp &b){return cp(a.x+b.x,a.y+b.y);}
		friend cp operator-(cs cp &a,cs cp &b){return cp(a.x-b.x,a.y-b.y);}
		friend cp operator*(cs cp &a,cs cp &b){return cp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
		void operator+=(cs cp &b){x+=b.x,y+=b.y;}
		void operator-=(cs cp &b){x-=b.x,y-=b.y;}
		void operator*=(cs cp &b){*this=*this*b;}
		void operator*=(double b){x*=b,y*=b;}
		void operator/=(double b){x/=b,y/=b;}
		cp conj()cs{return cp(x,-y);}
	};cs double PI=acos(-1),PI2=2*PI;
	cp omega(int i,int k){return cp(cos(PI2*i/k),sin(PI2*i/k));}
	
	cs int bit=17,SIZE=1<<bit|7;
	int r[SIZE];cp *w[bit+1];
	void init_omega(){
		for(int re i=1;i<=bit;++i)w[i]=new cp[1<<(i-1)];
		for(int re d=1;d<=bit;++d){
			cp wn=omega(1,1<<d);
			for(int re i=0;i<(1<<(d-1));++i)
			w[d][i]=(i&31)?w[d][i-1]*wn:omega(i,1<<d);
		}
	}
	void DFT(cp *A,int l){
		for(int re i=0;i<l;++i)
		if(i<r[i])std::swap(A[i],A[r[i]]);
		for(int re i=1,d=1;i<l;i<<=1,++d)
		for(int re j=0;j<l;j+=i<<1)
		for(int re k=0;k<i;++k){
			cp &t1=A[j+k],&t2=A[i+j+k];
			cp t=t2*w[d][k];t2=t1-t,t1+=t;
		}
	}
	void init_rev(int len){
		for(int re i=1;i<len;++i)r[i]=r[i>>1]>>1|((i&1)?len>>1:0);
	}
	void mul(int *a,int *b,int len,int *c){
		static cp A[SIZE],B[SIZE],C[SIZE],D[SIZE];
		for(int re i=0;i<len;++i){
			A[i]=cp(a[i]&0x7fff,a[i]>>15);
			B[i]=cp(b[i]&0x7fff,b[i]>>15);
		}init_rev(len);DFT(A,len),DFT(B,len);
		for(int re i=0;i<len;++i){
			int u=(len-i)&(len-1);
			C[i]=(A[i].conj()+A[u])*cp(0.5,0)*B[u];
			D[i]=(A[i].conj()-A[u])*cp(0,0.5)*B[u];
		}DFT(C,len),DFT(D,len);
		for(int re i=0;i<len;++i){
			ll x=C[i].x/len+0.5,y=(C[i].y+D[i].x)/len+0.5,z=D[i].y/len+0.5;
			x%=mod,y%=mod,z%=mod;c[i]=(x+(y<<15)%mod+(z<<30)%mod)%mod;
		}
	}
}

cs int N=1<<18|7;
int ifc[N];
void init_fac(int n){
	int fac=1;for(int re i=2;i<=n;++i)Mul(fac,i);
	ifc[n]=inv(fac);
	for(int re i=n;i;--i)ifc[i-1]=mul(ifc[i],i);
}

void calc(int *a,int *b,int n,int k){
	static int f[N],g[N],h[N],p[N],ip[N];
	int len=1;while(len<=n+n)len<<=1;int t=dec(k,n);
	for(int re i=0;i<=n;++i)f[i]=mul(a[i],mul(ifc[i],ifc[n-i]));
	for(int re i=n-1;i>=0;i-=2)f[i]=mod-f[i];
	for(int re i=0;i<=n+n;++i)g[i]=add(i,t);
	p[0]=g[0];
	for(int re i=1;i<=n+n;++i)p[i]=mul(p[i-1],g[i]);
	ip[n+n]=inv(p[n+n]);
	for(int re i=n+n;i;--i)ip[i-1]=mul(ip[i],g[i]);
	g[0]=ip[0];
	for(int re i=1;i<=n+n;++i)g[i]=mul(ip[i],p[i-1]);
	for(int re i=n+1;i<len;++i)f[i]=0;
	for(int re i=n+n+1;i<len;++i)g[i]=0;
	MTT::mul(f,g,len,h);
	int res=1,p1=dec(k,n),p2=k;
	for(int re i=0;i<=n;++i)Mul(res,add(t,i));
	for(int re i=0;i<=n;++i)g[i]=add(p1,i);
	p[0]=g[0];
	for(int re i=1;i<=n;++i)p[i]=mul(p[i-1],g[i]);
	ip[n]=inv(p[n]);
	for(int re i=n;i;--i)ip[i-1]=mul(ip[i],g[i]);
	g[0]=ip[0];
	for(int re i=1;i<=n;++i)g[i]=mul(ip[i],p[i-1]);
	for(int re i=0;i<=n;Inc(p2,1),++i)
	b[i]=mul(h[i+n],res),Mul(res,mul(g[i],p2+1));
}

int solve(int bl){
	static int a[N],b[N];
	int s=0,iv=inv(bl),res=1;init_fac(bl);
	for(int re p=bl;p;p>>=1)++s;a[0]=1,--s;
	for(int re p=0;s>=0;--s){
		if(p){
			calc(a,b,p,p+1);
			for(int re i=0;i<=p;++i)a[p+1+i]=b[i];
			a[p<<1|1]=0;calc(a,b,p<<1,mul(p,iv));
			p<<=1;for(int re i=0;i<=p;++i)Mul(a[i],b[i]);
		}
		if(bl>>s&1){
			for(int re i=0;i<=p;++i)Mul(a[i],add(mul(bl,i),p+1));
			p|=1,a[p]=1;
			for(int re i=1;i<=p;++i)Mul(a[p],add(mul(bl,p),i));
		}
	}
	for(int re i=0;i<bl;++i)Mul(res,a[i]);
	return res;
}
int calc_fac(int n){
	int bl=sqrt(n),res=solve(bl);
	for(int re i=bl*bl+1;i<=n;++i)
	Mul(res,i);return res;
}
int fac(int n){
	if(n>mod-1-n){
		int res=inv(calc_fac(mod-1-n));
		return (mod-1-n)&1?res:mod-res;
	}return calc_fac(n);
}

int n;
void work(){
	MTT::init_omega();
	int T;scanf("%d",&T);
	while(T--){
		scanf("%d%d",&n,&mod);
		cout<<fac(n)<<"\n";
	}
}

void file(){
#ifdef zxyoi
	freopen("factorial.in","r",stdin);
#endif
}
signed main(){file();work();return 0;}