拓展Lucas定理求C(n, m) % p

#include<bits/stdc++.h>
using namespace std;
const int maxn = 100000 + 10;
typedef long long LL;

LL Pow(LL n, LL m, LL mod) {
    LL ans = 1;
    while(m > 0) {
        if(m & 1) ans = (LL)ans * n % mod;
        n = (LL)n * n % mod; m >>= 1;
    }
    return ans;
}
LL Pow(LL n,LL m) {
    LL ans = 1;
    while(m > 0) {
        if(m & 1) ans = ans * n;
        n = n * n; m >>= 1;
    }
    return ans;
}
LL x, y;
LL exgcd(LL a, LL b) {
    if(a == 0) {
        x = 0, y = 1;
        return b;
    }LL r = exgcd(b%a, a);
    LL t = x; x = y - (b/a)*x; y = t;
    return r;
}
LL rev(LL a, LL b) { exgcd(a, b); return ((x % b) + b) % b; }
LL Calc(LL n, LL p, LL t) {
    if(n == 0) return 1;

    LL s = Pow(p, t), k = n / s, tmp = 1;
    for(LL i=1; i<=s; i++) if(i % p) tmp = (LL)tmp * i % s;

    LL ans = Pow(tmp, k, s);
    for(LL i=s*k + 1; i<=n; i++) if(i % p) ans = (LL)ans * i % s;

    return (LL)ans * Calc(n / p, p, t) % s;
}
LL C(LL n, LL m, LL p, LL t) {
    LL s = Pow(p, t), q = 0;
    for(LL i=n; i; i/=p) q += i / p;
    for(LL i=m; i; i/=p) q -= i / p;
    for(LL i=n-m; i; i/=p) q -= i / p;

    LL ans = Pow(p, q);
    LL a = Calc(n, p, t), b = Calc(m, p, t), c = Calc(n-m, p, t);
    return (LL)(ans * a * rev(b, s) * rev(c, s)) % s;
}
LL China(LL A[], LL M[], LL cnt) {
    LL ans = 0, m, n = 1;
    for(LL i=1; i<=cnt; i++) n *= M[i];
    for(LL i=1; i<=cnt; i++) {
        m = n / M[i];
        exgcd(M[i], m);
        ans = (ans + (LL)y * m * A[i]) % n;
    }
    return (ans + n) % n;
}
LL A[maxn], M[maxn], cnt;
LL Lucas(LL n, LL m, LL mod) {
    for(LL i=2; i*i <= mod; i++) if(mod % i == 0) {
        LL t = 0;
        while(mod % i == 0) t++, mod /= i;
        M[++cnt] = Pow(i, t);
        A[cnt] = C(n, m, i, t);
    }if(mod > 1) {
        M[++cnt] = mod;
        A[cnt] = C(n, m, mod, 1);
    }
    return China(A, M, cnt);
}
LL n, k, p;
int main() {
#ifndef ONLINE_JUDGE
    freopen("input.in", "r", stdin);
    freopen("output.out", "w", stdout);
#endif
    cin >> n >> k >> p;
    cout << Lucas(n, k, p) << endl;
    return 0;
}