【知識總結】擴展盧卡斯定理（exLucas）

擴展盧卡斯定理用於求如下式子（其中\(p\)不一定是質數）：

\[C_n^m\ mod\ p \]

我們將這個問題由總體到局部地分為三個層次解決。

層次一：原問題

首先對\(p\)進行質因數分解：

\[p=\prod_i p_i^{k_i} \]

顯然\(p_i^{k_i}\)是兩兩互質的，所以如果分別求出\(C_n^m\ mod\ p_i^{k_i}\)，就可以構造出若干個形如\(C_n^m=a_i\ mod\ p_i^{k_i}\)的方程，然后用中國剩余定理即可求解。

層次二：組合數模質數冪

現在的問題就轉化成了求如下式子（其中\(p\)是質數）：

\[C_n^m\ mod\ p^k \]

腦補一下組合數公式\(C_n^m=\frac{n!}{m!\times (n-m)!}\)，發現由於\(m!\)和\((n-m)!\)可能包含質因子\(p\)，所以不能直接求他們對於\(p^k\)的逆元。此時我們可以將\(n!\)、\(m!\)、\((n-m)!\)中的質因子\(p\)全部提出來，最后再乘回去即可。即變為下式（\(k1\)為\(n!\)中質因子\(p\)的次數，\(k2\)、\(k3\)同理）：

\[\frac{\frac{n!} {p^{k1}}}{\frac{m!}{p^{k2}}\times \frac{(n-m)!}{p^{k3}}}\times p^{k1-k2-k3} \]

\(\frac{m!}{p^{k2}}\)和\(\frac{(n-m)!}{p^{k3}}\)和\(p^k\)是互質的，可以直接求逆元。

層次三：階乘除去質因子后模質數冪

現在看看如何計算形如下式的式子。

\[\frac{n!}{p^{a}}\ mod\ p^k \]

先考慮如何計算\(n!\ mod\ p^k\)
舉個例子：\(n=22\)，\(p=3\)，\(k=2\)
把這個寫出來：
\(22!=1\times 2\times 3\times 4\times 5\times 6\times 7\times 8\times 9\times 10 \times 11\times 12\times 13\times 14\times 15\times 16\times 17\times 18\times 19 \times 20 \times 21 \times 22\)

把其中所有\(p\)（也就是\(3\)）的倍數提取出來，得到：

\(22!=3^7 \times (1\times 2\times 3\times 4\times 5\times 6\times 7)\times(1\times 2\times 4\times 5\times 7\times 8\times 10 \times 11\times 13\times 14\times 16\times 17\times 19 \times 20 \times 22 )\)

可以看出上式分為三個部分：第一個部分是\(3\)的冪，次數是小於等於\(22\)的\(3\)的倍數的個數，即\(\lfloor\frac{n}{p}\rfloor\)

第二個部分是一個階乘\(7!\)，即\(\lfloor\frac{n}{p}\rfloor!\)，可以遞歸解決

第三個部分是\(n!\)中與\(p\)互質的部分的乘積，這一部分具有如下性質：
\(1\times 2\times 4\times 5\times 7\times 8\equiv10 \times 11\times 13\times 14\times 16\times 17\ mod\ p^k\)

在模\(3^2\)的意義下\(10\)和\(1\)同余，\(11\)和\(2\)同余……寫成下式就比較顯然
（\(t\)是任意正整數）

\[\prod_{i,(i,p)=1}^{p^k}i\equiv\prod_{i,(i,p)=1}^{p^k}(i+tp^k)\ mod\ p^k \]

\(\prod_{i,(i,p)=1}^{p^k}i\)一共循環了\(\lfloor\frac{n}{p^k}\rfloor\)次，暴力求出\(\prod_{i,(i,p)=1}^{p^k}i\)然后用快速冪求它的\(\lfloor\frac{n}{p^k}\rfloor\)次冪。

最后還要乘上\(19\times 20 \times 22\)（即\(\prod_{i,(i,p)=1}^{n\ mod\ p^k}i\)），顯然這一段的長度一定小於\(p^k\)，暴力乘上去即可。

