【luogu P4720】【模板】扩展卢卡斯定理/exLucas（数论）（CRT）

【模板】扩展卢卡斯定理/exLucas

题目链接：luogu P4720

题目大意

要你求一个组合数在模一个数的结果。
模数不一定是质数。

思路

我们考虑质因数分解模数。

考虑分成了 (p=p_1^{a_1},p_2^{a_2},...)

然后你考虑求出：
(egin{cases}C_{n}^mmod p_1^{a_1}\ C_{n}^mmod p_2^{a_2}\ ……end{cases})

不难看出求出来的结果我们可以直接用 CRT（中国剩余定理）合并。

那接着就变成如何求 (C_n^mmod p^{k})。
首先我们用普通的式子：(dfrac{n!}{m!(n-m)!}mod p^k)

显然是不能直接搞，因为我们不能求出下面的两个的逆元。
然后又逆元的条件是这个数和模数互质，所以我们可以换一下式子：
(dfrac{frac{n!}{p^x}}{frac{m!}{p^y}frac{(n-m)!}{p^z}}p^{x-y-z}mod p^k)

然后 (x) 是 (n!) 中包含 (p) 因子的个数。（(y,z) 同理）

那问题就是求 (dfrac{n!}{p^x}mod p^k)
考虑处理一下 (n!)：
(=(pcdot 2pcdot3pcdot...)(1cdot 2cdot 3cdot...))
（左边是 (1sim n) 中 (p) 的倍数，右边就是不是倍数的）

然后提出左边的 (p)：
(=p^{leftlfloorfrac{n}{p} ight floor}(leftlfloordfrac{n}{p} ight floor)!prodlimits_{i=1,i otequiv0pmod p}^ni)

然后后面的显然是有循环节：
(=p^{leftlfloorfrac{n}{p} ight floor}(leftlfloordfrac{n}{p} ight floor)!(prodlimits_{i=1,i otequiv0pmod p}^{p^k}i)^{leftlfloorfrac{n}{p^k} ight floor}(prodlimits_{i=p^kleftlfloorfrac{n}{p^k} ight floor,i otequiv0pmod p}^{n}i))

两个分别代表循环节和余项。

然后因为取模，(p^{leftlfloorfrac{n}{p} ight floor}) 会没掉，但是 ((leftlfloordfrac{n}{p} ight floor)!) 中可能还有 (p)。
那就是一个递归下去的过程，设 (F(n)=dfrac{n!}{p^x})

然后递推就是：
(F(n)=F(leftlfloordfrac{n}{p} ight floor)(prodlimits_{i=1,i otequiv0pmod p}^{p^k}i)^{leftlfloorfrac{n}{p^k} ight floor}(prodlimits_{i=p^kleftlfloorfrac{n}{p^k} ight floor,i otequiv0pmod p}^{n}i))
（这个的复杂度是 (log_pn)，边界是 (f(0)=1)）

然后在原来的式子，就变成了 (dfrac{F(n)}{F(m)F(n-m)}p^{x-y-z}mod p^k)
这下分母的两个就是一定跟 (p^k) 互质，就可以用扩欧求啦。

然后你会发现还有一个 (p^{x-y-z})。
那我们就是要求 (F(n)=dfrac{n!}{p^x}) 中的 (x)。

想想上面的式子，我们要的就是被除掉的部分，就是 (p^{leftlfloorfrac{n}{p} ight floor}) 中的 (leftlfloordfrac{n}{p} ight floor)。
每个 (leftlfloordfrac{n}{p} ight floor) 加起来就是结果，所以：
(g(n)=leftlfloordfrac{n}{p} ight floor+g(leftlfloordfrac{n}{p} ight floor))
（复杂度也是 (log_pn)，边界是 (g(n)=0(n<p))）

然后再带入到原来的式子就是：
(dfrac{F(n)}{F(m)F(n-m)}p^{G(n)-G(m)-G(n-m)}mod p^k)

然后就好啦！！！

代码

#include<cstdio>
#define ll long long

using namespace std;

ll n, m, p;
ll a[1001], b[1001];

ll ksm(ll x, ll y, ll p) {
	ll re = 1;
	while (y) {
		if (y & 1) re = re * x % p;
		x = x * x % p;
		y >>= 1;
	} 
	return re;
}

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

ll inv(ll a, ll p) {//别用快速幂求，要用扩欧
	ll x, y;
	exgcd(a, p, x, y);
	return (x % p + p) % p;
}

ll F(ll n, ll p, ll pk) {
	if (!n) return 1;
	ll rou = 1, el = 1;
	for (ll i = 1; i <= pk; i++) {//循环节部分
		if (i % p) rou = rou * i % pk;
	}
	rou = ksm(rou, n / pk, pk);
	for (ll i = pk * (n / pk); i <= n; i++) {//余项部分
		if (i % p) el = el * (i % pk) % pk;
	}
	return F(n / p, p, pk) * rou % pk * el % pk;
}

ll G(ll n, ll p) {
	if (n < p) return 0;
	return n / p + G(n / p, p);
}

ll work_C(ll n, ll m, ll p, ll pk) {
	ll fz = F(n, p, pk), fm1 = inv(F(m, p, pk), pk), fm2 = inv(F(n - m, p, pk), pk);
	ll mi = ksm(p, G(n, p) - G(m, p) - G(n - m, p), pk);
	return fz * fm1 % pk * fm2 % pk * mi % pk;
}

ll exLucas(ll n, ll m, ll p) {
	ll tmpp = p, tot = 0, di;
	for (ll i = 2; i * i <= tmpp; i++) {
		if (tmpp % i == 0) {
			tot++; di = 1;
			while (tmpp % i == 0) {
				di *= i;
				tmpp /= i;
			}
			a[tot] = di;
			b[tot] = work_C(n, m, i, di);
		}
	}
	if (tmpp > 1) a[++tot] = tmpp, b[tot] = work_C(n, m, tmpp, tmpp);
	
	ll ans = 0;//CRT
	for (int i = 1; i <= tot; i++) {
		ll M = p / a[i], V = inv(M, a[i]);
		ans = (ans + b[i] * M % p * V % p) % p;
	}
	return ans;
}

int main() {
	scanf("%lld %lld %lld", &n, &m, &p);
	printf("%lld", exLucas(n, m, p));
	
	return 0;
}