exLucas

前置知识

CRT，Lucas。

exLucas

求 \(\binom nm\bmod p\)，\(1\le n,m\le 10^{18}\)，\(1\le p\le 10^6\)，不保证 \(p\) 为质数。

先将 \(p\) 分解质因数，设 \(p=\prod p_i^{\alpha_i}\)。

考虑这样一个方程组：

\[\begin{cases} \binom nm\equiv ans_1\pmod{p_1^{\alpha_1}}\\ \binom nm\equiv ans_2\pmod{p_2^{\alpha_2}}\\ \cdots\\ \binom nm\equiv ans_3\pmod{p_3^{\alpha_3}}\\ \end{cases} \]

如果我们能求出每一个 \(ans_i\)，那么我们就可以通过 CRT 合并出答案。

如果 \(\alpha_i=1\)

此时方程组就变成了 \(\binom nm\equiv ans_i\pmod{p_i}\)

\(ans_i\) 可通过 Lucas 求出。

一般情况

考察这个式子：

\[\binom nm\bmod p^{\alpha} \]

根据组合数定义，原式可变为：

\[\frac{n!}{m!(n-m)!}\bmod p^{\alpha} \]

然而 \(m!\) 和 \((n-m)!\) 的逆元没法求。

设 \(g(x)\) 表示将 \(x!\) 分解质因数后 \(p\) 的次数（即 \(g(x)=\left\lfloor\frac{x}{p}\right\rfloor+\left\lfloor\frac{x}{p^2}\right\rfloor+\left\lfloor\frac{x}{p^3}\right\rfloor+\cdots\)，可以 \(\mathcal O(\log_p x)\) 求出），则原式变为：

\[\left(\frac{\frac{n!}{p^{g(n)}}}{\frac{m!}{p^{g(m)}}\cdot\frac{(n-m)!}{p^{g(n-m)}}}\cdot p^{g(n)-g(m)-g(n-m)}\right)\bmod p^{\alpha} \]

可以发现 \(\frac{m!}{p^{g(m)}}\) 和 \(\frac{(n-m)!}{p^{g(n-m)}}\) 的逆元是可以求的，因为它们没有 \(p\) 这个因子。

于是我们的目标就变成了求这个式子：

\[f(n)=\frac{n!}{p^k}\bmod p^{\alpha} \]

把 \(n!\) 拆开：

\[n!=(p\cdot 2p\cdots)\cdot(1\cdot 2\cdots) \]

就是说把 \(p\) 的倍数分一组，其它的分一组，然后：

\[n!=p^{\left\lfloor\frac np\right\rfloor}\left(\left\lfloor\frac np\right\rfloor\right)!\prod_{1\le i\le n,p\nmid i}i \]

这个式子分为三个部分。第二部分可以递归求解（边界：\(f(0)=1\)）。递归到下一层的时候，第一部分的次数是 \(\left\lfloor\frac{n}{p^2}\right\rfloor\)，递归到再下一层的时候，第一部分的次数是 \(\left\lfloor\frac{n}{p^3}\right\rfloor\)，于是累积下来第一部分的次数就变成了 \(\left\lfloor\frac{n}{p}\right\rfloor+\left\lfloor\frac{n}{p^2}\right\rfloor+\left\lfloor\frac{n}{p^3}\right\rfloor+\cdots\)，它刚好是 \(g(n)\)，也就是 \(k\)，所以最终第一部分被消掉了。现在我们有：

\[f(n)=\left(f\left(\left\lfloor\frac np\right\rfloor\right)\prod_{1\le i\le n,p\nmid i}i\right)\bmod p^{\alpha} \]

下面我们来看第三部分。可以发现它模 \(p^{\alpha}\) 是有循环节的：

\[\prod_{1\le i\le n,p\nmid i}i=\left(\prod_{1\le i\le p^{\alpha},p\nmid i}i\right)^{\left\lfloor\frac{n}{p^{\alpha}}\right\rfloor}\left(\prod_{p^{\alpha}\left\lfloor\frac{n}{p^{\alpha}}\right\rfloor\le i\le n,p\nmid i}i\right) \]

到这里就可以暴力算了。

Code

#include<cstdio>
#include<algorithm>
#include<cmath>

#define int long long

const int N=1e6;

int A[N+10],M[N+10];

int gcd(int a,int b){return !b?a:gcd(b,a%b);}
inline int lcm(int a,int b){return a/gcd(a,b)*b;}

int exgcd(int a,int b,int &x,int &y){
	if(b==0)return x=1,y=0,a;
	int tx,ty,res=exgcd(b,a%b,tx,ty);
	x=ty,y=tx-a/b*ty;return res;
}
inline int inv(int a,int p){
	int x,y,tmp=exgcd(a,p,x,y);
	if(tmp!=1)return -1;else return (x%p+p)%p;
}

inline int qpow(int a,int b,int p){
	int res=1;
	while(b){if(b&1)res=1LL*res*a%p;a=1LL*a*a%p;b>>=1;}
	return res;
}

void merge(int &a1,int &b1,int a2,int b2){
	if(b2>b1)std::swap(a1,a2),std::swap(b1,b2);
	while(a1%b2!=a2)a1+=b1;
	b1=lcm(b1,b2);
}

std::pair<int,int> CRT(int *a,int *b,int num){
	int x=a[1],y=b[1];
	for(int i=2;i<=num;i++)merge(x,y,a[i],b[i]);
	return std::make_pair(x,y);
}

int f(int a,int p,int pk){
	if(a==0)return 1;
	int tmp1=1,tmp2=1,res=1;
	for(int i=1;i<=pk;i++)
		if(i%p!=0)tmp1=tmp1*i%pk;
	for(int i=pk*(a/pk);i<=a;i++)
		if(i%p!=0)tmp2=tmp2*(i%pk)%pk;
	res=qpow(tmp1,a/pk,pk)*tmp2%pk;
	res=res*f(a/p,p,pk)%pk;
	return res;
}

inline int g(int a,int p){
	int res=0;
	for(int i=p;i<=a;i*=p)res+=a/i;
	return res;
}

inline int comb(int a,int b,int p,int pk){
	int x=g(a,p),y=g(b,p),z=g(a-b,p);
	int xx=f(a,p,pk),yy=f(b,p,pk),zz=f(a-b,p,pk);
	int res=xx*inv(yy,pk)%pk*inv(zz,pk)%pk;
	res=res*qpow(p,x-y-z,pk)%pk;
	return res;
}

inline int C(int a,int b,int p){
	int tot=0;
	for(int i=2;i<=(int)(sqrt(p));i++)
		if(p%i==0){
			int palpha=1;
			while(p%i==0)palpha*=i,p/=i;
			A[++tot]=comb(a,b,i,palpha),M[tot]=palpha;
		}
	if(p>1)A[++tot]=comb(a,b,p,p),M[tot]=p;
	return CRT(A,M,tot).first;
}

signed main(){
	int n,m,p;scanf("%lld%lld%lld",&n,&m,&p);
	printf("%lld\n",C(n,m,p));
	return 0;
}