【数学】扩展卢卡斯定理

Description

求
$$\left(\genfrac{}{}{0ex}{}{n}{m}\right)\mathrm{m}\mathrm{o}\mathrm{d}p$$
其中 $p$ 较小且 不保证 $p$ 是质数。

Method

前置芝士：

• 中国剩余定理

因为 $p$ 不为质数，所以得使用扩展卢卡斯定理（$\mathrm{E}\mathrm{x}\mathrm{t}\mathrm{e}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{d}\mathrm{L}\mathrm{u}\mathrm{c}\mathrm{a}\mathrm{s},\mathrm{e}\mathrm{x}\mathrm{L}\mathrm{u}\mathrm{c}\mathrm{a}\mathrm{s}$）。

Part 1：中国剩余定理

首先按照 唯一分解定理 将 $p$ 分解质因数：
$$p=\prod _{i=1}^k{q}_i^{{\alpha }_i}$$
如果我们能够求出
$$\{\begin{array}{l}\left(\genfrac{}{}{0ex}{}{n}{m}\right)\mathrm{m}\mathrm{o}\mathrm{d}{q}_1^{{\alpha }_1}\\ \left(\genfrac{}{}{0ex}{}{n}{m}\right)\mathrm{m}\mathrm{o}\mathrm{d}{q}_2^{{\alpha }_2}\\ \cdots \\ \left(\genfrac{}{}{0ex}{}{n}{m}\right)\mathrm{m}\mathrm{o}\mathrm{d}{q}_k^{{\alpha }_k}\end{array}$$
那么因为 $\forall i\ne j,\gcd (q_i^{{\alpha }_i},q_j^{{\alpha }_j})=1$，所以可以用 CRT 合并。

那么问题就变成了求出 $\left(\genfrac{}{}{0ex}{}{n}{m}\right)\mathrm{m}\mathrm{o}\mathrm{d}{q}_i^{{\alpha }_i}$。

Part 2：移除分子分母中的质因数

$$\begin{array}{rl}\left(\genfrac{}{}{0ex}{}{n}{m}\right)& \equiv \frac{n!}{m!(n-m)!}\\ & \equiv n!\cdot \mathrm{inv}(m!)\cdot \mathrm{inv}[(n-m)!](\mathrm{m}\mathrm{o}\mathrm{d}{q}^{\alpha })\end{array}$$

但 $m!,(n-m)!$ 与 $q^{\alpha }$ 不一定互质，所以不一定存在逆元。

我们将分子和分母中的 $q$ 的倍数全部提出来，设 $n!,m!,(n-m)!$ 分别含 $x,y,z$ 个质因数 $q$。

$$\begin{array}{rl}\frac{n!}{m!(n-m)!}& =\frac{\frac{n!}{q^x}}{\frac{m!}{q^y}\cdot \frac{(n-m)!}{q^z}}\cdot q^{x-y-z}\\ & \equiv \frac{n!}{q^x}\cdot \mathrm{inv}\left(\frac{m!}{q^y}\right)\cdot \mathrm{inv}\left[\frac{(n-m)!}{q^z}\right](\mathrm{m}\mathrm{o}\mathrm{d}{q}^{\alpha })\end{array}$$

这时 $\frac{m!}{q^y},\frac{(n-m)!}{q^z}$ 都与 $q^{\alpha }$ 互质。

Part 3：计算

以 $n!$ 为例：

$$\begin{array}{rl}n!& =(q\times 2q\times \cdots \times \lfloor \frac{n}{q}\rfloor q)\times \prod _{q\nmid i}^{n}i\\ & =q^{\lfloor \frac{n}{q}\rfloor }\times \lfloor \frac{n}{q}\rfloor !\times \prod _{q\nmid i}^{n}i\end{array}$$
其中 $q^{\lfloor \frac{n}{q}\rfloor }$ 可用快速幂直接算，$\lfloor \frac{n}{q}\rfloor !$ 可递归求解，后面 $\prod _{q\nmid i}^{n}i$ 就不好算了。

我们发现因为模数是 $q^{\alpha }$，又有 $x+q^{\alpha }\equiv x(\mathrm{m}\mathrm{o}\mathrm{d}{q}^{\alpha })$，所以可以将 $q^{\alpha }$ 作为一个循环周期，暴力算出 $\prod _{q\nmid i}^{q^{\alpha }}i$，然后循环实际上有 $\lfloor \frac{n}{q^{\alpha }}\rfloor$ 个，直接用快速幂。

对于多出来的部分，长度一定小于 $q^{\alpha }$，也可以暴力算。

以 $n=19,p=3,\alpha =2$ 为例：
$$\begin{array}{rl}19!& =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\\ & =(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 (3\times 6\times 9\times 12\times 15\times 18)\\ & =(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 3^6\times (1\times 2\times 3\times 4\times 5\times 6)\\ & \equiv (1\times 2\times 4\times 5\times 7\times 8{)}^2\times 19\times 3^6\times 6!(\mathrm{m}\mathrm{o}\mathrm{d}{3}^2)\\ & =\prod _{3\nmid i}^{19}i\times 3^{\lfloor \frac{19}{3}\rfloor }\times \lfloor \frac{19}{3}\rfloor !\end{array}$$
再递归计算 $6!\mathrm{m}\mathrm{o}\mathrm{d}{3}^2$ 即可。

递归部分的代码如下：

ll cal(ll n, ll p, ll pa)
{
	if (!n)
	{
		return 1;
	}
	ll ans = 1;
	for (int i = 1; i <= pa; i++) //循环部分
	{
		if (i % p)
		{
			ans = ans * i % pa;
		}
	}
	ans = qpow(ans, n / pa, pa);
	for (int i = 1; i <= n % pa; i++) //多出部分
	{
		if (i % p)
		{
			ans = ans * i % pa;
		}
	}
	return ans * cal(n / p, p, pa) % pa;
}

注意到这部分其实求的就是 $\frac{n!}{q^x}$。

对于 $m!,(n-m)!$ 的处理同理。

Part 4：合并

通过 Part 3 可以求出 $\frac{n!}{q^x},\frac{m!}{q^y},\frac{(n-m)!}{q^z}$，现在只需要求出 $x,y,z$ 就可以得到 $q^{x-y-z}$ 了。

怎么算想必小奥都讲过了吧，比如说
$$x=\lfloor \frac{n}{q}\rfloor +\lfloor \frac{n}{q^2}\rfloor +\lfloor \frac{n}{q^3}\rfloor +\cdots$$
至此可求出 $\left(\genfrac{}{}{0ex}{}{n}{m}\right)\mathrm{m}\mathrm{o}\mathrm{d}{q}_i^{{\alpha }_i}$，再结合 Part 1 使用 CRT 合并就完成了求解 $\left(\genfrac{}{}{0ex}{}{n}{m}\right)\mathrm{m}\mathrm{o}\mathrm{d}p$ 的过程。

Code

//18 = 9 + 9 = 18.
#include <iostream>
#include <cstdio>
#define Debug(x) cout << #x << "=" << x << endl
typedef long long ll;
using namespace std;

ll qpow(ll a, ll b, ll p)
{
	ll base = a, ans = 1;
	while (b)
	{
		if (b & 1)
		{
			ans = ans * base % p;
		}
		base = base * base % p;
		b >>= 1;
	}
	return ans;
}

ll cal(ll n, ll p, ll pa)
{
	if (!n)
	{
		return 1;
	}
	ll ans = 1;
	for (int i = 1; i <= pa; i++)
	{
		if (i % p)
		{
			ans = ans * i % pa;
		}
	}
	ans = qpow(ans, n / pa, pa);
	for (int i = 1; i <= n % pa; i++)
	{
		if (i % p)
		{
			ans = ans * i % pa;
		}
	}
	return ans * cal(n / p, p, pa) % pa;
}

ll cnt_p(ll n, ll m, ll p)
{
	ll cnt = 0;
	for (ll i = p; i <= n; i *= p)
	{
		cnt += n / i;
	}
	for (ll i = p; i <= m; i *= p)
	{
		cnt -= m / i;
	}
	for (ll i = p; i <= n - m; i *= p)
	{
		cnt -= (n - m) / i;
	}
	return cnt;
}

ll x, y;

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

ll inv(ll a, ll p)
{
	exgcd(a, p);
	x = (x % p + p) % p;
	return x;
}

ll C(ll n, ll m, ll p, ll pa)
{
	ll a = cal(n, p, pa), b = cal(m, p, pa), c = cal(n - m, p, pa), cnt = cnt_p(n, m, p);
	return a * inv(b, pa) % pa * inv(c, pa) % pa * qpow(p, cnt, pa) % pa;
}

ll a[10], b[10];

ll CRT(int n)
{
	ll m = 1;
	for (int i = 1; i <= n; i++)
	{
		m *= a[i];
	}
	ll ans = 0;
	for (int i = 1; i <= n; i++)
	{
		ll mi = m / a[i];
		ll Mi = inv(mi, a[i]);
		ans = (ans + b[i] * mi % m * Mi % m) % m;
	}
	return ans;
}

ll exLucas(ll n, ll m, ll p)
{
	int k = 0;
	for (ll i = 2; i * i <= p; i++)
	{
		if (p % i == 0)
		{
			a[++k] = 1;
			while (p % i == 0)
			{
				a[k] *= i;
				p /= i;
			}
			b[k] = C(n, m, i, a[k]);
		}
	}
	if (p > 1)
	{
		a[++k] = p;
		b[k] = C(n, m, p, p);
	}
	return CRT(k);
}

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

Complexity

计算 $\frac{n!}{q^x}\mathrm{m}\mathrm{o}\mathrm{d}{q}^{\alpha }$ 的时间为 $O(q^{\alpha }\log n)$。

计算 $x$ 的时间为 $O(\log n)$。

所以计算 $\left(\genfrac{}{}{0ex}{}{n}{m}\right)\mathrm{m}\mathrm{o}\mathrm{d}{q}^{\alpha }$ 的时间为 $O(q^{\alpha }\log n)$。

整个预处理部分就是 $O(\sqrt{p}+\sum _{i=1}^k{q}_i^{{\alpha }_i}\log n)$。

CRT 部分是 $O(k\log \prod _{i=1}^k{q}_i^{{\alpha }_i})=O(k\log p)$ 的。

综上，exLucas 的时间复杂度为 $O(\sqrt{p}+\sum _{i=1}^k{q}_i^{{\alpha }_i}\log n+k\log p)\sim O(p\log n)$。

Problems

模板

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

模板 + 乘法原理

P2183 [国家集训队]礼物

答案为 $\left(\genfrac{}{}{0ex}{}{n}{w_1}\right)\times \left(\genfrac{}{}{0ex}{}{n-w_1}{w_2}\right)\times \left(\genfrac{}{}{0ex}{}{n-w_1-w_2}{w_3}\right)\times \cdots \mathrm{m}\mathrm{o}\mathrm{d}P$。

References

• [1] JokerJim：【算法】扩展卢卡斯详解
• [2] YangTY：Lucas 定理
• [3] OI Wiki：卢卡斯定理