$\therefore p^{e(n)} = p^{\lfloor \frac{n}{p} \rfloor} \times p^{e(\lfloor \frac{n}{p} \rfloor)}$。
$\therefore f(n) = \frac{S(n) \times p^{\lfloor \frac{n}{p} \rfloor} \times \lfloor \frac{n}{p} \rfloor !}{p^{\lfloor \frac{n}{p} \rfloor} \times p^{e(\lfloor \frac{n}{p} \rfloor)}} = \frac{S(n) \times \lfloor \frac{n}{p} \rfloor !}{p^{e(\lfloor \frac{n}{p} \rfloor)}}$。
仔细后面这一部分： $\frac{\lfloor \frac{n}{p} \rfloor !}{p^{e(\lfloor \frac{n}{p} \rfloor)}}$，根据 $f(n) = \frac{n!}{p^{e(n)}}$，则 $\frac{\lfloor \frac{n}{p} \rfloor !}{p^{e(\lfloor \frac{n}{p} \rfloor)}} = f(\lfloor \frac{n}{p} \rfloor)$。
$\therefore f(n) = S(n) \times f(\lfloor \frac{n}{p} \rfloor)$。
至此，我们已经推导了怎么快速求出 $f(n)$，接下来就是套公式 $C_n^m = \sum_i^q \frac{f(n)}{f(m) \times f(n-m)} \times {p_i}^{e(n) - e(m) - e(n-m)}$。
---
### 代码详解（各部分需要的函数怎么求？）
① 求乘法逆元部分
使用拓展欧几里得定理求逆元，此处不过多解释。
```cpp
#define pii pair<ll, ll>
pii exgcd (ll a, ll b) {
if (b == 0) return {1, 0};
pii res = exgcd(b, a % b);
ll x = res.first;
ll y = res.second;
return {y, x - y * (a / b)};
}
ll get_inv (ll x, ll _mod) {
pii res = exgcd(x, _mod);
ll inv = res.first;
return (inv % _mod + _mod) % _mod;
}
```
② 快速幂部分
此处需要注意！
由于我们需要求 ${p_i}^{k_i}$，然而给出的模数 $P$ 有可能是某一个 ${p_i}^{k_i}$ 的因数，此时 ${p_i}^{k_i} \bmod P = 0$，后续的操作中会出现对 ${p_i}^{k_i}$ 取模的情况。若 ${p_i}^{k_i}$ 为0，则会RE。所以我们要写两个快速幂（或者写一个带不同状态的快速幂）。用于求带取模的快速幂和不带取模的快速幂。
```cpp
ll power (ll a, ll b, ll _mod) {
ll res = 1 % _mod;
ll base = a % _mod;
while (b != 0) {
if (b & 1) res = res * base % _mod;
base = base * base % _mod;
b >>= 1;
}
return res;
}
ll power2 (ll a, ll b) { //不带模
ll res = 1;
ll base = a;
while (b != 0) {
if (b & 1) res = res * base;
base = base * base;
b >>= 1;
}
return res;
}
```
③ 求 $e(n)$ 函数部分
直接按照定义式 $e(n) = \sum_{i=1}^{p^i \le n} \lfloor \frac{n}{p^i} \rfloor$ 求即可。
```cpp
ll e (ll n, ll p) { //模p意义下的e(n)
ll res = 0;
while (n != 0) {
res += n / p;
n /= p;
}
return res;
}
```
④ 求 $f(n)$ 函数部分
如果直接求式子 $f(n) = S(n) \times f(\lfloor \frac{n}{p} \rfloor)$ 的值，时间复杂度是 $O(n)$，可能 TLE。
考虑优化求 $S(n)$ 值的过程。
观察连续的整数序列，发现它们在模 $p^k$ 时是存在循环节的，周期正好就是 $p^k$，即 $i \equiv i + p^k \pmod{p^k}
根据这个性质，可以用带余除法的思想，把长度为 n 的连乘序列拆成两部分：若干个长度为 p^k 的完整块，以及最后剩下的一个不完整块。
对于任意一个完整块，其中不能被 p 整除的数的乘积在模 p^k 意义下都是完全相等的。即 \prod_{i=1, p \nmid i}^{p^k} i \equiv \prod_{i=1, p \nmid i}^{p^k} (i + p^k) \equiv \prod_{i=1, p \nmid i}^{p^k} (i + 2p^k) \dots \pmod {p^k}。
所以只需要算出第一个完整块的乘积即可。令 \text{sum} = \prod_{i=1, p \nmid i}^{p^k} i (\bmod \ {p^k})。
ll f (ll n, ll p, ll k) {
if (n == 0) return 1;
ll pk = power2(p, k);
ll sum = 1;
for (ll i = 1; i <= pk; i++) {
if (i % p != 0) sum = sum * i % pk;
}
ll sn = power(sum, n / pk, pk);
for (ll i = 1; i <= n % pk; i++) {
if (i % p) sn = sn * i % pk;
}
return sn * f(n / p, p, k) % pk;
}
注意：当枚举质因数的 for 循环结束后，若剩余的数（代码中的 t）仍大于 1，说明 P 中包含一个大于 \sqrt{P} 的大质因子。由于该质因子的幂次必然为1，需将其作为底数、指数设为1，再单独执行一次局部解的求解与 CRT 合并。
ll exLucas (ll n, ll m) {
int t = mod;
ll res = 0;
for (ll i = 2; i * i <= t; i++) {
if (t % i == 0) { //分解质因数
ll k = 0;//指数
while (t % i == 0) {
t /= i;
k++;
}
ll c = C(n, m, i, k);
res = (res + crt(c, i, k)) % mod;
}
}
if (t > 1) {
ll c = C(n, m, t, 1);
res = (res + crt(c, t, 1)) % mod;
}
return res;
}
完结撒花！
附完整代码
#include<bits/stdc++.h>
using namespace std ;
#define ll long long
#define pii pair<ll, ll>
ll mod;
pii exgcd (ll a, ll b) {
if (b == 0) return {1, 0};
pii res = exgcd(b, a % b);
ll x = res.first;
ll y = res.second;
return {y, x - y * (a / b)};
}
ll get_inv (ll x, ll _mod) {
pii res = exgcd(x, _mod);
ll inv = res.first;
return (inv % _mod + _mod) % _mod;
}
ll power (ll a, ll b, ll _mod) {
ll res = 1 % _mod;
ll base = a % _mod;
while (b != 0) {
if (b & 1) res = res * base % _mod;
base = base * base % _mod;
b >>= 1;
}
return res;
}
ll power2 (ll a, ll b) { //不带模
ll res = 1;
ll base = a;
while (b != 0) {
if (b & 1) res = res * base;
base = base * base;
b >>= 1;
}
return res;
}
ll e (ll n, ll p) { //模p意义下的e(n)
ll res = 0;
while (n != 0) {
res += n / p;
n /= p;
}
return res;
}
ll f (ll n, ll p, ll k) { //模p^k意义下的f(n)
if (n == 0) return 1;
ll pk = power2(p, k);
ll sum = 1;
for (ll i = 1; i <= pk; i++) {
if (i % p != 0) sum = sum * i % pk;
}
ll sn = power(sum, n / pk, pk);
for (ll i = 1; i <= n % pk; i++) {
if (i % p) sn = sn * i % pk;
}
return sn * f(n / p, p, k) % pk;
}
ll C (ll n, ll m, ll p, ll k) { //求C(n, m) mod p^k
ll pk = power2(p, k);
if (m > n) return 0;
if (m == n) return 1;
if (m == 1) return n;
ll E = e(n, p) - e(m, p) - e(n - m, p);
if (E >= k) return 0;
return f(n, p, k) % pk * get_inv(f(m, p, k), pk) % pk * get_inv(f(n - m, p, k), pk) % pk * power(p, E, pk) % pk;
}
ll crt (ll a, ll p, ll k) {
ll pk = power2(p, k);
ll M_i = mod / pk;
return a * get_inv(M_i, pk) % mod * M_i % mod;
}
ll exLucas (ll n, ll m) {
ll t = mod;
ll res = 0;
for (ll i = 2; i * i <= t; i++) {
if (t % i == 0) {
ll k = 0;
while (t % i == 0) {
t /= i;
k++;
}
ll c = C(n, m, i, k);
res = (res + crt(c, i, k)) % mod;
}
}
if (t > 1) {
ll c = C(n, m, t, 1);
res = (res + crt(c, t, 1)) % mod;
}
return res;
}
main () {
ll n, m;
cin >> n >> m >> mod;
cout << exLucas(n, m);
return 0;
}