卢卡斯定理学习笔记

卢卡斯定理

常用于求解\(\binom{n}{m} (mod \ \ P)\) , 且\(n\), \(m\)值域特别大, 而 \(P\) 的值域比较小

\(P\)为素数

直接记定理就好了 : \[ \binom{n}{m} (mod \ \ P)= \binom{n / P}{m / P} \cdot \binom {n \% P}{m \% P} (mod\ \ P) \] 然后我们预处理出\(1 \sim P\)的阶乘和逆元, 那么求这个组合数的复杂度就是\(O(P + log_P^{n})\)

代码 :

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
#include <bits/stdc++.h>

using namespace std;

const int maxn = 1e6 + 5;

int n, m, fac[maxn], P;

inline int read(){
char ch = getchar(); int u = 0, f = 1;
while(!isdigit(ch)){if(ch == '-')f = -1; ch = getchar();}
while(isdigit(ch)){u = u * 10 + ch - 48; ch = getchar();} return u * f;
}

inline int ksm(int x, int k){
int cnt = 1;
while(k){
if(k & 1)cnt = 1ll * cnt * x % P;
x = 1ll * x * x % P; k >>= 1;
} return cnt;
}

inline int C(int n, int m, int P){
if(n < m)return 0;
return (int)(1ll * fac[n] * ksm(fac[n - m], P - 2) % P * ksm(fac[m], P - 2) % P);
}

inline Lucas(int n, int m, int P){
if(n == 0 || m == 0)return 1;
return (int)(1ll * Lucas(n / P, m / P, P) * C(n % P, m % P, P) % P);
}

int main(){
n = read(); m = read(); P = read();
fac[0] = 1;
for(register int i = 1; i <= P; i++)
fac[i] = 1ll * fac[i - 1] * i % P;
cout << Lucas(n, m, P);
return 0;
}

\(P\) 不为素数

其实这个东西和卢卡斯定理并没有什么关系

因为\(P\)不为素数, 所以 : \[ P = \prod P_i^{a_i} \] 我们考虑求出 \[ \binom {n}{m} (mod \ \ P_i^{a_i}) \] 然后直接用\(CRT\)合并就是答案了

那么现在问题转化为了如何快速求出\(\binom{n}{m}(mod\ \ P_i^{a_i})\)

因为实在膜 \(P_i^{a_i}\) 的意义下, 所以阶乘其实是有循环节的, 我们可以考虑直接暴力处理循环节和不完整的循环

难道这样就好了吗, 事实上并没有这么简单

因为\(P_i^{a_i}\)\(P_i\)的倍数, 所以很有可能再膜意义下就直接变为零了, 但是事实上并不是这样的, 我们观察组合数的公式 : \[ \binom {n}{m} = \frac{n!}{(n - m)!m!} \] 所以分母和分子的\(P_i\)是可已被约掉的, 所以我们不能直接取膜, 我们再考虑把是\(P_i\)的倍数的数给单独提出来搞一下

那么 \[ fac(n) = \big(fac(\lfloor \frac{n}{P_i^{a_i}} \rfloor) \cdot P_i^{\lfloor \frac{n}{P_i^{a_i}} \rfloor}\big) \cdot \big(\prod_{i = 1}^{P_i^{a_i}} i\big) \cdot \big(\prod_{i = 1}^{n \% P_i^{a_i}} i\big) (mod \ \ P_i^{a_i}) \] 我们发现这个是一个递归式, 那么直接递归就好了, 特别的如果\(n = 0\)就返回\(1\)

需要注意的是\(P_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
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
#include <bits/stdc++.h>

typedef long long LL;

using namespace std;

LL n, m, P;

inline LL read(){
char ch = getchar(); LL u = 0, f = 1;
while(!isdigit(ch)){if(ch == '-')f = -1; ch = getchar();}
while(isdigit(ch)){u = u * 10 + ch - 48; ch = getchar();} return u * f;
}

inline LL ksm(LL x, LL k, LL P){
LL cnt = 1;
while(k){
if(k & 1)cnt = cnt * x % P;
x = x * x % P; k >>= 1;
} return cnt;
}

inline LL Exgcd(LL a, LL b, LL &x, LL &y){
if(b == 0){x = 1; y = 0; return a;}
LL x0 = 0, y0 = 0;
LL d = Exgcd(b, a % b, y0, x0);
x = x0; y = y0 - ((a / b) * x0);
return d;
}

inline LL fac(LL n, LL Pi, LL Pk){
if(n == 0)return 1;
LL rnt = 1;
for(register int i = 2; i <= Pk; i++)
if(i % Pi)rnt = rnt * i % Pk;
rnt = ksm(rnt, n / Pk, Pk);
for(register int i = 2; i <= n % Pk; i++)
if(i % Pi)rnt = rnt * i % Pk;
return rnt * fac(n / Pi, Pi, Pk) % Pk;
}

inline LL inv(LL x, LL P){
LL inv, tmp;
Exgcd(x, P, inv, tmp);
return (inv + P) % P;
}

inline LL ExLucas(LL n, LL m, LL Pi, LL Pk){
LL C1 = fac(n, Pi, Pk), C2 = fac(n - m, Pi, Pk), C3 = fac(m, Pi, Pk);
LL tim = 0;
for(register LL i = n; i; i /= Pi)tim += i / Pi;
for(register LL i = n - m; i; i /= Pi)tim -= i / Pi;
for(register LL i = m; i; i /= Pi)tim -= i / Pi;
return C1 * inv(C2, Pk) % Pk * inv(C3, Pk) % Pk * ksm(Pi, tim, Pk) % Pk;
}

inline LL China(LL n, LL m, LL P){
LL A = 0, M = P, Pk = 1;
for(register LL i = 2; i <= P; i++)
if(P % i == 0){
Pk = 1;
while(P % i == 0)Pk *= i, P /= i;
A = (A + ExLucas(n, m, i, Pk) * inv(M / Pk, Pk) % M * (M / Pk) % M) % M;
}
return (A + M) % M;
}

int main(){
n = read(); m = read(); P = read();
cout << China(n, m, P);
return 0;
}