扩展Lucas定理 学习笔记

扩展Lucas定理学习笔记

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

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

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

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

\[p=\prod p_i^{c_i} \]

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

考虑求

\[\binom{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!\)为

\[\begin{aligned} n!&=(p\times2p\times3p\times\dots(\lfloor\frac np\rfloor)p)(1\times2\times3\dots\times n)\\ &=p^{\lfloor\frac np\rfloor}(\lfloor\frac np\rfloor)!\prod_{i=1}^{n}[p\nmid i]i \end{aligned} \]

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

\[g(n)=\lfloor\frac np\rfloor+g(\lfloor\frac np\rfloor) \]

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

\[f(n)=f(\lfloor\frac np\rfloor)\prod_{i=1}^{n}[p\nmid i]i \]

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

\[(\prod_{i=1}^{p^k}[p\nmid i]i)^{\lfloor\frac {n}{p^k}\rfloor}(\prod_{i=p^k}^n [p\nmid 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];
	}
}

最后，

\[ans\equiv\frac{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)\)，模板题的\(p\le 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\n",exLucas(n,m,mod));	
	return 0;
}