Lucas定理及其扩展

Lucas定理

定义

对于质数 \(p\)，有：$$\dbinom{n}{m} \mod p=\dbinom{n \mod p}{m \mod p} \dbinom{\lfloor \frac{n}{p} \rfloor}{\lfloor \frac{m}{p} \rfloor} \mod p$$

由于 \(n \mod p\) 和 \(m \mod p\) 都比模数 \(p\) 小，可以预处理，而 \(\tbinom{\lfloor \frac{n}{p} \rfloor}{\lfloor \frac{m}{p} \rfloor}\) 则可以再次调用 \(Lucas\) 定理求解。

时间复杂度：\(O(\log n)\)

意义

对于一般的组合数，我们有预处理 \(O(n)\) 和单次查询 \(O(1)\) 的阶乘算法，但是当 \(n\) 和 \(m\) 都特别大的时候，我们无法用数组来存那么多数的阶乘，于是使用牺牲时间换空间的方法。

\(Lucas\) 定理就是对于 \(n\) 和 \(m\) 特别大的时候，而模数 \(p\) 不是很大的时候，来求解组合数的算法。在求解过程中，发现只需要存 \(0 \sim p\) 的阶乘即可。但是单次查询的时间复杂度也由原来的 \(O(1)\) 变为 \(O(\log n)\)。

当然，当 \(n\) 和 \(m\) 以及 \(p\) 都非常大的时候，目前并没有更优的解法，只能用单次查询 \(O(n)\) 的暴力了。

Code

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N=2e4;
const int mod=1e4+7;
LL f[N+5];
void init(){
    f[0]=1;
    for(int i=1;i<=N;i++){
        f[i]=f[i-1]*i%mod;
    }
}
LL qpow(LL a,LL b){
    LL ans=1%mod;a%=mod;
    while(b){
        if(b&1)ans=ans*a%mod;
        a=a*a%mod;b>>=1;
    }
    return ans;
}
LL C(LL n,LL m){
    if(n<m)return 0;
    return f[n]*qpow(f[m],mod-2)%mod*qpow(f[n-m],mod-2)%mod;
}
LL Lucas(LL n,LL m,LL p){
    if(!m)return 1;
    return C(n%p,m%p)*lucas(n/p,m/p,p)%p;
}
int main(){
    init();
    int t;scanf("%d",&t);
    while(t--){
        LL n,m;scanf("%lld%lld",&n,&m);
        printf("%lld\n",Lucas(n,m,mod));
    }
    return 0;
}

扩展Lucas定理/exLucas定理

内容

\(exLucas\) 用来求解形如以下式子：$$\dbinom{n}{m} \mod p$$

其中，\(p\) 不一定是质数。

solution1

首先我们对 \(p\) 进行质因数分解：

\[p=\prod_{i=1}^k p_i^{c_i} \]

显然，每个 \(p_i^{c_i}\) 之间都是互质的。

solution2

那么现在的问题就变成了如果知道 \(a_i\)，求解：

\[\begin{cases} \dbinom{n}{m} \equiv a_1 (\mod p_1^{c_1}) \\ \dbinom{n}{m} \equiv a_2 (\mod p_2^{c_2}) \\ \cdots \\ \dbinom{n}{m} \equiv a_k (\mod p_k^{c_k}) \end{cases}\]

由于互质，因此可以使用中国剩余定理来求解。

solution3

那么现在我们来看 \(a_i=\tbinom{n}{m} \mod p_i^{c_i}\) 怎么求。

\[\dbinom{n}{m} \mod p_i^{c_i}=\frac{n!}{m!(n-m)!} \mod p_i^{c_i} \]

由于 \(m!\) 和 \((n-m)!\) 中可能包含因子 \(p_i\)，因此不能直接求在模 \(p_i^{c_i}\) 意义下的逆元，我们要先将 \(m!\) 和 \((n-m)!\) 中 的质因子 \(p_i\) 全部提出来，再乘回去：

\[\frac{\frac{n!}{p^{k_1}}}{\frac{m!}{p^{k_2}} \times \frac{(n-m)!}{p^{k_3}}} \times p^{k_1-k_2-k_3} \]

那么现在的 \(\frac{m!}{p^{k_2}}\) 和 \(\frac{(n-m)!}{p^{k_3}}\) 和 \(p^k\) 是互质的，可以直接求逆元。

solution4

那么我们再来看如何计算： $$\frac{n!}{p^a} \mod p^k$$

先考虑如何计算 \(n! \mod p^k\)。

举个例子，\(n=22,p=3,k=2\)

