Lucas 定理

前置知识：CRT，exgcd。

Lucas 定理

洛谷 P3807

证明见 oi-wiki

Lucas 定理：\(C_{n}^{m} = C_{n / p}^{m / p} \cdot C_{n \% p}^{m \% p} \pmod p\)，其中 \(p\) 为质数。一般适用于 \(p\) 较小的情况。

其中 \(C_{n \% p}^{m \% p} \pmod p\) 用逆进行求解。 当 \(n \% p < m \% p\) 时，返回 \(0\)。

\(C_{n / p}^{m / p}\) 继续用 \(Lucas\) 进行展开，递归求解。

int getC(int n, int m, int p) { // 返回 C[n][m] % p
  if (n < m) return 0; 
  return 1ll * fac[n] * facinv[n - m] * facinv[m] % p;
}

int Lucas(int n, int m, int p) { 
  if (!m) return 1;
  return 1ll * Lucas(n / p, m / p, p) * getC(n % p, m % p, p) % p; // 分两个部分处理
}

扩展卢卡斯定理 (exLucas)

洛谷 P4720

Lucas 定理中的 \(p\) 为质数，而 exLucas 解决 \(p\) 不为质数的情况。

模数为 \(p^k\)

先假设模数为 \(p^k\)，设 \(A(i)\) 为 \(i\) 质因数分解中 \(p\) 的指数，\(B(i) = \frac{i}{p ^ {A(i)}}\)，\(c = A(n!) -A(m!) - A((n - m)!)\) 则答案为 \(C_{n}^{m} = p^{c} \cdot \frac{B(n!)}{B(m!)B((n - m)!) }\)。可以发现若 \(c \geqslant k\)，则直接答案为 \(0\)。、

\(A(n!) = \sum\limits_{i = 1} \lfloor \frac{n}{p^i} \rfloor\)。

对于 \(B(n!)\)，以 \(n = 22, p^2 = 3^2 = 9\) 为例。最后一列等价于 \(B(\lfloor \frac{n}{p} \rfloor!) = B(7)\)，前两列每 \(\frac{p^2}{p} = p\) 行为 \(1\) 个循环节。

所以我们可以预处理 \(S(1 \sim p^k)\)，其中 \(S(0) = 1, S(i) = S(i - 1) \times \begin{cases} 1, p \mid i \\ i, p \nmid i \end{cases}\)，则 \(B(n!) = B(\lfloor \frac{n}{p} \rfloor!)S(p^k)^{\lfloor \frac{n}{p^k} \rfloor}S(n \% p^k)\)。

1 -> 1	2 -> 2	3 -> 1
4 -> 4	5 -> 5	6 -> 2
7- > 7	8 -> 8	9 -> 3
10 -> 1	11 -> 2	12 ->4
13 -> 4	14 -> 5	15 -> 5
16 -> 7	17 ->8	18 -> 6
19 -> 1	20 -> 2	21 -> 7
22 -> 4

int get(ll n, int x) { // A
  int ret = 0;
  i128 o = x;
  for (; o <= n; o *= x) {
    ret += n / o;
  }
  return ret;
}

for (int i = 1; i <= Mod; i++) { // Mod = p^k
  s[i] = s[i - 1] * ll(i % a[u] == 0 ? 1 : i) % Mod;
}

int work(ll n, int x) { // B
  ll ret = 1;
  for (; n; n /= x) {
    (ret *= qpow(s[Mod], n / Mod)) %= Mod;
    (ret *= s[n % Mod]) %= Mod;
  }
  return ret;
}

因为 \(B\) 中去除了 \(p\)，可以使用 exgcd 求出 \(B(m!), B((n - m)!)\) 的逆元，然后就可以求出，\(C_n^m\) 在模 \(p^k\) 意义下的逆元。

一般模数

对于一个一般模数，可以将 \(p\) 质因数分解得到 \(p = p_1^{k_1}p_2^{k_2}\dots p_n^{k_1}\)。

因为 \(p_i^{k_i}\) 两两互质，可以求出 \(C_{n}^m \equiv a_i \pmod {p_i^{k_i}}\)，然后使用中国剩余定理得到 \(C_{n}^m \% p\)。

时间复杂度懒得算了，参考 oi-wiki

#include <bits/stdc++.h>

using namespace std;
using ll = long long;
using i128 = __int128;

const int MAXP = 1e6 + 3;

ll n, m;
int p, Mod, k, a[25], b[25], s[MAXP] = {1}, x, y;

int get(ll n, int x) {
  int ret = 0;
  i128 o = x;
  for (; o <= n; o *= x) {
    ret += n / o;
  }
  return ret;
}

int qpow(int x, ll k) {
  ll ret = 1;
  if (!Mod) Mod = p + 1;
  for (; k; k >>= 1, x = 1ll * x * x % Mod) {
    if (k & 1) ret *= x;
  }
  return ret;
}

int work(ll n, int x) {
  ll ret = 1;
  for (; n; n /= x) {
    (ret *= qpow(s[Mod], n / Mod)) %= Mod;
    (ret *= s[n % Mod]) %= Mod;
  }
  return ret;
}

int exgcd(int a, int b) {
  if (!b) {
    x = 1, y = 0;
    return a;
  }
  int g = exgcd(b, a % b), t = y;
  y = x - a / b * y, x = t;
  return g;
}

int solve(int u) { // C(n, m) % (a[u] ^ b[u])
  int c = get(n, a[u]) - get(m, a[u]) - get(n - m, a[u]);
  if (c >= b[u]) return 0;
  for (int i = 1; i <= Mod; i++) {
    s[i] = s[i - 1] * ll(i % a[u] == 0 ? 1 : i) % Mod;
  }
  i128 ret = 1ll * qpow(a[u], c) * work(n, a[u]);
  exgcd(work(m, a[u]), Mod), ret *= x, exgcd(work(n - m, a[u]), Mod);
  return ret * x % Mod;
}

int main() {
  ios::sync_with_stdio(0), cin.tie(0);
  cin >> n >> m >> p;
  int t = p;
  for (int i = 2; i * i <= t; i++) { // 质因数分解
    if (t % i == 0) {
      a[++k] = i;
      while (t % i == 0) {
        t /= i, b[k]++;
      }
    }
  }
  if (t > 1) a[++k] = t, b[k] = 1;
  for (int i = 1; i <= k; i++) {
    Mod = 0, Mod = qpow(a[i], b[i]);
    b[i] = (solve(i) + Mod) % Mod, a[i] = Mod;
  }
  ll ans = 0;
  for (int i = 1; i <= k; i++) { // 中国剩余定理
    exgcd(p / a[i], a[i]);
    (ans += 1ll * b[i] * p / a[i] * x) %= p;
  }
  cout << (ans + p) % p;
  return 0;
}