[模板] 扩展Lucas定理

扩展Lucas定理

为什么需要扩展Lucas定理

我们都知道，Lucas定理是用来求组合数取模一类问题的定理，但是当模数 p 为合数时，我们就无法使用定理来递归求解。

---

扩展Lucas定理的核心思想

唯一分解定理 + CRT(中国剩余定理)

说人话便是求出每个C(n,m)%p^k,在根据中国剩余定理求解同余方程组的解 x 。

---

具体求法及推倒过程

C(m,n)mod p^k (p是质数)
=(n!/P^x)/((m!/P^y)/((n-m)!/P^z))*P^(x-y-z) mod P^k
=......(不好敲了)

我们需要把x,y,z求解出来便是本问题的核心

不妨先来举个例子:

22!
=(1 * 2 * 4 * 5 * 7 * 8)(10 * 11 * 13 * 14 * 16 * 17)(19 * 20 * 22)(3 * 6 * 9 * 12 * 15 * 18 * 21) (mod 3^2)
=3^7 * 7!(1 * 2 * 4 * 5 * 7 * 8)^2(19 * 20 * 22)

它由三部分:3的项，一个阶乘，一个循环节和余项

之所以说是三项是因为7！还可以递归，还有3的项可以提取出来，并且这个个数是（n/p)。

我们令 f[ n ] = n!/P = f[ (n/p) ] * 当前循环节 * 余项

g[ n ] = g[ (n/p) ] + (n/p),意义是当前阶乘可以拆得的含P的项的个数

#include <iostream>
#include <cstdio>
#include <cmath>
#define ll long long
#define re register
using namespace std;
inline void exgcd(ll a,ll b,ll &x,ll &y){
	if(!b){x=1;y=0;return;}
	exgcd(b,a%b,x,y);
	ll tmp=x;x=y;y=tmp-(a/b)*y;
}
inline ll gcd(ll a,ll b){
	if(!b)return a;
	return gcd(b,a%b);
}
inline ll inv(ll a,ll p){
	ll x,y;
	exgcd(a,p,x,y);
	return (x+p)%p;
}
inline ll lcm(ll a,ll b){
	return a/gcd(a,b)*b;
}
inline ll qmul(ll a,ll b,ll mod){
	ll t=0;a%=mod;b%=mod;
	while(b){
		if(b&1LL)t=(t+a)%mod;
		b>>=1LL;a=(a+a)%mod;
	}
	return t;
}
inline ll power(ll a,ll b,ll mod){
	if(!b)return 1; 
	ll k=power(a,b/2,mod);
	k=k*k%mod;
	if(b&1)k=k*a%mod;
	return k;
}
inline ll F(ll n,ll P,ll PK){
	if(n==0)return 1;
	ll rou=1;//循环节 
	ll rem=1;//余项
	for(re ll i=1;i<=PK;i++)if(i%P)rou=rou*i%PK;
	rou=power(rou,n/PK,PK);
	
	for(re ll i=PK*(n/PK);i<=n;i++)if(i%P)rem=rem*(i%PK)%PK;
	return F(n/P,P,PK)*rou%PK*rem%PK; 
}
inline ll G(ll n,ll P){
	if(n<P)return 0;
	return G(n/P,P)+(n/P);
}
inline ll C_PK(ll n,ll m,ll P,ll PK){
	ll fz=F(n,P,PK),fm1=inv(F(m,P,PK),PK),fm2=inv(F(n-m,P,PK),PK);
	ll mi=power(P,G(n,P)-G(m,P)-G(n-m,P),PK);
	return fz*fm1%PK*fm2%PK*mi%PK; 
} 
ll A[1005],B[1005];//A存 
//x=B(mod)
//main algorithm
inline ll exLucas(ll n,ll m,ll P){//并不在意p，而是在意p^k 
	ll ljc=P,cnt=0;
	for(re ll i=2;i*i<=P;i++){
		if(ljc%i==0){ 
			ll PK=1;
			while(ljc%i==0)PK*=i,ljc/=i;
			A[++cnt]=PK;B[cnt]=C_PK(n,m,i,PK);
		}
	}
	if(ljc>1){
		A[++cnt]=ljc;B[cnt]=C_PK(n,m,ljc,ljc);
	} 
	//CRT->last answer
	ll ans = 0;
	for(re ll i=1;i<=cnt;i++){
		ll M=P/A[i],T=inv(M,A[i]);
		ans=(ans+B[i]*M%P*T%P)%P;
	}
	return ans;
}
int main(){

	ll n,m,P;
	scanf("%lld%lld%lld",&n,&m,&P);
	printf("%lld",exLucas(n,m,P));
	return 0;
}