扩展卢卡斯定理学习笔记
扩展卢卡斯定理学习笔记
前置知识
- \(exgcd\)
- 中国剩余定理
一定数学能力
卢卡斯定理是用来求\(C_n^m \bmod{p}\)的工具,若\(p\)为质数,则\(C_n^m\equiv C_{\lfloor{\frac{n}{p}}\rfloor}^{\lfloor{\frac{m}{p}}\rfloor}*C_{n\%p}^{m\%p}(\bmod \ p)\)
但是当\(p\)不是质数的时候,就需要用到扩展卢卡斯定理了.下面简单说一下过程:
1
将模数\(P\)分解成\(\prod p_i^{k_i}\)的形式,也就是将\(P\)表示成若干个质数的指数的乘积形式.
2
设答案为\(ans\),则有$$ans\equiv C^{m}_{n}(\mod \ p_i ^{k_i})$$,也就是说,我们只需要求解\(C_n^m\)在模某个质数的几次幂下的结果,然后将这些结果用中国剩余定理合并就可以求出答案.
3
下面用\(p\)表示\(p_i\),\(k\)表示\(k_i\).
将\(n!,m!,(n-m)!\)中含有\(p_i\)的因子提出来,将上面的式子进行变形,则有:
那么此时式子的左半边都是与\(p^k\)互质的,可以直接用\(exgcd\)求逆元,右边可以快速幂解决.
4
现在我们需要解决$$\frac{n!}{p^{a1}}$$,根据别人得出的式子,我们可以知道:
其中\(num_i\)表示不含\(p\)因子的数字.
其中\((\prod_{i=1}^{p^k}num_i)^{\lfloor \frac{n}{p^k} \rfloor}\)可以直接枚举\(1\)到\(p^k\)的数字乘起来,然后次方做一遍快速幂.
\((\prod_{i=1}^{n\%p^k}num_i)\)直接枚举乘起来,\(\lfloor \frac{n}{p} \rfloor!\)递归求解,边界是当\(n=0\)时,返回\(1\).
因为我们要求\(\frac{n!}{p^{a1}}\),所以在乘的时候可以直接把式子第一部分的\(p^{\lfloor \frac{n}{p} \rfloor}\)直接忽略掉.
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^m\equiv \frac{\frac{n!}{p{a1}}}{\frac{m!}{p{a2}}*\frac{(n-m)!}{p{a3}}}*p(\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
再根据$$ans\equiv 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;
}
