1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
| vector<pair<int, int>> factorize(int n) {
vector<pair<int, int>> factors;
for (int i = 2; static_cast<long long>(i) * i <= n; i++) {
if (n % i == 0) {
int t = 0;
for (; n % i == 0; n /= i)
++t;
factors.emplace_back(i, t);
}
}
if (n > 1)
factors.emplace_back(n, 1);
return factors;
}
constexpr int power(int base, i64 exp) {
int res = 1;
for (; exp > 0; base *= base, exp /= 2) {
if (exp % 2 == 1) {
res *= base;
}
}
return res;
}
constexpr int power(int base, i64 exp, int mod) {
int res = 1 % mod;
for (; exp > 0; base = 1LL * base * base % mod, exp /= 2) {
if (exp % 2 == 1) {
res = 1LL * res * base % mod;
}
}
return res;
}
int inverse(int a, int m) {
int g = m, r = a, x = 0, y = 1;
while (r != 0) {
int q = g / r;
g %= r;
swap(g, r);
x -= q * y;
swap(x, y);
}
return x < 0 ? x + m : x;
}
int solveModuloEquations(const vector<pair<int, int>> &e) {
int m = 1;
for (size_t i = 0; i < e.size(); i++) {
m *= e[i].first;
}
int res = 0;
for (size_t i = 0; i < e.size(); i++) {
int p = e[i].first;
res = (res + 1LL * e[i].second * (m / p) * inverse(m / p, p)) % m;
}
return res;
}
constexpr int N = 1E5;
class Binomial {
const int mod;
private:
const vector<pair<int, int>> factors;
vector<int> pk;
vector<vector<int>> prod;
static constexpr i64 exponent(i64 n, int p) {
i64 res = 0;
for (n /= p; n > 0; n /= p) {
res += n;
}
return res;
}
int product(i64 n, size_t i) {
int res = 1;
int p = factors[i].first;
for (; n > 0; n /= p) {
res = 1LL * res * power(prod[i].back(), n / pk[i], pk[i]) % pk[i] * prod[i][n % pk[i]] % pk[i];
}
return res;
}
public:
Binomial(int mod) : mod(mod), factors(factorize(mod)) {
pk.resize(factors.size());
prod.resize(factors.size());
for (size_t i = 0; i < factors.size(); i++) {
int p = factors[i].first;
int k = factors[i].second;
pk[i] = power(p, k);
prod[i].resize(min(N + 1, pk[i]));
prod[i][0] = 1;
for (int j = 1; j < prod[i].size(); j++) {
if (j % p == 0) {
prod[i][j] = prod[i][j - 1];
} else {
prod[i][j] = 1LL * prod[i][j - 1] * j % pk[i];
}
}
}
}
int operator()(i64 n, i64 m) {
if (n < m || m < 0) {
return 0;
}
vector<pair<int, int>> ans(factors.size());
for (int i = 0; i < factors.size(); i++) {
int p = factors[i].first;
int k = factors[i].second;
int e = exponent(n, p) - exponent(m, p) - exponent(n - m, p);
if (e >= k) {
ans[i] = make_pair(pk[i], 0);
} else {
int pn = product(n, i);
int pm = product(m, i);
int pd = product(n - m, i);
int res = 1LL * pn * inverse(pm, pk[i]) % pk[i] * inverse(pd, pk[i]) % pk[i] * power(p, e) % pk[i];
ans[i] = make_pair(pk[i], res);
}
}
return solveModuloEquations(ans);
}
};
|
