数论补完计划 Part5 组合计数

1
#include <bits/stdc++.h>
2
using namespace std;
3

4
typedef long long ll;
5
typedef pair<ll, ll> Equation;
6

7
ll exgcd(ll a, ll b, ll &x, ll &y) {
8
    if (!b) {
9
        x = 1, y = 0;
10
        return a;
11
    }
12
    ll d = exgcd(b, a % b, y, x);
13
    y -= a / b * x;
14
    return d;
15
}
16

17
ll inv(ll a, ll p) {
18
    ll ans, tmp;
19
    exgcd(a, p, ans, tmp);
20
    return (ans % p + p) % p;
21
}
22

23
ll qmul(ll a, ll b, ll p) {
24
    ll res = 0;
25
    while (b) {
26
        if (b & 1) res = (res + a) % p;
27
        a = (a + a) % p;
28
        b >>= 1;
29
    }
30
    return res % p;
31
}
32

33
ll qpow(ll a, ll b, ll p) {
34
    ll res = 1;
35
    while (b) {
36
        if (b & 1) res = (res % p * a % p) % p;
37
        a = (a * a) % p;
38
        b >>= 1;
39
    }
40
    return res % p;
41
}
42

43
ll mul(initializer_list<ll> args, ll p) {
44
    // 很方便的变长参数列表求连乘，int128 防止爆范围
45
    __int128 res = 1;
46
    for (ll i: args) {
47
        res = (res % p * i % p) % p;
48
    }
49
    return (ll) (res % p);
50
}
51

52
ll CRT(vector<Equation> &v) {
53
    // 中国剩余定理部分
54
    ll M = 1;
55
    ll res = 0;
56
    for (Equation eq: v) M *= eq.second;
57
    for (Equation i: v) {
58
        ll Mi = M / i.second;
59
        res = (res + qmul(qmul(i.first, Mi, M), inv(Mi, i.second), M)) % M;
60
    }
61
    return res % M;
62
}
63

64
ll resolve(ll a, ll p, ll ppow) {
65
    // 递归过程求解，即实现例3.1
66
    if (!a) return 1; // 1的阶乘模质数当然无论如何都是1
67
    ll s = 1;
68
    for (ll i = 1; i <= ppow; i++) {
69
        if (i % p) s = mul({s, i}, ppow); // 完整组循环节连乘
70
    }
71
    s = qpow(s, a / ppow, ppow); // 对完整组的积求幂
72
    for (ll i = a / ppow * ppow + 1; i <= a; i++) {
73
        if (i % p) s = mul({s, i}, ppow); // 对不完整组求积
74
    }
75
    return mul({s, resolve(a / p, p, ppow)}, ppow); // 递归过程
76
}
77

78
ll getA(ll n, ll m, ll p, ll ppow) {
79
    // 获取 CRT 方程组中的 a
80
    ll x = 0, y = 0, z = 0;
81
    // 统计 n, m, n - m 分别是 p 的多少次幂
82
    for (ll i = n; i; i /= p) x += i / p;
83
    for (ll i = m; i; i /= p) y += i / p;
84
    for (ll i = n - m; i; i /= p) z += i / p;
85
    // 套公式
86
    return mul({
87
        qpow(p, x - y - z, ppow),
88
        resolve(n, p, ppow),
89
        inv(resolve(m, p, ppow), ppow),
90
        inv(resolve(n - m, p, ppow), ppow)
91
    }, ppow);
92
}
93

94
ll exLucas(ll n, ll m, ll p) {
95
    // 分解质因数
96
    vector<Equation> E;
97
    for (ll i = 2; i * i <= p; i++) {
98
        if (p % i == 0) {
99
            // 统计对应质数以及它的幂
100
            ll c = 1;
101
            while (p % i == 0) {
102
                c *= i;
103
                p /= i;
104
            }
105
            E.emplace_back(getA(n, m, i, c), c);
106
        }
107
    }
108
    if (p > 1) E.emplace_back(getA(n, m, p, p), p); // 最后一个质数
109
    return CRT(E);
110
}
111

112
int main() {
113
    ios::sync_with_stdio(false);
114
    cin.tie(nullptr);
115
    cout.tie(nullptr);
116

117
    ll n, m, p;
118
    cin >> n >> m >> p;
119
    cout << exLucas(n, m, p) << endl;
120

121
    return 0;
122
}