卢卡斯定理学习笔记

卢卡斯定理

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

P P P为素数

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

代码 :

#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 不为素数

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

因为 P P P不为素数, 所以 :
P = ∏ P i a i P = \prod P_i^{a_i} P=Piai
我们考虑求出
( n m ) ( m o d    P i a i ) \binom {n}{m} (mod \ \ P_i^{a_i}) (mn)(mod  Piai)
然后直接用 C R T CRT CRT合并就是答案了

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

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

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

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

f a c ( n ) = ( f a c ( ⌊ n P i ⌋ ) ⋅ P i ⌊ n P i ⌋ ) ⋅ ( ∏ i = 1 P k i ) ⌊ n P k ⌋ ⋅ ( ∏ i = 1 n   m o d   P k i ) ( m o d P i a i ) fac(n) = \left(fac(\lfloor \frac{n}{P_i} \rfloor) \cdot P_i^{\lfloor \frac{n}{P_i} \rfloor}\right) \cdot \left(\prod_{i = 1}^{P_k} i\right)^{\lfloor\frac{n}{P_k}\rfloor} \cdot \left(\prod_{i = 1}^{n ~mod~ P_k} i\right) \pmod {P_i^{a_i}} fac(n)=(fac(Pin)PiPin)(i=1Pki)Pkn(i=1n mod Pki)(modPiai)
我们发现这个是一个递归式, 那么直接递归就好了, 特别的如果 n = 0 n = 0 n=0就返回 1 1 1

需要注意的是 P i P_i Pi这个东西不能在递归的时候直接乘上, 我们考虑, 在递归完之后, 或者递归之前直接处理出来就好了

代码 :

#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;
}
  • 0
    点赞
  • 2
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值