组合意义天地灭。
Lucas 定理
问题 \(1\)：给定 \(n, m \in \mathbb{N}\) 与 \(p \in \mathbb{P}\)，其中 \(n\) 与 \(m\) 相当大，而 \(p\) 则相对较小，要求计算 \(\binom{n}{m} \bmod p\) 的值。
一般的预处理逆元以及递推的方法在 \(n, m\) 充分大时均会失效，我们需要新的工具来解决这类问题。
对于 \(n, m \in \mathbb{N}\) 与 \(p \in \mathbb{P}\)，有 \(\binom{n}{m} \equiv \prod{\binom{n_i}{m_i}} \pmod{p}\)，其中 \(n = \sum_{i \in \mathbb{N}} n_i p^i\)，\(m = \sum_{i \in \mathbb{N}} m_i p^i\)，也即 \(n\) 与 \(m\) 的 \(p\) 进制表示。
上边的这一坨就是 Lucas 定理了。不过这貌似跟我们平常使用的递归形式不太一样？
对于 \(n, m \in \mathbb{N}\) 与 \(p \in \mathbb{P}\)，有 \(\binom{n}{m} \equiv \binom{n \bmod p}{m \bmod p} \binom{\lfloor \frac{n}{p} \rfloor}{\lfloor \frac{m}{p} \rfloor} \pmod{p}\)。
这是我们平常见到的形式。不过，稍作观察就可以发现，这两种形式是等价的，在此不做证明。
下文中，如无特殊说明，均有 \(n, m \in \mathbb{N}\) 且 \(p \in \mathbb{P}\)。
下面我们尝试证明 Lucas 定理。
引理 \(1\)：对于任意 \(x \in \mathbb{N}\) 且 \(x < p\)，\(\binom{p}{x} \equiv 0 \pmod{p}\)。
这是好证明的，我们直接从定义下手即可。由于 \(\binom{p}{x} = \frac{p!}{x!(p - x)!}\)，又 \(x < p\)，因此 \(x!\) 中一定没有质因子 \(p\)，也即 \(x!\) 在 \(\bmod p\) 意义下存在逆元 \((x!)^{-1}\)。而 \((p - x)!\) 同理，因此有：
\[\binom{p}{x} \equiv p!(x!)^{-1}((p - x)!)^{-1} \pmod{p} \]右边显然为 \(p\) 的倍数。\(\square\)
引理 \(2\)：对于任意 \(x \in \mathbb{N}\)，有 \((1 + x)^p \equiv 1 + x^p \pmod{p}\)。
通过二项式定理，我们得到 \((1 + x)^p = \sum_{0 \le i \le p} \binom{p}{i} x^i\)，而根据引理 \(1\)，我们立即得到上式中除去首项 \(\binom{p}{0} x^0 = 1\) 与末项 \(\binom{p}{p} x^p = x^p\) 外，其余所有项在 \(\bmod p\) 意义下均为 \(0\)，因此引理得证。\(\square\)
现在来做一些推导：
\[\begin{aligned} (1 + x)^m &= (1 + x)^{{\lfloor\frac{m}{p}\rfloor} p + m \bmod p} \\ &= ((1 + x)^p)^{\lfloor\frac{m}{p}\rfloor} (1 + x)^{m \bmod p} \\ &\equiv (1 + x^p)^{\lfloor\frac{m}{p}\rfloor} (1 + x)^{m \bmod p} \\ &= \sum_{0 \le i \le {\lfloor\frac{m}{p}\rfloor}} \binom{\lfloor\frac{m}{p}\rfloor}{i} x^{ip} \sum_{0 \le i \le m \bmod p} \binom{m \bmod p}{i} x^i \\ &= \sum_{0 \le i \le {\lfloor\frac{m}{p}\rfloor}} \sum_{0 \le j \le m \bmod p} \binom{\lfloor\frac{m}{p}\rfloor}{i} \binom{m \bmod p}{j} x^{ip + j} \pmod{p} \end{aligned} \]第 \(3\) 步的转化使用了引理 \(2\)。容易看出 \(ip + j\) 取遍了所有 \([0, m]\) 的整数，于是考虑枚举 \(k = ip + j\)，则有：
\[(1 + x)^m \equiv \sum_{0 \le k \le m} \binom{\lfloor\frac{m}{p}\rfloor}{\lfloor\frac{k}{p}\rfloor} \binom{m \bmod p}{k \bmod p} x^k \pmod{p} \]又由二项式定理，我们立即得到：
\[\sum_{0 \le k \le m} \binom{m}{k} x^k \equiv \sum_{0 \le k \le m} \binom{\lfloor\frac{m}{p}\rfloor}{\lfloor\frac{k}{p}\rfloor} \binom{m \bmod p}{k \bmod p} x^k \pmod{p} \]那么 Lucas 定理是显而易见的。\(\square\)
实际使用中，由于给出的 \(p\) 一般不会超过 \(10^6\)，我们通常考虑直接计算 \(\binom{n \bmod p}{m \bmod p}\) 的值，而 \(\binom{\lfloor \frac{n}{p} \rfloor}{\lfloor \frac{m}{p} \rfloor}\) 的值则递归的使用 Lucas 定理进行计算，其边界条件为 \(m = 0\) 时返回 \(1\)。如果预处理所有 \(i \in [0, p)\) 的阶乘与其逆元的话，我们便可以 \(\mathcal{O}(1)\) 计算前边一部分的值，那么容易得到一个 \(\mathcal{O}(p \log p) + \mathcal{O}(\log_p n)\) 的算法。
扩展 Lucas 定理
这玩意虽然叫定理，但事实上只是个算法。
发现 Lucas 定理只能解决 \(p \in \mathbb{P}\) 的问题，而某些时候题目给定的模数不一定是质数，这个时候我们该怎么办呢？
问题 \(2\)：给定 \(n, m, p \in \mathbb{N}\)，其中 \(n\) 与 \(m\) 相当大，要求计算 \(\binom{n}{m} \bmod p\) 的值。
由于未保证 \(p\) 为质数，普通的 Lucas 定理对此也无能为力。我们考虑将 \(p\) 质因数分解，对每一个因子求解后得到一组同余方程，再用 CRT 来合并答案。
对于 \(p \in \mathbb{P}\)，有：
\[\begin{aligned} \binom{n}{m} \bmod p^k &= \frac{n!}{m!(n - m)!} \bmod p^k \\ &= \frac{\frac{n!}{p^{g(n)}}}{\frac{m!}{p^{g(m)}} \frac{(n - m)!}{p^{g(n - m)}}} p^{g(n) - g(m) - g(n - m)} \bmod p^k \\ \end{aligned} \]其中 \(g(x)\) 为 \(x!\) 中质因子 \(p\) 的次数。由于 \(\frac{x!}{p^{g(x)}} \perp p^k\)，那么 \(\frac{x!}{p^{g(x)}}\) 存在 \(\bmod p^k\) 的逆元，于是我们考虑求出这个式子的值：
\[\frac{n!}{p^{g(n)}} \bmod p^k \]考虑记上边的式子为 \(f(n)\)，然后再来做一些推导，我们考虑将 \(n!\) 中所有的 \(p\) 的倍数提出来，容易看出一共有 \(\lfloor\frac{n}{m}\rfloor\) 个这样的数：
\[\begin{aligned} f(n) &= \frac{n!}{p^{g(n)}} \\ &= \frac{(p^{\lfloor\frac{n}{p}\rfloor} \prod _{1 \le i \le {\lfloor\frac{n}{p}\rfloor}} i)(\prod_{1 \le i \le n, i \not \equiv 0 \pmod{p}} i)}{p^{g(n)}} \\ &= \frac{(p^{\lfloor\frac{n}{p}\rfloor} ({\lfloor\frac{n}{p}\rfloor})!) (\prod_{1 \le i \le n, i \not \equiv 0 \pmod{p}} i)}{p^{g(n)}} \\ \end{aligned} \]由于 \(p\) 的倍数在 \(1 \sim n\) 中分布均匀，并且 \(1 \sim p^k\) 在 \(\bmod p^k\) 下构成一个完全剩余系，那么后面的积式是有循环节的，那么有：
\[\begin{aligned} f(n) &= \frac{(p^{\lfloor\frac{n}{p}\rfloor} ({\lfloor\frac{n}{p}\rfloor})!) (\prod_{1 \le i \le n, i \not \equiv 0 \pmod{p}} i)}{p^{g(n)}} \\ &\equiv \frac{(p^{\lfloor\frac{n}{p}\rfloor} ({\lfloor\frac{n}{p}\rfloor})!)(\prod_{1 \le i < p^k, i \not \equiv 0 \pmod{p}} i)^{\lfloor\frac{n}{p^k}\rfloor} (\prod_{\lfloor\frac{n}{p^k}\rfloor p^k \le i \le n, i \not \equiv 0 \pmod{p}} i)}{p^{g(n)}} \pmod{p^k} \end{aligned} \]这步变换不难理解。\(\prod_{1 \le i \le n, i \not \equiv 0 \pmod{p}} i\) 根据 \(\bmod p^k\) 的值可以分为 \(\lfloor\frac{n}{p^k}\rfloor\) 块，其中前边的块在没有去掉 \(p\) 的倍数前均构成 \(p^k\) 的一个完系，可以用快速幂求得，而最后一块是散块直接暴力累加就行。
下一步的化简就需要用到 \(g\) 的性质了。容易注意到 \(g\) 的一种递归式：
\[g(n) = \lfloor \frac{n}{p} \rfloor + g(\lfloor \frac{n}{p} \rfloor) \]这非常显而易见，不予证明。那么我们可以将式子作进一步化简：
\[\begin{aligned} f(n) &\equiv \frac{(p^{\lfloor\frac{n}{p}\rfloor} ({\lfloor\frac{n}{p}\rfloor})!)(\prod_{1 \le i < p^k, i \not \equiv 0 \pmod{p}} i)^{\lfloor\frac{n}{p^k}\rfloor} (\prod_{\lfloor\frac{n}{p^k}\rfloor p^k \le i \le n, i \not \equiv 0 \pmod{p}} i)}{p^{g(n)}} \\ &= \frac{p^{\lfloor\frac{n}{p}\rfloor} ({\lfloor\frac{n}{p}\rfloor})!}{p^{g(n)}}(\prod_{1 \le i < p^k, i \not \equiv 0 \pmod{p}} i)^{\lfloor\frac{n}{p^k}\rfloor} (\prod_{\lfloor\frac{n}{p^k}\rfloor p^k \le i \le n, i \not \equiv 0 \pmod{p}} i) \\ &= \frac{p^{\lfloor\frac{n}{p}\rfloor} p^{g({\lfloor\frac{n}{p}\rfloor})} f({\lfloor\frac{n}{p}\rfloor})}{p^{g(n)}}(\prod_{1 \le i < p^k, i \not \equiv 0 \pmod{p}} i)^{\lfloor\frac{n}{p^k}\rfloor} (\prod_{\lfloor\frac{n}{p^k}\rfloor p^k \le i \le n, i \not \equiv 0 \pmod{p}} i) \\ &= f({\lfloor\frac{n}{p}\rfloor}) (\prod_{1 \le i < p^k, i \not \equiv 0 \pmod{p}} i)^{\lfloor\frac{n}{p^k}\rfloor} (\prod_{\lfloor\frac{n}{p^k}\rfloor p^k \le i \le n, i \not \equiv 0 \pmod{p}} i) \pmod{p^k} \end{aligned} \]于是就做完了。最后再用 CRT 合并一下就可以求得答案。
通过分析可以得到，时间复杂度为 \(\mathcal{O}(p^k \log_p n)\)。
constexpr int N = 1050;
i64 n, m, p;
i64 a[N], b[N];
i64 fastPow(i64 base, i64 exp, i64 mod) {
i64 res = 1;
for (; exp; exp >>= 1) {
if (exp & 1) (res *= base) %= mod;
(base *= base) %= mod;
}
return res;
}
i64 exgcd(i64 a, i64 b, i64 &x, i64 &y) {
if (!b) return x = 1, y = 0, a;
i64 d = exgcd(b, a % b, y, x);
y -= a / b * x;
return d;
}
i64 inv(i64 num, i64 mod) {
i64 x, y;
exgcd(num, mod, x, y);
return (x + mod) % mod;
}
i64 f(i64 n, i64 p, i64 pk) {
if (!n) return 1;
i64 cir = 1, rst = 1;
for (i64 i = 1; i <= pk; i++) {
if (i % p) (cir *= i % pk) %= pk;
}
cir = fastPow(cir, n / pk, pk);
for (i64 i = pk * (n / pk); i <= n; i++) {
if (i % p) (rst *= i % pk) %= pk;
}
return f(n / p, p, pk) * cir % pk * rst % pk;
}
i64 g(i64 n, i64 p) {
i64 res = 0, cpy = p;
for (; p <= n; p *= cpy) res += n / p;
return res;
}
i64 C(i64 n, i64 m, i64 p, i64 pk) {
i64 a = f(n, p, pk), b = inv(f(m, p, pk), pk), c = inv(f(n - m, p, pk), pk);
return a * b % pk * c % pk *
fastPow(p, g(n, p) - g(m, p) - g(n - m, p), pk) % pk;
}
i64 exLucas(i64 n, i64 m, i64 mod) {
i64 cpy = mod, res = 0;
int cnt = 0;
for (i64 i = 2; i * i <= cpy; i++) {
if (cpy % i == 0) {
i64 pk = 1;
while (cpy % i == 0) cpy /= i, pk *= i;
a[++cnt] = pk, b[cnt] = C(n, m, i, pk);
}
}
if (cpy != 1) a[++cnt] = cpy, b[cnt] = C(n, m, cpy, cpy);
for (int i = 1; i <= cnt; i++) {
i64 k = mod / a[i], t = inv(k, a[i]);
(res += b[i] * k % mod * t % mod) %= mod;
}
return res;
}