【P4720】 扩展卢卡斯

要注意的点还是很多的比如求n！ mod (p^k)时， 对剩余不足一周期的部分（即计算remain时）i可能很大， remain * i很容易爆long long，注意要取模。

#include <algorithm>
#include <iostream>

using namespace std;

using ll = long long;

ll gcd(ll a, ll b){
    return b == 0? a: gcd(b, a % b);
}
ll exgcd(ll a, ll b, ll & x, ll & y){
    if (!b){
        x = 1, y = 0;
        return a;
    }
    ll r = exgcd(b, a % b, y, x);
    y -= (a / b) * x;
    return r;
}

// (a, mod) = 1
ll inv(int a, int mod){
    ll x, y;
    exgcd(a, mod, x, y);
    return (x + mod) % mod;
}

ll qpow(ll a, ll b, ll mod){
    ll res = 1;
    a %= mod;
    while (b){
        if (b & 1ll) res = (res * a) % mod;
        a = (a * a) % mod;
        b >>= 1;
    }
    return res;
}

// cal  n!/p^x mod(p^k)
// get x separately later
ll F(ll n, ll p, ll pk){
    if (n == 0) return 1;
    ll cy = 1, remain = 1;
    for (ll i = 1; i <= pk; ++i) {
        if (i % p) cy = (cy * i) % pk;
    }
    cy = qpow(cy, n / pk, pk);
    for (ll i = (pk * (n / pk)); i <= n; ++i) {
        // i % pk
        if (i % p) remain = remain * (i % pk) % pk;
    }
    return F(n / p, p, pk) * cy % pk * remain % pk;
}

// count x in n! / p^x
ll count(ll n, ll p){
    if (n < p) return 0;
    return count(n / p, p) + n / p;
}

ll C_pk(ll n, ll m, ll p, ll pk){
    ll u = F(n, p, pk);
    ll v = inv(F(m, p, pk), pk);
    ll w = inv(F(n - m, p, pk), pk);
    ll exp = qpow(p, count(n, p) - count(m, p) - count(n - m, p), pk);
    return u * v % pk * w % pk * exp % pk;
}

ll A[1001], B[1001];
// x = b(mod a)

ll exLucas(ll n, ll m, ll p){
    int tot = 0;
    ll p2 = p;
    for (ll div = 2; div * div <= p; ++div) {
        if (p2 % div == 0){
            ll pk = 1;
            while (p2 % div == 0){
                p2 /= div;
                pk *= div;
            }
            A[++tot] = pk;
            B[tot] = C_pk(n, m, div, pk);
        }
    }
    if (p2 > 1){
        A[++tot] = p2;
        B[tot] = C_pk(n, m, p2, p2);
    }
    //A[i] is co-prime
    // use CRT to find the minial x
    ll res = 0;
    for (int i = 1; i <= tot; ++i) {
        ll M = p / A[i];
        res = (res + B[i] * M % p * inv(M, A[i]) % p) % p;
    }
    return res;
}


#define debug

int main(){
#ifdef debug
    freopen("in.txt", "r", stdin);
#endif
    ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);
    ll n, m, p;
    cin >> n >> m >> p;
    cout << exLucas(n, m, p) << '\n';
    return 0;
}