隐藏
Lucas定理学习笔记 | Bill Yang's Blog

路终会有尽头,但视野总能看到更远的地方。

0%

Lucas定理学习笔记

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

Lucas定理

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

Lucas定理

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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$个同余方程:

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

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

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

  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}$。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
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)=O(p^a)$。

姥爷们赏瓶冰阔落吧~