扩展卢卡斯定理学习笔记

扩展卢卡斯定理学习笔记

【模板】扩展卢卡斯

前置知识

1. (exgcd)
2. 中国剩余定理
3. 一定数学能力

卢卡斯定理是用来求(C_n^m mod{p})的工具,若(p)为质数,则(C_n^mequiv C_{lfloor{frac{n}{p}} floor}^{lfloor{frac{m}{p}} floor}*C_{n\%p}^{m\%p}(mod p))

但是当(p)不是质数的时候,就需要用到扩展卢卡斯定理了.下面简单说一下过程:

1

将模数(P)分解成(prod p_i^{k_i})的形式,也就是将(P)表示成若干个质数的指数的乘积形式.

2

设答案为(ans),则有$$ansequiv C^{m}_{n}(mod p_i ^{k_i})$$,也就是说,我们只需要求解(C_n^m)在模某个质数的几次幂下的结果,然后将这些结果用中国剩余定理合并就可以求出答案.

3

下面用(p)表示(p_i),(k)表示(k_i).

[C_n^m=frac{n!}{m!*(n-m)!} ]

将(n!,m!,(n-m)!)中含有(p_i)的因子提出来,将上面的式子进行变形,则有:

[C_n^mequiv frac{frac{n!}{p^{a1}}}{frac{m!}{p^{a2}}*frac{(n-m)!}{p^{a3}}}*p^{a1-a2-a3}(mod p^k) ]

那么此时式子的左半边都是与(p^k)互质的,可以直接用(exgcd)求逆元,右边可以快速幂解决.

4

现在我们需要解决$$frac{n!}{p^{a1}}$$,根据别人得出的式子,我们可以知道:

[n!equiv p^{lfloor frac{n}{p} floor}*lfloor frac{n}{p} floor!*(prod_{i=1}^{p^k}num_i)^{lfloor frac{n}{p^k} floor}*(prod_{i=1}^{n\%p^k}num_i)(mod p^k) ]

其中(num_i)表示不含(p)因子的数字.
其中((prod_{i=1}^{p^k}num_i)^{lfloor frac{n}{p^k} floor})可以直接枚举(1)到(p^k)的数字乘起来,然后次方做一遍快速幂.
((prod_{i=1}^{n\%p^k}num_i))直接枚举乘起来,(lfloor frac{n}{p} floor!)递归求解,边界是当(n=0)时,返回(1).
因为我们要求(frac{n!}{p^{a1}}),所以在乘的时候可以直接把式子第一部分的(p^{lfloor frac{n}{p} floor})直接忽略掉.

int fac(int n, int pi, int pk){
	if(!n) return 1;
	int res = 1;
	for(int i = 2; i < pk; i++)
		if(i % pi) (res *= i) %= pk;
	res = qpow(res, n/pk, pk);
	for(int i = 2; i <= n%pk; i++)
		if(i % pi) (res *= i) %= pk;
	return fac(n/pi, pi, pk)*res%pk;
}

5

根据$$C_n^mequiv frac{frac{n!}{p^{{a1}}}{frac{m!}{p}{a2}}*frac{(n-m)!}{p^{a3}}}*p{a1-a2-a3}(mod p^k)$$,计算出了左边式子的上下部分,就可以直接用(exgcd)求出下面两个结果在模(p^k)意义下的逆元,然后乘起来.

inline int C(int n, int m, int pi, int pk){
	if(n < m) return 0;
	int r1 = fac(n, pi, pk), r2 = fac(m, pi, pk);
	int r3 = fac(n-m, pi, pk), cnt = 0;
	for(int i = n; i; i /= pi) cnt += i/pi;
	for(int i = m; i; i /= pi) cnt -= i/pi;
	for(int i = n-m; i; i /= pi) cnt -= i/pi;
	return r1*inv(r2, pk)%pk*inv(r3, pk)%pk*qpow(pi, cnt, pk)%pk;
}

