扩展Lucas定理 学习笔记

扩展Lucas定理学习笔记

前置知识：Lucas定理，中国剩余定理

一般的(Lucas)定理求解(inom{n}{m}pmod p)时，要求(p)是质数，事实上如果(p)不是质数，用阶乘逆元求解组合数也无法实现，因为当(p)不是质数时，很可能不会存在逆元((a)关于(p)存在逆元要求(gcd(a,p)=1)，至于为什么考虑逆元定义的同余方程在(gcd(a,p) ot=1)时无解)

而扩展(Lucas)，就是解决(p)不为质数时的问题。

首先，对(p)质因数分解：

[p=prod p_i^{c_i} ]

那么接下来我们可以分别求出所有的答案关于(p_i^{c_i})取模的结果，然后用中国剩余定理合并即可。

考虑求

[inom{n}{m}pmod {p^k}=frac{n!}{m!(n-m)!}pmod {p^k} ]

我们最大的问题还是没有逆元，但我们可以对于(x)，将所有(p)质因子除掉，那么剩下的部分就有逆元了！即：

[ans=frac{frac{n!}{p^x}}{frac{m!}{p^y}frac{(n-m)!}{p^z}}p^{x-y-z}pmod {p^k} ]

我们设(g(n))表示(n!)中因子(p)的次数，(f(n))表示

[frac{n!}{p^{g(n)}}pmod {p^k} ]

既然要除(p)的次幂，那我们就把(p)的倍数与其他数分开呗，于是我们拆(n!)为

[egin{aligned} n!&=(p imes2p imes3p imesdots(lfloorfrac np floor)p)(1 imes2 imes3dots imes n)\ &=p^{lfloorfrac np floor}(lfloorfrac np floor)!prod_{i=1}^{n}[p mid i]i end{aligned} ]

前面的(p)的次幂，直接计入(g)中，但((lfloorfrac np floor)!)中可能还有(p)的倍数，因此我们递归下去处理：

[g(n)=lfloorfrac np floor+g(lfloorfrac np floor) ]

同时这些(p)次幂在(f)中也会被除掉，但((lfloorfrac np floor)!)中还有一些也要除掉，我们同样递归处理：

[f(n)=f(lfloorfrac np floor)prod_{i=1}^{n}[p mid i]i ]

最后这个连乘看起来比较棘手，但显然在(mod p^k)意义下它是有循环节的，于是将它分解为若干个长度为(p^k)的段以及最后剩下的一段，即：

[(prod_{i=1}^{p^k}[p mid i]i)^{lfloorfrac {n}{p^k} floor}(prod_{i=p^k}^n [p mid i]i) ]

这部分暴力做即可，单次复杂度是(mathcal O(p^k))的，因为一共会递归(log_p n)次，所以总复杂度是(mathcal O(p^klog_pn))。

代码如下：

inline int F(ll n,int p,int pk){
	if(!n) return 1;
	int rec=1,rep=1;
	for(ll i=1;i<pk;++i)
		if(i%p!=0) rep=1ll*rep*i%pk;
	for(ll i=1ll*pk*(n/pk);i<=n;++i)
		if(i%p!=0) rec=1ll*rec*(i%pk)%pk;
	return 1ll*F(n/p,p,pk)*ksm(rep,n/pk)%pk*rec%pk;
}
inline int G(ll n,int p,int pk){
	if(n<p) return 0;
	return n/p+G(n/p,p,pk);
}

你会发现(rep)与(rec)的求解每次都是类似的，可以通过在计算前预处理达到(mathcal O(p^k))，这个尤其在模数确定多组数据时极其重要。

inline int F(ll n,int p,int pk){
	if(!n) return 1;
	return 1ll*F(n/p,p,pk)*ksm(rep[pk],n/pk)%pk*rep[n%pk]%pk;
}
inline int G(ll n,int p,int pk){
	if(n<p) return 0;
	return n/p+G(n/p,p,pk);
}
inline void init(int p,int pk){
	rep[0]=1;
	for(int i=1;i<=pk;++i){
		if(i%p!=0) rep[i]=1ll*rep[i-1]*i%pk;
		else rep[i]=rep[i-1];
	}
}

最后，

[ansequivfrac{f(n)}{f(m)f(n-m)}p^{g(n)-g(m)-g(n-m)}pmod {p^k} ]

(CRT)合并即可，总复杂度为(mathcal O(sum{p_i^{k_i}})le mathcal O(p))，模板题的(ple 10^6)稳过。也能作模数比较大但十分特殊的情况。

模板题代码：

#include<bits/stdc++.h>
using namespace std;
#define pb push_back
#define mp make_pair
typedef long long ll;
typedef pair<int,int> pii;
const int N=1e6+10;
ll n,m;
int mod;
vector<pii> ve;
inline int ksm(int x,ll y){
	int ret=1;
	for(;y;x=1ll*x*x%mod,y>>=1) (y&1)&&(ret=1ll*ret*x%mod);
	return ret;
}
inline void exgcd(int a,int b,int &x,int &y){
	if(b==0){
		x=1;y=0;
		return ;
	}
	exgcd(b,a%b,x,y);
	int t=x;x=y;y=t-a/b*y;
}
int rep[N];
inline int F(ll n,int p,int pk){
	if(!n) return 1;
	return 1ll*F(n/p,p,pk)*ksm(rep[pk],n/pk)%pk*rep[n%pk]%pk;
}
inline int G(ll n,int p,int pk){
	if(n<p) return 0;
	return n/p+G(n/p,p,pk);
}
inline void init(int p,int pk){
	rep[0]=1;
	for(int i=1;i<=pk;++i){
		if(i%p!=0) rep[i]=1ll*rep[i-1]*i%pk;
		else rep[i]=rep[i-1];
	}
}
inline int inv(int a,int mod){
	int x,y;
	exgcd(a,mod,x,y);
	return (x%mod+mod)%mod;
}
inline int binom_pk(ll n,ll m,int p,int pk){
	int ans=1ll*F(n,p,pk)*inv(F(m,p,pk),pk)%pk*inv(F(n-m,p,pk),pk)%pk;
	int c=G(n,p,pk)-G(m,p,pk)-G(n-m,p,pk);
	ans=1ll*ans*ksm(p,c)%pk;
	return ans;
}
inline int exLucas(ll n,ll m,int mod){
	int t=mod;
	for(int i=2;i*i<=mod;++i){
		if(t%i==0){
			int an=1;
			while(t%i==0) t/=i,an*=i;
			ve.pb(mp(i,an));
		}
	}
	if(t>1) ve.pb(mp(t,t));
	int ans=0;
	for(pii u:ve){
		int p=u.first,pk=u.second;
		init(p,pk);
		int ret=binom_pk(n,m,p,pk);
		int M=mod/pk;
		ret=1ll*ret*M%mod*inv(M%pk,pk)%mod;
		ans=(ans+ret)%mod;
	}
	return ans;
}
int main(){
//	freopen("P4720_10.in","r",stdin);
	scanf("%lld%lld%d",&n,&m,&mod);
	printf("%d
",exLucas(n,m,mod));	
	return 0;
}