ex_Lucas定理

Lucas定理(p为质数)：

\(C_n^m=C_{n/p}^{m/p}*C_{n\ mod\ p}^{m\ mod\ p}\)

可是p不为质数怎么办呢？

ex_Lucas定理 (p不为质数)

• 思路
  因为Lucas定理只能解决质数的情况，于是我们把P分解质因数, \(P=mul(p^k)\)
  然后对于每个\(p^k\)求出对应的\(md\ =C_n^m\ mod\ pk\),然后用CRT(中剩)合并出最后的答案,是不是非常有趣
  \(C_n^m=\frac{n!} {(n-m)!m!}\) (mod \(p^k\)意义下)
  因为要求逆元但是\((n-m)!m!\)可能不与\(p^k\)互质,所以我们提出，质因数\(p\):
  \(C_n^m=\frac{\frac{n!}{p^a}} {\frac{(n-m)!}{p^b}\ \ \frac{m!}{p^c}}*p^{a-b-c}\)

所以求:\(n!\ mod\ p^k\)
乘积显然会以\(p^k\)为循环节,注意乘的时候要记得跳过\(p\)的倍数,然后最后一段不完整的组暴力处理
注意代码处理中我们会直接跳过p的倍数但是,对于 \(k\)属于[\(1\),\(n/p\)] , \(p*k\)中\(k\)会被我们跳过,因此我们采用递归的手段，因此多处理\((n/p)!\ mod\ p^k\)

• 代码:

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=1e7+5;
ll r[N],md[N],jc[N];
ll ksm(ll t,ll b,ll mod) {
	ll mul=1;
	for(ll i=b;i;i>>=1,t=t*t%mod) if(i&1) mul=mul*t%mod;
	return mul;
}
ll ex_gcd(ll a,ll b,ll &x,ll &y) {
	if(b==0) {x=1; y=0; return a;}
	ll d=ex_gcd(b,a%b,x,y);
	ll xx=x; x=y; y=xx-a/b*y;
	return d;
}
ll inv(ll u,ll p) {
	ll v,k,d=ex_gcd(u,p,v,k);
	v%=(p/d);
	if(v<0) v+=(p/d);
	return v;
}
ll fac(ll n,ll pk,ll p) {
	if(n==0||n==1) return 1;
	return fac(n/p,pk,p)*ksm(jc[pk],n/pk,pk)%pk*jc[n%pk]%pk;
}
ll Cnt(ll n,ll p) {
	ll sum=0;
	for(ll i=n;i;i/=p) sum+=i/p;
	return sum;
}
ll C_modpk(ll n,ll m,ll p,ll pk)	{
	jc[0]=1; 
	for(ll i=1;i<=pk;i++) 	//加速技巧
		if(i%p) jc[i]=jc[i-1]*i%pk;
		else jc[i]=jc[i-1];
	return fac(n,pk,p)*inv(fac(m,pk,p),pk)%pk*inv(fac(n-m,pk,p),pk)%pk*ksm(p,(Cnt(n,p)-Cnt(m,p)-Cnt(n-m,p)),pk)%pk;
}
ll CRT(int tot) {
	ll mul=1,sum=0,x,y;
	for(int i=1;i<=tot;i++) mul*=md[i];
	for(int i=1;i<=tot;i++) {
		ll d=ex_gcd(mul/md[i],md[i],x,y);
		sum+=(mul/md[i])*x%mul*r[i]%mul%mul,sum=(sum%mul+mul)%mul;
	}
	return sum;
}
ll ex_Lucus(ll n,ll m,ll P) {
	int tot=0;
	for(ll p=2;p*p<=P;p++) {
		ll pk=1;
		if(!(P%p)) {
			while(!(P%p)) {P/=p; pk*=p; }
			md[++tot]=pk; r[tot]=C_modpk(n,m,p,pk);
		}
	} 
	if(P>1) {
		md[++tot]=P; r[tot]=C_modpk(n,m,P,P);
	}
	ll ans=CRT(tot);
	return ans;
}
int main() {
	ll n,m,p;
	scanf("%lld%lld%lld",&n,&m,&p);
	if(n<m) swap(n,m);
	printf("%lld",ex_Lucus(n,m,p));
	return 0;
}