6

再根据$$ansequiv C^{m}_{n}(mod p ^{k})$$,我们可以用中国剩余定理合并这些结果(这里我用的是扩展中国剩余定理)

inline int exlucas(int n, int m, int p){
	int pi, pk, res = 0, x, y, gcd, c, M = 1;
	for(int i = 2; i*i <= p; i++){
		if(p % i == 0){
			pi = i, pk = 1;
			while(p % i == 0) p /= i, pk *= i;
			gcd = exgcd(M, pk, x, y), c = Mod((C(n, m, pi, pk)-res)%pk+pk, pk);
			x = Mod(x*(c/gcd)%pk+pk, pk);
			res += M*x, M *= pk/gcd, res %= M;
		}
	}
	if(p > 1){ // 最后分解质因数后可能存在一个大于sqrt(p)的大质数
		pi = pk = p; gcd = exgcd(M, pk, x, y);
		c = Mod((C(n, m, pi, pk)-res)%pk+pk, pk);
		x = Mod(x*(c/gcd)%pk+pk, pk);
		res += M*x, M *= pk/gcd, res %= M;
	}
	return res;
}

下面放一下完整代码

#include<bits/stdc++.h>
using namespace std;
typedef int _int;
#define int long long

int n, m, p;

int exgcd(int a, int b, int &x, int &y){
	if(!b){ x = 1, y = 0; return a; }
	int gcd = exgcd(b, a%b, x, y), tmp;
	tmp = x, x = y, y = tmp-a/b*y;
	return gcd;
}

inline int Mod(int x, int mod){ return x > mod ? x-mod : x; }

inline int inv(int a, int mod){
	int x, y; exgcd(a, mod, x, y);
	return (x%mod+mod)%mod;
}

inline int qpow(int x, int n, int mod){
	int res = 1;
	for(; n; x = x*x%mod, n >>= 1)
		if(n & 1) (res *= x) %= mod;
	return res;
}

int fac(int n, int pi, int pk){
	if(!n) return 1;
	int res = 1;
	for(int i = 2; i < pk; i++)
		if(i % pi) (res *= i) %= pk;
	res = qpow(res, n/pk, pk);
	for(int i = 2; i <= n%pk; i++)
		if(i % pi) (res *= i) %= pk;
	return fac(n/pi, pi, pk)*res%pk;
}

inline int C(int n, int m, int pi, int pk){
	if(n < m) return 0;
	int r1 = fac(n, pi, pk), r2 = fac(m, pi, pk);
	int r3 = fac(n-m, pi, pk), cnt = 0;
	for(int i = n; i; i /= pi) cnt += i/pi;
	for(int i = m; i; i /= pi) cnt -= i/pi;
	for(int i = n-m; i; i /= pi) cnt -= i/pi;
	return r1*inv(r2, pk)%pk*inv(r3, pk)%pk*qpow(pi, cnt, pk)%pk;
}

inline int exlucas(int n, int m, int p){
	int pi, pk, res = 0, x, y, gcd, c, M = 1;
	for(int i = 2; i*i <= p; i++){
		if(p % i == 0){
			pi = i, pk = 1;
			while(p % i == 0) p /= i, pk *= i;
			gcd = exgcd(M, pk, x, y), c = Mod((C(n, m, pi, pk)-res)%pk+pk, pk);
			x = Mod(x*(c/gcd)%pk+pk, pk);
			res += M*x, M *= pk/gcd, res %= M;
		}
	}
	if(p > 1){
		pi = pk = p; gcd = exgcd(M, pk, x, y);
		c = Mod((C(n, m, pi, pk)-res)%pk+pk, pk);
		x = Mod(x*(c/gcd)%pk+pk, pk);
		res += M*x, M *= pk/gcd, res %= M;
	}
	return res;
}

_int main(){
	cin >> n >> m >> p;
	cout << exlucas(n, m, p) << endl;
	return 0;
}