扩展卢卡斯定理

10 / 07 / 2020 | 最后修改于 05 / 24 / 2023

OI Math

求：

\binom n m \bmod p

其中 $p$ 不保证为质数。

前置知识

中国剩余定理（CRT）

ExGCD求逆元

考虑令 $p=p_1^{k_1}p_2^{k_2}p_3^{k_3} \dots$

原问题转化为求：

\begin{cases} C_n^m \bmod p_1^{k_1} \\ C_n^m \bmod p_2^{k_2} \\ C_n^m \bmod p_3^{k_3} \\ \dotsm \end{cases}

的一个解，使用中国剩余定理合并即可。

针对其中一个式子 $C_n^m \bmod p^{k}$ 即 $\frac {n!} {m!(n-m)!} \bmod p^k$ ，使用求逆元的方法是不行的，因为 $p$ 不保证为质数， $m!(n-m)!$ 可能不存在逆元。

设 $f(n)=n! \bmod p^k$ ，进一步考虑，既然不互质，那就把他们的公因子都约掉，我们把阶乘展开看一看：

\begin{align} n! &= 1 \times 2 \times 3 \cdots n \\ &= (p \cdot 2p \cdot 3p \cdots) \cdot [1 \times 2 \times 3 \cdots (p-1) \cdot (p + 1) \cdot (p + 2) \cdots] \\ &= p^{\lfloor \frac{n}{p} \rfloor} (1 \times 2 \times 3 \cdots) \cdot [1 \times 2 \times 3 \cdots (p-1)] \cdot [(p + 1) \cdot (p + 2) \cdots)] \\ &= \large p^{\lfloor \frac{n}{p} \rfloor} (\frac{n}{p})! \prod_{p \nmid i}^n i \\ &= \large p^{\lfloor \frac{n}{p} \rfloor} (\frac{n}{p})! (\prod_{p \nmid i}^{p^k} i)^{\frac {n}{p^k}} \text{rest} \\ \end{align}

其中 $\text{rest}$ 是指最后不够一整段 $p^k$ 的剩下的项， $\displaystyle \large (\prod_{p \nmid i}^{p^k} i)$ 是在 $\bmod p^k$ 意义下的一个循环节，这样的循环节有 $\frac n {p^k}$ 个。

每这样一次，我们都可以从阶乘提取出 $p^{\lfloor \frac n p \rfloor}$ 这个因子，并把原问题变为一个规模更小的问题 $f(\frac n p)$ ，可以递归进行，并在过程中累加 $\lfloor \frac n p \rfloor$ ，记为 $h(n)$ ，最后的 $h(n)$ 就是 $n!$ 中 $p$ 这个因子的数量。

至此，我们成功把 $n!$ 里的 $p$ 的因子分离出来。回归原来的问题 $C_n^m \bmod p^{k}$ ：

\large \frac {\frac {n!}{p^{h(n)}}} {\frac {m!}{p^{h(m)}}\frac {(n-m)!}{p^{h(n-m)}}} p^{h(n)-h(m)-h(n-m)} \bmod p^k

即

\frac{f(n)}{f(m)f(n-m)} p^{h(n)-h(m)-h(n-m)} \bmod p^k

由于 $f(n)$ 已经是被约掉 $p$ 因子的，一定有逆元，于是CRT的一个式子就解决了，最后合并就好了。

Talk is cheap, show me your code

#include <cstdio>

typedef long long ll;
const int N = 100000;
int p, cnt, pri[N], k[N], cpoy_p;
ll n, m;

inline ll qpow(ll base, ll exp, const int &mod) {
  ll res = 1;
  while (exp) {
    if (exp & 1) res = res * base % mod;
    base = base * base % mod;
    exp >>= 1;
  }
  return res;
}

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

namespace CRT {

int a[N], ans = 0, pk;
ll inv, y;

inline int solve() {
  for (int i = 1; i <= cnt; ++i) {
    pk = qpow(pri[i], k[i], p + 1);
    exgcd(p / pk, pk, inv, y);
    ans += a[i] * inv % p * (p / pk) % p;
    if (ans >= p) ans -= p;
    if (ans < 0) ans += p;
  }
  return ans;
}

}  // namespace CRT

namespace Exlucas {

int prime, mod;  // mod = p^k
ll h;            //指数也要开long long

inline void MOD(ll &x, const int &mod) {
  if (x >= mod) x %= mod;
}

ll f(ll n) {
  if (n == 0) return 1;
  h += n / prime;
  ll res = 1, tmp = 1;
  res *= f(n / prime);
  MOD(res, mod);
  for (int i = 1; i <= mod; ++i) {
    if (i % prime) (tmp *= i) %= mod;
  }
  res *= qpow(tmp, n / mod, mod);
  MOD(res, mod);
  for (ll i = n / mod * mod; i <= n; ++i) {
    if (i % prime) {
      res *= i % mod;
      MOD(res, mod);
    }
  }
  return res;
}

inline ll solve(int i) {
  prime = pri[i], mod = qpow(prime, k[i], p + 1);
  int exp = 0;
  ll res = 1, inv, y;
  h = 0, res *= f(n), exp += h;

  h = 0, exgcd(f(m), mod, inv, y);
  inv = (inv % mod + mod) % mod;
  res *= inv;
  MOD(res, mod);
  exp -= h;

  h = 0, exgcd(f(n - m), mod, inv, y);
  inv = (inv % mod + mod) % mod;
  res *= inv;
  MOD(res, mod);
  exp -= h;

  return res * qpow(prime, exp, mod) % mod;
}

}  // namespace Exlucas

signed main() {
#ifndef ONLINE_JUDGE
#ifdef LOCAL
  freopen("testdata.in", "r", stdin);
  freopen("testdata.out", "w", stdout);
#else
  freopen("ExLucas.in", "r", stdin);
  freopen("ExLucas.out", "w", stdout);
#endif
#endif

  scanf("%lld%lld%d", &n, &m, &p);
  cpoy_p = p;
  for (int i = 2; i * i <= p; ++i) {
    if (cpoy_p % i == 0) {
      pri[++cnt] = i;
      while (cpoy_p % i == 0) { cpoy_p /= i, k[cnt]++; }
    }
  }
  if (cpoy_p != 1) pri[++cnt] = cpoy_p, k[cnt] = 1;
  for (int i = 1; i <= cnt; ++i) { CRT::a[i] = Exlucas::solve(i); }
  printf("%d", CRT::solve());
  return 0;
}