Skip to content

浅谈扩展卢卡斯定理

约 1345 字大约 4 分钟

数学卢卡斯定理

2018-02-15

扩展卢卡斯定理,嗯~~~ 首先要需要学会扩展欧几里得算法中国剩余定理。至于卢卡斯定理,它和扩展卢卡斯定理并没什么联系吧……

扯dan

(建议跳过) 首先讲一下扩展卢卡斯用于什么地方。 对于组合数取模C(n,m)%p当n和m范围极大但模数p较小时我们会考虑用卢卡斯定理来求解,但使用卢卡斯定理要求p是一个素数。所以当p不是素数是怎么办呢?这时就要用到扩展卢卡斯定理了,

浅谈扩展Lucas定理

我们知道:p=pikip=\prod p_i^{k_i} 那么显然,如果我们知道C(n,m)C(n,m)modpiki\mod p_i^{k_i}意义下的答案,那么显然可以用CRT(中国剩余定理)进行合并。所以现在我们只需考虑C(n,m)modpikiC(n,m) \mod p_i^{k_i}即可。 我们回归到组合数取模的最简单最暴力的方法。那就是:

n!m!(nm)!(modpiki)\frac{n!}{m!(n-m)!} \pmod{p_i^{k_i}}

我们将n!n!m!m!(nm)!(n-m)!三者分开来看,暴力算是O(n)O(n)的,不但优雅到起飞,而且还会各种出锅。那么开始考虑优化。分成三部分没有错,我们要分别处理n!(modpiki)n!\pmod{p_i^{k_i}}以及inv(m!,piki)inv(m!, p_i^{k_i})inv((nm)!,piki)inv((n-m)!, p_i^{k_i})inv(a,b)inv(a,b)表示a在模b意义下的逆元)。
我们先考虑n!modpikin! \mod p_i^{k_i}n!n!modpiki\mod p_i^{k_i}的意义下,可以看做一个循环节连乘乘积的若干次幂乘上**pip_i的若干次幂再乘上一些奇奇怪怪的东西**(其实也不奇怪)。嗯,这样讲的确十分令人费解,举个栗子(网上讲扩展Lucas的都是这个栗子):
n=19pi=3ki=2n=19 \quad p_i=3 \quad k_i=2
n!=12341819n!=1*2*3*4*\cdots*18*19
n!=(124578)(101113141617)36(123456)19n!=(1*2*4*5*7*8)*(10*11*13*14*16*17)*3^6*(1*2*3*4*5*6)*19
因为 124578101113141617 (mod32)1*2*4*5*7*8 \equiv 10*11*13*14*16*17 \ (\mod 3^2)
所以 n!=(124578)2366!19mod32n!=(1*2*4*5*7*8)^2*3^6*6!*19 \mod 3^2(不要问:这一串不是应该等于0吗?这样的问题,后面会解释)
然后看到 6!6! 可以直接把它递归处理掉,而余下的19就是上面所说的奇怪的部分,它就是凑不齐一个循环节的部分,根据栗子不难看出循环节的长度是小于pikip_i^{k_i}的,而凑不齐的就跟少了,所以我们就可以直接暴力求了,嗯,优雅。 然后我们看剩下的两个逆元的部分,因为m!m!(nm)!(n-m)!中可能含有pip_i因此不保证互质,所以可能求不出逆元,所以我们要将m!m!(nm)!(n-m)!中的所有pip_i都先提出来用n!n!中所有的pip_i除以它,再把所得结果乘回去。这就是上文栗子里的363^6(实际上不是363^6而是383^8,要将6!6!先递归处理掉)。也就是说实际的计算方法是这样的:

n!inv(m!,piki)inv((nm)!,piki)piq1q2q3modpikin!*inv(m!,p_i^{k_i})*inv((n-m)!,p_i^{k_i})*p_i^{q1-q2-q3} \mod p_i^{k_i}

这里的阶乘都不是完整的阶乘而是完整的阶乘中除去了完整的阶乘中所包含的所有的pip_i,并把所有pip_i都提到后面去统一计算。 不懂的可以看看代码:

#include<iostream>
#include<algorithm>
#include<cstdio>
#include<cstring>
#define ll long long
using namespace std;
ll n,m,mod;
inline ll pow(ll base,ll x,ll p){
	ll ans=1;
	while(x){
		if(x&1) ans=ans*base%p;
		base=base*base%p;
		x>>=1;
	}
	return ans;
}
inline ll exgcd(ll a,ll b,ll &x,ll &y){
	if(!b){
		x=1;y=0;
		return a;
	}
	ll r=exgcd(b,a%b,x,y);
	ll t=x;
	x=y;
	y=t-a/b*y;
	return r;
}
inline ll inv(ll a,ll b){
	ll x,y;
	return exgcd(a,b,x,y)==1?(x+b)%b:-1;
}
inline ll CRT(ll x,ll p){return x*inv(mod/p,p)%mod*mod/p%mod;}
inline ll fac(ll x,ll p,ll pp){
	if(!x) return 1;
	ll res=1;
	for(int i=2;i<=p;i++) if(i%pp) res=res*i%p;
	res=pow(res,x/p,p);
	for(int i=2;i<=x%p;i++) if(i%pp) res=res*i%p;
	return res*fac(x/pp,p,pp)%p;//递归处理n/pi
}
inline ll C(ll n,ll m,ll p,ll pp){
	ll fn=fac(n,p,pp),fm=fac(m,p,pp),fnm=fac(n-m,p,pp);
	ll k=0;
	for(ll i=n;i;i/=pp) k+=i/pp;//计算n!中pi的个数
	for(ll i=m;i;i/=pp) k-=i/pp;
	for(ll i=n-m;i;i/=pp) k-=i/pp;
	return fn*inv(fm,p)%p*inv(fnm,p)%p*pow(pp,k,p)%p;
}
inline ll exlucas(ll n,ll m){
	ll t=mod,ans=0;
	for(int i=2;i*i<=mod;i++)
		if(t%i==0){
			ll p=1;
			while(t%i==0) p*=i,t/=i;
			ans=(ans+CRT(C(n,m,p,i),p))%mod;
		}
	if(t>1) ans=(ans+CRT(C(n,m,t,t),t)%mod);
	return ans%mod;
}
int main(){
	scanf("%lld%lld%lld",&n,&m,&mod);
	printf("%lld\n",exlucas(n,m));
	return 0;
}

接着扯dan

扩展Lucas适用在模数在1e7左右吧…… 其实如果模数比较大的话,理论上可以先将模数分解成若干个互质的因子,对每个因子在分别扩展Lucas,最后CRT和并(其实可以分解成不同的因子也许也行,可以exCRT和并)。 但是如果模数是个很大的质数怎么办?嗯,这个问题留给读者自己思考。

水平有限,讲的不好,还请指教,感激不尽。