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

参考链接：https://oi.wiki/math/number-theory/lucas/#exlucas-算法

解题思路：

扩展卢卡斯定理（Extended Lucas Theorem, exLucas），用于在模数 \(p\) 不一定是质数 的情况下，求组合数 \(C_n^m \pmod p\) 的值。

素数幂模的情形

首先考虑模数是质数的幂（\(p^a\)）的形式。

由 Lucas定理 得：

\[n! = p^{v_p(n!)} (n!)_p \]

其中，\(v_p(n!)\) 是 \(n!\) 的质因子分解中 \(p\) 的幂次。\((n!)_p\) 是 \(n!\) 中除去所有因数 \(p\) 之后的结果。

显然，\(p^{v_p(n!)}\) 和 \((n!)_p\) 互质。

因此，组合数 \(\binom{n}{k}\)（即 \(C_n^m\)）可以写做

\[\binom{n}{m} = \frac{ n! }{ (n-m)! \cdot m! } \]

\[= \frac{ p^{v_p(n!)} (n!)_p }{ ( p^{v_p((n-m)!)} ((n-m)!)_p ) \cdot ( p^{v_p(m!)} ) (m!)_p } \]

\[= p^{ v_p(n!) - v_p((n-m)!) - v_p(m!) } \cdot \frac{ (n!)_p }{ ((n-m)!)_p \cdot (m!)_p } \]

式子中的 \(v_p(n!)\)、\(v_p((n-m)!)\)、\(v_p(m!)\) 可以通过 勒让德公式 计算；
至于 \((n!)_p \pmod {p^a}\) 的求法，可以参考 这篇博客 ，使用递归分块法进行求解。

一般模数的情形

对于 \(m\) 是一般合数的情形，只需要首先对它做质因数分解：

\[m = p_1^{\alpha_1} p_2^{\alpha_2} \cdot p_s^{\alpha_s} \]

然后，分别计算出模 \(p_i^{\alpha_i}\) 下组合数 \(\binom{n}{m}\) 的余数，就得到 \(s\) 个同余方程：

\[\begin{cases} \binom{n}{m} \equiv r_1 \pmod { p_1^{\alpha_1} } \\ \binom{n}{m} \equiv r_2 \pmod { p_2^{\alpha_2} } \\ \binom{n}{m} \equiv r_3 \pmod { p_3^{\alpha_3} } \\ \binom{n}{m} \equiv r_4 \pmod { p_4^{\alpha_4} } \end{cases} \]

最后，利用 中国剩余定理 求出模 \(m\) 的余数即可。

示例程序

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
const int maxn = 1e6 + 5;

ll qpow(ll a, ll b, ll p) {
    ll res = 1;
    for (ll t = a%p; b; b >>= 1, t = t * t % p) {
        if (b & 1ll)
            (res *= t) %= p;
    }
    return res;
}

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

ll inv(ll a, ll n) {
    ll x, y;
    ll d = exgcd(a, n, x, y);
    assert(d == 1);
    return (x + n) % n;
}

ll n, m, p, fac[maxn] = {1};

void init(ll p, ll q) {
    for (int i = 1; i < q; i++) {
        fac[i] = fac[i-1];
        if (i % p)
            (fac[i] *= i) %= q;
    }
}

ll legendre(ll n, ll p) {
    ll res = n;
    for (; n; n /= p)
        res -= n % p;
    return res / (p - 1);
}

ll f(ll n, ll p, ll q) {
    if (n < q) return fac[n];
    ll c = n / q, res = fac[n % q];
    if (c % 2 && !(p == 2 && q >= 8))
        res *= -1;
    res = (res + q) % q;
    return res;
}

ll cal(ll n, ll p, ll q) {
    ll res = 1;
    for (ll i = n; i; i /= p)
        (res *= f(i, p, q)) %= q;
    return res;
}

// C(n, m) p^a=q
ll get_r(ll n, ll m, ll p, ll a, ll q) {
    ll cnt = legendre(n, p) - legendre(n-m, p) - legendre(m, p);
    if (cnt >= a)
        return 0;
    init(p, q);
    ll res = qpow(p, cnt, q) * cal(n, p, q) % q * inv(cal(n-m, p, q), q) % q * inv(cal(m, p, q), q) % q;
    return (res + q) % q;
}

struct Node {
    ll b, r; // b除数, r余数
};
vector<Node> crts;

ll crt() {
    ll n = 1, x = 0;
    for (auto &[b, r] : crts) {
        n *= b;
    }
    for (auto &[b, r] : crts) {
        ll m = n / b;
        ll c = m * inv(m, b);
        x += r * c;
    }
    return (x % n + n) % n;
}

int main() {
    cin >> n >> m >> p;
    for (ll i = 2; i * i <= p; i++) {
        if (p % i == 0) {
            ll a = 0, q = 1;
            for (; p % i == 0; a++, q *= i, p /= i);
            ll r = get_r(n, m, i, a, q);
            crts.push_back({q, r});
        }
    }
    if (p > 1) {
        ll r = get_r(n, m, p, 1, p);
        crts.push_back({p, r});
    }

    ll ans = crt();
    printf("%lld\n", ans);

    return 0;
}