如上三部分的乘積就是\(n!\)。最終要求的是\(\frac{n!}{p^{a}}\ mod\ p^k\)，分母全部由上述第一部分和第二部分貢獻（第三部分和\(p\)互質）。而遞歸計算第二部分的時候已經除去了第二部分中的因子\(p\)，所以最終的答案就是上述第二部分遞歸返回的結果和第三部分的乘積（與第一部分無關）。

結合代碼方便理解：

ll fac(const ll n, const ll p, const ll pk)
{
	if (!n)
		return 1;
	ll ans = 1;
	for (int i = 1; i < pk; i++)
		if (i % p)
			ans = ans * i % pk;
	ans = power(ans, n / pk, pk);
	for (int i = 1; i <= n % pk; i++)
		if (i % p)
			ans = ans * i % pk;
	return ans * fac(n / p, p, pk) % pk;
}

層次二：組合數模質數冪

回到這個式子

\[\frac{\frac{n!} {p^{k1}}}{\frac{m!}{p^{k2}}\times \frac{(n-m)!}{p^{k3}}}\times p^{k1-k2-k3} \]

可以很容易地把它轉換成代碼（注意i要開long long）：

ll C(const ll n, const ll m, const ll p, const ll pk)
{
	if (n < m)
		return 0;
	ll f1 = fac(n, p, pk), f2 = fac(m, p, pk), f3 = fac(n - m, p, pk), 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 f1 * inv(f2, pk) % pk * inv(f3, pk) % pk * power(p, cnt, pk) % pk;
}

層次一：原問題

完整代碼（題目：洛谷4720【模板】擴展盧卡斯）：

#include <cstdio>
#include <algorithm>
#include <cstring>
#include <iostream>
#include <climits>
#include <cmath>
using namespace std;
namespace zyt
{
	const int N = 1e6;
	typedef long long ll;
	ll n, m, p;
	inline ll power(ll a, ll b, const ll p = LLONG_MAX)
	{
		ll ans = 1;
		while (b)
		{
			if (b & 1)
				ans = ans * a % p;
			a = a * a % p;
			b >>= 1;
		}
		return ans;
	}
	ll fac(const ll n, const ll p, const ll pk)
	{
		if (!n)
			return 1;
		ll ans = 1;
		for (int i = 1; i < pk; i++)
			if (i % p)
				ans = ans * i % pk;
		ans = power(ans, n / pk, pk);
		for (int i = 1; i <= n % pk; i++)
			if (i % p)
				ans = ans * i % pk;
		return ans * fac(n / p, p, pk) % pk;
	}
	ll exgcd(const ll a, const ll b, ll &x, ll &y)
	{
		if (!b)
		{
			x = 1, y = 0;
			return a;
		}
		ll xx, yy, g = exgcd(b, a % b, xx, yy);
		x = yy;
		y = xx - a / b * yy;
		return g;
	}
	ll inv(const ll a, const ll p)
	{
		ll x, y;
		exgcd(a, p, x, y);
		return (x % p + p) % p;
	}
	ll C(const ll n, const ll m, const ll p, const ll pk)
	{
		if (n < m)
			return 0;
		ll f1 = fac(n, p, pk), f2 = fac(m, p, pk), f3 = fac(n - m, p, pk), 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 f1 * inv(f2, pk) % pk * inv(f3, pk) % pk * power(p, cnt, pk) % pk;
	}
	ll a[N], c[N];
	int cnt;
	inline ll CRT()
	{
		ll M = 1, ans = 0;
		for (int i = 0; i < cnt; i++)
			M *= c[i];
		for (int i = 0; i < cnt; i++)
			ans = (ans + a[i] * (M / c[i]) % M * inv(M / c[i], c[i]) % M) % M;
		return ans;
	}
	ll exlucas(const ll n, const ll m, ll p)
	{
		ll tmp = sqrt(p);
		for (int i = 2; p > 1 && i <= tmp; i++)
		{
			ll tmp = 1;
			while (p % i == 0)
				p /= i, tmp *= i;
			if (tmp > 1)
				a[cnt] = C(n, m, i, tmp), c[cnt++] = tmp;
		}
		if (p > 1)
			a[cnt] = C(n, m, p, p), c[cnt++] = p;
		return CRT();
	}
	int work()
	{
		ios::sync_with_stdio(false);
		cin >> n >> m >> p;
		cout << exlucas(n, m, p);
		return 0;
	}
}
int main()
{
	return zyt::work();
}