扩展卢卡斯定理 学习笔记

个人感觉这玩意属于是省选数论中的boss了，很恶心
考虑求模\(P\)意义下的组合数\(\binom{n}{m}\).
卢卡斯定理可以十分简单地解决\(P\)为素数时的情况
如果\(P\)不为质数呢?扩展卢卡斯定理出现了.
套路地，令\(P=\prod_{i=1}^{r} {p_i}^{\alpha_i}\)
那么我们只用求出对于每个\(p_i\)地答案，然后使用\(\text{CRT}\)合并就行了
现在问题转化为求$$\binom{n}{m}=\frac{n!}{m!(n-m)!} \mod {p}^k$$
由于我们没有办法取逆元，所以把含\(p\)的项全部提出来,求剩下的部分
怎么求提出来的\(p\)的次数呢?递归即可，当然，可以写成一个\(\text{for}\),这部分复杂度\(O(\log_{p}{n})\)
剩下的部分就没有\(p\)了，可以直接求逆元
那么我们现在要求的就是\(n!\)去掉\(p\)的这个东西
考虑$$n!=(n/p)！* p^{n/p}*S$$
其中 $$S=\prod_{i=1}^{n} ([i \perp p] i)$$
容易发现 \(S\)中的因子是会以\(p^k\)为循环节循环的.
所以可以\(O(p^k)\)暴力
当然，这样子的打法十分简陋，复杂度是单次询问\(O(p\log p)\)
看上去很弱，但是实际应用中若模数一定通常可以预处理出和模数有关的一些式子，可优化为预处理\(O(p)\),单次询问\(O(\log p)\)
还有个优化，就是当你提出来的\(p\)的次数比模数中\(p\)的次数还要大时，可以直接返回0
这个优化效果十分明显

#include <cstdio>
#include <iostream>
#define LL long long
using namespace std;
inline LL read() {
	LL res = 0, flag = 0; char ch = getchar();
	for(; !isdigit(ch); ch = getchar()) if(ch == '-') flag = 1;
	for(; isdigit(ch); ch = getchar()) res = (res << 1) + (res << 3) + (ch ^ 48);
	if(flag) res = ~res + 1;
	return res;
}
inline LL fpow(LL base, LL pnt, LL P) {
	LL res = 1; 
	for(; pnt; pnt >>= 1, base = base * base % P) 
		if(pnt & 1) res = res * base % P;
	return res;
}
inline void exgcd(LL a, LL b, LL &x, LL &y) {
	if(!b) x = 1, y = 0;
	else exgcd(b, a % b, y, x), y -= a / b * x;
}
inline LL inv(LL x, LL P) {
	LL y, z;
	exgcd(x, P, y, z);
	y %= P; if(y < 0) y += P;
	return y;
}
inline LL CRT (LL cnt, LL *p, LL *a, LL P, LL *c) {
	LL res = 0, M;
	for(int i = 1; i <= cnt; ++i) 
		M = P / p[i], res += M * inv(M, p[i]) % P * a[i], res %= P;
	return res;
}
LL calc(LL n, LL p, LL P) {
	if(!n) return 1;
	LL s = 1;
	for(LL i = 1; i <= P; ++i) if(i % p) s = s * i % P;
	s = fpow(s, n / P, P);
	for(LL i = n / P * P + 1; i <= n; ++i) if(i % p) s = i % P * s % P;
	return s * calc(n / p, p, P) % P;
}
LL multiLucas(LL n, LL m, LL p, LL P) {
	LL cnt = 0;
	for(LL i = n; i; i /= p) cnt += i / p;
	for(LL i = m; i; i /= p) cnt -= i / p;
	for(LL i = n - m; i; i /= p) cnt -= i / p; 
	return fpow(p, cnt, P) * calc(n, p, P) % P * inv(calc(m, p, P) * calc(n - m, p, P) % P, P) % P;
}
LL exLucas(LL n, LL m, LL P) {
	LL cnt = 0, MOD = P;
	LL p[200], a[200], c[200];
	for(LL i = 2; i * i <= P; ++i) {
		if(!(P % i)) {
			p[++cnt] = 1, c[cnt] = i;
			while(!(P % i)) P /= i, p[cnt] *= i;
			a[cnt] = multiLucas(n, m, i, p[cnt]);
		}
	}
	if(P > 1) p[++cnt] = P, a[cnt] = multiLucas(n, m, P, P);
	return CRT(cnt, p, a, MOD, c);
}
LL n, m, p;
int main() {
	n = read(), m = read(), p = read(),
	printf("%lld\n",exLucas(n, m, p));
	return 0;
}