Lucas定理学习笔记

Lucas定理用于计算组合数取模$C_n^x\bmod p$，适用于$n$很大而$p$较小的情况。

Lucas定理

若要计算$C_n^m\bmod p$，其中$p$是素数。

Lucas定理

令$n=ap+b,m=cp+d$。
故：

$$C_n^m=C_{ap+b}^{cp+d}\equiv C_a^cC_b^d\pmod p$$

对于$C_b^d$，范围在$p$以内，可以通过阶乘算出来。
对于$C_a^c$，利用Lucas定理递归计算。

LL Quick_Pow(LL a,LL b) {
	LL ans=1;
	for(; b; b>>=1,a=a*a%p)if(b&1)ans=ans*a%p;
	return ans;
}

LL C(LL n,LL m) {
	if(m>n)return 0;
	return f[n]*inv[m]%p*inv[n-m]%p;
}

LL Lucas(LL n,LL m) {
	if(m==0)return 1;
	return (C(n%p,m%p)*Lucas(n/p,m/p))%p;
}

证明

Lucas定理我也不会证，贴个大佬的博客。

时间复杂度

若忽略阶乘（对阶乘进行预处理）与求逆元的时间（线性求逆元），那么一次计算的时间复杂度为$O(\log_pn)$。

扩展Lucas定理

Lucas定理只能用于求解组合数模素数，若$p$是一个非素数怎么办呢？
当$p=p_1^1\cdot p_2^1\cdots p_k^1$时，可以直接使用CRT合并。
但当$p=p_1^{a_1}p_2^{a_2}\cdots p_k^{a_k}$时，不能直接使用CRT合并，此时就需要用到扩展Lucas定理。

扩展Lucas定理

首先对分解后的$p$写出$k$个同余方程：

$$x\equiv b_1\pmod{p_1^{a_1}}$$ $$x\equiv b_2\pmod{p_2^{a_2}}$$ $$\cdots$$ $$x\equiv b_k\pmod{p_k^{a_k}}$$

若能够求出$b_i$，就可以使用CRT合并了。
因此现在考虑如何求出$x\bmod p_i^{a_i}$。
将组合数写成阶乘形式：

$$C_n^m\bmod p_i^{a_i}=\frac{n!}{m!(n-m)!}\bmod p_i^{a_i}$$

接下来我们考虑将$x!$进行拆分，下面我们用$16!\bmod 3^2$举例。

将$x!$的所有$p_i$的倍数单独拿出，提取公因数得到三部分：

$$16!=(1\times2\times4\times5\times7\times8\times10\times11\times13\times14\times16)\times3^5\times(1\times2\times3\times4\times5)$$

1. 第三部分，$1\times2\times3\times4\times5$，发现是子问题$\frac n{p_i}!$，递归处理。
2. 第一部分，$1\times2\times4\times5\times7\times8\times10\times11\times13\times14\times16$，这部分不含有$p_i$，分块计算，按照$p_i^{a_i}$分块，显然每块得到的答案相同，因此可以使用$O(p_i^{a_i})$的时间暴力计算，然后使用快速幂得到所有答案。
3. 第二部分，$3^5$，暂时无视，等会儿再说。

• 当我们完成上述过程后，我们就求出了$x!$中所有除开$p_i$的数的乘积，记为$f(x)$，这个乘积与$p_i^{a_i}$互素，故一定存在逆元，可以使用$\frac{f(n)}{f(m)f(n-m)}$计算出$C_n^m\bmod p_i^{a_i}$的一部分。
• 对于$C_n^m\bmod p_i^{a_i}$的另一部分，也就是我们在计算$x!$中被无视的所有$p_i$，记其个数为$g(x)$，显然我们可以在$O(\log x)$的时间内得到$g(x)$，那么$p^{g(n)-g(m)-g(n-m)}$就是剩下的部分。

将上述两部分乘起来即可得到$C_n^m\bmod p_i^{a_i}$。

void Exgcd(LL a,LL b,LL &gcd,LL &x,LL &y) {
	if(!b)gcd=a,x=1,y=0;
	else Exgcd(b,a%b,gcd,y,x),y-=x*(a/b);
}
 
LL inv(LL a,LL mod) {
	LL gcd,x,y;
	Exgcd(a,mod,gcd,x,y);
	return (x+mod)%mod;
}
 
LL Quick_Pow(LL a,LL b,LL mod) {
	LL ans=1;
	for(; b; b>>=1,a=a*a%mod)if(b&1)ans=ans*a%mod;
	return ans;
}
 
LL Fac(LL n,LL p,LL num) {
	if(n==0)return 1;
	LL sum=1;
	for(LL i=2; i<=num; i++)if(i%p)sum=sum*i%num;
	sum=Quick_Pow(sum,n/num,num);
	for(int i=2; i<=n%num; i++)if(i%p)sum=sum*i%num;
	return sum*Fac(n/p,p,num)%num;
}
 
LL Cal(LL n,LL p) {
	LL sum=0;
	for(LL i=n; i; i/=p)sum+=i/p;
	return sum;
}
 
LL C(LL n,LL m,LL p,LL num) {
	if(n<m)return 0;
	LL x=Fac(n,p,num),y=Fac(m,p,num),z=Fac(n-m,p,num);
	LL sum=Cal(n,p)-Cal(m,p)-Cal(n-m,p);
	return x*inv(y,num)%num*inv(z,num)%num*Quick_Pow(p,sum,num)%num;
}
 
LL ExLucas(LL n,LL m,LL mod) {
	LL x=mod,ans=0;
	for(LL i=2; i<=sqrt(mod); i++)
		if(x%i==0) {
			LL num=1;
			while(x%i==0)x/=i,num*=i;
			LL Mi=mod/num;
			ans=(ans+Mi*inv(Mi,num)%mod*C(n,m,i,num)%mod)%mod;
		}
	if(x>1)ans=(ans+mod/x*inv(mod/x,x)%mod*C(n,m,x,x)%mod)%mod;
	return ans;
}

时间复杂度

$$T(n)=T(\frac np)+O(p^a)$$

根据主定理可以得到$T(n)=O(p^a)$。