\(22!=1\times2\times3\times4\times5\times6\times7\times8\times9\times10\times11\times12\times13\times14\times15\times16\times17\times18\times19\times20\times21\times22\)

将 \(3\) 的倍数提出来。

\(22!=3^7\times(1\times2\times3\times4\times5\times6\times7)\times(1\times2\times4\times5\times7\times8\times10\times11\times13\times14\times16\times17\times19\times20\times22)\)

显然上式分为三个部分。

第一个部分：\(3^{\lfloor \frac{n}{p} \rfloor}\)

第二个部分：\(\lfloor \frac{n}{p} \rfloor !\)

第三个部分：\(n!\) 中与 \(p\) 互质的数的乘积。

对于第一个部分，先不管。

对于第二个部分，发现形式与 \(n!\) 一样，可以递归求解。

对于第三个部分，有如下性质：$$1\times2\times4\times5\times7\times8 \equiv 10\times11\times13\times14\times16\times17 (\mod p^k)$$

即 \(\prod\limits_{i,(i,p)=1}^{p^k}\) \(i \equiv \prod\limits_{i,(i,p)=1}^{p^k}\) \((i+t p^k) \mod p^k\)

\(\prod\limits_{i,(i,p)=1}^{p^k}\) \(i\) 一共重复了 \(\lfloor \frac{n}{p^k} \rfloor\) 次，暴力求出 \(\prod\limits_{i,(i,p)=1}^{p^k}\) \(i\)，然后快速幂即可。

而最后还要乘上 \(\prod\limits_{i,(i,p)=1}^{n \mod p^k}\) \(i\)，这一段直接暴力乘即可。

以上三个部分乘在一起就是 \(n!\)，但是最终求的是 \(\frac{n!}{p^a} \mod p^k\)，因此第一部分会和分母消掉。

Code

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N=1110000;
LL qpow(LL a,LL b,LL p){
	LL ans=1%p;a%=p;
	while(b){
		if(b&1)ans=ans*a%p;
		a=a*a%p;b>>=1;
	}
	return ans;
}
LL fac(LL n,LL p,LL pk){
	if(!n)return 1;
	LL ans=1;
	for(LL i=1;i<pk;i++){
		if(i%p!=0){
			ans=ans*i%pk;
		}
	}
	ans=qpow(ans,n/pk,pk);
	for(LL i=1;i<=n%pk;i++){
		if(i%p!=0){
			ans=ans*i%pk;
		}
	}
	return ans*fac(n/p,p,pk)%pk;
}
void exgcd(LL a,LL b,LL &d,LL &x,LL &y){
	if(!b){d=a;x=1;y=0;}
	else{
		exgcd(b,a%b,d,y,x);
		y-=(a/b)*x;
	}
}
LL inv(LL a,LL p){
	LL A,B,d,X,Y,K;
	A=a,B=p,K=1;
	exgcd(A,B,d,X,Y);
	X=X*(K/d);
	LL dx=(B/d);
	X=(X%dx+dx)%dx;
	return X;
}
LL C(LL n,LL m,LL p,LL pk){
	if(n<m)return 0;
	LL f1=fac(n,p,pk),f2=fac(m,p,pk),f3=fac(n-m,p,pk),cnt=0;
	for(LL i=n;i;i/=p){
		cnt=cnt+i/p;
	}
	for(LL i=m;i;i/=p){
		cnt=cnt-i/p;
	}
	for(LL i=n-m;i;i/=p){
		cnt=cnt-i/p;
	}
	return f1*inv(f2,pk)%pk*inv(f3,pk)%pk*qpow(p,cnt,pk)%pk;
}
LL a[N],c[N],cnt;
LL CRT(){
	LL m=1,ans=0;
	for(int i=1;i<=cnt;i++){
		m*=c[i];
	}
	for(int i=1;i<=cnt;i++){
		ans=(ans+a[i]*(m/c[i])%m*inv(m/c[i],c[i])%m)%m;
	}
	return ans;
}
LL exLucas(LL n,LL m,LL p){
	cnt=0;
	for(LL i=2;i*i<=p;i++){
		LL tmp=1;
		while(p%i==0){
			tmp*=i;
			p/=i;
		}
		if(tmp>1){
			cnt++;
			a[cnt]=C(n,m,i,tmp);
			c[cnt]=tmp;
		}
	}
	if(p>1){
        cnt++;
		a[cnt]=C(n,m,p,p);
		c[cnt]=p;
	}
	return CRT();
}
int main(){
	LL n,m,p;
	scanf("%lld%lld%lld",&n,&m,&p);
	printf("%lld",exLucas(n,m,p));
	return 0;
}